[R] hessian in solnp

[CAMARDA Carlo Giovanni]

In the commented snippet below, a behavior of the solnp() I cannot fully explain: when a constraint is added, hessian matrix is obviously changing, but in a way I don't understand. I know that the model below could be estimated differently and in a much efficient way. However I wanted to present the simplest example I could think. 

In particular, I encountered this issue because I aim to analytically compute variance-covariance matrix in a (non-linear) model with equality constraints. Any practical hint to find a solution to this more general issue would be also welcome. 

Thanks, GC 




## simulating a simple linear model 
n <- 100 
x <- sort(runif(n)) 
X <- cbind(1,x) 
k <- ncol(X) 
y <- X%*%c(3,4) + rnorm(n, sd=0.5) 
## estimating with lm() 
fitlm <- lm(y~x) 


## using solnp 
library(Rsolnp) 
## objecting function 
fn1 <- function(par){ 
y.hat <- X%*%par 
res <- y-y.hat 
RSS <- t(res)%*%res/2 
return(c(RSS)) 
} 
## using solnp without any constraint 
solUR <- solnp(c(3.2, 3.5), fun = fn1) 
## getting the same fitted objects as in lm() 
## coefficients 
cbind(solUR$pars, fitlm$coef) 
## and standard errors by variance-covariance matrix computed 
## with hessian internally estimated 
(H.UR <- solUR$hessian) 
y.UR <- X%*%solUR$pars 
RSS.UR <- t(y-y.UR)%*%(y-y.UR) 
sigma2.UR <- RSS.UR/(n-k) 
V.UR <- c(sigma2.UR)*solve(H.UR) 
se.UR <- sqrt(diag(V.UR)) 
cbind(se.UR,summary(fitlm)$coef[,2]) 


## adding equality constraint 
eqn1 <- function(par) sum(par) 
c <- sum(solUR$pars) 
## I use the sum of the previously estimated coefficients 
solR <- solnp(c(3.2, 3.5), fun = fn1, eqfun = eqn1, eqB = c) 
## as expected we get the same coefficients 
cbind(solR$pars,fitlm$coef) 
## but the hessian is rather different 
(H.R <- solR$hessian) 
## which leads to completely different (much larger here) 
## standard errors 
y.R <- X%*%solR$pars 
RSS.R <- t(y-y.R)%*%(y-y.R) 
sigma2.R <- RSS.R/(n-k) 
V.R <- c(sigma2.R)*solve(H.R) 
se.R <- sqrt(diag(V.R)) 
cbind(se.R,summary(fitlm)$coef[,2]) 
## I noticed that differences in the hessian 
## is ~constant over both dimensions 
## and it depends (also) on sample size 
H.R-H.UR 




	[[alternative HTML version deleted]]