#################################################################################
# This is a set of test routines for the sechol.s function.  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 Gauss format.							#
#################################################################################

# FIRST EXAMPLE
X <- matrix(c(5,2,3,2.95,2,1,2,1,5,2,3,3),ncol=3)
Y <- c(9,11,-5,-2)
b.hat <- solve(t(X)%*%X)%*%t(X)%*%Y
Y.hat <- X%*%b.hat
e <- Y - Y.hat
sigma <- as.numeric(t(e)%*%e)
b.var <- sigma/(nrow(X)-ncol(X)) * solve(t(X)%*%X)
b.cor <- b.var
for (i in 1:nrow(b.cor))  {
    for (j in 1:ncol(b.cor))  {
	if (i != j)
	    b.cor[i,j] <- b.cor[i,j]/sqrt(b.var[i,i]*b.var[j,j])
	else
	    b.cor[i,j] <- sqrt(b.var[i,j])
    }
}

# SECOND EXAMPLE
X <- matrix(c(5,2,3,2.99,2,1,2,1,5,2,3,3),ncol=3)
Y <- c(9,11,-5,-2)
b.hat <- solve(t(X)%*%X)%*%t(X)%*%Y
Y.hat <- X%*%b.hat
e <- Y - Y.hat
sigma <- as.numeric(t(e)%*%e)
b.var <- sigma/(nrow(X)-ncol(X)) * solve(t(X)%*%X)
b.cor <- b.var
for (i in 1:nrow(b.cor))  {
    for (j in 1:ncol(b.cor))  {
	if (i != j)
	    b.cor[i,j] <- b.cor[i,j]/sqrt(b.var[i,i]*b.var[j,j])
	else
	    b.cor[i,j] <- sqrt(b.var[i,j])
    }
}

ginv <- function(X, tol=sqrt(.Machine$double.eps))  {
	s <- svd(X)
	nz <- s$d > tol * s$d[1]
	if(any(nz)) s$v[,nz] %*% (t(s$u[, nz])/s$d[nz]) else X*0
}
mpinv <- function(X, tol = sqrt(.Machine$double.eps)) {
	s <- svd(X)
	e <- s$d
	e[e > tol] <- 1/e[e > tol]
	s$v %*% diag(e) %*% t(s$u))
}

# THIRD EXAMPLE
X <- matrix(c(5,2,3,2.999,2,1,2,1,5,2,3,3),ncol=3)
Y <- c(9,11,-5,-2)
b.hat <- mpinv(t(X)%*%X)%*%t(X)%*%Y
Y.hat <- X%*%b.hat
e <- Y - Y.hat
sigma <- as.numeric(t(e)%*%e)
b.var <- sigma * mpinv(t(X)%*%X)
b.cor <- b.var
for (i in 1:nrow(b.cor))  {
    for (j in 1:ncol(b.cor))  
            b.cor[i,j] <- b.cor[i,j]/sqrt(b.var[i,i]*b.var[j,j])
}
diag(b.cor) <- sqrt(diag(b.var))

# FIRST CHOLESKY EXAMPLE
S <- matrix(c(2,0,2.4,0,2,0,2.4,0,3),ncol=3)
chol(S)
T <- sechol(S)
t(T) %*% T

# SECOND CHOLESKY EXAMPLE
S <- matrix(c(2,0,2.5,0,2,0,2.5,0,3),ncol=3)
sechol(S)

# THIRD CHOLESKY EXAMPLE
S <- matrix(c(2,0,10,0,2,0,10,0,3),ncol=3)
T <- sechol(S)
t(T) %*% T