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