## Examples for Applied Time Series Analysis, He Shuyuan.
## Chapter 2, AR model

## simulate AR(p) model.
## a is vector a_1, \dots, a_p
## Model is X_t = a_1 X_{t-1} + \dots + a_p X_{t-p} + \epsilon_t
## \Var(\epsilon_t) = sigma^2
## This version use manual iteration, is obsolete
ar.gen.old <- function(n, a, sigma=1.0, plot.it=FALSE, n0=1000){
  ##set.seed(1)
  n2 <- n0 + n
  p <- length(a)
  arev <- rev(a)
  x2 <- numeric(n2)
  eps <- rnorm(n2, 0, sigma)
  for(tt in seq(p+1, n2)){
    x2[tt] <- eps[tt] + sum(x2[(tt-p):(tt-1)] * arev)
  }
  x <- x2[(n0+1):n2]
  x <- ts(x)
  attr(x, "model") <- "AR"
  attr(x, "a") <- a
  attr(x, "sigma") <- sigma
  if(plot.it) plot(x)
  x
}

## simulate AR(p) model.
## a is vector a_1, \dots, a_p
## Model is X_t = a_1 X_{t-1} + \dots + a_p X_{t-p} + \epsilon_t
## \Var(\epsilon_t) = sigma^2
## This version use the ``filter'' function.
ar.gen <- function(n, a, sigma=1.0, by.roots=FALSE,
                   plot.it=FALSE, n0=1000,
                   x0=numeric(length(a))){
  if(by.roots){
    require(polynom)
    cf <- Re(coef(poly.calc(a)))
    cf <- cf / cf[1]
    a <- -cf[-1]
  }
  ##set.seed(1)
  n2 <- n0 + n
  p <- length(a)
  eps <- rnorm(n2, 0, sigma)
  x2 <- filter(eps, a, method="recursive", side=1, init=x0)
  x <- x2[(n0+1):n2]
  x <- ts(x)
  attr(x, "model") <- "AR"
  attr(x, "a") <- a
  attr(x, "sigma") <- sigma
  if(plot.it) plot(x)
  x
}

## AR theoretical spectrum given AR coefficients
ar.true.spectrum <- function(a, ngrid=256, sigma=1, plot.it=TRUE){
  p <- length(a)
  freqs <- seq(from=0, to=pi, length=ngrid)
  spec <- numeric(ngrid)
  for(ii in seq(ngrid)){
    spec[ii] <- 1 - sum(complex(mod=a, arg=freqs[ii]*seq(p)))
  }
  spec = sigma^2 / (2*pi) / abs(spec)^2
  if(plot.it){
    plot(freqs, spec, type='l',
         main="AR True Spectral Density",
         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)))
  }
  list(frequencies=freqs, spectrum=spec,
       ar.coefficients=a, sigma=sigma)
}

demo.ar1 <- function(){
  as <- c(0.1, 0.3, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99, 1.0, 1.001)
  x0 <- 10
  ##ns <- c(rep(100, 3*2), rep(500, 3*2), rep(1000, 2),
  ##        rep(2000,2), rep(5000,6))
  ns <- rep(500, 2*length(as))
  ii <- 0
  for(a in as){
    for(a1 in c(a, -a)){
      ii <- ii+1
      n <- ns[ii]
      x <- ar.gen(n, a1, sigma=0.2, n0=0, x0=x0)
      plot(x, main=paste("AR(1): x0=10, a=", a1, sep=""),
           ylim=c(-12, 12))
      abline(h=0, col="red")
      locator(1)
    }
  }
}

