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.
