Decomposition forecasting

The decomposition_model function allows for a separate approach to the seasonal and seasonally adjusted components. The latter is often combined with a model covered later including Holt’s method or ARIMA models.

Two Examples

us_retail_employment <- us_employment %>%
    filter(year(Month) >= 1990, Title == "Retail Trade")
fit_dcmp <- us_retail_employment %>%
    model(stlf = decomposition_model(STL(Employed ~ trend(window = 7), robust = TRUE),
        NAIVE(season_adjust)))
fit_dcmp %>%
    forecast() %>%
    autoplot(us_retail_employment) + labs(y = "Number of people", title = "US retail employment")

fit_dcmp %>%
    gg_tsresiduals()

fit_dcmp <- us_retail_employment %>%
    model(stlf = decomposition_model(STL(Employed ~ trend(window = 7), robust = TRUE),
        RW(season_adjust)))
fit_dcmp %>%
    forecast() %>%
    autoplot(us_retail_employment) + labs(y = "Number of people", title = "US retail employment")

fit_dcmp %>%
    gg_tsresiduals()

Evaluating Forecast Accuracy of Points

• On training and test [usually 20%]
  • filter and slice
• Scale dependent errors
  • MAE: mean absolute error
  • RMSE: root mean squared error
• Percentage errors
  • MAPE: mean absolute percentage error
  • sMAPE: symmetric MAPE
• Scaled errors [def in 5.8]
  • MASE: mean absolute scaled error
  • RMSSE: root mean square scaled error

accuracy

Evaluating Accuracy of Distributions

accuracy(., list(...))

• Quantile scores quantile_score
• Winkler scores winkler_score
• Continuous ranked probability scores CRPS
• Skill scores skill_score() with a statistic for comparison.

Key Principles

• Set the forecasting task clearly and concisely
• Implement a systematic approach
• Document and justify
• Systematically evaluate forecasts
• Segregate forecasters and users

Delphi

Delphi

Analogies

Structured Analogies

Adjustments

• Use them sparingly.
• Make them structured.

A Simplified Example

load("USEmployment.RData")
us_employment %>%
    data.frame() %>%
    group_by(Series_ID) %>%
    summarise(Title = first(Title)) %>%
    mutate(series_id = Series_ID) %>%
    ungroup() %>%
    select(-Series_ID) -> Names.List
US.Employment.T <- left_join(US.Employment, Names.List, by = c(series_id = "series_id")) %>%
    mutate(YM = yearmonth(date)) %>%
    rename(Employed = value) %>%
    as_tsibble(., index = YM, key = Title)
USET <- US.Employment.T %>%
    filter(YM > yearmonth("1990-01"), Title %in% c("Retail Trade")) %>%
    as_tsibble(., index = YM, key = Title)

A Seasonal Linear Model

y_{t} = Month_{t} + \epsilon The fitted/predicted value is \hat{y} = \hat{\alpha} + \hat{\beta}_{M}*Month_{t} and the residual is:

\hat{e} = y - \hat{y} or expanding and distributing: \hat{e} = y - \hat{\alpha} - \hat{\beta}_{M}*Month_{t}

The residual

USRTM %>%
    gg_tsresiduals()

Residuals

fit_USET %>%
    gg_tsresiduals()

How Does it Look?

augment(fit_USET) %>%
    select(YM, Employed, .fitted) %>%
    pivot_longer(c(Employed, .fitted)) %>%
    autoplot()

Knots

fit_USET2 <- USET %>%
    mutate(Month = as.character(month(YM, label = TRUE))) %>%
    model(TSLM(Employed ~ trend(knots = yearmonth(c("2001-09", "2008-10", "2020-03"))) +
        Month))

This Is Cheating

fit_USET2 <- USET %>%
    mutate(Year = as.character(year(YM)), Month = as.character(month(YM, label = TRUE))) %>%
    model(TSLM(Employed ~ Year + Month))

Evaluation

1. fit %>% gg_tsresiduals() and augment(fit) %>% features(.innov, ljung_box, lag = ??, dof = ??) for time series information in residuals.
   
2. plot of .innov or .resid against the predictors. This many require joining the residuals to the original tsibble, e.g. data %>% left_join(residuals(fit), by="Index").
   
3. plot of residuals [.resid] and fitted.values [.fitted].
   
4. outliers and influence are always worthwhile to consider. Regular regression diagnostics are relevant here. It is just a tibble with a index: timetk.

Handy predictors

1. trend()
2. season() the dummy/binary variable month(label=TRUE), weekday and the trap
3. interventions, trading days, holidays, and distributed lags.
4. the Fourier series.

Forecasting and Scenarios

Take the model of quarterly consumption as a function of income.

fit_cons <- us_change %>%
    model(TSLM(Consumption ~ Income))
new_cons <- scenarios(`Average increase` = new_data(us_change, 4) %>%
    mutate(Income = mean(us_change$Income)), `Extreme increase` = new_data(us_change,
    4) %>%
    mutate(Income = c(mean(us_change$Income), mean(us_change$Income) + 3, mean(us_change$Income) +
        6, mean(us_change$Income) + 9)), names_to = "Scenario")
fcast <- forecast(fit_cons, new_cons)
us_change %>%
    autoplot(Consumption) + autolayer(fcast) + labs(y = "% change in US consumption")

Knots and their Importance

This section points out splines and ties.

Correlation and Causation

ARE NOT THE SAME.

A Bit More on Linear Models

Nonlinearity; they argue against nonlinear trends. So now what? Knots. Trends are not always monotonic. Witness the example of population in Afghanistan. NB: But this probably gets unrealistic pretty fast.

