Chapter 3: ARIMA Models

1. Autoregressive Moving Average Models

• AR, MA and ARMA models
• Backshift operator, mean square convergence
• Causal and explosive, stochastically equivalent
• Nonuniqueness and invertibility
• Matching, Taylor expansion for functions of complex variables
• Causal and invertible conditions for ARMA, parameter redundancy

2. Difference Equations

• General homogeneous difference Equations
• General solution
• ACF of AR(2), \(\psi\) weights for ARMA

3. Autocorrelation and Partial Autocorrelation

• Definition and motivation
• Basics of regression
• Partial correlation and conditional correlation
• ACF and PACF for AR and MA
• Interpret the ACF and PACF plots

4. Forecasting

• Mean square loss function and conditional mean
• Linear predictor and best linear predictor
• Prediction equations: derivation and matrix form, mean square prediction error
• The Durbin-Levinson algorithm: induction proof and projection proof
• The innovation algorithm: idea and comparison with the Durbin-Levinson
• Forecasting ARMA: forecasts based on finite and infinite past, truncated prediction

5. Estimation

• Method of moments: Yule-Walker equations, the Durbin-Levinson algorithm
• MLE and least square estimation: unconditional and conditional versions
• Large sample properties

6. Integrated Models for Nonstationary Data

• Trend model: deterministic and stochastic trends
• ARIMA models
• Forecasting

Code

Example 3.2 AR(p)

par(mfrow=c(2,1))                         
# in the expressions below, ~ is a space and == is equal
tsplot(arima.sim(list(order=c(1,0,0), ar=.9), n=100), ylab="x", main=(expression(AR(1)~~~phi==+.9))) 
tsplot(arima.sim(list(order=c(1,0,0), ar=-.9), n=100), ylab="x", main=(expression(AR(1)~~~phi==-.9))) 

Example 3.5 MA(1)

par(mfrow=c(2,1))                                   
tsplot(arima.sim(list(order=c(0,0,1), ma=.9), n=100), ylab="x", main=(expression(MA(1)~~~theta==+.9)))    
tsplot(arima.sim(list(order=c(0,0,1), ma=-.9), n=100), ylab="x", main=(expression(MA(1)~~~theta==-.9)))    

Example 3.11 AR(2) with complex roots

par(mfrow=c(1,3))
ACF = ARMAacf(ar=c(1,-.25), ma=0, 50)
plot(ACF, type="h", xlab="lag", main="equal real roots")
abline(h=0)

ACF = ARMAacf(ar=c(.75,-.125), ma=0, 50)
plot(ACF, type="h", xlab="lag", main="distinct real roots")
abline(h=0)

ACF = ARMAacf(ar=c(1.5,-.75), ma=0, 50)
plot(ACF, type="h", xlab="lag", main="complex roots")
abline(h=0)

Example 3.16 ACF and PACF of AR(p)

ar2.acf = ARMAacf(ar=c(1.5,-.75), ma=0, 24)[-1]
ar2.pacf = ARMAacf(ar=c(1.5,-.75), ma=0, 24, pacf=TRUE)
par(mfrow=c(1,2))
plot(ar2.acf, type="h", xlab="lag")
abline(h=0)
plot(ar2.pacf, type="h", xlab="lag")
abline(h=0)

Example 3.18 Preliminary analysis of Rec

tsplot(rec, ylab="", main="Recruitment") 

acf2(rec, 48)     # will produce values and a graphic 

##         ACF  PACF
##  [1,]  0.92  0.92
##  [2,]  0.78 -0.44
##  [3,]  0.63 -0.05
##  [4,]  0.48 -0.02
##  [5,]  0.36  0.07
##  [6,]  0.26 -0.03
##  [7,]  0.18 -0.03
##  [8,]  0.13  0.04
##  [9,]  0.09  0.05
## [10,]  0.07 -0.02
## [11,]  0.06 -0.05
## [12,]  0.02 -0.14
## [13,] -0.04 -0.15
## [14,] -0.12 -0.05
## [15,] -0.19  0.05
## [16,] -0.24  0.01
## [17,] -0.27  0.01
## [18,] -0.27  0.02
## [19,] -0.24  0.09
## [20,] -0.19  0.11
## [21,] -0.11  0.03
## [22,] -0.03 -0.03
## [23,]  0.03 -0.01
## [24,]  0.06 -0.07
## [25,]  0.06 -0.12
## [26,]  0.02 -0.03
## [27,] -0.02  0.05
## [28,] -0.06 -0.08
## [29,] -0.09 -0.04
## [30,] -0.12 -0.03
## [31,] -0.13  0.06
## [32,] -0.11  0.05
## [33,] -0.05  0.15
## [34,]  0.02  0.09
## [35,]  0.08 -0.04
## [36,]  0.12 -0.10
## [37,]  0.10 -0.09
## [38,]  0.06 -0.02
## [39,]  0.01  0.05
## [40,] -0.02  0.08
## [41,] -0.03 -0.02
## [42,] -0.03 -0.01
## [43,] -0.02 -0.02
## [44,]  0.01  0.05
## [45,]  0.06  0.01
## [46,]  0.12  0.05
## [47,]  0.17  0.08
## [48,]  0.20 -0.04

(regr = ar.ols(rec, order=2, demean=F, intercept=TRUE))  # regression

