tag:blogger.com,1999:blog-6288862798546085706.post4495689915931072556..comments2020-08-07T15:25:35.217-04:00Comments on Econometrics By Simulation: Estimating Random Coefficients on XFrancishttp://www.blogger.com/profile/16658586705916884436noreply@blogger.comBlogger5125tag:blogger.com,1999:blog-6288862798546085706.post-59042265226693606122012-11-15T23:47:43.727-05:002012-11-15T23:47:43.727-05:00On second glance, you are right. I have not inclu...On second glance, you are right. I have not included the x2 term and the x term properly in myNormalReg2Francishttps://www.blogger.com/profile/16658586705916884436noreply@blogger.comtag:blogger.com,1999:blog-6288862798546085706.post-84545984289260105942012-11-15T13:34:53.708-05:002012-11-15T13:34:53.708-05:00I have a few more thoughts on the matter.
I don&#...I have a few more thoughts on the matter.<br /><br />I don't think IFGLS is more efficient (I know you did not say that but I suggested it) and as for the MLE being GLS, that is not the case. <br /><br />From Wooldridge's MIT book Chapter 7 he says that GLS is used when the variance-covariance matrix is "known". When it is not known, it can be estimated via methods as you purpose, FGLS. MLE is related. It just states that the variance is unknown and the coefficients are unknown and we are going to use the normality assumption to estimate them together.Francishttps://www.blogger.com/profile/16658586705916884436noreply@blogger.comtag:blogger.com,1999:blog-6288862798546085706.post-29334210691368039182012-11-15T12:34:48.806-05:002012-11-15T12:34:48.806-05:00Hi Szabolcs Lorincz,
Thanks for this thorough com...Hi Szabolcs Lorincz,<br /><br />Thanks for this thorough comment. <br /><br />1. I have run your code and it definitely seems to be working in a superior fashion to the MLE. For some reason I thought that the MLE estimator was supposedly more efficient than the FGLS estimator. However, I seem to be wrong. I will follow up this post with another post simulating the results of both methods to see which one is less biased and has lower variance. It seems clear to me that the FGLS estimator definitely is faster.<br /><br />2. I think that is a good point about the correlation between the random coefficients of variables.<br /><br />3. I don't think I am mistaken about including the x and x2 term in myNormalReg2. I included it in the MLE ly syntax in order to maximize compatibility of the code.<br /><br />ml model lf myNormalReg2 (y = x) (var_u: ) (var_v: x2, noconstant) (cov_uv: x, noconstant)<br /><br />By specifying directly x2 and x into the above equations it tells Stata to include them in the estimation. The above must be specified or else the MLE would not be able to converge.Francishttps://www.blogger.com/profile/16658586705916884436noreply@blogger.comtag:blogger.com,1999:blog-6288862798546085706.post-29691551098408125492012-11-13T08:41:00.209-05:002012-11-13T08:41:00.209-05:00Second part of Stata code (ML vs. IFGLS)
/* itera...Second part of Stata code (ML vs. IFGLS)<br /><br />/* iterated feasible GLS estimation */<br />matrix Chat_IFGLS=J(1,1,0)<br />matrix dif_IFGLS=J(1,1,0)<br />timer clear 1<br />timer on 1<br />local quiet quietly <br />capture drop weight<br />generate weight=1<br />local x=0<br />local _cons=0<br />local varv=0<br />local varu=0<br />local cov_uv=0<br />local difference=1<br />local tolerance=0.0000001<br />local maxit=1000<br />local i=0<br />while (`difference'>`tolerance' & `i'<`maxit') {<br /> local i=`i'+1<br /> `quiet' reg y x [aweight=weight]<br /> local dx=abs(`x'-_b["x"])<br /> local d_cons=abs(`_cons'-_b["_cons"])<br /> local x=_b["x"]<br /> local _cons=_b["_cons"]<br /> capture drop res<br /> predict res, res<br /> capture drop res2<br /> generate res2=res^2<br /> `quiet' reg res2 x2 x varu, noconstant<br /> local dvarv=abs(`varv'-_b["x2"])<br /> local dvaru=abs(`varu'-_b["varu"])<br /> local dcov_uv=abs(`cov_uv'-_b["x"])<br /> local varv=_b["x2"]<br /> local varu=_b["varu"]<br /> local cov_uv=_b["x"]<br /> local difference=max(`dx', `d_cons', `dvarv', `dvaru', `dcov_uv')<br /> capture drop fit<br /> predict fit, xb<br /> `quiet' replace weight=1/(fit^0.5) if fit>0 & fit!=.<br /> quietly count if fit<=0 | fit==.<br /> if (r(N)==0) {<br /> display "FGLS iteration `i', max. abs. difference in parameters: " %7.5f = `difference'<br /> }<br /> if (r(N)>0) {<br /> display "FGLS iteration `i', max. abs. difference in parameters: " %7.5f = `difference' ", # of obs. with inconsistent weight: " r(N) " (" 100*r(N)/_N "%)"<br /> }<br />}<br />reg y x [aweight=weight]<br />local cov_uv=`cov_uv'/2<br />matrix Chat_IFGLS = (`varv', `cov_uv' \ `cov_uv', `varu')<br />matrix dif_IFGLS=Chat_IFGLS-C<br />timer off 1<br />quietly timer list 1<br />local t1=r(t1)<br />local timerm=floor(`t1'/60)<br />local timers=round(`t1'-floor(`t1'/60)*60)<br />display "`timerm'm `timers's"<br /><br />display "True variance-covariance matrix"<br />matrix list C<br />display "ML estimate"<br />matrix list Chat_ML<br />display "IFGLS estimate"<br />matrix list Chat_IFGLS<br />display "Difference of ML estimate from true VC matrix"<br />matrix list dif_ML<br />display "Difference of IFGLS estimate from true VC matrix"<br />matrix list dif_IFGLSAnonymousnoreply@blogger.comtag:blogger.com,1999:blog-6288862798546085706.post-84662454979561392632012-11-13T08:38:36.812-05:002012-11-13T08:38:36.812-05:00Hi Francis,
congratulations to your blog - it'...Hi Francis,<br /><br />congratulations to your blog - it's very interesting and useful. I'd like to comment on this post.<br /><br />1. I think what you have formulated is a Generalized Least Squares model of heteroscedasticity. GLS can be estimated not only by Maximum Likelihood but also by Feasible GLS which starts by OLS and then estimates the error variance from the OLS residuals and estimates a weighted least squares model. When iterated until convergence, FGLS is ML. I copied below a code that compares your ML estimator with the IFGLS method (only the first part can go into this post, I'll copy the second part into the next, just put them into one Stata do file). The advantage of the latter is that it's a (sequence of) linear estimator(s) and it's probably faster. Note that I tried to deal with the possible non-PSD issue by not updating the weighting vector for the problematic observations. (Though this issue never arises in the example - probably a further advantage of the IFGLS.)<br /><br />2. I think in terms of cross-variances it would be more interesting to think of correlations between random terms on coefficients of different observed regressors. Similarly to, say, discrete choice models like BLP, it might be of interest to see if "preferences" for different "characteristics" (different columns of x when there are more regressors) are correlated or not. This is probably more relevant than the correlation between the random terms and the unobserved error term, as it can tell something about the relationship between the observed variables. Of course, I guess identification depends on whether one has enough variation in those variables.<br /><br />3. It seems to me that there is a problem in your second estimator (myNormalReg2). In the line you define the variance (gen `sigma2' = (exp(`var_u') + exp(`var_v') + 2*`cov_uv')), the expression is not a function of x and x-square.<br /><br />I hope you find my comments useful.<br /><br />Cheers,<br /><br />Szabolcs Lorincz<br /><br />************************************************<br /><br />clear<br />clear matrix<br />clear mata<br />set more off<br />set mem 3g<br />set obs 100000<br />set seed 121<br /><br />gen x = runiform()*5<br />matrix C = (9, 2.5 \ 2.5, 23)<br />drawnorm v u, cov(C)<br />gen y = 4*x + v*x + u<br />gen x2 = x^2<br />gen varu = 1<br /><br />/* ML estimation */<br />matrix Chat_ML=J(1,1,0)<br />matrix dif_ML=J(1,1,0)<br />timer clear 1<br />timer on 1<br />cap program drop myNormalReg<br />program define myNormalReg<br /> args lnlk xb sigma2<br /> qui replace `lnlk' = -ln(sqrt(`sigma2') * sqrt(2*_pi)) - ($ML_y-`xb')^2/(2*`sigma2')<br />end<br />ml model lf myNormalReg (y = x) (sigma2: varu x2 x, noconstant)<br />ml maximize<br />local varu = [sigma2]_b[varu]<br />local varv = [sigma2]_b[x2]<br />local cov_uv = [sigma2]_b[x]/2<br />matrix Chat_ML = (`varv', `cov_uv' \ `cov_uv', `varu')<br />matrix dif_ML=Chat_ML-C<br />timer off 1<br />quietly timer list 1<br />local t1=r(t1)<br />local timerm=floor(`t1'/60)<br />local timers=round(`t1'-floor(`t1'/60)*60)<br />display "`timerm'm `timers's"Anonymousnoreply@blogger.com