Time Series and Forecasting
Building SARIMA Models
Maria Eduarda Silva
[Link]
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 1 / 38
1 General Procedure (Box-Jenkins approach)
2 Characterization of the time series
3 Identify unit / seasonal unit roots
4 Identify the dependence orders of the model
5 Diagnostics
Statistical Significance of the Model
Residual Analysis
Model Selection
6 Overfitting
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 2 / 38
General Procedure (Box-Jenkins approach)
Basic steps to fitting SARIMA models to time series data
Plot and identify important characteristics of the data
Consider transforming the data if necessary: Box-Cox to stabilize
variance
If the data is non stationary: take first and/or seasonal differences
until data appears stationary
Examine the ACF/PACF to check stationarity and model order
Parameter estimation
Diagnostics
I Adequacy of the model- analysis of the residuals
I Statistical significance of the model
Use AICC to choose among models
Compute forecasts
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 3 / 38
Box-Jenkins approach
Step 1 Identify the model
Step 2 Estimate the model
Step 3 Diagnostics check
Adequacy Are the residuals uncorrelated?
Statistical significance Are all the parameters statistically
significant?
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 4 / 38
Box-Jenkins approach
Step 1 Identify the model
Step 2 Estimate the model
Step 3 Diagnostics check
Adequacy Are the residuals uncorrelated?
Statistical significance Are all the parameters statistically
significant?
If both Adequacy and Statistical significance are true you found an
adequate model
If at least one of Adequacy or Statistical significance is false return to
Step1
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 4 / 38
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 5 / 38
Characterization of the time series
Plot the data in a chronogram.
Check for:
I discontinuities such as level changes
I unusual observations- outliers
I changes in variance
I seasonality
I trend
I cycles
Some deterministic componentes due to physical phenomena may be
present and may be removed by deterministic functions: yearly
periodicities
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 6 / 38
Exemple: a series with a discontinuity
220
200
180
160
140
120
100
80
60
0 50 100 150 200 250 300 350 400 450
Figure: Australian beer production Jan 1956 - Abril de 1990.
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 7 / 38
Exemple: outlier
(a) (b)
3000 3000
2500 2500
2000 2000
1500 1500
1000 1000
500 500
0 0
0 50 100 150 0 50 100 150
Figure: Monthly sales (kl) of red wine in Australia Jan 1980 - Out 1991: the
outliers was an input error at t = 75 (a) original series (b).
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 8 / 38
Exemple: trend, seasonality and heterocedasticity
700
600
500
400
300
200
100
0 50 100 150
Figure: Number of airline passengers (×103 ) Jan 49 -Dec 60.
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 9 / 38
Box-Cox transforms
Stabilize the variance
Xtλ −1
(
λ se λ 6= 0
Ut =
log Xt se λ = 0
Choose λ that minimizes variance of data
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 10 / 38
Example
6.5
5.5
4.5
0 50 100 150
Figure: Log airline passengers
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 11 / 38
Transformations to stabilize the variance (Hyndman,
[Link], page 101 )
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 12 / 38
Back-transformations (Hyndman,
[Link], page 104)
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 13 / 38
Other transforms
Length of the month: since the different months of the year have
different number of days and also because of leap year, one may
adjust to the length of the month as follows:
365.25/12
Wt = Xt ×
no days in month t
Number of working days: after adjusting for the length of the month
mean number of working in a month
Wt = Xt ×
number of working days in month t
Adjust for moving holidays and interventions in general.
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 14 / 38
Exemple: monthly milk production per cow
(a) (b)
1000 1000
950 950
900 900
850 850
800 800
Libras
Libras
750 750
700 700
650 650
600 600
550 550
0 20 40 60 80 100 120 140 160 180 0 20 40 60 80 100 120 140 160
Figure: monthly milk production per cow (a) adjusted for length of month (b).
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 15 / 38
Statistical tests to determine the required order of
differencing
Augmented Dickey Fuller test: null hypothesis is that the data are
non-stationary and non-seasonal.
Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test: null hypothesis is
that the data are stationary and non-seasonal.
Other tests available for seasonal data.
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 16 / 38
Example:Federal Reserve Board Production Index
data
library(tseries)
[Link](prodn,alternative="s")
Augmented Dickey-Fuller Test
data: prodn
Dickey-Fuller = -2.9333, Lag order = 7, p-value = 0.183
alternative hypothesis: stationary
p-value> 0.05 indicates the need for first difference
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 17 / 38
KPSS: Test for unit root
Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test: reverses the
hypotheses, so the null-hypothesis is that the data are stationary in
level or trend
In this case, small p-values (e.g., less than 0.05) suggest that
differencing is required.
[Link](prodn)
KPSS Test for Level Stationarity
data: prodn
KPSS Level = 7.3711, Truncation lag parameter = 4, p-value = 0.01
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 18 / 38
ndiffs, nsdiffs from package forecast
Determine the lowest number of non-seasonal and seasonal necessary for
the series to become stationary
library(forecast)
ndiffs(WWWusage)
[1] 1
nsdiffs(log(AirPassengers))
[1] 1
ndiffs(diff(log(AirPassengers),12))
[1] 1
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 19 / 38
Identify the dependence orders of the model
Compute and plot Sample ACF, SACF, e Sample PACF, SPACF
Try to identify tentative orders for AR and/or MA components
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 20 / 38
Example: identify orders of the model: recruitment
data
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 21 / 38
library(astsa)
str(rec)
plot(rec)
acf2(rec)
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 22 / 38
The parameters must be significantly different from zero: at a 5% level
parameter θ estimated by θ̂ with standard error se is significantly different
from zero if 0 ∈
/ θ̂ ± 2 se.
Coefficients: ar1 ar2 intercept 1.3512 -0.4612 61.8585 s.e. 0.0416 0.0417
4.0039
both parameters and mean are statistically different from 0
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 23 / 38
Testing the residuals
The residuals must be UNCORRELATED
Bartlett test: if the residuals are approximately iid then the sample
acf of the residuals is N(0, 1/n).
Ljung-Box test: under the hypothesis of iid residuals
QLB = n(n + 2) hj=1 ρ̂2 (j)/(n − j) ∼ χ2h and graph the p-values
P
Normal probability plot to check for departures from Gaussianity
The sarima function produces the necessary plots
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 24 / 38
rec data
sarima(rec,2,0,0)
Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
Q), period = S), xreg = xmean, [Link] = FALSE, [Link] = li
REPORT = 1, reltol = tol))
Coefficients:
ar1 ar2 xmean
1.3512 -0.4612 61.8585
s.e. 0.0416 0.0417 4.0039
sigma^2 estimated as 89.33: log likelihood = -1661.51, aic = 3331.02
$AIC#$
[1] 5.505631
$AICc#$
[1] 5.510243
$BIC#$
[1] 4.532889
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 25 / 38
Checking the residuals for rec data
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 26 / 38
rec data
sarima(rec,3,0,0)
Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
Q), period = S), xreg = xmean, [Link] = FALSE, [Link] = li
REPORT = 1, reltol = tol))
Coefficients:
ar1 ar2 ar3 xmean
1.3318 -0.4043 -0.0421 61.9256
s.e. 0.0469 0.0759 0.0469 3.8411
sigma^2 estimated as 89.17: log likelihood = -1661.11, aic = 3332.22
The coefficiente for the AR(3) is not significant and the residuals checks
worsened.
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 27 / 38
Checking the residuals for rec data
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 28 / 38
rec data
sarima(rec,p=2,d=0,q=0,P=2,S=12)
Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
Q), period = S), xreg = xmean, [Link] = FALSE, [Link] = li
REPORT = 1, reltol = tol))
Coefficients:
ar1 ar2 sar1 sar2 xmean
1.3256 -0.4217 0.1116 0.1641 61.6089
s.e. 0.0431 0.0437 0.0460 0.0481 6.0800
sigma^2 estimated as 85.54: log likelihood = -1652.14, aic = 3316.28
$AIC#$
[1] 5.471106
$AICc#$
[1] 5.475937
$BIC#$
[1] 4.516536
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 29 / 38
Checking the residuals for rec data
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 30 / 38
Information Criteria
Akaike (1969, 1973, 1974) suggested measuring the goodness of a model
by balancing the error of the fit against the number of parameters in the
model. Thus Akaike Information Criteria was born. Later developed into
AICc and BIC (Bayesian Information Criteria)
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 31 / 38
Model selection for rec data
2 adequate models for rec data. Choose the model with minimum AIC
(BIC): SARIMA(2,0,0,2,0,0)12
Model
Coefficients AR(2) SARIMA(2,0,0,2,0,0)12
Mean 61.86 61.61
(4.0) (6.08)
AR1 1.354 1.33
(0.040) (0.040)
AR2 -0.46 -0.42
(0.040) (0.040)
SAR1 0.11
(0.05)
SAR2 0.16
(0.05)
AIC 5.50 5.47
BIC 4.53 4.51
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 32 / 38
Forecasting for rec data
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 33 / 38
Forecasting for rec data
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 34 / 38
Forecasting for milk data
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 35 / 38
Overfitting
Be aware of overfitting
More is not always synonym of better
Overfitting leads to less precise estimators
Adding more parameters may fit the data better but may also lead to
bad forecasts
Example: The fit dor the U.S. population by official census from 1910
to 1990 is perfect but the forecasts are terrible: negative population
sometime in 2002! The fit is obtained from polinomial of degree 8!
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 36 / 38
Example:US population
xpop=USPop$population[13:21]
tt=USPop$year[13:21]
tt1=tt-mean(tt)
[Link]=[Link]([Link],new)
plot(USPop$year[13:22],USPop$population[13:22],pch=19,xlim=c(1900,2020),yli
lines(tt,[Link]$[Link])
lines(c(1990,seq(1991,2010,1)),c([Link]$[Link][21],[Link]))
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 37 / 38
Hands on
Now try with the following data sets:
Quarterly U.S. GNP, gnp from the astsa package
AirPassengers
Maria Eduarda Silva ([Link]) Time Series and Forecasting Building SARIMA Models 38 / 38