## Examples of Applied Time Series Analysis, By He Shuyuan
## Chpater 5, Predictions


demo.pred.arma42 <- function(n=40){
  a <- c(-0.9, -1.4, -0.7, -0.6)
  b <- c(0.5, -0.4)
  ng <- 30

  x <- arma.gen(n, a, b, plot.it=F)
  gams <- arma.gamma(ng, a, b, sigma=1)

  n.use <- 14  ## use n.use points in prediction
  m.pred <- 7  ## predit m.pred steps
  n.start <- n - m.pred - n.use + 1  ## which x are affected
  ## predict usging true ACV
  Ga <- matrix(0, nrow=n.use, ncol=n.use)
  for(ii in seq(n.use)) for(jj in seq(n.use)) {
    ind <- abs(ii-jj)+1
    Ga[ii,jj] <- gams[ind]
  }
  Gar <- solve(Ga)

  x.use <- x[n.start:(n.start+n.use-1)]
  y.pred <- x
  y.pred[] <- NA
  errs <- numeric(m.pred)
  for(k in seq(m.pred)){
    g <- gams[(n.use+k):(k+1)]
    a <- Gar %*% g
    y.pred[n.start+n.use-1+k] <- c(crossprod(a, x.use))
    errs[k] <- gams[1] - g %*% Gar %*% g
  }
  lb <- y.pred[(n-m.pred+1):n]-1.96*errs
  ub <- y.pred[(n-m.pred+1):n]+1.96*errs
  yl <- range(c(lb,ub,x,y.pred), na.rm=T)
  plot(1:n, x, type="n",
       main="Prediction of ARMA(4,2)",
       xlab="t", ylab="y",
       ylim=yl)
  lines(1:(n.start), x[1:n.start],
        lwd=2)  ## unused
  lines(n.start:(n.start+n.use-1), x[n.start:(n.start+n.use-1)],
        col="green", lwd=2) ## used for predition
  lines((n-m.pred):n, x[(n-m.pred):n],
        type="b", col="green", lty=3, lwd=2)  ## true values
  lines((n-m.pred+1):n, y.pred[(n-m.pred+1):n],
        type="b", col="red", lwd=2)  ## predictons
  lines((n-m.pred+1):n, y.pred[(n-m.pred+1):n]+1.96*errs,
        type="l", col="red", lty=2)  ## upper bounds
  lines((n-m.pred+1):n, y.pred[(n-m.pred+1):n]-1.96*errs,
        type="l", col="red", lty=2)  ## lower bounds

  invisible()
}