Time Series II: AR, MA, ARMA, ARIMA

rm(list=ls())
# white noise
set.seed(1)
w <- rnorm(100)
plot(w, type = "l")
x <- seq(-3,3, length = 1000)
hist(rnorm(100), prob = T); points(x, dnorm(x), type = "l")
set.seed(2)
acf(rnorm(100))
# random walk
set.seed(2)
x <- w <- rnorm(1000)
for (t in 2:1000) x[t] <- x[t - 1] + w[t]
plot(x, type = "l")
acf(x)
acf(diff(x))
www <- "http://elena.aut.ac.nz/~pcowpert/ts/HP.txt"
HP.dat <- read.table(www, header = T) ; attach(HP.dat)
plot (as.ts(Price))
DP <- diff(Price) ; plot (as.ts(DP)) ; acf (DP)
mean(DP) + c(-2, 2) * sd(DP)/sqrt(length(DP))
# AR(1)
rho <- function(k, alpha) alpha^k
layout(1:2)
plot(0:10, rho(0:10, 0.7), type = "b")
plot(0:10, rho(0:10, -0.7), type = "b")
layout(1:1)
set.seed(1)
x <- w <- rnorm(100)
for (t in 2:100) x[t] <- 0.7 * x[t - 1] + w[t]
plot(x, type = "l")
acf(x)
pacf(x)
x.ar <- ar(x, method = "mle")
x.ar$order
x.ar$ar
x.ar$ar + c(-2, 2) * sqrt(x.ar$asy.var)
www = "http://elena.aut.ac.nz/~pcowpert/ts/pounds_nz.dat"
Z = read.table(www, header = T)
Z.ts = ts(Z, st = 1991, fr = 4)
Z.ar <- ar(Z.ts)
mean(Z.ts)
Z.ar$order
Z.ar$ar
Z.ar$ar + c(-2, 2) * sqrt(Z.ar$asy.var)
acf(Z.ar$res[-1])
www = "http://elena.aut.ac.nz/~pcowpert/ts/global.dat"
Global = scan(www)
Global.ts = ts(Global, st = c(1856, 1), end = c(2005, 12), fr = 12)
Global.ar <- ar(aggregate(Global.ts, FUN = mean), method = "mle")
mean(aggregate(Global.ts, FUN = mean))
Global.ar$order
Global.ar$ar
acf(Global.ar$res[-(1:Global.ar$order)], lag = 50)
# MA: correlogram and simulation
# correlogram
rho = function(k, beta){
  q = length(beta) - 1
  if (k > q) ACF = 0
  else{
    s1 = 0
    s2 = 0
    for (i in 1:(q - k + 1)){
      s1 = s1 + beta[i] * beta[i + k]
    }
    for (i in (1:(q + 1))){
      s2 = s2 + beta[i]^2
    }
    ACF = s1/s2
  }
  ACF 
}
 
beta <- c(1, 0.7, 0.5, 0.2) 
rho.k <- rep(1, 10)
for (k in 1:10) rho.k[k] <- rho(k, beta)
plot(0:10, c(1, rho.k), pch = 4, ylab = expression(rho[k])) 
abline(0, 0) 
 
# simulation of MA(3) 
set.seed(1) 
b = c(0.8, 0.6, 0.4) 
x = w = rnorm(1000)
for (t in c(4:1000)){
  for (j in c(1:3)){
    x[t] = x[t] + b[j] * w[t - j]
  }
}
plot(x, type = "l")
acf(x)
x.ma <- arima(x, order = c(0, 0, 3))
x.ma
x.ma <- arima(x, order = c(0, 0, 3), include.mean=FALSE)
x.ma
# fit MA model
www = "http://elena.aut.ac.nz/~pcowpert/ts/pounds_nz.dat"
x <- read.table(www, header = T)
x.ts <- ts(x, st = 1991, fr = 4)
x.ma <- arima(x.ts, order = c(0, 0, 1))
x.ma
# ARMA simulation
set.seed(1)
x = arima.sim(n = 10000, list(ar = -0.6, ma = 0.5))
coef(arima(x, order = c(1, 0, 1)))
AIC(arima(x, order = c(1, 0, 1)))
acf(resid(arima(x, order = c(1, 0, 1))))
www = "http://elena.aut.ac.nz/~pcowpert/ts/cbe.dat"
CBE = read.table(www, header = T)
Elec.ts = ts(CBE[,3], start = 1958, freq = 12)
Time = 1:length(Elec.ts)
Imth = cycle(Elec.ts)
Elec.lm = lm(log(Elec.ts) ~ Time + I(Time^2) + factor(Imth))
summary(Elec.lm)
acf(resid(Elec.lm))
# train ARMA model
best.order = c(0, 0, 0)
best.aic = Inf
for (i in c(0:2)){
  for (j in c(0:2)){
    fit.aic = AIC(arima(resid(Elec.lm), order = c(i, 0, j)))
    if (fit.aic < best.aic){
      best.order = c(i, 0, j)
      best.arma = arima(resid(Elec.lm), order = best.order)
      best.aic = fit.aic           
    }
  }
}
best.order
acf(resid(best.arma))
# prediction
new.time <- seq(length(Elec.ts), length = 36)
new.data <- data.frame(Time = new.time, Imth = rep(1:12, 3))
predict.lm <- predict(Elec.lm, new.data)
predict.arma <- predict(best.arma, n.ahead = 36)
elec.pred <- ts(exp(predict.lm + predict.arma$pred), start = 1991, freq = 12)
ts.plot(cbind(Elec.ts, elec.pred), lty = 1:2)
www = "http://elena.aut.ac.nz/~pcowpert/ts/wave.dat"
wave.dat <- read.table(www, header = T)
attach(wave.dat)
layout(1:3)
plot(as.ts(waveht), ylab = 'Wave height')
acf(waveht)
pacf (waveht)
wave.arma <- arima(waveht, order = c(4,0,4))
acf(wave.arma$res[-(1:4)])
pacf(wave.arma$res[-(1:4)])
hist(wave.arma$res[-(1:4)], xlab='height / mm', main='')
# train Non-seasonal ARIMA model
www = "http://elena.aut.ac.nz/~pcowpert/ts/cbe.dat"
CBE = read.table(www, header = T)
Elec.ts = ts(CBE[, 3], start = 1958, freq = 12)
layout(c(1, 1, 2, 3))
plot(Elec.ts)
plot(diff(Elec.ts))
plot(diff(log(Elec.ts)))
layout(c(1, 1))
# ARIMA(1, 1, 1) simulation and inference
set.seed(1)
x = w = rnorm(1000)
for (i in 3:1000){
  x[i] = 0.5 * x[i - 1] + x[i - 1] - 0.5 * x[i - 2] + w[i] + 0.3 * w[i - 1]
}
plot(x, type = "l")
arima(x = x, order  = c(1, 1, 1))
x = arima.sim(model = list(order = c(1, 1, 1), ar = 0.5, ma = 0.3), n = 1000)
arima(x = x, order  = c(1, 1, 1))
# IMA(1, 1)
Beer.ts = ts(CBE[, 2], start = 1958, freq = 12)
Beer.ima = arima(Beer.ts, order = c(0, 1, 1))
Beer.ima
acf(resid(Beer.ima))
# prediction
Beer.1991 = predict(Beer.ima, n.ahead = 12)
plot(Beer.1991$pred)
plot(Beer.1991$se)