Srinivasan Parthasarathy
2008-Sep-16 18:40 UTC
[R] Maximum likelihood estimation of a truncated regression model
Hi, I have a quick question regarding estimation of a truncation regression model (truncated above at 1) using MLE in R. I will be most grateful to you if you can help me out. The model is linear and the relationship is "dhat = bhat0+Z*bhat+e", where dhat is the dependent variable >0 and upper truncated at 1; bhat0 is the intercept; Z is the independent variable and is a uniform random variate between 0 and 1; and e is the error term. I realised that R doesn't have a built-in function to estimate truncated regression models as does STATA, LIMDEP etc. I tried the survival and FEAR packages and couldn't fit it for my case. So I wrote the log likelihood function of the truncated regression model and reparametrised it using Olsen (1978) so that the function is globally concave and has an unique maximiser. I used a quasi-Newton method (BFGS) to maximise my function. I also used Newton-Raphson method (maxNR) to maximise my function. The (naive) code can be seen below. sw1<-function(theta,dhat,z) { gamma<-theta[1:2] eta<-theta[3] d1<-dhat*eta-z%*%gamma d2<- eta-z%*%gamma logl<- (-n*0.5*log(2*pi))+(n*log(eta))-(0.5*sum(d1*d1))-(sum(log(pnorm(d2)))) return(-logl) } try(j<-optim(c(a,b,c),sw1,method="BFGS",dhat=dhat,z=z),silent=TRUE) q[i,]=abs(j$par[2]/j$par[3]) I am trying to iterate the above procedure 2000 times and compute the RMSE of the estimated bhat. In the above code, a, b, and c are the intercept, bhat and 1/sigma estimates from OLS of the linear model [I get the standard error estimate from OLS and obtain the sigma by multiplying it by sqrt(n) where n=200]. I realised that in many of the iterations, the above code doesn't optimize as the bhat estimate I get from the above code is worser than the one I get from STATA although optim (or maxNR) says the function has converged successfully. I realise that the starting parameter values could be a reason for these sub-optimal solutions. 1. Is there an optimization routine that doesn't need starting parameter values? 2. Can you please tell me what could be the possible error in the above code? 3. Is there a routine built in R that could do truncated regression for my situation? Thank you very much for your help. Best Regards, Srini