\*********************************************************************************/
\* 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;