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
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
athttp://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
[[alternative HTML version deleted]]