SARIMA Models

Time Series analysis using Seasonal AutoRegressive Integrated Moving Average models. The SARIMA(p,d,q)x(P,D,Q)s is given by:

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

Chicken Data

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

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

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

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

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)    236.56         0 *
McLeod-Li Q                Q ~ chisq(20)     69.81         0 *
Turning points T    (T-110)/5.4 ~ N(0,1)        85         0 *
Diff signs S         (S-83)/3.7 ~ N(0,1)        89    0.1088
Rank P         (P-6930.5)/361.3 ~ N(0,1)      6323    0.0927

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 -5.00    0.01
[2,]   1 -5.54    0.01
[3,]   2 -5.42    0.01
[4,]   3 -4.45    0.01
[5,]   4 -4.15    0.01
Type 2: with drift no trend 
     lag   ADF p.value
[1,]   0 -4.98    0.01
[2,]   1 -5.53    0.01
[3,]   2 -5.41    0.01
[4,]   3 -4.43    0.01
[5,]   4 -4.14    0.01
Type 3: with drift and trend 
     lag   ADF p.value
[1,]   0 -5.03    0.01
[2,]   1 -5.58    0.01
[3,]   2 -5.46    0.01
[4,]   3 -4.47    0.01
[5,]   4 -4.17    0.01
---- 
Note: in fact, p.value = 0.01 means p.value <= 0.01 

# Results in stationary data. No differencing is needed.

SARIMA models

# Potential Values
# based on the PACF
# p = 0, 1, 2
# P = 0, 1, 2
# based on the ACF
# q = 0, 1, 2, 3, 9, 10, 11
# Q = 0, 1, 2,
# Eliminate Trend and Season
# d = 1
# D = 1

# SARIMA (1,1,1)x(1,1,1)12
fit1 = arima(chicken,
            order = c(1,1,1),
            seasonal = list(order=c(1,1,1),period=12))