demo.ar.roots <- function(){
  n <- 1024
  rtlis <- list(c(complex(mod=1.09, arg=pi/3*c(1,-1)),
                  complex(mod=1.098, arg=pi*2/3*c(1,-1))),
                c(complex(mod=1.264, arg=pi/3*c(1,-1)),
                  complex(mod=1.273, arg=pi*2/3*c(1,-1))),
                c(complex(mod=1.635, arg=pi/3*c(1,-1)),
                  complex(mod=1.647, arg=pi*2/3*c(1,-1))),
                complex(mod=1.02, arg=pi/6*c(1,-1)),
                complex(mod=1.02, arg=pi/2*c(1,-1)),
                c(complex(mod=1.05, arg=pi/6*c(1,-1)),
                  complex(mod=1.05, arg=pi/2*c(1,-1))))
  tits <- c("AR(4): 1.09exp(+-i pi/3), 1.098exp(+-i 2pi/3)",
            "AR(4): 1.264exp(+-i pi/3), 1.273exp(+-i 2pi/3)",
            "AR(4): 1.635exp(+-i pi/3), 1.647exp(+-i 2pi/3)",
            "AR(2): mod=1.02 arg=+-pi/6",
            "AR(2): mod=1.02 arg=+-pi/2",
            "AR(4): mod=1.05 arg=+-pi/6,+-pi/2")
  oldpar <- par(mfrow=c(3,1)); on.exit(par(oldpar))
  for(ii in seq(along=rtlis)){
    rt = rtlis[[ii]]
    y <- ar.gen(n, rt, sigma=1.0, by.roots=TRUE,
                plot.it=FALSE)
    plot(window(y, 1, 60), main=tits[ii])
    acf(y)
    ##spectrum(y, taper=0.2)
    ar.true.spectrum(attr(y, "a"))
    locator(1)
  }
}

## example 2.5.1
demo.ar1.example <- function(){
  oldpar <- par(mfrow=c(2,2)); on.exit(par(oldpar))
  n <- 100
  a1 <- 0.85

  y1 <- ar.gen(n, a1, sigma=1.0, plot.it=FALSE)
  y2 <- ar.gen(n, -a1, sigma=1.0, plot.it=FALSE)
  plot(y1, main="AR(1): a = 0.85")
  acf1 <- exp(seq(from=0, by=log(a1), length=21))
  plot(0:20, acf1, type="l",
       xlab="k", ylab="ACF",
       main="ACF of AR(1): a = 0.85")
  ar.true.spectrum(attr(y1, "a"))
  plot(y1, type="n", axes=F, xlab="", ylab="")
  locator(1)

  plot(y2, main="AR(1): a = -0.85")
  acf2 <- (-a1)^seq(from=0, length=21)
  plot(0:20, acf2, type="l",
       xlab="k", ylab="ACF",
       main="ACF of AR(1): a = -0.85")
  abline(h=0)
  ar.true.spectrum(attr(y2, "a"))
  plot(y1, type="n", axes=F, xlab="", ylab="")
  invisible()
}

ar2.stable.graph <- function(){
  oldpar <- par(mfrow=c(1,1), mar=c(1,1,0.5,0.5));
  on.exit(par(oldpar))
  plot(c(-2.2, 2.2), c(-1.2, 1.2),
       type="n", axes=F,
       xlab=expression(a[1]),
       ylab=expression(a[2]))
  nn <- 100
  x1 <- seq(-2, 2, length=nn)
  y1 <- -1/4 * x1^2
  polygon(x1, y1, col="red") # imaginary roots
  polygon(c(0,x1), c(1,y1), col="blue") # real roots
  lines(c(-2,0), c(-1, 1))
  lines(c(2,0), c(-1,1))
  lines(c(-2.2, 2.2), c(0,0), lwd=1.5)
  lines(c(0,0), c(-1.5, 1.5), lwd=1.5)
  text(0.15, 1, "1")
  text(0.15, -1.1, "-1")
  lines(c(2,2), c(-0.05,0.05))
  text(2, 0.15, "2")
  lines(c(-2,-2), c(-0.05,0.05))
  text(-2, 0.15, "-2")

  ## dev.off()
  invisible()
}

f <- function(){
  postscript(file="ar2cond.eps", paper="special",
             width=8/2.54, height=4/2.54, horizontal=F, onefile=F)
  ar2.stable.graph()
  dev.off()

  pdf(file="ar2cond.pdf", paper="special",
      width=8/2.54, height=4/2.54, horizontal=F, onefile=F)
  ar2.stable.graph()
  dev.off()
}


