编辑代码

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))
  ## 用Rarima函数估计模型参数。最大似然估计法。
  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()