global_economy %>%
    filter(Country == "Afghanistan") %>%
    autoplot(Population/1e+06) + labs(y = "Population (millions)") + geom_ribbon(aes(xmin = 1979.98,
    xmax = 1989.13), fill = "pink", alpha = 0.4) + annotate("text", x = 1984.5,
    y = 10, label = "Soviet-Afghan war", col = "red", size = 3)

Compare

fit <- global_economy %>%
    filter(Country == "Afghanistan") %>%
    model(linear = TSLM(Population ~ trend()), piecewise = TSLM(Population ~
        trend(knots = c(1980, 1989))))
fit %>%
    report()

Compare

# A tibble: 2 × 16
  Country     .model r_squ…¹ adj_r…²  sigma2 stati…³  p_value    df log_lik
  <fct>       <chr>    <dbl>   <dbl>   <dbl>   <dbl>    <dbl> <int>   <dbl>
1 Afghanistan linear   0.838   0.835 1.02e13    290. 8.37e-24     2   -950.
2 Afghanistan piece…   0.999   0.999 9.05e10  12929. 4.34e-77     4   -812.
# … with 7 more variables: AIC <dbl>, AICc <dbl>, BIC <dbl>, CV <dbl>,
#   deviance <dbl>, df.residual <int>, rank <int>, and abbreviated
#   variable names ¹​r_squared, ²​adj_r_squared, ³​statistic

Handling Multiple Models

mable stores models. The rows are given by the key and the columns are given by the model names. For diagnostics on a particular model/subset, we will need a filter/select combination.

global_economy %>%
  model(
    linear = TSLM(Population ~ trend()),
    piecewise = TSLM(Population ~ trend(knots = c(1980, 1989)))
  ) %>%
  filter(Country=="Afghanistan") %>%
  select(piecewise) %>%
  gg_tsresiduals()

Residual Assessment

fit1 <- global_economy %>%
    filter(Country == "Afghanistan") %>%
    model(linear = TSLM(Population ~ trend()))
gg_tsresiduals(fit1)

Residual Assessment

Residual Assessment

fit2 <- global_economy %>%
    filter(Country == "Afghanistan") %>%
    model(piecewise = TSLM(Population ~ trend(knots = c(1980, 1989))))
gg_tsresiduals(fit2)

Residual Assessment

Fit Assessment

Forecasts?

fc <- fit %>%
    forecast(h = "5 years")
autoplot(fc) + autolayer(filter(global_economy %>%
    filter(Country == "Afghanistan")), Population)

Simple Exponential Smoothing

Appropriate for absence of seasons and trends. We could use naive (most recent observation receives all the weight, \alpha=1) and average (all values equally weighted, \alpha=0) methods for such data. Simple exponential smoothing occupies the space between with parameter \alpha such that

\hat{y}_{T+1|T} = \alpha y_{T} + \alpha(1 - \alpha)y_{T-1} + \alpha(1-\alpha)^{2}y_{T-2} + \ldots + \alpha(1-\alpha)^{k}y_{T-k}

Weighted Averages and the Component Form

Describe the outcome as two separate but combineable forms: a Forecast equation \hat{y}_{t+h|t} = S_{t} and a Smoothing Equation S_{t} = \alpha y_{t} + (1-\alpha) S_{t-1}.

If we combine them, we get:

\hat{y}_{t+h|t} = \alpha y_{t} + (1-\alpha) S_{t-1}

The Model: ETS

the error component (either M multiplicative or A additive)

the trend

the season

Some comments on options – the optimization method.

Damping

Flatten the trend into the future. Expand Holt’s trend model with an additional damping parameter – \phi.

\hat{y}_{t+h|t} = S_{t} + (\phi + \phi^2 + \ldots + \phi^{h})b_{t}

S_{t} = \alpha y_{t} + (1-\alpha)(S_{t-1} + \phi b_{t-1})

b_{t} = \beta^{*}(S_{t} - S_{t-1}) + (1-\beta^{*})\phi b_{t-1}

Damping

Add the damping terms in the trend to the previous:

HW-Damp

A Taxonomy [and the code]

A glorious taxonomy for the models that all fits in one simple code expression.

model(ETS(outcome ~ error("A"/"M") 
+ trend(method=c("N","A","Ad")) 
+ season(method = c("N", "A", "M"))))

Model Selection and IC [and Caveats]

NB: Multiplicative errors and negative values [also, next slide]

Information criterion [IC] provide a measure of goodness of fit for models that may differ in number of parameters or underlying form. Minimize it.

ETS will do this automagically if not given guidance.

aus_busines <- tourism %>%
    filter(Purpose == "Business") %>%
    summarise(Trips = sum(Trips)/1000)
fit <- aus_busines %>%
    model(ETS(Trips))
report(fit)

Result

Series: Trips 
Model: ETS(M,N,A) 
  Smoothing parameters:
    alpha = 0.473 
    gamma = 1e-04 

  Initial states:
 l[0]   s[0] s[-1]  s[-2]  s[-3]
 3.96 0.0727  0.39 0.0983 -0.561

  sigma^2:  0.0057

 AIC AICc  BIC 
 164  165  180 

What does it look like?

components(fit) %>%
    autoplot() + labs(title = "ETS(M,N,A) components")

Innovations and Residuals [Multiplicative]

fit %>%
    augment() %>%
    select(.innov, .resid) %>%
    pivot_longer(c(.innov, .resid)) %>%
    autoplot()

Forecast Variances

Not all are known analytically, but we can always use distributions of step ahead paths to simulate and summarise them.

A Simplified Example

