[R] Error in maximum likelihood estimation.

[Dong-hyun Oh]

Dear UseRs,

I wrote the following function to use MLE.

---------------------------------------------
mlog <- function(theta, nx = 1, nz = 1, dt){
   beta <- matrix(theta[1:(nx+1)], ncol = 1)
   delta <- matrix(theta[(nx+2):(nx+nz+1)], ncol = 1)
   sigma2 <- theta[nx+nz+2]
   gamma <- theta[nx+nz+3]
   y <- as.matrix(dt[, 1], ncol = 1)
   x <- as.matrix(data.frame(1, as.matrix(dt[, 2:(nx+1)], ncol = 2)))
   z <- as.matrix(dt[, (nx+2):(nx+nz+1)], ncol = nz)

   d <- z %*% delta / (gamma * sigma2)^.5
   mustar <- (1-gamma) * z %*% delta - gamma * ( y - x %*% beta)
   sigmastar <- (gamma * (1-gamma) * sigma2)^.5
   dstar <- mustar / sigmastar

   loglik <- (-0.5 * nrow(x) *(log(2*pi) + log(sigma2))
              -0.5 * sum(( y - x %*% beta + z %*% delta)^2/sigma2)
              -sum(log(pnorm(d))) + sum(log(pnorm(dstar))))
   return(-loglik)
}
-----------------------------------------------

Loglikelihood function is from page 21of Battese and Coelli (1993).  
(You can download this article at http://www.une.edu.au/economics/publications/econometrics/emwp69.PDF 
  )

To test the above function with an artificial data set, I created the  
following data.frame.

-----------------------------------------------
x1 <- abs(rnorm(100))*100
x2 <- abs(rnorm(100))*10
z1 <- abs(rnorm(100))*5
z2 <- abs(rnorm(100))*7
y <- abs(0.3 + 0.3* log(x1) + 0.7* log(x2))
dat <- data.frame(log(y), log(x1), log(x2), z1, z2)
------------------------------------------------

The starting value I set up is as follows:
-----------------------------------------------
theta.start <- c(0.09, 0.008, 0.008, 0.023, 0.0008, 0.008, 0.008)
----------------------------------------------

Applying nlm() function to the above function with the starting values  
gives the following error message.

------------------------------------------------
out <- nlm(mlog, theta.start, nx = 2 , nz = 2, dt = dat, print.level =  
2, hessian\
  = T, iterlim = 500)

iteration = 0
Step:
[1] 0 0 0 0 0 0 0
Parameter:
[1] 0.0900 0.0080 0.0080 0.0230 0.0008 0.0080 0.0080
Function Value
[1] 6116.2
Gradient:
[1]  -10704.8  -46465.7  -22536.4   47168.2   54542.9 -776512.5      
579.9

Error in nlm(mlog, theta.start, nx = 2, nz = 2, dt = dat, print.level  
= 2,  :
   (converted from warning) NA/Inf replaced by maximum positive value

------------------------------------------------

My questions are
1. Is my loglikelihood function set up appropriately?
2. How to set up starting values efficiently? I read a thread shown at https://stat.ethz.ch/pipermail/r-help/2005-September/079617.html 
  , but I cannot figure out how to set up starting values with  
expand.grid() of parameters I used (beta, delta, gamma, sigma2).

Thank you in advance.

Sincerely,
Dong-hyun Oh