demo.ar2.example <- function(){
  oldpar <- par(mfrow=c(3,1)); on.exit(par(oldpar))
  n <- 100
  acomplex <- cbind(2/1.1*cos(c(pi/6, pi/2, 2*pi/3)),
                    -1/1.1^2)
  alist <- rbind(c(0.1, 0.5),
                 c(-0.1, 0.5),
                 c(0, 0.8),
                 c(1, -0.1),
                 c(-1, -0.1),
                 acomplex)
  tlist.y <- paste("AR(2): a1=", round(alist[,1],3),
                   " a2=", round(alist[,2],3), sep="")
  tlist.rho <- c("反号实根, 正根近单位圆, 取震荡的正值衰减",
                 "反号实根, 负根近单位圆, 正负交替衰减",
                 "根为相反数, 非负震荡衰减",
                 "两正根, 单调正值衰减",
                 "两负根, 正负交替衰减",
                 "共轭复根, 以周期12震荡衰减",
                 "共轭复根, 以周期4震荡衰减",
                 "共轭复根, 以周期3震荡衰减")

  m <- 20
  acfv <- numeric(m+1)
  acfv[1] <- 1
  for(ii in seq(nrow(alist))){
    y <- ar.gen(n, alist[ii,], sigma=1.0, plot.it=FALSE)
    a1 <- alist[ii,1];
    a2 <- alist[ii,2];
    acfv[2] <- a1 / (1 - a2)
    for(jj in seq(3, m+1)){
      acfv[jj] <- a1*acfv[jj-1] + a2*acfv[jj-2]
    }
    plot(y, main=tlist.y[ii], type="b")
    plot(0:m, acfv, type="l",
         xlab="k", ylab="ACF",
         main=tlist.rho[ii])
    abline(h=0)
    ar.true.spectrum(attr(y, "a"))
    ##plot(y, type="n", axes=F, xlab="", ylab="")
    locator(1)
  }
  invisible()
}

##demo.ar2.example()

## 随机点作出AR(2)稳定域和容许域。
ar2.allow1 <- function(nrep=1000){
  oldpar <- par(mfrow=c(2,1)); on.exit(par(oldpar))
  ip <- 0
  amat <- matrix(0, nrep, 2)
  rhomat <- amat
  for(ii in seq(nrep)){
    a1 <- runif(1,-2,2)
    a2 <- runif(1,-1,1)
    if(abs(a2)<1){
      if(a2+a1<1 && a2-a1<1){
        ip <- ip+1
        amat[ip,1] <- a1
        amat[ip,2] <- a2

        rho1 <- a1/(1-a2)
        rho2 <- a1*rho1 + a2
        rhomat[ip,1] <- rho1
        rhomat[ip,2] <- rho2
      }
    }
  }

  a1 <- amat[1:ip,1]
  a2 <- amat[1:ip,2]
  rho1 <- rhomat[1:ip,1]
  rho2 <- rhomat[1:ip,2]
  plot(a1, a2, xlim=c(-2,2), ylim=c(-1,1),
       main="AR(2): a2 by a1")
  plot(rho1, rho2, xlim=c(-1,1), ylim=c(-1,1),
       main="AR(2): rho2 by rho1")
}

## 随机点作出AR(2)容许域和稳定域。
ar2.allow <- function(nrep=10000){
  oldpar <- par(mfrow=c(2,1)); on.exit(par(oldpar))
  ip <- 0
  amat <- matrix(0, nrep, 2)
  rhomat <- amat
  for(ii in seq(nrep)){
    rho1 <- runif(1,-1,1)
    rho2 <- runif(1,-1,1)
    if(abs(rho1)<1 && abs(rho2)<1){
      if(rho1^2 < (1 + rho2)/2){
        ip <- ip+1
        rhomat[ip,1] <- rho1
        rhomat[ip,2] <- rho2

        a1 <- (rho1 - rho1*rho2)/(1-rho1^2)
        a2 <- (rho2 - rho1^2)/(1-rho1^2)
        amat[ip,1] <- a1
        amat[ip,2] <- a2
      }
    }
  }

  a1 <- amat[1:ip,1]
  a2 <- amat[1:ip,2]
  rho1 <- rhomat[1:ip,1]
  rho2 <- rhomat[1:ip,2]
  plot(rho1, rho2, 
       main="AR(2): rho2 by rho1")
  plot(a1, a2, 
       main="AR(2): a2 by a1")
}

demo02 <- function(){
  demos <- list("AR(1) demo" = demo.ar1,
                "AR roots demo" = demo.ar.roots,
                "AR(1) example" = demo.ar1.example,
                "AR(2) example" = demo.ar2.example
                )
  ndemos <- length(demos)
  prompts <- paste("   ", seq(ndemos), ". ", names(demos), "\n", sep="")
  cat("Demos of Chapter 2:\n",
      prompts)
  sel <- scan(n=1)
  demos[[sel]]()
}



##tmp.x <- ar.gen(50, 0.5)
##demo.ar1()