load(url("https://github.com/robertwwalker/rww-science/raw/master/static/xaringan/CH7HA/USEmployment.RData"))
us_employment %>%
    data.frame() %>%
    group_by(Series_ID) %>%
    summarise(Title = first(Title)) %>%
    mutate(series_id = Series_ID) %>%
    ungroup() %>%
    select(-Series_ID) -> Names.List
US.Employment.T <- left_join(US.Employment, Names.List, by = c(series_id = "series_id")) %>%
    mutate(YM = yearmonth(date)) %>%
    rename(Employed = value) %>%
    as_tsibble(., index = YM, key = Title)
USET <- US.Employment.T %>%
    filter(YM > yearmonth("1990-01"), Title %in% c("Retail Trade")) %>%
    as_tsibble(., index = YM, key = Title)
USEM <- US.Employment.T %>%
    filter(YM > yearmonth("1990-01"), Title %in% c("Manufacturing")) %>%
    as_tsibble(., index = YM, key = Title)

A Model

USETS <- USET %>%
    model(ETS(Employed))
USETS %>%
    report()

Series: Employed 
Model: ETS(M,N,A) 
  Smoothing parameters:
    alpha = 0.967 
    gamma = 1e-04 

  Initial states:
  l[0] s[0] s[-1] s[-2] s[-3] s[-4] s[-5] s[-6] s[-7] s[-8] s[-9] s[-10]
 13428  584   407  41.5 -89.7 -7.18 -12.9 -19.3  -128  -225  -231   -253
 s[-11]
  -65.9

  sigma^2:  1e-04

 AIC AICc  BIC 
5872 5874 5931 

The residual

USETS %>%
    gg_tsresiduals()

They don’t always end up MNA

USETSM <- USEM %>%
    model(ETS(Employed))
USETSM %>%
    report()

Series: Employed 
Model: ETS(A,Ad,A) 
  Smoothing parameters:
    alpha = 0.836 
    beta  = 0.128 
    gamma = 0.164 
    phi   = 0.929 

  Initial states:
  l[0]  b[0] s[0] s[-1] s[-2] s[-3] s[-4] s[-5] s[-6] s[-7] s[-8] s[-9]
 17951 -41.1 24.7  33.6  50.4  80.6    99  17.2  81.3 -37.8 -85.8   -66
 s[-10] s[-11]
  -89.2   -108

  sigma^2:  7864

 AIC AICc  BIC 
5573 5575 5644 

The flow:

1. Stationarity and Differencing
   
2. Non-seasonal ARIMA models
   
3. Estimation and order selection
   
4. ARIMA in R
   
5. Forecasting ARIMA models
   
6. Seasonal ARIMA models
   
7. Comparison: ETS vs. ARIMA

Stationarity and Differencing

Definition:

If Y_{t} is a stationary time series, then for all s, the distribution y_{t}, ..., y_{t+s} doesn’t depend on t.

The implies that - the mean (roughly horizontal), - the variance (constant variance), - and the covariance (no predictable long-term patterns) between elements of Y are constants.

Stationary

library(hrbrthemes)
gafa_stock %>%
    filter(Symbol == "GOOG", year(Date) == 2018) %>%
    autoplot(Close) + ylab("Google closing stock price") + xlab("Day") + theme_ipsum_rc()

Stationary?

gafa_stock %>%
    filter(Symbol == "GOOG", year(Date) == 2018) %>%
    autoplot(difference(Close)) + ylab("Google closing stock price") + xlab("Day") +
    theme_ipsum_rc()

Stationary?

as_tsibble(fma::strikes) %>%
    autoplot(value) + ylab("Number of strikes") + xlab("Year") + theme_ipsum_rc()

Stationary?

as_tsibble(fma::hsales) %>%
    autoplot(value) + labs(x = "Year", y = "Total sales", title = "Sales of new one-family houses, USA") +
    theme_ipsum_rc()

Stationary?

as_tsibble(fma::eggs) %>%
    autoplot(value) + labs(x = "Year", y = "$", title = "Price of a dozen eggs in 1993 dollars") +
    theme_ipsum_rc()

Stationary?

aus_livestock %>%
    filter(Animal == "Pigs", State == "Victoria", year(Month) >= 2010) %>%
    autoplot(Count/1000) + labs(x = "Year", y = "thousands", title = "Number of pigs slaughtered in Victoria") +
    theme_ipsum_rc()

Stationary?

pelt %>%
    autoplot(Lynx) + labs(x = "Year", y = "Number trapped", title = "Annual Canadian Lynx Trappings") +
    theme_ipsum_rc()

Definition

If \{y_t\} is a stationary time series, then for all s, the distribution of (y_t,\dots,y_{t+s}) does not depend on t.

Transformations help to stabilize the variance.

For ARIMA modelling, we also need to stabilize the mean.

Non-stationarity in the mean

Identifying non-stationary series

• time plot.
• The ACF of stationary data drops to zero relatively quickly
• The ACF of non-stationary data decreases slowly.
• For non-stationary data, the value of r_1 is often large and positive.

Google Price

google_2018 <- gafa_stock %>%
    filter(Symbol == "GOOG", year(Date) == 2018) %>%
    mutate(trading_day = row_number()) %>%
    update_tsibble(index = trading_day, regular = TRUE)

google_2018 %>%
    autoplot(Close) + theme_ipsum_rc()

google_2018 %>%
    ACF(Close) %>%
    autoplot() + theme_ipsum_rc()

Google change

gafa_stock %>%
    filter(Symbol == "GOOG", year(Date) == 2018) %>%
    autoplot(difference(Close)) + labs(y = "Change in Google closing stock price ( $USD)",
    x = "Day") + theme_ipsum_rc()