test(fit1$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     26.41    0.1526
McLeod-Li Q                Q ~ chisq(20)     26.03     0.165
Turning points T  (T-118.7)/5.6 ~ N(0,1)       113     0.314
Diff signs S       (S-89.5)/3.9 ~ N(0,1)        90    0.8976
Rank P           (P-8055)/404.2 ~ N(0,1)      7885     0.674

confint(fit1)

           2.5 %     97.5 %
ar1   0.46878966  0.7763441
ma1   0.06514496  0.4192910
sar1 -0.16334192  0.2108877
sma1 -1.09510965 -0.6989791

fit1

Call:
arima(x = chicken, order = c(1, 1, 1), seasonal = list(order = c(1, 1, 1), period = 12))

Coefficients:
         ar1     ma1    sar1     sma1
      0.6226  0.2422  0.0238  -0.8970
s.e.  0.0785  0.0903  0.0955   0.1011

sigma^2 estimated as 0.3296:  log likelihood = -154,  aic = 318

# SARIMA (1,1,1)x(2,1,1)12
fit2 = arima(chicken,
            order = c(1,1,1),
            seasonal = list(order=c(2,1,1),period=12))
test(fit2$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     24.24    0.2322
McLeod-Li Q                Q ~ chisq(20)     25.54    0.1816
Turning points T  (T-118.7)/5.6 ~ N(0,1)       115    0.5147
Diff signs S       (S-89.5)/3.9 ~ N(0,1)        92    0.5198
Rank P           (P-8055)/404.2 ~ N(0,1)      7864    0.6365

confint(fit2)

           2.5 %      97.5 %
ar1   0.45832127  0.77089296
ma1   0.04960726  0.40221312
sar1 -0.28700500  0.10441229
sar2 -0.41013311 -0.05367363
sma1 -0.93736485 -0.59308042

fit2

Call:
arima(x = chicken, order = c(1, 1, 1), seasonal = list(order = c(2, 1, 1), period = 12))

Coefficients:
         ar1     ma1     sar1     sar2     sma1
      0.6146  0.2259  -0.0913  -0.2319  -0.7652
s.e.  0.0797  0.0900   0.0999   0.0909   0.0878

sigma^2 estimated as 0.3227:  log likelihood = -151.14,  aic = 314.28

# Auto fit
bestmodel = auto.arima(chicken, trace= TRUE)

 Fitting models using approximations to speed things up...

 ARIMA(2,1,2)(1,0,1)[12] with drift         : Inf
 ARIMA(0,1,0)            with drift         : 512.1327
 ARIMA(1,1,0)(1,0,0)[12] with drift         : 357.8553
 ARIMA(0,1,1)(0,0,1)[12] with drift         : 399.0111
 ARIMA(0,1,0)                               : 521.4674
 ARIMA(1,1,0)            with drift         : 380.5252
 ARIMA(1,1,0)(2,0,0)[12] with drift         : 361.1736
 ARIMA(1,1,0)(1,0,1)[12] with drift         : Inf
 ARIMA(1,1,0)(0,0,1)[12] with drift         : 361.7262
 ARIMA(1,1,0)(2,0,1)[12] with drift         : Inf
 ARIMA(0,1,0)(1,0,0)[12] with drift         : 496.251
 ARIMA(2,1,0)(1,0,0)[12] with drift         : 349.6098
 ARIMA(2,1,0)            with drift         : 361.2204
 ARIMA(2,1,0)(2,0,0)[12] with drift         : 352.9452
 ARIMA(2,1,0)(1,0,1)[12] with drift         : Inf
 ARIMA(2,1,0)(0,0,1)[12] with drift         : 347.5637
 ARIMA(2,1,0)(0,0,2)[12] with drift         : 348.8669
 ARIMA(2,1,0)(1,0,2)[12] with drift         : Inf
 ARIMA(3,1,0)(0,0,1)[12] with drift         : 346.8076
 ARIMA(3,1,0)            with drift         : 359.6763
 ARIMA(3,1,0)(1,0,1)[12] with drift         : Inf
 ARIMA(3,1,0)(0,0,2)[12] with drift         : 348.5706
 ARIMA(3,1,0)(1,0,0)[12] with drift         : 350.3946
 ARIMA(3,1,0)(1,0,2)[12] with drift         : Inf
 ARIMA(4,1,0)(0,0,1)[12] with drift         : 349.2894
 ARIMA(3,1,1)(0,0,1)[12] with drift         : 347.6645
 ARIMA(2,1,1)(0,0,1)[12] with drift         : 346.754
 ARIMA(2,1,1)            with drift         : 359.2506
 ARIMA(2,1,1)(1,0,1)[12] with drift         : Inf
 ARIMA(2,1,1)(0,0,2)[12] with drift         : 348.4092
 ARIMA(2,1,1)(1,0,0)[12] with drift         : 350.2436
 ARIMA(2,1,1)(1,0,2)[12] with drift         : Inf
 ARIMA(1,1,1)(0,0,1)[12] with drift         : 354.7731
 ARIMA(2,1,2)(0,0,1)[12] with drift         : 348.0198
 ARIMA(1,1,2)(0,0,1)[12] with drift         : 348.7836
 ARIMA(3,1,2)(0,0,1)[12] with drift         : 349.5839
 ARIMA(2,1,1)(0,0,1)[12]                    : 347.2775

 Now re-fitting the best model(s) without approximations...

 ARIMA(2,1,1)(0,0,1)[12] with drift         : 351.5017

 Best model: ARIMA(2,1,1)(0,0,1)[12] with drift         

test(bestmodel$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     20.39    0.4339
McLeod-Li Q                Q ~ chisq(20)     25.62    0.1789
Turning points T  (T-118.7)/5.6 ~ N(0,1)       126    0.1926
Diff signs S       (S-89.5)/3.9 ~ N(0,1)        96    0.0942
Rank P           (P-8055)/404.2 ~ N(0,1)      8124    0.8644

# SARIMA (2,1,1)x(0,1,1)12
fit3 = arima(chicken,
            order = c(2,1,1),
            seasonal = list(order=c(0,1,1),period=12))
test(fit3$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     24.19    0.2342
McLeod-Li Q                Q ~ chisq(20)     24.05    0.2403
Turning points T  (T-118.7)/5.6 ~ N(0,1)       109    0.0859
Diff signs S       (S-89.5)/3.9 ~ N(0,1)        90    0.8976
Rank P           (P-8055)/404.2 ~ N(0,1)      7886    0.6758

confint(fit3)

          2.5 %     97.5 %
ar1   0.3663525  1.3562175
ar2  -0.5707739  0.1810233
ma1  -0.4705880  0.5126463
sma1 -1.0460511 -0.7262906

fit3

Call:
arima(x = chicken, order = c(2, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))

Coefficients:
         ar1      ar2     ma1     sma1
      0.8613  -0.1949  0.0210  -0.8862
s.e.  0.2525   0.1918  0.2508   0.0816

sigma^2 estimated as 0.329:  log likelihood = -153.62,  aic = 317.24

# SARIMA (1,1,1)x(0,1,1)12
fit4 = arima(chicken,
            order = c(1,1,1),
            seasonal = list(order=c(0,1,1),period=12))
test(fit4$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     26.57    0.1478
McLeod-Li Q                Q ~ chisq(20)     25.56    0.1808
Turning points T  (T-118.7)/5.6 ~ N(0,1)       113     0.314
Diff signs S       (S-89.5)/3.9 ~ N(0,1)        90    0.8976
Rank P           (P-8055)/404.2 ~ N(0,1)      7871    0.6489

confint(fit4)

           2.5 %     97.5 %
ar1   0.46810473  0.7760082
ma1   0.06587317  0.4202834
sma1 -1.04237314 -0.7267588

fit4

Call:
arima(x = chicken, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))

Coefficients:
         ar1     ma1     sma1
      0.6221  0.2431  -0.8846
s.e.  0.0785  0.0904   0.0805

sigma^2 estimated as 0.3309:  log likelihood = -154.03,  aic = 316.06

# SARIMA (1,1,0)x(0,1,1)12
fit5 = arima(chicken,
            order = c(1,1,0),
            seasonal = list(order=c(0,1,1),period=12))
test(fit5$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)      45.3     0.001 *
McLeod-Li Q                Q ~ chisq(20)     27.08     0.133
Turning points T  (T-118.7)/5.6 ~ N(0,1)       105    0.0152 *
Diff signs S       (S-89.5)/3.9 ~ N(0,1)        93    0.3675
Rank P           (P-8055)/404.2 ~ N(0,1)      7762    0.4685

confint(fit5)

          2.5 %     97.5 %
ar1   0.6273770  0.8369338
sma1 -0.9925754 -0.7215922

#fit5 is not adequate
fit5

Call:
arima(x = chicken, order = c(1, 1, 0), seasonal = list(order = c(0, 1, 1), period = 12))

Coefficients:
         ar1     sma1
      0.7322  -0.8571
s.e.  0.0535   0.0691

sigma^2 estimated as 0.3486:  log likelihood = -157.22,  aic = 320.44

## Compare
AIC.values = c(
      AutoModel = bestmodel$aic,
       fit1 = fit1$aic,
       fit2 = fit2$aic,
       fit3 = fit3$aic,
       fit4 = fit4$aic
      )
AIC.values

AutoModel      fit1      fit2      fit3      fit4 
 351.0134  317.9969  314.2806  317.2427  316.0594 

The equation of fit3: SARIMA(2,1,1)x(0,1,1)12 is:

\[ (1-0.861B+0.195B^2)(1-B)(1-B^{12})X_t = (1 + 0.021B)(1-0.886B^{12})Z_t, \qquad Z_t \sim WN(0,0.329) \]

The equation of fit2: SARIMA (1,1,1)x(2,1,1)12 is:

\[ (1-0.614B)(1+0.091B^{12}+0.232B^{24})(1-B)(1-B^{12})X_t = (1 + 0.226B)(1-0.765B^{12})Z_t, \qquad Z_t \sim WN(0,0.3227) \]

Forecasting with Time Series

chicken_train = chicken[1:178]


# SARIMA (1,1,1)x(1,1,1)12
fit1 = arima(chicken_train,
            order = c(1,1,1),
            seasonal = list(order=c(1,1,1),period=12))

# SARIMA (1,1,1)x(2,1,1)12
fit2 = arima(chicken_train,
            order = c(1,1,1),
            seasonal = list(order=c(2,1,1),period=12))

# Auto fit
bestmodel = auto.arima(chicken_train, trace= TRUE)

 Fitting models using approximations to speed things up...

 ARIMA(2,1,2) with drift         : 358.0721
 ARIMA(0,1,0) with drift         : 507.8724
 ARIMA(1,1,0) with drift         : 377.5175
 ARIMA(0,1,1) with drift         : 412.1149
 ARIMA(0,1,0)                    : 517.4894
 ARIMA(1,1,2) with drift         : 361.6613
 ARIMA(2,1,1) with drift         : 355.8188
 ARIMA(1,1,1) with drift         : 368.1617
 ARIMA(2,1,0) with drift         : 358.1683
 ARIMA(3,1,1) with drift         : 356.2727
 ARIMA(3,1,0) with drift         : 356.2762
 ARIMA(3,1,2) with drift         : Inf
 ARIMA(2,1,1)                    : 359.48

 Now re-fitting the best model(s) without approximations...

 ARIMA(2,1,1) with drift         : 360.3365

 Best model: ARIMA(2,1,1) with drift         

bestmodel = arima(chicken_train, order = c(2,1,1))
# SARIMA (2,1,1)x(0,1,1)12
fit3 = arima(chicken_train,
            order = c(2,1,1),
            seasonal = list(order=c(0,1,1),period=12))

# SARIMA (1,1,1)x(0,1,1)12
fit4 = arima(chicken_train,
            order = c(1,1,1),
            seasonal = list(order=c(0,1,1),period=12))

# SARIMA (1,1,0)x(0,1,1)12
fit5 = arima(chicken_train,
            order = c(1,1,0),
            seasonal = list(order=c(0,1,1),period=12))

Forecasting using the forecast function.

f1 = aTSA::forecast(object = fit1,lead = 2,output = FALSE)
f2 = aTSA::forecast(object = fit2,lead = 2,output = FALSE)
f3 = aTSA::forecast(object = fit3,lead = 2,output = FALSE)
f4 = aTSA::forecast(object = fit4,lead = 2,output = FALSE)
f5 = aTSA::forecast(object = fit5,lead = 2,output = FALSE)
f.auto = aTSA::forecast(object = bestmodel,lead = 2,output = FALSE)

Evaluate the performance of the forecasting:

MAPE1 = mean(abs(chicken[179:180]-f1[,2])/chicken[179:180])*100
MAPE2 = mean(abs(chicken[179:180]-f2[,2])/chicken[179:180])*100
MAPE3 = mean(abs(chicken[179:180]-f3[,2])/chicken[179:180])*100
MAPE4 = mean(abs(chicken[179:180]-f4[,2])/chicken[179:180])*100
MAPE5 = mean(abs(chicken[179:180]-f5[,2])/chicken[179:180])*100
MAPEauto = mean(abs(chicken[179:180]-f.auto[,2])/chicken[179:180])*100


MAPEs = c(fit1_mape = MAPE1,
          fit2_mape = MAPE2,
          fit3_mape = MAPE3,
          fit4_mape = MAPE4,
          fit5_mape = MAPE5,
          fitauto_mape = MAPEauto)
MAPEs

   fit1_mape    fit2_mape    fit3_mape    fit4_mape    fit5_mape fitauto_mape 
   1.0718671    0.8199072    1.1577641    1.0802866    1.0175245    0.8646009 

Wine Data: using Fixed Argument in arima()

plot(wine,type="o")

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

# data model
M=c("diff",1,"diff",12)
r = Resid(sqrt(wine),M) # sqrt based on Box-Cox Transformation
test(r)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     70.12         0 *
McLeod-Li Q                Q ~ chisq(20)     23.68    0.2566
Turning points T   (T-84.7)/4.8 ~ N(0,1)        86    0.7792
Diff signs S         (S-64)/3.3 ~ N(0,1)        64         1
Rank P           (P-4128)/245.6 ~ N(0,1)      4112    0.9481

adf.test(r)

Augmented Dickey-Fuller Test 
alternative: stationary 
 
Type 1: no drift no trend 
     lag    ADF p.value
[1,]   0 -16.36    0.01
[2,]   1 -11.79    0.01
[3,]   2  -9.98    0.01
[4,]   3 -10.56    0.01
[5,]   4  -7.66    0.01
Type 2: with drift no trend 
     lag    ADF p.value
[1,]   0 -16.29    0.01
[2,]   1 -11.74    0.01
[3,]   2  -9.94    0.01
[4,]   3 -10.52    0.01
[5,]   4  -7.63    0.01
Type 3: with drift and trend 
     lag    ADF p.value
[1,]   0 -16.23    0.01
[2,]   1 -11.70    0.01
[3,]   2  -9.91    0.01
[4,]   3 -10.48    0.01
[5,]   4  -7.61    0.01
---- 
Note: in fact, p.value = 0.01 means p.value <= 0.01 

# d=1 and D=1
# p = 0 1 2 3 4 7 11
# P = 0 1
# q = 0 1 5 10
# Q = 0 1 

automodel = auto.arima(y = sqrt(wine),trace = TRUE, seasonal = TRUE)

 ARIMA(2,1,2) with drift         : Inf
 ARIMA(0,1,0) with drift         : 856.1598
 ARIMA(1,1,0) with drift         : 849.2294
 ARIMA(0,1,1) with drift         : 846.0082
 ARIMA(0,1,0)                    : 854.2631
 ARIMA(1,1,1) with drift         : Inf
 ARIMA(0,1,2) with drift         : Inf
 ARIMA(1,1,2) with drift         : Inf
 ARIMA(0,1,1)                    : 844.3247
 ARIMA(1,1,1)                    : 829.9837
 ARIMA(1,1,0)                    : 847.4173
 ARIMA(2,1,1)                    : 832.0936
 ARIMA(1,1,2)                    : 832.0946
 ARIMA(0,1,2)                    : 836.6607
 ARIMA(2,1,0)                    : 845.4098
 ARIMA(2,1,2)                    : 834.0203

 Best model: ARIMA(1,1,1)                    

fit1 = arima(wine, order = c(2,1,5) ,
            seasonal = list(order=c(0,1,1),period=12), fixed = c(NA,0,0,0,NA,0,0,NA))
test(fit1$residuals)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     33.31    0.0312 *
McLeod-Li Q                Q ~ chisq(20)     29.91    0.0713
Turning points T     (T-93.3)/5 ~ N(0,1)        91    0.6402
Diff signs S       (S-70.5)/3.5 ~ N(0,1)        70    0.8848
Rank P         (P-5005.5)/283.5 ~ N(0,1)      5101    0.7362

confint(fit1)

          2.5 %      97.5 %
ar1  -0.6237578 -0.31360292
ar2          NA          NA
ma1          NA          NA
ma2          NA          NA
ma3  -0.3440913  0.01895912
ma4          NA          NA
ma5          NA          NA
sma1 -0.7568829 -0.41608966

US employment data

library(fpp3)

us_empl = us_employment %>% 
  filter(year(Month)>2009, Title=="All Employees, Total Nonfarm") %>%
  select(-Series_ID)
 
us_empl %>% 
  autoplot(Employed)

auto.arima(us_empl$Employed)

Series: us_empl$Employed 
ARIMA(3,1,1) with drift 

Coefficients:
         ar1      ar2      ar3      ma1     drift
      0.5173  -0.2112  -0.3101  -0.9059  198.7101
s.e.  0.0944   0.0974   0.0926   0.0499    8.0855

sigma^2 = 715029:  log likelihood = -945.12
AIC=1902.24   AICc=1903.01   BIC=1918.76

empl_data = us_empl$Employed

acf(empl_data,lag.max = 100)

M = c("diff",1,"diff",12)
r = Resid(empl_data,M)
test(r)

Null hypothesis: Residuals are iid noise.
Test                        Distribution Statistic   p-value
Ljung-Box Q                Q ~ chisq(20)     36.26    0.0143 *
McLeod-Li Q                Q ~ chisq(20)     16.13    0.7087
Turning points T     (T-68)/4.3 ~ N(0,1)        73    0.2408
Diff signs S         (S-51.5)/3 ~ N(0,1)        54     0.398
Rank P             (P-2678)/178 ~ N(0,1)      2339    0.0569