ARIMA Models

Time Series analysis using AutoRegressive Integrated Moving Average models is appropriate for stationary data. The ARIMA(p,d,q) is given by:

\[ \phi(B) (1-B)^d X_t = \theta(B) Z_t, \qquad Z_t \sim WN(0,\sigma^2) \]

where \(\Phi(B) = 1-\phi_1X_{t-1} - \ldots - \phi_pX_{t-p}\) and \(\Theta(B) = 1+\theta_1Z_{t-1} + \ldots + \theta_qZ_{t-q}\)

Quarterly U.S. GNP Gross National Product

library(car)
library(ltsa)
library(itsmr)
library(astsa)
library(itsmr)
library(forecast)

# plot
plot(gnp,type="o")

# ACF plot
acf(gnp,lag.max = 100)

# data model
M=c("diff",1)
r = Resid(gnp,M)
test(r)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     91.98         0 *
McLeod-Li Q                Q ~ chisq(20)     51.48     1e-04 *
Turning points T  (T-146.7)/6.3 ~ N(0,1)       145    0.7899
Diff signs S      (S-110.5)/4.3 ~ N(0,1)       115    0.2965
Rank P        (P-12265.5)/553.1 ~ N(0,1)     15350         0 *

Test stationary

# Augmented Dickey-Fuller Test
# H0: X_t is a non-stationary time series
# H1: X_t is a stationary time series

library(aTSA)
adf.test(r)

Augmented Dickey-Fuller Test 
alternative: stationary 
 
Type 1: no drift no trend 
     lag   ADF p.value
[1,]   0 -9.76    0.01
[2,]   1 -6.70    0.01
[3,]   2 -5.98    0.01
[4,]   3 -5.52    0.01
[5,]   4 -5.52    0.01
Type 2: with drift no trend 
     lag   ADF p.value
[1,]   0 -9.74    0.01
[2,]   1 -6.69    0.01
[3,]   2 -5.97    0.01
[4,]   3 -5.50    0.01
[5,]   4 -5.51    0.01
Type 3: with drift and trend 
     lag    ADF p.value
[1,]   0 -10.64    0.01
[2,]   1  -7.43    0.01
[3,]   2  -6.79    0.01
[4,]   3  -6.40    0.01
[5,]   4  -6.59    0.01
---- 
Note: in fact, p.value = 0.01 means p.value <= 0.01 

# Results in stationary data. No differencing is needed, then d=0.

ARIMA models

