Hello,
I'm using optim to program a set of mle regression procedures for non-normal
disturbances. This is for teaching and expository purposes only. I've
successfully programmed the normal, generalized gamma, gamma, weibull,
exponential, and lognormal regression functions. And optim returns
reasonable answers for all of these compared with the identical optimization
problems in STATA and LIMDEP. However, I'm having a problem with beta
regression. The parameterization of beta regression that I am using solves
for the second distributional parameter in terms of the first and xbeta.
Therefore, there is only a single distribution parameter. This
parameterization works fine in STATA, LIMDEP, and SHAZAM. However, the
results returned by optim in R are not reasonable in terms of the value of
the log likelihood and parameter estimates.
Here is my code. Does anyone see a problem with what I'm doing here? Any
advice would be appreciated.
Thanks.
B. Dan Wood
# Define the function to be optimized
llik.beta <- function(par,X,YP) {
X <- as.matrix(x)
YP <- as.vector(y)
xbeta <- X %*% par[1:K]
p <- par[K1:K1]
sum(
-lgamma(p)
+lgamma(p+(p/xbeta-p))
-lgamma(p/xbeta-p)
+(p-1)*log(YP)
+log(1-YP)*(p/xbeta-p-1)
)
}
llik.beta
# Now use the above function to estimate the model. First, create a set of
reasonable start values.
startvalues <- c(coefficients(ols.model),74)
startvalues
# Now call optim
mod.beta <- optim(startvalues,llik.beta, method = "BFGS", control
list(trace,
maxit=1000,fnscale = -1), hessian = TRUE)
mod.beta
B. Dan Wood, Professor and University Faculty Fellow
Texas A&M University
Department of Political Science
4348 TAMU
College Station, TX 77843-4348
Office: (979) 845-1610
Home: (979) 690-0390
Voice Mail: (979) 492-7599
FAX: (979) 847-8924
<http://www-polisci.tamu.edu/bdanwood
<http://www-polisci.tamu.edu/bdanwood>>
[[alternative HTML version deleted]]