google_2018 %>%
    ACF(difference(Close)) %>%
    autoplot() + theme_ipsum_rc()

Conclusion

The price is unlikely stationary, the change appears to be.

Differencing

Definition:
- The differenced series is the change between each observation in the original series: y'_t = y_t - y_{t-1}.

• Differencing helps to stabilize the mean.

• The differenced series will have only T-1 values since it is not possible to calculate a difference y_1' for the first observation.

Sometimes one isn’t enough

Instead of the change, we need model the change in change – a second difference.

y^{''}_{t} = y^{'}_{t} - y^{'}_{t - 1}

y^{''}_{t} = (y_t - y_{t-1}) - (y_{t-1}-y_{t-2})

y^{''}_{t} = y_t - 2y_{t-1} +y_{t-2}

• y_t'' will have T-2 values.
• In practice, it is almost never necessary to go beyond second-order differences.

Seasonal differencing

A seasonal difference is the difference between an observation and the corresponding observation from the previous year. y^{'seasonal}_t = y_t - y_{t-m} where m= number of seasons.

• For monthly data m=12.
• For quarterly data m=4.

Electricity Generation

The oscillation appears more variable in the latter half than the former half.

usmelec <- fpp2::usmelec %>%
    as_tsibble() %>%
    mutate(Generation = value)
usmelec %>%
    autoplot() + theme_ipsum_rc()

Electricity Generation (log)

usmelec %>%
    autoplot(log(Generation)) + theme_ipsum_rc()

Electricity Generation (seasonal difference and log)

usmelec %>%
    autoplot(log(Generation) %>%
        difference(12)) + theme_ipsum_rc()

Electricity Generation (seasonal difference and log)

Patterns seem to remain.

usmelec %>%
    ACF(log(Generation) %>%
        difference(12)) %>%
    autoplot() + theme_ipsum_rc()

Electricity Generation (seasonal difference and log with difference)

usmelec %>%
    autoplot(log(Generation) %>%
        difference(12) %>%
        difference()) + theme_ipsum_rc()

Electricity Generation (seasonal difference and log with difference)

usmelec %>%
    ACF(log(Generation) %>%
        difference(12) %>%
        difference()) %>%
    autoplot() + theme_ipsum_rc()

Electricity production

• Seasonally differenced series is closer to being stationary.
• Remaining non-stationarity can be removed with further first difference.

If y'_t = y_t - y_{t-12} denotes seasonally differenced series, then twice-differenced series is

y^*_t = y'_t - y'_{t-1}

y^*_t = (y_t - y_{t-12}) - (y_{t-1} - y_{t-13})

y^*_t = y_t - y_{t-1} - y_{t-12} + y_{t-13}

Seasonal differencing

When both seasonal and first differences are applied \dots

• it makes no difference which is done first – the result will be the same.
• If seasonality is strong, we recommend that seasonal differencing be done first because sometimes the resulting series will be stationary and there will be no need for further first difference.

It is important that if differencing is used, the differences are interpretable.

Interpretation of differencing

• first differences are the change between one observation and the next;
• seasonal differences are the change between one year to the next.

But taking lag 3 differences for yearly data, for example, results in a model which cannot be sensibly interpreted.

Unit root tests

Statistical tests can suggest the required order of differencing.

1. Augmented Dickey Fuller test: null hypothesis is that the data are non-stationary and non-seasonal.
2. Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test: null hypothesis is that the data are stationary and non-seasonal.
3. Other tests available for seasonal data.

You May Recall from Time Series features

KPSS test

Google prices won’t work; differences are indicated.

google_2018 %>%
    features(Close, features = list(unitroot_kpss, unitroot_ndiffs))

# A tibble: 1 × 4
  Symbol kpss_stat kpss_pvalue ndiffs
  <chr>      <dbl>       <dbl>  <int>
1 GOOG       0.573      0.0252      1

Automatically selecting differences

STL decomposition: y_t = T_t+S_t+R_t

Seasonal strength F_s = \max\big(0, 1-\frac{\text{Var}(R_t)}{\text{Var}(S_t+R_t)}\big)

If F_s > 0.64, do one seasonal difference.

usmelec %>%
    mutate(log_gen = log(Generation)) %>%
    features(log_gen, list(unitroot_nsdiffs, feat_stl))

# A tibble: 1 × 10
  nsdiffs trend_…¹ seaso…² seaso…³ seaso…⁴ spikin…⁵ linea…⁶ curva…⁷ stl_e…⁸
    <int>    <dbl>   <dbl>   <dbl>   <dbl>    <dbl>   <dbl>   <dbl>   <dbl>
1       1    0.994   0.943       7       4 1.40e-12    5.65  -0.911   0.207
# … with 1 more variable: stl_e_acf10 <dbl>, and abbreviated variable
#   names ¹​trend_strength, ²​seasonal_strength_year, ³​seasonal_peak_year,
#   ⁴​seasonal_trough_year, ⁵​spikiness, ⁶​linearity, ⁷​curvature, ⁸​stl_e_acf1

Automatically selecting differences

usmelec %>%
    mutate(log_gen = log(Generation)) %>%
    features(log_gen, unitroot_nsdiffs)

# A tibble: 1 × 1
  nsdiffs
    <int>
1       1

usmelec %>%
    mutate(d_log_gen = difference(log(Generation), 12)) %>%
    features(d_log_gen, unitroot_ndiffs)

# A tibble: 1 × 1
  ndiffs
   <int>
1      1

Try it out with tourism

For the tourism dataset, compute the total number of trips and find an appropriate differencing (after transformation if necessary) to obtain stationary data.

tourism %>%
    summarise(Trips = sum(Trips)) %>%
    autoplot() + labs(title = "Total trips in tourism", x = "Quarter", y = "Trips") +
    theme_ipsum_rc()

The Backshift Operator

A very useful notational device is the backward shift operator, B, which is used as follows: B y_{t} = y_{t - 1}

In other words, B, operating on y_{t}, has the effect of shifting the data back one period. Two applications of B to y_{t} shifts the data back two periods: B(By_{t}) = B^{2}y_{t} = y_{t-2}

For monthly data, if we wish to shift attention to “the same month last year”, then B^{12} is used, and the notation is B^{12}y_{t} = y_{t-12}.

Backshift notation

The backward shift operator is convenient for describing the process of differencing.

A first difference can be written as y'_{t} = y_{t} - y_{t-1} = y_t - By_{t} = (1 - B)y_{t}

Note that a first difference is represented by (1 - B).

Similarly, if second-order differences (i.e., first differences of first differences) have to be computed, then:

y''_{t} = y_{t} - 2y_{t - 1} + y_{t - 2} = (1 - B)^{2} y_{t}

Backshift notation

• Second-order difference is denoted (1- B)^{2}.
• Second-order difference is not the same as a second difference, which would be denoted 1- B^{2};
• In general, a dth-order difference can be written as (1 - B)^{d} y_{t}
  
• A seasonal difference followed by a first difference can be written as (1-B)(1-B^m)y_t

Backshift notation

The “backshift” notation is convenient because the terms can be multiplied together to see the combined effect.

(1-B)(1-B^m)y_t = (1 - B - B^m + B^{m+1})y_t

(1-B)(1-B^m)y_t = y_t-y_{t-1}-y_{t-m}+y_{t-m-1}

For monthly data, m=12 and we obtain the same result as earlier.

Autoregressive models

Autoregressive (AR) models:

y_{t} = c + \phi_{1}y_{t - 1} + \phi_{2}y_{t - 2} + \cdots + \phi_{p}y_{t - p} + \varepsilon_{t}

where \varepsilon_t is white noise. This is a multiple regression with lagged values of y_t as predictors.

AR(1) model

y_{t} = 18 -0.8 y_{t - 1} + \varepsilon_{t}

\varepsilon_t\sim N(0,1), T=100.

AR(1) model

y_{t} = c + \phi_1 y_{t - 1} + \varepsilon_{t}

• When \phi_1=0, y_t is equivalent to WN
• When \phi_1=1 and c=0, y_t is equivalent to a RW
• When \phi_1=1 and c\ne0, y_t is equivalent to a RW with drift
• When \phi_1<0, y_t tends to oscillate between positive and negative values.

AR(2) model

y_t = 8 + 1.3y_{t-1} - 0.7 y_{t-2} + \varepsilon_t

\varepsilon_t\sim N(0,1), \qquad, T=100.

p2

Stationarity conditions

We normally restrict autoregressive models to stationary data, and then some constraints on the values of the parameters are required.

General condition for stationarity

Complex roots of 1-\phi_1 z - \phi_2 z^2 - \dots - \phi_pz^p lie outside the unit circle on the complex plane.

• For p=1: -1<\phi_1<1.
• For p=2:-1<\phi_2<1\qquad \phi_2+\phi_1 < 1 \qquad \phi_2 -\phi_1 < 1.
• More complicated conditions hold for p\ge3.
• Estimation software takes care of this.

Moving Average (MA) models

Moving Average (MA) models: y_{t} = c + \varepsilon_t + \theta_{1}\varepsilon_{t - 1} + \theta_{2}\varepsilon_{t - 2} + \cdots + \theta_{q}\varepsilon_{t - q}, where \varepsilon_t is white noise.

This is a multiple regression with past errors as predictors. Don’t confuse this with moving average smoothing!

MA(1) model

y_t = 20 + \varepsilon_t + 0.8 \varepsilon_{t-1}

\varepsilon_t\sim N(0,1) and T=100.

p1

MA(2) model

y_t = \varepsilon_t -\varepsilon_{t-1} + 0.8 \varepsilon_{t-2}

\varepsilon_t\sim N(0,1) and T=100.

p2

Invertibility

• Any MA( q ) process can be written as an AR( \infty ) process if we impose some constraints on the MA parameters.
• Then the MA model is called invertible.
• Invertible models have some mathematical properties that make them easier to use in practice.
• Invertibility of an ARIMA model is equivalent to forecastability of an ETS model.

Invertibility

General condition for invertibility: Complex roots of 1+\theta_1 z + \theta_2 z^2 + \dots + \theta_qz^q lie outside the unit circle on the complex plane.

• For q=1: -1<\theta_1<1.
• For q=2: -1<\theta_2<1 \qquad \theta_2+\theta_1 >-1 \qquad \theta_1 -\theta_2 < 1
• More complicated conditions hold for q\ge3.
• Estimation software takes care of this.

ARMA models

Autoregressive Moving Average models:

y_{t} = c + \phi_{1}y_{t - 1} + \cdots + \phi_{p}y_{t - p} + \theta_{1}\varepsilon_{t - 1} + \cdots + \theta_{q}\varepsilon_{t - q} + \varepsilon_{t}.

• Predictors include both lagged values of y_t and lagged errors.
• Conditions on coefficients ensure stationarity.
• Conditions on coefficients ensure invertibility.

Autoregressive Integrated Moving Average models

• Combine ARMA model with differencing.
• (1-B)^d y_t follows an ARMA model.

ARIMA models

Autoregressive Integrated Moving Average models

ARIMA( p, d, q ) model

AR: p = order of the autoregressive part
I: d = degree of first differencing involved
MA: q = order of the moving average part.

• White noise model: ARIMA(0,0,0)
  
• Random walk: ARIMA(0,1,0) with no constant
  
• Random walk with drift: ARIMA(0,1,0) with \rlap{const.}
  
• AR( p ): ARIMA( p,0,0)
  
• MA( q ): ARIMA(0,0,q )

Backshift notation for ARIMA

• ARMA model:

y_{t} = c + \phi_{1}By_{t} + \cdots + \phi_pB^py_{t} + \varepsilon_{t} + \theta_{1}B\varepsilon_{t} + \cdots + \theta_qB^q\varepsilon_{t}

or

y_{t} = (1-\phi_1B - \cdots - \phi_p B^p) y_t = c + (1 + \theta_1 B + \cdots + \theta_q B^q)\varepsilon_t

ARIMA(1,1,1) model

\underbrace{(1 - \phi_{1} B)}_{\text{AR(1)}} \underbrace{(1 - B)}_{\text{First Difference}} y_{t} = c + \underbrace{(1 + \theta_{1} B) \varepsilon_{t}}_{\text{MA(1)}}

Written out: y_t = c + y_{t-1} + \phi_1 y_{t-1}- \phi_1 y_{t-2} + \theta_1\varepsilon_{t-1} + \varepsilon_t

US consumption expenditure

us_change %>%
    autoplot(Consumption) + xlab("Year") + ylab("Quarterly percentage change") +
    ggtitle("US consumption") + theme_ipsum_rc()

US personal consumption

fit <- us_change %>%
    model(arima = ARIMA(Consumption ~ PDQ(0, 0, 0)))
report(fit)

Series: Consumption 
Model: ARIMA(1,0,3) w/ mean 

Coefficients:
        ar1     ma1     ma2     ma3  constant
      0.573  -0.362  0.0925  0.1934    0.3160
s.e.  0.150   0.161  0.0787  0.0824    0.0371

sigma^2 estimated as 0.3334:  log likelihood=-170
AIC=352   AICc=352   BIC=372

ARIMA(1,0,3) model:

y_t = 0.316 + 0.573y_{t-1} -0.362 \varepsilon_{t-1} + 0.0925 \varepsilon_{t-2} + 0.193 \varepsilon_{t-3} + \varepsilon_{t}

where \varepsilon_t is white noise with a standard deviation of 0.577 = \sqrt{0.333}.

US personal consumption

fit %>%
    forecast(h = 10) %>%
    autoplot(tail(us_change, 80)) + theme_ipsum_rc()

Understanding ARIMA models

• If c=0 and d=0, the long-term forecasts will go to zero.
• If c=0 and d=1, the long-term forecasts will go to a non-zero constant.
• If c=0 and d=2, the long-term forecasts will follow a straight line.
• If c\ne0 and d=0, the long-term forecasts will go to the mean of the data.
• If c\ne0 and d=1, the long-term forecasts will follow a straight line.
• If c\ne0 and d=2, the long-term forecasts will follow a quadratic trend.

Forecast variance and d

• The higher the value of d, the more rapidly the prediction intervals increase in size.
• For d=0, the long-term forecast standard deviation will go to the standard deviation of the historical data.

Cyclic behaviour

• For cyclic forecasts, p\ge2 and some restrictions on coefficients are required.
• If p=2, we need \phi_1^2+4\phi_2<0. Then average cycle of length

(2\pi)/\left[\text{arccos}(-\phi_1(1-\phi_2)/(4\phi_2))\right].

ML Estimation

Having identified the model order, we need to estimate the parameters c, \phi_1,\dots,\phi_p, \theta_1,\dots,\theta_q.

• MLE is very similar to least squares estimation obtained by minimizing \sum_{t-1}^T e_t^2
• The ARIMA() function allows CLS or MLE estimation.
• Non-linear optimization must be used in either case.
• Different software will give different estimates.

Partial autocorrelations

Partial autocorrelations measure relationship between y_{t} and y_{t - k}, when the effects of other time lags 1, 2, 3, \dots, k - 1 are removed.

\alpha_k = \text{kth partial autocorrelation coefficient}

\alpha_k = \text{equal to the estimate of } \phi_k \text{ in regression }

y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \dots + \phi_k y_{t-k}.

• Varying number of terms on RHS gives \alpha_k for different values of k.
• \alpha_1=\rho_1
  • same critical values of \pm 1.96/\sqrt{T} as for ACF.
• Last significant \alpha_k indicates the order of an AR model.

Mink trapping

mink <- as_tsibble(fma::mink)
mink %>%
    autoplot(value) + xlab("Year") + ylab("Minks trapped (thousands)") + ggtitle("Annual number of minks trapped") +
    theme_ipsum_rc()

Mink trapping

p1 <- mink %>%
    ACF(value) %>%
    autoplot()
p2 <- mink %>%
    PACF(value) %>%
    autoplot()
p1 | p2

Mink trapping

mink %>%
    gg_tsdisplay(value, plot_type = "partial")

ACF and PACF interpretation

AR(1)

\rho_k = \phi_1^k \qquad \text{ for } k=1,2,\dots

\alpha_1 = \phi_1 \qquad\alpha_k = 0\qquad\text{for} k=2,3,\dots

So we have an AR(1) model when

• autocorrelations exponentially decay
• there is a single significant partial autocorrelation.

ACF and PACF interpretation

AR(p)

• ACF dies out in an exponential or damped sine-wave manner
• PACF has all zero spikes beyond the pth spike

So we have an AR(p) model when

• the ACF is exponentially decaying or sinusoidal
• there is a significant spike at lag p in PACF, but none beyond p

ACF and PACF interpretation

MA(1)

\rho_1 = \theta_1\qquad \rho_k = 0 \qquad \text{ for } k=2,3,\dots

\alpha_k = -(-\theta_1)^k

So we have an MA(1) model when

• the PACF is exponentially decaying and
• there is a single significant spike in ACF

ACF and PACF interpretation

MA( q )

• PACF dies out in an exponential or damped sine-wave manner
• ACF has all zero spikes beyond the qth spike

So we have an MA( q ) model when

• the PACF is exponentially decaying or sinusoidal
• there is a significant spike at lag q in ACF, but none beyond q

Choosing your own model

web_usage <- as_tsibble(WWWusage)
web_usage %>%
    gg_tsdisplay(value, plot_type = "partial")

Choosing your own model

web_usage %>%
    mutate(diff = difference(value)) %>%
    gg_tsdisplay(diff, plot_type = "partial")

Choosing your own model

fit <- web_usage %>%
    model(arima = ARIMA(value ~ pdq(3, 1, 0)))
report(fit)

Series: value 
Model: ARIMA(3,1,0) 

Coefficients:
        ar1     ar2     ar3
      1.151  -0.661  0.3407
s.e.  0.095   0.135  0.0941

sigma^2 estimated as 9.656:  log likelihood=-252
AIC=512   AICc=512   BIC=522

Choosing your own model

web_usage %>%
    model(ARIMA(value ~ pdq(d = 1))) %>%
    report()

Series: value 
Model: ARIMA(1,1,1) 

Coefficients:
         ar1     ma1
      0.6504  0.5256
s.e.  0.0842  0.0896

sigma^2 estimated as 9.995:  log likelihood=-254
AIC=514   AICc=515   BIC=522

Choosing your own model

web_usage %>%
    model(ARIMA(value ~ pdq(d = 1), stepwise = FALSE, approximation = FALSE)) %>%
    report()

Series: value 
Model: ARIMA(3,1,0) 

Coefficients:
        ar1     ar2     ar3
      1.151  -0.661  0.3407
s.e.  0.095   0.135  0.0941

sigma^2 estimated as 9.656:  log likelihood=-252
AIC=512   AICc=512   BIC=522

Choosing your own model

gg_tsresiduals(fit)

Choosing your own model

augment(fit) %>%
    features(.resid, ljung_box, lag = 10, dof = 3)

# A tibble: 1 × 3
  .model lb_stat lb_pvalue
  <chr>    <dbl>     <dbl>
1 arima     4.49     0.722

Choosing your own model

fit %>%
    forecast(h = 10) %>%
    autoplot(web_usage)

Modelling procedure with ARIMA()

1. Plot the data. Identify any unusual observations.
   
2. If necessary, transform the data (using a Box-Cox transformation) to stabilize the variance.
   
3. If the data are non-stationary: take first differences of the data until the data are stationary. 4. Examine the ACF/PACF: Is an AR(p) or MA(q) model appropriate?
   
4. Try your chosen model(s), and use the AICc to search for a better model.
   
5. Check the residuals from your chosen model by plotting the ACF of the residuals, and doing a portmanteau test of the residuals. If they do not look like white noise, try a modified model.
   
6. Once the residuals look like white noise, calculate forecasts.

Automatic modelling procedure with ARIMA()

1. Plot the data. Identify any unusual observations.
   
2. If necessary, transform the data (using a Box-Cox transformation) to stabilize the variance.
   
3. Use ARIMA to automatically select a model.
   
4. Check the residuals from your chosen model by plotting the ACF of the residuals, and doing a portmanteau test of the residuals. If they do not look like white noise, try a modified model.
   
5. Once the residuals look like white noise, calculate forecasts.

Modelling procedure

link here

Seasonally adjusted electrical equipment

elecequip <- as_tsibble(fpp2::elecequip)
dcmp <- elecequip %>%
    model(STL(value ~ season(window = "periodic"))) %>%
    components() %>%
    select(-.model)
dcmp %>%
    as_tsibble %>%
    autoplot(season_adjust) + xlab("Year") + ylab("Seasonally adjusted new orders index")

Seasonally adjusted electrical equipment

1. Time plot shows sudden changes, particularly big drop in 2008/2009 due to global economic environment. Otherwise nothing unusual and no need for data adjustments.
2. No evidence of changing variance, so no Box-Cox transformation.
3. Data are clearly non-stationary, so we take first differences.

Seasonally adjusted electrical equipment

dcmp %>%
    mutate(diff = difference(season_adjust)) %>%
    gg_tsdisplay(diff, plot_type = "partial")

Seasonally adjusted electrical equipment

4. PACF is suggestive of AR(3). So initial candidate model is ARIMA(3,1,0). No other obvious candidates.
5. Fit ARIMA(3,1,0) model along with variations: ARIMA(4,1,0), ARIMA(2,1,0), ARIMA(3,1,1), etc. ARIMA(3,1,1) has smallest AICc value.

Seasonally adjusted electrical equipment

fit <- dcmp %>%
    model(arima = ARIMA(season_adjust))
report(fit)

Series: season_adjust 
Model: ARIMA(3,1,0) 

Coefficients:
          ar1      ar2     ar3
      -0.3418  -0.0426  0.3185
s.e.   0.0681   0.0725  0.0682

sigma^2 estimated as 9.639:  log likelihood=-494
AIC=996   AICc=996   BIC=1009

Seasonally adjusted electrical equipment

fit <- dcmp %>%
    model(arima = ARIMA(season_adjust, approximation = FALSE))
report(fit)

Series: season_adjust 
Model: ARIMA(3,1,1) 

Coefficients:
         ar1     ar2     ar3     ma1
      0.0044  0.0916  0.3698  -0.392
s.e.  0.2201  0.0984  0.0669   0.243

sigma^2 estimated as 9.577:  log likelihood=-493
AIC=995   AICc=996   BIC=1012

Seasonally adjusted electrical equipment

6. ACF plot of residuals from ARIMA(3,1,1) model look like white noise.

Seasonally adjusted electrical equipment

augment(fit) %>%
    features(.resid, ljung_box, lag = 24, dof = 4)

# A tibble: 1 × 3
  .model lb_stat lb_pvalue
  <chr>    <dbl>     <dbl>
1 arima     24.0     0.241

Seasonally adjusted electrical equipment

fit %>%
    forecast %>%
    autoplot(dcmp)

Point forecasts

1. Rearrange ARIMA equation so y_t is on LHS.
   
2. Rewrite equation by replacing t by T+h.
   
3. On RHS, replace future observations by their forecasts, future errors by zero, and past errors by corresponding residuals.

Start with h=1. Repeat for h=2,3,\dots.

Intervals

Intervals

Let’s Try….

For the GDP data (from global_economy):

• fit a suitable ARIMA model to the logged data for all countries
• check the residual diagnostics for Austria;
• produce forecasts of your fitted model for Austria.

Seasonal ARIMA models

where m = number of observations per year.

Common ARIMA models

The US Census Bureau uses the following models most often:

ARIMA(0,1,1)(0,1,1)_m& with log transformation
ARIMA(0,1,2)(0,1,1)_m& with log transformation
ARIMA(2,1,0)(0,1,1)_m& with log transformation
ARIMA(0,2,2)(0,1,1)_m& with log transformation
ARIMA(2,1,2)(0,1,1)_m& with no transformation

Seasonal ARIMA models

The seasonal part of an AR or MA model will be seen in the seasonal lags of the PACF and ACF.

ARIMA(0,0,0)(0,0,1)_{12} will show:

• a spike at lag 12 in the ACF but no other significant spikes.
• The PACF will show exponential decay in the seasonal lags; that is, at lags 12, 24, 36, .

ARIMA(0,0,0)(1,0,0)_{12} will show:

• exponential decay in the seasonal lags of the ACF
• a single significant spike at lag 12 in the PACF.

US electricity production

usmelec %>%
    autoplot(log(Generation)) + theme_ipsum_rc()

US electricity production

usmelec %>%
    autoplot(log(Generation) %>%
        difference(12)) + theme_ipsum_rc()

US electricity production

usmelec %>%
    autoplot(log(Generation) %>%
        difference(12) %>%
        difference()) + theme_ipsum_rc()

US electricity production

usmelec %>%
    gg_tsdisplay(log(Generation) %>%
        difference(12) %>%
        difference(), plot_type = "partial")

US electricity production

• d=1 and D=1 seems necessary
• P=0 and Q=1 suggested by seasonal lags
• p=0 and q=3 suggested by non-seasonal lags.

US electricity production

usmelec %>%
    model(arima = ARIMA(log(Generation) ~ pdq(0, 1, 3) + PDQ(0, 1, 1))) %>%
    report()

Series: Generation 
Model: ARIMA(0,1,3)(0,1,1)[12] 
Transformation: log(Generation) 

Coefficients:
          ma1      ma2      ma3     sma1
      -0.4266  -0.2496  -0.0439  -0.8358
s.e.   0.0462   0.0516   0.0430   0.0262

sigma^2 estimated as 0.0006904:  log likelihood=1045
AIC=-2080   AICc=-2080   BIC=-2059

US electricity production

usmelec %>%
    model(arima = ARIMA(log(Generation))) %>%
    report()

Series: Generation 
Model: ARIMA(1,1,1)(2,1,1)[12] 
Transformation: log(Generation) 

Coefficients:
         ar1      ma1    sar1     sar2     sma1
      0.4116  -0.8483  0.0100  -0.1017  -0.8204
s.e.  0.0617   0.0348  0.0561   0.0529   0.0357

sigma^2 estimated as 0.0006841:  log likelihood=1047
AIC=-2082   AICc=-2082   BIC=-2057

US electricity production

fit <- usmelec %>%
    model(arima = ARIMA(log(Generation)))
gg_tsresiduals(fit)

US electricity production

augment(fit) %>%
    features(.resid, ljung_box, lag = 24, dof = 5)

# A tibble: 1 × 3
  .model lb_stat lb_pvalue
  <chr>    <dbl>     <dbl>
1 arima     44.5  0.000799

US electricity production

usmelec %>%
    model(arima = ARIMA(log(Generation))) %>%
    forecast(h = "3 years") %>%
    autoplot(usmelec) + theme_ipsum_rc()

US electricity production

usmelec %>%
    model(arima = ARIMA(log(Generation))) %>%
    forecast(h = "3 years") %>%
    autoplot({
        usmelec %>%
            filter(year(index) > 2010)
    }) + theme_ipsum_rc()

ARIMA vs ETS

• Myth that ARIMA models are more general than exponential smoothing.

• Linear exponential smoothing models all special cases of ARIMA models.

• Non-linear exponential smoothing models have no equivalent ARIMA counterparts.

• Many ARIMA models have no exponential smoothing counterparts.

• ETS models all non-stationary. Models with seasonality or non-damped trend (or both) have two unit roots; all other models have one unit \rlap{root.}

ARIMA vs ETS

Save Stuff

save.image("./data/FullWorkspace.RData")