Hmm, I tried replacing the x's in the model with their principal component scores, and suddenly everything converges as a greased lightning:
> Z <- princomp(cbind(x1,x2,x3))$scores > Z <- as.data.frame(princomp(cbind(x1,x2,x3))$scores) > names(Z)<- letters[1:3] > optimx(rep(0,8),lnl, hessian=TRUE, y1=y1, y2=y2, + x1=Z$a, x2=Z$b, x3=Z$c, method="BFGS") par 1 0.59107682, 0.01043568, -0.20639153, 0.25902746, 0.24675162, -0.01426477, 0.18045177, -0.23657327 fvalues method fns grs itns conv KKT1 KKT2 xtimes 1 -76.73878 BFGS 157 37 NULL 0 TRUE TRUE 0.055 Warning messages: 1: In max(logpar) : no non-missing arguments to max; returning -Inf 2: In min(logpar) : no non-missing arguments to min; returning Inf The likelihood appears to be spot on, but, obviously, the parameter estimates need to be back-transformed to the original x1,x2,x3 system. I'll leave that as the proverbial "exercise for the reader"... However, the whole thing leads me to a suspicion that maybe our numerical gradients are misbehaving in cases with high local collinearity? -pd On May 9, 2011, at 13:40 , Alex Olssen wrote: > Hi Mike, > > Mike said > "is this it, page 1559?" > > That is the front page yes, page 15*6*9 has the table, of which the > model labelled 18s is the one I replicated. > > "did you actually cut/paste code anywhere and is your first > coefficient -.19 or -.019? > Presumably typos would be one possible problem." > > -0.19 is not a typo, it is pi10 in the paper, and a1 in my Stata > estimation - as far as I can tell cutting and pasting is not the > problem. > > "generally it helps if we could at least see the equations to check your > code against typos ( note page number ?) in lnl that may fix part of the > mystery. Is full text available > on author's site, doesn't come up on citeseer AFAICT," > > Unfortunately I do not know of any place to get the full version of > this paper that doesn't require access to a database such as JSTOR. > The fact that this likelihood function reproduces the published > results in Stata makes me confident that it is correct - I also have > read a lot on systems estimation and it is a pretty standard > likelihood function. > > "I guess one question would be " what is beta" in lnl supposed to be - > it isn't used anywhere but I will also mentioned I'm not that familiar > with the R code ( I'm trying to work through this to learn R and the > optimizers). > > maybe some words would help, is sigma supposed to be 2x2 or 8x8 and what are > e1 and e2 supposed to be?" > > The code above is certainly not as elegant as it could be - in this > case the beta is rather superfluous. It is just a vector to hold > parameters - but since I called all the parameters theta anyway there > is no need for it. e1 and e2 are the residuals from the first and > second equations of the system. Sigma is a 2x2 matrix which is the > outer product of the two vectors of residuals. > > Kind regards, > > Alex > > > > On 9 May 2011 23:12, Mike Marchywka <marchy...@hotmail.com> wrote: >> >> >> >> >> >> >> ---------------------------------------- >>> Date: Mon, 9 May 2011 22:06:38 +1200 >>> From: alex.ols...@gmail.com >>> To: pda...@gmail.com >>> CC: r-help@r-project.org; da...@otter-rsch.com >>> Subject: Re: [R] maximum likelihood convergence reproducing Anderson >>> Blundell 1982 Econometrica R vs Stata >>> >>> Peter said >>> >>> "Ahem! You might get us interested in your problem, but not to the >>> level that we are going to install Stata and Tsp and actually dig out >>> and study the scientific paper you are talking about. Please cite the >>> results and explain the differences." >>> >>> Apologies Peter, will do, >>> >>> The results which I can emulate in Stata but not (yet) in R are reported >>> below. >> >> did you actually cut/paste code anywhere and is your first coefficient -.19 >> or -.019? >> Presumably typos would be one possible problem. >> >>> They come from Econometrica Vol. 50, No. 6 (Nov., 1982), pp. 1569 >> >> >> is this it, page 1559? >> >> http://www.jstor.org/pss/1913396 >> >> generally it helps if we could at least see the equations to check your >> code against typos ( note page number ?) in lnl that may fix part of the >> mystery. Is full text available >> on author's site, doesn't come up on citeseer AFAICT, >> >> >> http://citeseerx.ist.psu.edu/search?q=blundell+1982&sort=ascdate >> >> I guess one question would be " what is beta" in lnl supposed to be - >> it isn't used anywhere but I will also mentioned I'm not that familiar >> with the R code ( I'm trying to work through this to learn R and the >> optimizers). >> >> maybe some words would help, is sigma supposed to be 2x2 or 8x8 and what are >> e1 and e2 supposed to be? >> >> >> >> >>> >>> TABLE II - model 18s >>> >>> coef std err >>> p10 -0.19 0.078 >>> p11 0.220 0.019 >>> p12 -0.148 0.021 >>> p13 -0.072 >>> p20 0.893 0.072 >>> p21 -0.148 >>> p22 0.050 0.035 >>> p23 0.098 >>> >>> The results which I produced in Stata are reported below. >>> I spent the last hour rewriting the code to reproduce this - since I >>> am now at home and not at work :( >>> My results are "identical" to those published. The estimates are for >>> a 3 equation symmetrical singular system. >>> I have not bothered to report symmetrical results and have backed out >>> an extra estimate using adding up constraints. >>> I have also backed out all standard errors using the delta method. >>> >>> . ereturn display >>> ------------------------------------------------------------------------------ >>> | Coef. Std. Err. z P>|z| [95% Conf. Interval] >>> -------------+---------------------------------------------------------------- >>> a | >>> a1 | -.0188115 .0767759 -0.25 0.806 -.1692895 .1316664 >>> a2 | .8926598 .0704068 12.68 0.000 .7546651 1.030655 >>> a3 | .1261517 .0590193 2.14 0.033 .010476 .2418275 >>> -------------+---------------------------------------------------------------- >>> g | >>> g11 | .2199442 .0184075 11.95 0.000 .183866 .2560223 >>> g12 | -.1476856 .0211982 -6.97 0.000 -.1892334 -.1061378 >>> g13 | -.0722586 .0145154 -4.98 0.000 -.1007082 -.0438089 >>> g22 | .0496865 .0348052 1.43 0.153 -.0185305 .1179034 >>> g23 | .0979991 .0174397 5.62 0.000 .0638179 .1321803 >>> g33 | -.0257405 .0113869 -2.26 0.024 -.0480584 -.0034226 >>> ------------------------------------------------------------------------------ >>> >>> In R I cannot get results like this - I think it is probably to do >>> with my inability at using the optimisers well. >>> Any pointers would be appreciated. >>> >>> Peter said "Are we maximizing over the same parameter space? You say >>> that the estimates from the paper gives a log-likelihood of 54.04, but >>> the exact solution clocked in at 76.74, which in my book is rather >>> larger." >>> >>> I meant +54.04 > -76.74. It is quite common to get positive >>> log-likelihoods in these system estimation. >>> >>> Kind regards, >>> >>> Alex >>> >>> On 9 May 2011 19:04, peter dalgaard wrote: >>>> >>>> On May 9, 2011, at 06:07 , Alex Olssen wrote: >>>> >>>>> Thank you all for your input. >>>>> >>>>> Unfortunately my problem is not yet resolved. Before I respond to >>>>> individual comments I make a clarification: >>>>> >>>>> In Stata, using the same likelihood function as above, I can reproduce >>>>> EXACTLY (to 3 decimal places or more, which is exactly considering I >>>>> am using different software) the results from model 8 of the paper. >>>>> >>>>> I take this as an indication that I am using the same likelihood >>>>> function as the authors, and that it does indeed work. >>>>> The reason I am trying to estimate the model in R is because while >>>>> Stata reproduces model 8 perfectly it has convergence >>>>> difficulties for some of the other models. >>>>> >>>>> Peter Dalgaard, >>>>> >>>>> "Better starting values would help. In this case, almost too good >>>>> values are available: >>>>> >>>>> start <- c(coef(lm(y1~x1+x2+x3)), coef(lm(y2~x1+x2+x3))) >>>>> >>>>> which appears to be the _exact_ solution." >>>>> >>>>> Thanks for the suggestion. Using these starting values produces the >>>>> exact estimate that Dave Fournier emailed me. >>>>> If these are the exact solution then why did the author publish >>>>> different answers which are completely reproducible in >>>>> Stata and Tsp? >>>> >>>> >>>> Ahem! You might get us interested in your problem, but not to the level >>>> that we are going to install Stata and Tsp and actually dig out and study >>>> the scientific paper you are talking about. Please cite the results and >>>> explain the differences. >>>> >>>> Are we maximizing over the same parameter space? You say that the >>>> estimates from the paper gives a log-likelihood of 54.04, but the exact >>>> solution clocked in at 76.74, which in my book is rather larger. >>>> >>>> Confused.... >>>> >>>> -p >>>> >>>> >>>>> >>>>> Ravi, >>>>> >>>>> Thanks for introducing optimx to me, I am new to R. I completely >>>>> agree that you can get higher log-likelihood values >>>>> than what those obtained with optim and the starting values suggested >>>>> by Peter. In fact, in Stata, when I reproduce >>>>> the results of model 8 to more than 3 dp I get a log-likelihood of >>>>> 54.039139. >>>>> >>>>> Furthermore if I estimate model 8 without symmetry imposed on the >>>>> system I reproduce the Likelihood Ratio reported >>>>> in the paper to 3 decimal places as well, suggesting that the >>>>> log-likelihoods I am reporting differ from those in the paper >>>>> only due to a constant. >>>>> >>>>> Thanks for your comments, >>>>> >>>>> I am still highly interested in knowing why the results of the >>>>> optimisation in R are so different to those in Stata? >>>>> >>>>> I might try making my convergence requirements more stringent. >>>>> >>>>> Kind regards, >>>>> >>>>> Alex >>>> >>>> -- >>>> Peter Dalgaard >>>> Center for Statistics, Copenhagen Business School >>>> Solbjerg Plads 3, 2000 Frederiksberg, Denmark >>>> Phone: (+45)38153501 >>>> Email: pd....@cbs.dk Priv: pda...@gmail.com >>>> >>>> >>> >>> ______________________________________________ >>> R-help@r-project.org mailing list >>> https://stat.ethz.ch/mailman/listinfo/r-help >>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html >>> and provide commented, minimal, self-contained, reproducible code. >> > > ______________________________________________ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. -- Peter Dalgaard Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd....@cbs.dk Priv: pda...@gmail.com ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.