# Autocorrelation and lag plots
# Functions: acf, lagplot

# Try to understand how properties of the time series plot are reflected in
# the estimated autocorrelation, and what additional information
# is in the lag plot. 

#1: lynx
par(mfrow=c(2,1))
plot(lynx)
plot(log(lynx))
acf(lynx)
acf(log(lynx))
pacf(lynx)
pacf(log(lynx))
par(mfrow=c(1,1))
lag.plot(log(lynx),lags=12,layout=c(3,4))

#2: Airlines
plot(log(AirPassengers))
acf(log(AirPassengers))
acf(diff(diff(log(AirPassengers)),lag=12))
pacf(diff(diff(log(AirPassengers)),lag=12))
airl.stl <- stl(log(AirPassengers),s.window=7,t.window=21)
acf(airl.stl$time.series[,3])

#3: Ocean waves : "The change in the level of ambient noise in the ocean"
# (4 measurements per second), from Percival& Walden, Spectral Analysis for
# Physical Applications
ocwave <- ts(scan("http://stat.ethz.ch/Teaching/Datasets/ocwave.dat"),
          start=1,deltat=0.25)
plot(ocwave)
acf(ocwave)
acf(ocwave,lag.max=200)
pacf(ocwave)

#4: Simulation of AR
phi <- 0.8
phi <- -0.8
x <- filter(ts(rnorm(200)),filter=0.8,method="recursive",init=rnorm(1))
plot(x)
acf(x)
pacf(x)

#4: Simulation of moving average
n <- 1000
phi <- 0.8
phi <- -0.8
x <- rnorm(n)
x <- ts(x[-1]+phi*x[-n])
plot(x)
acf(x)
pacf(x)

#5: Simulation of Arch
              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)
}
x <- f.arch.sim(n=2000,a=0.8)
acf(x)
acf(x^2)
par(mfrow=c(1,1))
lag.plot(x,lags=4,layout=c(2,2))
lag.plot(x^2,lags=4,layout=c(2,2))
lag.plot(log(abs(x)),lags=4,layout=c(2,2))

# 6: Financial time series (SMI)
SMI.ret <- diff(log(EuStockMarkets)[,2])
plot(SMI.ret)
acf(SMI.ret)
acf(SMI.ret^2)
acf(SMI.ret[1:1300])
acf(SMI.ret[1:1300]^2)

#6: Linear increase
acf(1:100,lag.max=50)
   # Why can an estimated autocorrelation become negative when x[t+h] > x[t]
   # for all t ?


#7: Simulation of a chaotic system
f.2 <- function(x) 3.8*x*(1-x) #logistic growth
n <- 500
sigma <- 0.01 # s.d. of the noise. See what happens if sigma=0
x <- ts(rnorm(500))
x[1] <- 0.5
for (i in (2:500)) {
    x[i] <- f.2(x[i-1]) + sigma*x[i]
    x[i] <- max(0,min(1,x[i])) #keep the process in [0,1]
  }
plot(x)
acf(x)
lag.plot(x,lags=6,layout=c(2,3))


#8: Multivariate time series
sales1 <- ts.union(BJsales, lead = BJsales.lead)
plot(sales1)
acf(diff(sales1))
sales2 <- ts.intersect(BJsales, lead3 = lag(BJsales.lead, -3))
acf(diff(sales2))
plot(EuStockMarkets)
plot(diff(EuStockMarkets))
plot(diff(log(EuStockMarkets)))
acf(diff(log(EuStockMarkets))) # no structure for lags different from zero
acf((diff(log(EuStockMarkets)))^2) # small structure for lags different from zero