Alex Olssen
2011-May-06 12:29 UTC
[R] maximum likelihood convergence reproducing Anderson Blundell 1982 Econometrica R vs Stata
Dear R-help,
I am trying to reproduce some results presented in a paper by Anderson
and Blundell in 1982 in Econometrica using R.
The estimation I want to reproduce concerns maximum likelihood
estimation of a singular equation system.
I can estimate the static model successfully in Stata but for the
dynamic models I have difficulty getting convergence.
My R program which uses the same likelihood function as in Stata has
convergence properties even for the static case.
I have copied my R program and the data below. I realise the code
could be made more elegant - but it is short enough.
Any ideas would be highly appreciated.
## model 18
lnl <- function(theta,y1, y2, x1, x2, x3) {
n <- length(y1)
beta <- theta[1:8]
e1 <- y1 - theta[1] - theta[2]*x1 - theta[3]*x2 - theta[4]*x3
e2 <- y2 - theta[5] - theta[6]*x1 - theta[7]*x2 - theta[8]*x3
e <- cbind(e1, e2)
sigma <- t(e)%*%e
logl <- -1*n/2*(2*(1+log(2*pi)) + log(det(sigma)))
return(-logl)
}
p <- optim(0*c(1:8), lnl, method="BFGS", hessian=TRUE, y1=y1,
y2=y2,
x1=x1, x2=x2, x3=x3)
"year","y1","y2","x1","x2","x3"
1929,0.554779,0.266051,9.87415,8.60371,3.75673
1930,0.516336,0.297473,9.68621,8.50492,3.80692
1931,0.508201,0.324199,9.4701,8.27596,3.80437
1932,0.500482,0.33958,9.24692,7.99221,3.76251
1933,0.501695,0.276974,9.35356,7.98968,3.69071
1934,0.591426,0.287008,9.42084,8.0362,3.63564
1935,0.565047,0.244096,9.53972,8.15803,3.59285
1936,0.605954,0.239187,9.6914,8.32009,3.56678
1937,0.620161,0.218232,9.76817,8.42001,3.57381
1938,0.592091,0.243161,9.51295,8.19771,3.6024
1939,0.613115,0.217042,9.68047,8.30987,3.58147
1940,0.632455,0.215269,9.78417,8.49624,3.57744
1941,0.663139,0.184409,10.0606,8.69868,3.6095
1942,0.698179,0.164348,10.2892,8.84523,3.66664
1943,0.70459,0.146865,10.4731,8.93024,3.65388
1944,0.694067,0.161722,10.4465,8.96044,3.62434
1945,0.674668,0.197231,10.279,8.82522,3.61489
1946,0.635916,0.204232,10.1536,8.77547,3.67562
1947,0.642855,0.187224,10.2053,8.77481,3.82632
1948,0.641063,0.186566,10.2227,8.83821,3.96038
1949,0.646317,0.203646,10.1127,8.82364,4.0447
1950,0.645476,0.187497,10.2067,8.84161,4.08128
1951,0.63803,0.197361,10.2773,8.9401,4.10951
1952,0.634626,0.209992,10.283,9.01603,4.1693
1953,0.631144,0.219287,10.3217,9.06317,4.21727
1954,0.593088,0.235335,10.2101,9.05664,4.2567
1955,0.60736,0.227035,10.272,9.07566,4.29193
1956,0.607204,0.246631,10.2743,9.12407,4.32252
1957,0.586994,0.256784,10.2396,9.1588,4.37792
1958,0.548281,0.271022,10.1248,9.14025,4.42641
1959,0.553401,0.261815,10.2012,9.1598,4.4346
1960,0.552105,0.275137,10.1846,9.19297,4.43173
1961,0.544133,0.280783,10.1479,9.19533,4.44407
1962,0.55382,0.281286,10.197,9.21544,4.45074
1963,0.549951,0.28303,10.2036,9.22841,4.46403
1964,0.547204,0.291287,10.2271,9.23954,4.48447
1965,0.55511,0.281313,10.2882,9.26531,4.52057
1966,0.558182,0.280151,10.353,9.31675,4.58156
1967,0.545735,0.294385,10.3351,9.35382,4.65983
1968,0.538964,0.294593,10.3525,9.38361,4.71804
1969,0.542764,0.299927,10.3676,9.40725,4.76329
1970,0.534595,0.315319,10.2968,9.39139,4.81136
1971,0.545591,0.315828,10.2592,9.34121,4.84082
peter dalgaard
2011-May-07 08:46 UTC
[R] maximum likelihood convergence reproducing Anderson Blundell 1982 Econometrica R vs Stata
On May 6, 2011, at 14:29 , Alex Olssen wrote:> Dear R-help, > > I am trying to reproduce some results presented in a paper by Anderson > and Blundell in 1982 in Econometrica using R. > The estimation I want to reproduce concerns maximum likelihood > estimation of a singular equation system. > I can estimate the static model successfully in Stata but for the > dynamic models I have difficulty getting convergence. > My R program which uses the same likelihood function as in Stata has > convergence properties even for the static case. > > I have copied my R program and the data below. I realise the code > could be made more elegant - but it is short enough. > > Any ideas would be highly appreciated.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. Apart from that, it seems that the conjugate gradient methods have difficulties with this likelihood, for some less than obvious reason. Increasing the maxit gets you closer but still not satisfactory. I would suggest trying out the experimental optimx package. Apparently, some of the algorithms in there are much better at handling this likelihood, notably "nlm" and "nlminb".> > ## model 18 > lnl <- function(theta,y1, y2, x1, x2, x3) { > n <- length(y1) > beta <- theta[1:8] > e1 <- y1 - theta[1] - theta[2]*x1 - theta[3]*x2 - theta[4]*x3 > e2 <- y2 - theta[5] - theta[6]*x1 - theta[7]*x2 - theta[8]*x3 > e <- cbind(e1, e2) > sigma <- t(e)%*%e > logl <- -1*n/2*(2*(1+log(2*pi)) + log(det(sigma))) > return(-logl) > } > p <- optim(0*c(1:8), lnl, method="BFGS", hessian=TRUE, y1=y1, y2=y2, > x1=x1, x2=x2, x3=x3) > > "year","y1","y2","x1","x2","x3" > 1929,0.554779,0.266051,9.87415,8.60371,3.75673 > 1930,0.516336,0.297473,9.68621,8.50492,3.80692 > 1931,0.508201,0.324199,9.4701,8.27596,3.80437 > 1932,0.500482,0.33958,9.24692,7.99221,3.76251 > 1933,0.501695,0.276974,9.35356,7.98968,3.69071 > 1934,0.591426,0.287008,9.42084,8.0362,3.63564 > 1935,0.565047,0.244096,9.53972,8.15803,3.59285 > 1936,0.605954,0.239187,9.6914,8.32009,3.56678 > 1937,0.620161,0.218232,9.76817,8.42001,3.57381 > 1938,0.592091,0.243161,9.51295,8.19771,3.6024 > 1939,0.613115,0.217042,9.68047,8.30987,3.58147 > 1940,0.632455,0.215269,9.78417,8.49624,3.57744 > 1941,0.663139,0.184409,10.0606,8.69868,3.6095 > 1942,0.698179,0.164348,10.2892,8.84523,3.66664 > 1943,0.70459,0.146865,10.4731,8.93024,3.65388 > 1944,0.694067,0.161722,10.4465,8.96044,3.62434 > 1945,0.674668,0.197231,10.279,8.82522,3.61489 > 1946,0.635916,0.204232,10.1536,8.77547,3.67562 > 1947,0.642855,0.187224,10.2053,8.77481,3.82632 > 1948,0.641063,0.186566,10.2227,8.83821,3.96038 > 1949,0.646317,0.203646,10.1127,8.82364,4.0447 > 1950,0.645476,0.187497,10.2067,8.84161,4.08128 > 1951,0.63803,0.197361,10.2773,8.9401,4.10951 > 1952,0.634626,0.209992,10.283,9.01603,4.1693 > 1953,0.631144,0.219287,10.3217,9.06317,4.21727 > 1954,0.593088,0.235335,10.2101,9.05664,4.2567 > 1955,0.60736,0.227035,10.272,9.07566,4.29193 > 1956,0.607204,0.246631,10.2743,9.12407,4.32252 > 1957,0.586994,0.256784,10.2396,9.1588,4.37792 > 1958,0.548281,0.271022,10.1248,9.14025,4.42641 > 1959,0.553401,0.261815,10.2012,9.1598,4.4346 > 1960,0.552105,0.275137,10.1846,9.19297,4.43173 > 1961,0.544133,0.280783,10.1479,9.19533,4.44407 > 1962,0.55382,0.281286,10.197,9.21544,4.45074 > 1963,0.549951,0.28303,10.2036,9.22841,4.46403 > 1964,0.547204,0.291287,10.2271,9.23954,4.48447 > 1965,0.55511,0.281313,10.2882,9.26531,4.52057 > 1966,0.558182,0.280151,10.353,9.31675,4.58156 > 1967,0.545735,0.294385,10.3351,9.35382,4.65983 > 1968,0.538964,0.294593,10.3525,9.38361,4.71804 > 1969,0.542764,0.299927,10.3676,9.40725,4.76329 > 1970,0.534595,0.315319,10.2968,9.39139,4.81136 > 1971,0.545591,0.315828,10.2592,9.34121,4.84082 > > ______________________________________________ > R-help at 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.mes at cbs.dk Priv: PDalgd at gmail.com
Ravi Varadhan
2011-May-07 15:51 UTC
[R] maximum likelihood convergence reproducing Anderson Blundell 1982 Econometrica R vs Stata
There is something strange in this problem. I think the log-likelihood is
incorrect. See the results below from "optimx". You can get much
larger log-likelihood values than for the exact solution that Peter provided.
## model 18
lnl <- function(theta,y1, y2, x1, x2, x3) {
n <- length(y1)
beta <- theta[1:8]
e1 <- y1 - theta[1] - theta[2]*x1 - theta[3]*x2 - theta[4]*x3
e2 <- y2 - theta[5] - theta[6]*x1 - theta[7]*x2 - theta[8]*x3
e <- cbind(e1, e2)
sigma <- t(e)%*%e
logl <- -1*n/2*(2*(1+log(2*pi)) + log(det(sigma))) # it looks like there
is something wrong here
return(-logl)
}
data <- read.table("e:/computing/optimx_example.dat", header=TRUE,
sep=",")
attach(data)
require(optimx)
start <- c(coef(lm(y1~x1+x2+x3)), coef(lm(y2~x1+x2+x3)))
# the warnings can be safely ignored in the "optimx" calls
p1 <- optimx(start, lnl, hessian=TRUE, y1=y1, y2=y2,
+ x1=x1, x2=x2, x3=x3, control=list(all.methods=TRUE, maxit=1500))
p2 <- optimx(rep(0,8), lnl, hessian=TRUE, y1=y1, y2=y2,
+ x1=x1, x2=x2, x3=x3, control=list(all.methods=TRUE, maxit=1500))
p3 <- optimx(rep(0.5,8), lnl, hessian=TRUE, y1=y1, y2=y2,
+ x1=x1, x2=x2, x3=x3, control=list(all.methods=TRUE, maxit=1500))
Ravi.
________________________________________
From: r-help-bounces at r-project.org [r-help-bounces at r-project.org] On
Behalf Of peter dalgaard [pdalgd at gmail.com]
Sent: Saturday, May 07, 2011 4:46 AM
To: Alex Olssen
Cc: r-help at r-project.org
Subject: Re: [R] maximum likelihood convergence reproducing Anderson Blundell
1982 Econometrica R vs Stata
On May 6, 2011, at 14:29 , Alex Olssen wrote:
> Dear R-help,
>
> I am trying to reproduce some results presented in a paper by Anderson
> and Blundell in 1982 in Econometrica using R.
> The estimation I want to reproduce concerns maximum likelihood
> estimation of a singular equation system.
> I can estimate the static model successfully in Stata but for the
> dynamic models I have difficulty getting convergence.
> My R program which uses the same likelihood function as in Stata has
> convergence properties even for the static case.
>
> I have copied my R program and the data below. I realise the code
> could be made more elegant - but it is short enough.
>
> Any ideas would be highly appreciated.
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.
Apart from that, it seems that the conjugate gradient methods have difficulties
with this likelihood, for some less than obvious reason. Increasing the maxit
gets you closer but still not satisfactory.
I would suggest trying out the experimental optimx package. Apparently, some of
the algorithms in there are much better at handling this likelihood, notably
"nlm" and "nlminb".
>
> ## model 18
> lnl <- function(theta,y1, y2, x1, x2, x3) {
> n <- length(y1)
> beta <- theta[1:8]
> e1 <- y1 - theta[1] - theta[2]*x1 - theta[3]*x2 - theta[4]*x3
> e2 <- y2 - theta[5] - theta[6]*x1 - theta[7]*x2 - theta[8]*x3
> e <- cbind(e1, e2)
> sigma <- t(e)%*%e
> logl <- -1*n/2*(2*(1+log(2*pi)) + log(det(sigma)))
> return(-logl)
> }
> p <- optim(0*c(1:8), lnl, method="BFGS", hessian=TRUE, y1=y1,
y2=y2,
> x1=x1, x2=x2, x3=x3)
>
>
"year","y1","y2","x1","x2","x3"
> 1929,0.554779,0.266051,9.87415,8.60371,3.75673
> 1930,0.516336,0.297473,9.68621,8.50492,3.80692
> 1931,0.508201,0.324199,9.4701,8.27596,3.80437
> 1932,0.500482,0.33958,9.24692,7.99221,3.76251
> 1933,0.501695,0.276974,9.35356,7.98968,3.69071
> 1934,0.591426,0.287008,9.42084,8.0362,3.63564
> 1935,0.565047,0.244096,9.53972,8.15803,3.59285
> 1936,0.605954,0.239187,9.6914,8.32009,3.56678
> 1937,0.620161,0.218232,9.76817,8.42001,3.57381
> 1938,0.592091,0.243161,9.51295,8.19771,3.6024
> 1939,0.613115,0.217042,9.68047,8.30987,3.58147
> 1940,0.632455,0.215269,9.78417,8.49624,3.57744
> 1941,0.663139,0.184409,10.0606,8.69868,3.6095
> 1942,0.698179,0.164348,10.2892,8.84523,3.66664
> 1943,0.70459,0.146865,10.4731,8.93024,3.65388
> 1944,0.694067,0.161722,10.4465,8.96044,3.62434
> 1945,0.674668,0.197231,10.279,8.82522,3.61489
> 1946,0.635916,0.204232,10.1536,8.77547,3.67562
> 1947,0.642855,0.187224,10.2053,8.77481,3.82632
> 1948,0.641063,0.186566,10.2227,8.83821,3.96038
> 1949,0.646317,0.203646,10.1127,8.82364,4.0447
> 1950,0.645476,0.187497,10.2067,8.84161,4.08128
> 1951,0.63803,0.197361,10.2773,8.9401,4.10951
> 1952,0.634626,0.209992,10.283,9.01603,4.1693
> 1953,0.631144,0.219287,10.3217,9.06317,4.21727
> 1954,0.593088,0.235335,10.2101,9.05664,4.2567
> 1955,0.60736,0.227035,10.272,9.07566,4.29193
> 1956,0.607204,0.246631,10.2743,9.12407,4.32252
> 1957,0.586994,0.256784,10.2396,9.1588,4.37792
> 1958,0.548281,0.271022,10.1248,9.14025,4.42641
> 1959,0.553401,0.261815,10.2012,9.1598,4.4346
> 1960,0.552105,0.275137,10.1846,9.19297,4.43173
> 1961,0.544133,0.280783,10.1479,9.19533,4.44407
> 1962,0.55382,0.281286,10.197,9.21544,4.45074
> 1963,0.549951,0.28303,10.2036,9.22841,4.46403
> 1964,0.547204,0.291287,10.2271,9.23954,4.48447
> 1965,0.55511,0.281313,10.2882,9.26531,4.52057
> 1966,0.558182,0.280151,10.353,9.31675,4.58156
> 1967,0.545735,0.294385,10.3351,9.35382,4.65983
> 1968,0.538964,0.294593,10.3525,9.38361,4.71804
> 1969,0.542764,0.299927,10.3676,9.40725,4.76329
> 1970,0.534595,0.315319,10.2968,9.39139,4.81136
> 1971,0.545591,0.315828,10.2592,9.34121,4.84082
>
> ______________________________________________
> R-help at 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.mes at cbs.dk Priv: PDalgd at gmail.com
______________________________________________
R-help at 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
2011-May-09 07:04 UTC
[R] maximum likelihood convergence reproducing Anderson Blundell 1982 Econometrica R vs Stata
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.mes at cbs.dk Priv: PDalgd at gmail.com
I wonder if someone with more experience than me on using R to summarise
by group wants to post a reply to this
http://www.analyticbridge.com/group/sasandstatisticalprogramming/forum/t
opics/why-still-use-sas-with-a-lot
To save everyone having to follow the link, the text is copied below
"SAS has some nice features, such as the SQL procedure or simple
"group
by" features. Try to compute correlations "by group" in R: say
you have
2,000 groups, 2 variables e.g. salary and education level, and 2 million
observations - you want to compute correlation between salary and
education within each group.
It is not obvious, your best bet is to use some R package (see sample
code on Analyticbridge to do it), and the solution is painful, you can
not return both correlation and stdev "by group", as the function can
return only one argument, not a vector. So if you want to return not
just two, but say 100 metrics, it becomes a nightmare."
________________________________________________________________________
The Numerical Algorithms Group Ltd is a company registered in England
and Wales with company number 1249803. The registered office is:
Wilkinson House, Jordan Hill Road, Oxford OX2 8DR, United Kingdom.
This e-mail has been scanned for all viruses by Star. Th...{{dropped:4}}
Alex Olssen
2011-May-09 10:06 UTC
[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.
They come from Econometrica Vol. 50, No. 6 (Nov., 1982), pp. 1569
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 <pdalgd at gmail.com>
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.mes at cbs.dk ?Priv: PDalgd at gmail.com
>
>