## 
## Call:
## ar.ols(x = rec, order.max = 2, demean = F, intercept = TRUE)
## 
## Coefficients:
##       1        2  
##  1.3541  -0.4632  
## 
## Intercept: 6.737 (1.111) 
## 
## Order selected 2  sigma^2 estimated as  89.72

regr$asy.se.coef  # standard errors                             

## $x.mean
## [1] 1.110599
## 
## $ar
## [1] 0.04178901 0.04187942

Example 3.25 Forecasting the Recruitment Series

regr = ar.ols(rec, order=2, demean=FALSE, intercept=TRUE)
fore = predict(regr, n.ahead=24)
ts.plot(rec, fore$pred, col=1:2, xlim=c(1980,1990), ylab="Recruitment")
lines(fore$pred, type="p", col=2)
lines(fore$pred+fore$se, lty="dashed", col=4)
lines(fore$pred-fore$se, lty="dashed", col=4)

Example 3.28 Yule-Walker Estimation of the Recruitment Series

rec.yw = ar.yw(rec, order=2)
rec.yw$x.mean  # = 62.26278 (mean estimate)

## [1] 62.26278

rec.yw$ar      # = 1.3315874, -.4445447  (parameter estimates)

## [1]  1.3315874 -0.4445447

sqrt(diag(rec.yw$asy.var.coef))  # = .04222637, .04222637  (standard errors)

## [1] 0.04222637 0.04222637

rec.yw$var.pred  # = 94.79912 (error variance estimate)

## [1] 94.79912

rec.pr = predict(rec.yw, n.ahead=24)
U = rec.pr$pred + rec.pr$se
L = rec.pr$pred - rec.pr$se
minx = min(rec,L); maxx = max(rec,U)
ts.plot(rec, rec.pr$pred, xlim=c(1980,1990), ylim=c(minx,maxx)) 
lines(rec.pr$pred, col="red", type="o")
lines(U, col="blue", lty="dashed")
lines(L, col="blue", lty="dashed")

Example 3.31 MLE for the Recruitment Series

rec.mle = ar.mle(rec, order=2)
rec.mle$x.mean

## [1] 62.26153

rec.mle$ar

## [1]  1.3512809 -0.4612736

sqrt(diag(rec.mle$asy.var.coef))

## [1] 0.04099159 0.04099159

rec.mle$var.pred

## [1] 89.33597

rec.pr = predict(rec.mle, n.ahead=24)
U = rec.pr$pred + rec.pr$se
L = rec.pr$pred - rec.pr$se
minx = min(rec,L); maxx = max(rec,U)
ts.plot(rec, rec.pr$pred, xlim=c(1980,1990), ylim=c(minx,maxx)) 
lines(rec.pr$pred, col="red", type="o")
lines(U, col="blue", lty="dashed")
lines(L, col="blue", lty="dashed")

Example 3.33 Fitting the Glacial Varve Series

x = diff(log(varve))
# Evaluate Sc on a Grid
c(0) -> w -> z
c() -> Sc -> Sz -> Szw
num = length(x)
th = seq(-.3,-.94,-.01)
for (p in 1:length(th)){
  for (i in 2:num){ w[i] = x[i]-th[p]*w[i-1] }
  Sc[p] = sum(w^2) }
plot(th, Sc, type="l", ylab=expression(S[c](theta)), xlab=expression(theta),
     lwd=2)
# Gauss-Newton Estimation
r = acf(x, lag=1, plot=FALSE)$acf[-1]
rstart = (1-sqrt(1-4*(r^2)))/(2*r) # from (3.105)
c(0) -> w -> z
c() -> Sc -> Sz -> Szw -> para
niter = 12
para[1] = rstart
for (p in 1:niter){
  for (i in 2:num){ w[i] = x[i]-para[p]*w[i-1]
  z[i] = w[i-1]-para[p]*z[i-1] }
  Sc[p] = sum(w^2)
  Sz[p] = sum(z^2)
  Szw[p] = sum(z*w)
  para[p+1] = para[p] + Szw[p]/Sz[p] }

round(cbind(iteration=0:(niter-1), thetahat=para[1:niter] , Sc , Sz ), 3)

##       iteration thetahat      Sc      Sz
##  [1,]         0   -0.495 158.739 171.240
##  [2,]         1   -0.668 150.747 235.266
##  [3,]         2   -0.733 149.264 300.562
##  [4,]         3   -0.756 149.031 336.823
##  [5,]         4   -0.766 148.990 354.173
##  [6,]         5   -0.769 148.982 362.167
##  [7,]         6   -0.771 148.980 365.801
##  [8,]         7   -0.772 148.980 367.446
##  [9,]         8   -0.772 148.980 368.188
## [10,]         9   -0.772 148.980 368.522
## [11,]        10   -0.773 148.980 368.673
## [12,]        11   -0.773 148.980 368.741

abline(v = para[1:12], lty=2)
points(para[1:12], Sc[1:12], pch=16)

Example 3.38 IMA(1,1) and EWMA

set.seed(666)    
x = arima.sim(list(order = c(0,1,1), ma = -0.8), n = 100)
(x.ima = HoltWinters(x, beta=FALSE, gamma=FALSE))  # α is 1-λ here

## Holt-Winters exponential smoothing without trend and without seasonal component.
## 
## Call:
## HoltWinters(x = x, beta = FALSE, gamma = FALSE)
## 
## Smoothing parameters:
##  alpha: 0.1663072
##  beta : FALSE
##  gamma: FALSE
## 
## Coefficients:
##        [,1]
## a -2.241533

plot(x.ima)