# Distribution of estimated acf

# R has two options to compute confidence limits for the acf
acf(log(lynx))
acf(log(lynx),ci.type="ma")
# For explanation use
help(acf)
help(plot.acf)

#simulations which show the distribution of estimated autocorrelations
par(mfrow=c(2,2))
resv <- matrix(0,nrow=200,ncol=10)
for (i in (1:200)) resv[i,] <- f.sim.acf(x.sim=f.ar.sim,beta=0.8)
for (i in (1:200)) resv[i,] <- f.sim.acf(x.sim=f.ma.sim,a=1)
for (i in (1:200)) resv[i,] <- f.sim.acf(x.sim=f.arch.sim,a=0.8)
for (i in (1:200)) resv[i,] <- f.sim.acf(x.sim=rnorm)
boxplot(resv,ylab="boxplots of acf",xlab="lag")
round(sqrt(apply(resv,2,var)),3) # compare with the value 1/sqrt(200)=0.71 

# Functions
f.sim.acf <- function(n=200,lag=10,x.sim,...)
{
  ## Purpose: Simulation of the autocorrelation function for various models
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 20 Nov 2000, 14:30
  x <- x.sim(n,...)
  acf(x,lag.max=lag,plot=F)$acf[-1]}

f.ar.sim <- function(n, beta, sigma = 1.0)
{
    x <- ts(rnorm(n+100, 0, sigma), start = -99)
    x <- filter(x, beta, method = "recursive")
    as.ts(x[-(1:100)])
}

f.ma.sim <- function(n, a)
{
    x <- rnorm(n+1)
    x <- x[-1] + a*x[-(n+1)]
    as.ts(x)
}

f.arch.sim <- function(n=500,init=200,a)
{
  ## Purpose: Simulating the Arch(1) process
  ## -------------------------------------------------------------------------
  ## Arguments: n=length of the series, init=length of the initial segment
  ## that is discarded, a= Parameter of the model
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date:  5 Feb 97, 16:49
  m <- n+init
  x <- rnorm(m)
  for (i in (2:m)) x[i] <- sqrt(1+a*x[i-1]^2)*x[i]
  ts(x[(init+1):m],start=1)
}