arma.spec <- function(a, b, sigma, ngrid=201){
p <- length(a)
q <- length(b)
freqs <- seq(from=0, to=pi, length=ngrid)
spec1 <- numeric(ngrid)
spec2 <- numeric(ngrid)
for(ii in seq(ngrid)){
spec1[ii] <- 1 + sum(complex(mod=b, arg=freqs[ii]*seq(q)))
spec2[ii] <- 1 - sum(complex(mod=a, arg=freqs[ii]*seq(p)))
}
spec = sigma^2 / (2*pi) * abs(spec1)^2 / abs(spec2)^2
list(frequency=freqs, spec=spec)
}
demo.arma.spec.sim <- function(n=300, ngrids=201){
a <- -c(0.9, 1.4, 0.7, 0.6)
b <- c(0.5, -0.4)
sig <- 1.0
x <- arima.sim(
model=list(ar=a, ma=b), n=n,
rand.gen=function(n, ...) rnorm(n, 0, sig))
## 用R的arima函数估计模型参数。最大似然估计法。
fit.mle <- arima(x, order=c(4,0,2),
include.mean=FALSE)
sres1 <- arma.spec(a=a, b=b, sigma=sig)
freqs <- sres1$frequency
spec <- sres1$spec
plot(freqs, spec, type='l', lwd=2,
main="ARMA(4,2) 谱密度估计",
xlab="frequency", ylab="spectrum",
axes=FALSE)
axis(2)
axis(1, at=(0:6)/6*pi,
labels=c(0, expression(pi/6),
expression(pi/3), expression(pi/2),
expression(2*pi/3), expression(5*pi/6), expression(pi)))
box()
## 添加估计值
sres2 <- arma.spec(
a=coef(fit.mle)[1:4],
b = coef(fit.mle)[5:6],
sigma=fit.mle$sigma2)
lines(freqs, sres2$spec, lty=2)
## 图例
legend("topleft", lty=c(1,2), lwd=c(2,1), legend=c("真实", "估计"))
invisible()
}
set.seed(1)
demo.arma.spec.sim()