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

[Spencer Graves]

	  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