[R] calculate true autocovariance from power spectrum
Bao, Wenkai
wbao at mail.smu.edu
Tue Apr 12 22:31:45 CEST 2011
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
More information about the R-help
mailing list