[R] R-help: gls with correlation=corARMA
Spencer Graves
spencer.graves at pdf.com
Sun Dec 11 22:44:52 CET 2005
The error message is misleading. It should say something like,
"Error in corARMA(q = 25, value = -ma.coefs, fixed = T) : The moving
average process specified is not invertible, having roots outside the
unit circle." Instead it says, "Error in corARMA(q = 25, value =
-ma.coefs, fixed = T) : All parameters must be less than 1 in absolute
value." I'm copying Doug Bates on this reply in case he wants to try to
fix this.
I got an answer just by shrinking your ma.coefs' by a factor of 0.8:
mod.gls=gls(obs~model,correlation=corARMA(q=25,value=0.8*ma.coefs,fixed=T),
method="ML")
This seemed to produce an answer for me; it least it did not give me
an error message.
In case you are interested in how I determined this, I will outline
the steps I took in analyzing this problem. First, I copied the web
address you gave for the data into a web browser to make sure it was
honest text and not something that might corrupt my computer. You are
to be commended for providing an example that allowed me to replicate
your problem. If the example had been smaller and simpler, it would
have made my job easier and might have gotten you an earlier reply from
someone else. Then I ran your code and got the error you reported:
...
> mod.gls=gls(obs~model,
+ correlation=corARMA(q=25,value=ma.coefs,fixed=T),
+ method="ML")
Error in corARMA(q = 25, value = ma.coefs, fixed = T) :
All parameters must be less than 1 in absolute value
Next, I considered ways to simplify this problem and still get the
same error message. I decided to try the "corARMA" part by itself:
> corARMA(q=25,value=ma.coefs,fixed=T)
Error in corARMA(q = 25, value = ma.coefs, fixed = T) :
All parameters must be less than 1 in absolute value
>
Progress. Then I typed "corARMA" at a command prompt and copied the
code into a scrit file. The I typed "debug(corARMA)" and repeated the
"corARMA(...)" command. After tracing through the corARMA code line by
line, I found that the error message is issued from
'.C("ARMA_unconstCoef", ...)'. I gave that up: This approach did not
help in this case, thoug it has in others.
Then I tried some simpler examples: 'corARMA(q=1,value=.5,fixed=T)'
and 'corARMA(q=1,value=-.5,fixed=T)' did NOT give me that error message,
but 'corARMA(q=2,value=c(.8, -.5),fixed=T)' did.
Then I checked a time series book for the conditions for
invertibility. I found that all the roots of the characteristic
equation must lie outside the unit circle. So I checked the following:
> round(Mod(polyroot(c(1, ma.coefs))), 3)
[1] 1.069 0.995 0.995 0.995 0.995 0.995 0.995 0.995 0.995 1.069 1.069
1.069
[13] 0.995 0.995 0.995 0.995 1.069 1.069 1.069 1.069 1.069 1.069 1.069 1.069
[25] 1.930
Then I shrunk the ma.coefs' by 0.999 and got larger roots but still
some inside the unit circle. So I tried 0.99 and 0.9 with the same
result. With 0.8, all the roots were outside the unit circle.
hope this helps.
spencer graves
gaffigan at sfos.uaf.edu wrote:
> Dear Madams/Sirs,
>
> Hello. I am using the gls function to specify an arma correlation during
> estimation in my model. The parameter values which I am sending the
> corARMA function are from a previous fit using arima. I have had some
> success with the method, however in other cases I get the following error
> from gls: "All parameters must be less than 1 in absolute value". None
> of the parameters (individually) are greater than or equal to 1.
> Please copy the code below into R to reproduce the error. Thanks.
>
> Is my logic incorrect? In the corARMA function, there's a call to
> pre-compiled C code with the name "ARMA_unconstCoef". Is the source
> code for such compiled code freely available for download?
> Thanks for your suggestions.
>
> Sincerely
>
> Steve Gaffigan
>
> data=read.table("http://ak.aoos.org/data/sample_070989.dat",header=T)
> attach(data)
> mod.ols=lm(obs~model)
> mod.sma=arima(residuals(mod.ols),order=c(0,0,1),seasonal=list(order=c(0,0,2),period=12))
> theta.1=mod.sma$coef[1]
> THETA.1=mod.sma$coef[2]
> THETA.2=mod.sma$coef[3]
> ma.coefs=c(-theta.1,double(10),-THETA.1,theta.1*THETA.1,double(10),-THETA.2,theta.1*THETA.2)
> library(nlme)
> mod.gls=gls(obs~model,correlation=corARMA(q=25,value=ma.coefs,fixed=T),method="ML")
> detach(data)
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
--
Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA
spencer.graves at pdf.com
www.pdf.com <http://www.pdf.com>
Tel: 408-938-4420
Fax: 408-280-7915
More information about the R-help
mailing list