Hi there,

I have been using the nlme::gls package created in R to fit a pretty

simple model (linear with AR error)

y(t) = beta*x(t) + e(t) where e(t) ~ rho*e(t-1) + Z(t)

and Z(t)~ N(0,sig^2)

I call the R routine

glsObj <- nlme::gls(y ~ x -1, data=data, correlation =

nlme::corAR1(form= ~x), method="ML")

All seems fine.

In addition, I have also coded the likelihood myself and maximized it

for beta, rho and sigma.

I get the exact same estimates of beta and rho, (as nlme::gls) but the

estimate of sigma is not the same and i can not figure out why.

The maximum likelihood estimator for sigma under this model is

sig^2 = (( 1-rho^2)u(1)^2 + sum((u(t)- rho*u(t-1))^2)/n

where the sum is t=2,...,n and

u(t) = y(t) - X(t)*beta

I have read the mixed-effects models in S and S-Plus book (nlme::gls

code is based directly on this) and this problem is specified on page

204 eq (5.5). I have also calculated sigma based on (5.7) -after the

transformation documented (5.2) -and i do not get the same value as

either the package or my implementation.

Any advice would be most welcomed. Is there a bug in the estimation of

sigma in this package?

Thanks

Andy

--

Andy Beet

Ecosystem Dynamics & Assessment Branch

Northeast Fisheries Science Center

NOAA Fisheries Service

166 Water Street

Woods Hole, MA 02543

tel: 508-495-2073

[[alternative HTML version deleted]]

______________________________________________

[hidden email] mailing list -- To UNSUBSCRIBE and more, see

https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide

http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.