Re: [R] R-help: gls with correlation=corARMA

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Mon 12 Dec 2005 - 08:44:52 EST

          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@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@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@pdf.com
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915

______________________________________________
R-help@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
Received on Mon Dec 12 08:54:40 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:41:35 EST