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

[gaffigan@sfos.uaf.edu]

Thank you for your extra efforts in pinpointing the source of my problem.
That is a smart workaround, to reduce the parameters until invertibility 
conditions are satisfied.  The concern in my case is what effect the 
constant will have on the covariance matrix and the resulting estimates 
for the slope in the model.  I might look at other cases where the 
invertibility condition was satisfied and compare the slope estimates 
and standard errors for the unscaled and scaled time series parameters.  

The example below was one of many series.  I should have been quicker to 
realize that someone might use the same form as my example to mess 
with someones machine.  I have used tryCatch over the weekend to skip the 
gls fit in cases where this error occurred (4282 out of 187283 series).  
There were no errors when I fit a (0,0,1)x(0,0,1) model (e.g. q=13), 
instead of q=25.  Thank you again for your time and explanation.

Sincerely,

Steve Gaffigan


On Sun, 11 Dec 2005, Spencer Graves wrote:

> 	  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
> 
> 

-- 
Alaska Ocean Observing System
School of Fisheries and Ocean Sciences : University of Alaska Fairbanks