[R] calculate true autocovariance from power spectrum

From: Bao, Wenkai <wbao_at_mail.smu.edu>
Date: Tue, 12 Apr 2011 20:31:45 +0000


I know using ARMAacf function can do the job for ARMA model, but it is not calculating from power spectrum. I have been trying to code with the following algorithm:

Since

|1-theta1*exp(2*pi*f*i)-...-thetaq*[exp(2*pi*f*i)]^q|^2

P(f)=sigma2*----------------------------------------------------------------------  ,

|1-phi1*exp(2*pi*f*i)-...-phip*[exp(2*pi*f*i)]^q|^2

autocovariances are the inverse Fourier transform of P(f), thus

           .5                                        .5
acv_k=int P(f)*exp(2*pi*f*i*k)df=2*int P(f)*cos(2*pi*f*k)df  , k is the lag 0,1,...
          -.5                                         0

So I make a R function trueacf accordingly:

####################################################
trueacf=function(phi,theta,sigma2=1,lags=20) {
 if(sum(abs(phi))!=0) {p=length(phi)} else p=0 ## decide AR order  if(sum(abs(theta))!=0) {q=length(theta)} else q=0 ## decide MA order

 acv=c()
 for (k in 1:(lags+1)) ## acv[1] is the variance of the process  {
  newfunc=function(f)
  {
   if(q==0) {comp1=1} else comp1=Mod(c(1,-theta)%*%exp(2*pi*(1i)*f)^c(0:q))^2    if(p==0) {comp2=1} else comp2=Mod(c(1,-phi)%*%exp(2*pi*(1i)*f)^c(0:p))^2    result.temp=sigma2*comp1/comp2*cos(2*pi*f*(k-1))    return(result.temp)
  }
  acv[k]=2*integrate(newfunc,0,.5)$value  }
 acf=acv/acv[1]
 list(acv=acv,acf=acf)
}

######################################################
Strangely, it keeps giving me this error msg: Error in c(1, -phi) %*% exp(2 * pi * (0+1i) * f)^c(0:p) :   non-conformable arguments
In addition: Warning message:
In exp(2 * pi * (0+1i) * f)^c(0:p) :
  longer object length is not a multiple of shorter object length

I checked the code many times, and I could not find any mistake. As a matter of fact, if I copy and paste the definitions of the parameters and newfunc, I can get the newfunc function run correctly. So I do not understand why it is error when the definitions are in the function trueacf.

Thank you very much for help in advance!

Wenkai Bao
Ph.D Student,
Department of Statistical Science,
Southern Methodist University
Dallas, TX 75275



R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Tue 12 Apr 2011 - 21:02:18 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Tue 12 Apr 2011 - 22:00:30 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive