\*********************************************************************************/
\* This is the Gill/Murray cholesky routine.  Reference:		         */  
\*										 */
\* Gill, Jeff and Gary King. ``What to do When Your Hessian is Not Invertible:   */
\* Alternatives to Model Respecification in Nonlinear Estimation,'' Sociological */
\* Methods and Research, Vol. 32, No. 1 (2004): Pp. 54--87. 		         */
\*										 */
\* Abstract									 */
\* 									         */
\* What should a researcher do when statistical analysis software terminates 	 */
\* before completion with a message that the Hessian is not invertable? The 	 */
\* standard textbook advice is to respecify the model, but this is another way 	 */
\* of saying that the researcher should change the question being asked. 	 */
\* Obviously, however, computer programs should not be in the business of 	 */
\* deciding what questions are worthy of study. Although noninvertable 		 */
\* Hessians are sometimes signals of poorly posed questions, nonsensical 	 */
\* models, or inappropriate estimators, they also frequently occur when 	 */
\* information about the quantities of interest exists in the data, through 	 */
\* the likelihood function. We explain the problem in some detail and lay out 	 */
\* two preliminary proposals for ways of dealing with noninvertable Hessians 	 */
\* without changing the question asked. 					 */
\* 										 */
\* Also available is the software to implement the procedure described in this 	 */
\* paper in R format.    							 */
\*********************************************************************************/

proc gmchol(A);
   /* calculates the Gill-Murray generalized choleski decomposition */
   /* input matrix A must be non-singular and symmetric */
   /* Author: Jeff Gill. Part of the Hessian Project. */
   local i,j,k,n,sum,R,theta_j,norm_A,phi_j,delta,xi_j,gamm,E,beta_j;
   n = rows(A);
   R = eye(n);
   E = zeros(n,n);
   norm_A = maxc(sumc(abs(A)));
   gamm = maxc(abs(diag(A))); 
   delta = maxc(maxc(__macheps*norm_A~__macheps));
   for j (1, n, 1); 
      theta_j = 0;
      for i (1, n, 1);
	 sum = 0;
	 for k (1, (i-1), 1);	
	    sum = sum + R[k,i]*R[k,j];
	 endfor;
	 R[i,j] = (A[i,j] - sum)/R[i,i];
	 if (A[i,j] -sum) > theta_j;
	    theta_j = A[i,j] - sum;
         endif;
	 if i > j;
	    R[i,j] = 0;
         endif;
      endfor;
      sum = 0;
      for k (1, (j-1), 1);	
         sum = sum + R[k,j]^2;
      endfor;
      phi_j = A[j,j] - sum;
      if (j+1) <= n;
         xi_j = maxc(abs(A[(j+1):n,j]));
      else;
         xi_j = maxc(abs(A[n,j]));
      endif;
      beta_j = sqrt(maxc(maxc(gamm~(xi_j/n)~__macheps)));
      if delta >= maxc(abs(phi_j)~((theta_j^2)/(beta_j^2)));
	 E[j,j] = delta - phi_j;
      elseif abs(phi_j) >= maxc(((delta^2)/(beta_j^2))~delta);
	 E[j,j] = abs(phi_j) - phi_j;
      elseif ((theta_j^2)/(beta_j^2)) >= maxc(delta~abs(phi_j));
	 E[j,j] = ((theta_j^2)/(beta_j^2)) - phi_j;
      endif;
      R[j,j] = sqrt(A[j,j] - sum + E[j,j]);
   endfor;
   retp(R'R);
endp;