# p = 1, 2
# q = 1, 2, 3
# ARIMA(2,0,1) == ARMA(2,1)
a21=arima(x = r,order = c(2,0,1))
test(a21$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     20.89    0.4038
McLeod-Li Q                Q ~ chisq(20)     72.35         0 *
Turning points T  (T-146.7)/6.3 ~ N(0,1)       146    0.9151
Diff signs S      (S-110.5)/4.3 ~ N(0,1)       117    0.1316
Rank P        (P-12265.5)/553.1 ~ N(0,1)     13768    0.0066 *

# ARIMA(1,0,1) == ARMA(1,1)
a11=arima(x = r,order = c(1,0,1))
test(a11$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)      24.2    0.2337
McLeod-Li Q                Q ~ chisq(20)     72.62         0 *
Turning points T  (T-146.7)/6.3 ~ N(0,1)       156    0.1358
Diff signs S      (S-110.5)/4.3 ~ N(0,1)       116     0.202
Rank P        (P-12265.5)/553.1 ~ N(0,1)     13745    0.0075 *

# ARIMA(2,0,0) == ARMA(2,0)
a2=arima(x = r,order = c(2,0,0))
test(a2$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     20.72     0.414
McLeod-Li Q                Q ~ chisq(20)     72.18         0 *
Turning points T  (T-146.7)/6.3 ~ N(0,1)       146    0.9151
Diff signs S      (S-110.5)/4.3 ~ N(0,1)       116     0.202
Rank P        (P-12265.5)/553.1 ~ N(0,1)     13778    0.0062 *

#abest=autofit(r,p = 0:2,q=0:2)
#test(abest)

We will try a transformation to reduce the effect of variance change.

l=powerTransform(gnp)
l$roundlam

gnp 
  0 

## we can use log() transformation

plot(log(gnp),type="o")

acf(log(gnp),lag.max = 200)

# data model
M=c("log","diff",1)
r = Resid(gnp,M)
test(r)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)      70.3         0 *
McLeod-Li Q                Q ~ chisq(20)     41.93    0.0028 *
Turning points T  (T-146.7)/6.3 ~ N(0,1)       145    0.7899
Diff signs S      (S-110.5)/4.3 ~ N(0,1)       114    0.4168
Rank P        (P-12265.5)/553.1 ~ N(0,1)     11482    0.1566

# Stationary
adf.test(r) # pvalue <0.05 => data is stationary

Augmented Dickey-Fuller Test 
alternative: stationary 
 
Type 1: no drift no trend 
     lag    ADF p.value
[1,]   0 -10.32    0.01
[2,]   1  -7.77    0.01
[3,]   2  -7.64    0.01
[4,]   3  -7.59    0.01
[5,]   4  -7.38    0.01
Type 2: with drift no trend 
     lag    ADF p.value
[1,]   0 -10.29    0.01
[2,]   1  -7.76    0.01
[3,]   2  -7.62    0.01
[4,]   3  -7.57    0.01
[5,]   4  -7.36    0.01
Type 3: with drift and trend 
     lag    ADF p.value
[1,]   0 -10.34    0.01
[2,]   1  -7.82    0.01
[3,]   2  -7.68    0.01
[4,]   3  -7.65    0.01
[5,]   4  -7.44    0.01
---- 
Note: in fact, p.value = 0.01 means p.value <= 0.01 

# ARIMA(1,0,1) == ARMA(1,1)
arma11=arima(x = r,order = c(1,0,1))
test(arma11$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     24.99    0.2019
McLeod-Li Q                Q ~ chisq(20)     29.54    0.0776
Turning points T  (T-146.7)/6.3 ~ N(0,1)       161     0.022 *
Diff signs S      (S-110.5)/4.3 ~ N(0,1)       117    0.1316
Rank P        (P-12265.5)/553.1 ~ N(0,1)     11709    0.3144

# ARIMA(2,0,1) == ARMA(2,1)
arma21=arima(x = r,order = c(2,0,1))
test(arma21$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     22.59    0.3095
McLeod-Li Q                Q ~ chisq(20)     30.72     0.059
Turning points T  (T-146.7)/6.3 ~ N(0,1)       155    0.1829
Diff signs S      (S-110.5)/4.3 ~ N(0,1)       116     0.202
Rank P        (P-12265.5)/553.1 ~ N(0,1)     11712     0.317

# ARIMA(2,0,0) == ARMA(2,0) == AR(2)
ar2=arima(x = r,order = c(2,0,0))
test(ar2$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     24.48    0.2219
McLeod-Li Q                Q ~ chisq(20)     30.72    0.0591
Turning points T  (T-146.7)/6.3 ~ N(0,1)       159    0.0487 *
Diff signs S      (S-110.5)/4.3 ~ N(0,1)       116     0.202
Rank P        (P-12265.5)/553.1 ~ N(0,1)     11728    0.3312

# ARIMA(1,0,2) == ARMA(1,2)
arma12=arima(x = r,order = c(1,0,2))
test(arma12$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     20.29      0.44
McLeod-Li Q                Q ~ chisq(20)     29.63     0.076
Turning points T  (T-146.7)/6.3 ~ N(0,1)       153    0.3114
Diff signs S      (S-110.5)/4.3 ~ N(0,1)       115    0.2965
Rank P        (P-12265.5)/553.1 ~ N(0,1)     11724    0.3276

# ARIMA(0,0,2) == MA(2)
ma2=arima(x = r,order = c(0,0,2))
test(ma2$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     21.58    0.3637
McLeod-Li Q                Q ~ chisq(20)     29.61    0.0764
Turning points T  (T-146.7)/6.3 ~ N(0,1)       153    0.3114
Diff signs S      (S-110.5)/4.3 ~ N(0,1)       115    0.2965
Rank P        (P-12265.5)/553.1 ~ N(0,1)     11672    0.2833

## Comparing AICs
c(MA2=ma2$aic, ARMA12=arma12$aic, ARMA21=arma21$aic)

      MA2    ARMA12    ARMA21 
-1431.929 -1430.948 -1429.673 

The equation of MA(2):

\[ (1-B)\log(X_t) = 0.3Z_{t-1} + 0.2Z_{t-2} + Z_t, \qquad Z_t \sim WN(0,8.919e-05) \]

Example Quarterly Time Series of the Number of Australian Residents

plot.ts(austres, type="o")

plota(austres)

M=c("diff",1)
e = Resid(austres,M)
test(e)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)    137.43         0 *
McLeod-Li Q                Q ~ chisq(20)     40.76     0.004 *
Turning points T   (T-57.3)/3.9 ~ N(0,1)        52     0.173
Diff signs S       (S-43.5)/2.7 ~ N(0,1)        45    0.5818
Rank P           (P-1914)/138.7 ~ N(0,1)      2321    0.0033 *

aTSA::adf.test(e)

Augmented Dickey-Fuller Test 
alternative: stationary 
 
Type 1: no drift no trend 
     lag   ADF p.value
[1,]   0 -4.63  0.0100
[2,]   1 -3.56  0.0100
[3,]   2 -2.49  0.0148
[4,]   3 -1.87  0.0625
Type 2: with drift no trend 
     lag   ADF p.value
[1,]   0 -4.60   0.010
[2,]   1 -3.53   0.010
[3,]   2 -2.47   0.145
[4,]   3 -1.85   0.386
Type 3: with drift and trend 
     lag   ADF p.value
[1,]   0 -4.85  0.0100
[2,]   1 -3.87  0.0195
[3,]   2 -2.71  0.2795
[4,]   3 -1.98  0.5786
---- 
Note: in fact, p.value = 0.01 means p.value <= 0.01 

# since the non-stationary was not rejected. We need to do another differencing operator.

aTSA::adf.test(diff(e))

Augmented Dickey-Fuller Test 
alternative: stationary 
 
Type 1: no drift no trend 
     lag    ADF p.value
[1,]   0 -12.80    0.01
[2,]   1 -10.54    0.01
[3,]   2  -8.92    0.01
[4,]   3  -6.35    0.01
Type 2: with drift no trend 
     lag    ADF p.value
[1,]   0 -12.73    0.01
[2,]   1 -10.49    0.01
[3,]   2  -8.87    0.01
[4,]   3  -6.32    0.01
Type 3: with drift and trend 
     lag    ADF p.value
[1,]   0 -12.65    0.01
[2,]   1 -10.42    0.01
[3,]   2  -8.81    0.01
[4,]   3  -6.28    0.01
---- 
Note: in fact, p.value = 0.01 means p.value <= 0.01 

test(diff(e))

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     49.28     3e-04 *
McLeod-Li Q                Q ~ chisq(20)        25    0.2015
Turning points T   (T-56.7)/3.9 ~ N(0,1)        56     0.864
Diff signs S         (S-43)/2.7 ~ N(0,1)        46    0.2679
Rank P         (P-1870.5)/136.4 ~ N(0,1)      1917    0.7331

# Fit ARIMA(p,d=1,q)
a1=arima(x = austres,order = c(2,1,2))
test(a1$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     26.56     0.148
McLeod-Li Q                Q ~ chisq(20)     15.09    0.7714
Turning points T     (T-58)/3.9 ~ N(0,1)        54    0.3096
Diff signs S         (S-44)/2.7 ~ N(0,1)        45     0.715
Rank P           (P-1958)/141.1 ~ N(0,1)      1966    0.9548

# Fit ARIMA(p,d=2,q)
a2=arima(x = austres,order = c(2,2,2))
test(a2$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     14.64    0.7968
McLeod-Li Q                Q ~ chisq(20)     17.71    0.6064
Turning points T     (T-58)/3.9 ~ N(0,1)        58         1
Diff signs S         (S-44)/2.7 ~ N(0,1)        47    0.2733
Rank P           (P-1958)/141.1 ~ N(0,1)      2052    0.5053

# Fit ARIMA(p,d=1,q) - not iid
a3=arima(x = austres,order = c(1,1,0))
test(a3$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     48.12     4e-04 *
McLeod-Li Q                Q ~ chisq(20)     25.49    0.1834
Turning points T     (T-58)/3.9 ~ N(0,1)        57    0.7995
Diff signs S         (S-44)/2.7 ~ N(0,1)        47    0.2733
Rank P           (P-1958)/141.1 ~ N(0,1)      1944     0.921

# Fit ARIMA(p,d=1,q) - not iid
a4=arima(x = austres,order = c(1,1,1))
test(a4$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)      32.2    0.0412 *
McLeod-Li Q                Q ~ chisq(20)     13.49    0.8554
Turning points T     (T-58)/3.9 ~ N(0,1)        52    0.1275
Diff signs S         (S-44)/2.7 ~ N(0,1)        45     0.715
Rank P           (P-1958)/141.1 ~ N(0,1)      2017    0.6758

# Fit ARIMA(p,d=2,q) - not iid
a5=arima(x = austres,order = c(1,2,0))
test(a5$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     52.83     1e-04 *
McLeod-Li Q                Q ~ chisq(20)     23.67    0.2571
Turning points T     (T-58)/3.9 ~ N(0,1)        52    0.1275
Diff signs S         (S-44)/2.7 ~ N(0,1)        45     0.715
Rank P           (P-1958)/141.1 ~ N(0,1)      2049    0.5189

# Fit ARIMA(p,d=2,q) 
a6=arima(x = austres,order = c(1,2,1))
test(a6$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     26.65    0.1453
McLeod-Li Q                Q ~ chisq(20)     13.52    0.8539
Turning points T     (T-58)/3.9 ~ N(0,1)        56    0.6115
Diff signs S         (S-44)/2.7 ~ N(0,1)        45     0.715
Rank P           (P-1958)/141.1 ~ N(0,1)      2095    0.3315

# Fit ARIMA(p,d=2,q) = IMA(2,1)
a7=arima(x = austres,order = c(0,2,1))
test(a7$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     28.67    0.0945
McLeod-Li Q                Q ~ chisq(20)     12.84    0.8843
Turning points T     (T-58)/3.9 ~ N(0,1)        54    0.3096
Diff signs S         (S-44)/2.7 ~ N(0,1)        46    0.4652
Rank P           (P-1958)/141.1 ~ N(0,1)      2108    0.2877

c(a1$aic,a2$aic,a6$aic,a7$aic)

[1] 668.2778 654.5503 654.4183 652.9911

The equation of ARIMA(2,1,2):

\[ (1-0.6843B-0.3121B^2)(1-B)(X_t) = (1 - 0.1961B -0.2975B^2)Z_t, \qquad Z_t \sim WN(0,99.47) \]

The equation of IMA(2,1):

\[ (1-B)^2 X_t = -0.59Z_{t-1} + Z_t; Z_t \sim WN(0,\sigma^2 = 101.17) \]