# 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)
}

# Constructing a moving average with a periodic behavior

# We construct a polynomial Q such that |Q(exp(2*pi*lambda))|^2
# has a maximum around 1/11. The moving average process with
# the spectral density equal to |Q(exp(-2*pi*lambda))|^2
# has then an approximately periodic behavior with period about 11

# Q is constructed by multiplying polynomials which have zeroes at
# +1, -1, exp(+-2.5i), exp(+-1.8i) and exp(+-1.1i)
nu <- seq(0,0.5,length=(1001))*2*pi
z <- exp(-nu*1i)
Q <- (1-z)^2*(1+z)^2*(1-2*cos(2.5)*z+z^2)^2*
  (1-cos(1.8)*z+z^2)^2*(1-cos(1.1)*z+z^2)^2
plot(nu/(2*pi), Mod(Q)^2,type="l")
2*pi/nu[Mod(Q) == max(Mod(Q))] # 1/frequency where spectral density is maximal
# Simulating the corresponding moving average process by applying successive
# moving averages for each factor in Q.
x <- rnorm(1000)
x <- filter(x,c(1,-1),"convolution",sides=1)
x <- filter(x,c(1,-1),"convolution",sides=1)
x <- filter(x,c(1,1),"convolution",sides=1)
x <- filter(x,c(1,1),"convolution",sides=1)
x <- filter(x,c(1,-2*cos(2.5),1),"convolution",sides=1)
x <- filter(x,c(1,-2*cos(2.5),1),"convolution",sides=1)
x <- filter(x,c(1,-2*cos(1.8),1),"convolution",sides=1)
x <- filter(x,c(1,-2*cos(1.8),1),"convolution",sides=1)
x <- filter(x,c(1,-2*cos(1.1),1),"convolution",sides=1)
x <- filter(x,c(1,-2*cos(1.1),1),"convolution",sides=1)
plot(x)
acf(x[17:1000])