Historical perspective

• Developed in the 1950s and 1960s as methods (algorithms) to produce point forecasts.
• Combine a “level”, “trend” (slope) and “seasonal” component to describe a time series.
• The rate of change of the components are controlled by “smoothing parameters”: \(\alpha\), \(\beta\) and \(\gamma\) respectively.
• Need to choose best values for the smoothing parameters (and initial states).
• Equivalent ETS state space models developed in the 1990s and 2000s.

Big idea: control the rate of change (1/2)

\(\alpha\) controls the flexibility of the level (Level - the average value for a specific time period)

• If \(\alpha = 0\), the level never updates (mean)
• If \(\alpha = 1\), the level updates completely (naive)

Big idea: control the rate of change (2/2)

\(\gamma\) controls the flexibility of the seasonality

• If \(\gamma = 0\), the seasonality is fixed (seasonal means)
• If \(\gamma = 1\), the seasonality updates completely (seasonal naive)

A model for levels, trends, and seasonalities

We want a model that captures the level (\(\ell_t\)), trend (\(b_t\)) and seasonality (\(s_t\)).

ETS models

General notation

Simple methods

Time series \(y_1, y_2, \dots, y_T\).

Naive (Random walk forecasts)

\[\hat{y}_{T+h|T} = y_T\]

Simple Exponential Smoothing

Forecast equation

\[ \hat{y}_{T+1|T} = \alpha y_T + \alpha(1-\alpha) y_{T-1} + \alpha(1-\alpha)^2 y_{T-2} + \cdots \]

where \(0 \le \alpha \le 1\).

Simple Exponential Smoothing

Component form

\[\begin{align*} \text{Forecast equation} && \hat{y}_{t+h|t} &= \ell_{t} \\ \text{Smoothing equation} && \ell_{t} &= \alpha y_{t} + (1 - \alpha)\ell_{t-1} \end{align*}\]

• \(\ell_t\) is the level (or the smoothed value) of the series at time t.
• \(\hat{y}_{t+1|t} = \alpha y_t + (1-\alpha) \hat{y}_{t|t-1}\)
• The process has to start somewhere, so we let the first fitted value at time 1 be denoted by \(\ell_0\) (which we will have to estimate) . . .

Simple Exponential Smoothing

Iterate to get exponentially weighted moving average form.

Weighted average form

\[ \hat{y}_{T+1|T} = \sum_{j=0}^{T-1} \alpha(1-\alpha)^j y_{T-j} + (1-\alpha)^T \ell_{0} \]

Optimizing smoothing parameters

• Need to choose best values for \(\alpha\) and \(\ell_0\).
• Similarly to regression, choose optimal parameters by minimizing SSE:

Models and methods

Methods

• Algorithms that return point forecasts.

ETS(A,N,N): SES with additive errors

Component form

\[\begin{align*} \text{Forecast equation} && \hat{y}_{t+h|t} &= \ell_{t} \\ \text{Smoothing equation} && \ell_{t} &= \alpha y_{t} + (1 - \alpha)\ell_{t-1} \end{align*}\]

Forecast error: \(e_t = y_t - \hat{y}_{t|t-1} = y_t - \ell_{t-1}\).

Error correction form

\[\begin{align*} y_t &= \ell_{t-1} + e_t \\ \ell_{t} &= \ell_{t-1} + \alpha (y_t - \ell_{t-1}) \\ &= \ell_{t-1} + \alpha e_t \end{align*}\]

Specify probability distribution for \(e_t\), we assume \(e_t = \varepsilon_t \sim \text{NID}(0, \sigma^2)\).

ETS(A,N,N): SES with additive errors

\[\begin{align*} \text{Measurement equation} && y_t &= \ell_{t-1} + \varepsilon_t \\ \text{State equation} && \ell_t &= \ell_{t-1} + \alpha \varepsilon_t \end{align*}\]

where \(\varepsilon_t \sim \text{NID}(0, \sigma^2)\).

• “Innovations” or “single source of error” because equations have the same error process, \(\varepsilon_t\).
• Measurement equation: relationship between observations and states.
• State equation(s): evolution of the state(s) through time.

ETS(M,N,N): SES with multiplicative errors

• Specify relative errors \(\varepsilon_t = \frac{y_t - \hat{y}_{t|t-1}}{\hat{y}_{t|t-1}} \sim \text{NID}(0, \sigma^2)\)
• Substituting \(\hat{y}_{t|t-1} = \ell_{t-1}\) gives:
  • \(y_t = \ell_{t-1} + \ell_{t-1} \varepsilon_t\)
  • \(e_t = y_t - \hat{y}_{t|t-1} = \ell_{t-1} \varepsilon_t\) . . .

ETS(M,N,N): SES with multiplicative errors

Measurement equation

\[\begin{align*} \text{Measurement equation} && y_t &= \ell_{t-1} (1 + \varepsilon_t) \\ \text{State equation} && \ell_t &= \ell_{t-1} (1 + \alpha \varepsilon_t) \end{align*}\]

• Models with additive and multiplicative errors with the same parameters generate the same point forecasts but different prediction intervals.

Example: Algerian Exports (2/4)

components(fit) |> autoplot()

Example: Algerian Exports (4/4)

fit |>
  forecast(h = 5) |>
  autoplot(algeria_economy) +
  labs(y = "% of GDP", title = "Exports: Algeria")

Holt’s linear trend

Component form

\[\begin{align*} \text{Forecast} && \hat{y}_{t+h|t} &= \ell_{t} + hb_{t} \\ \text{Level} && \ell_{t} &= \alpha y_{t} + (1 - \alpha)(\ell_{t-1} + b_{t-1}) \\ \text{Trend} && b_{t} &= \beta^*(\ell_{t} - \ell_{t-1}) + (1 - \beta^*)b_{t-1} \end{align*}\]

• Two smoothing parameters \(\alpha\) and \(\beta^*\) (\(0 \le \alpha, \beta^* \le 1\)).
• \(\ell_t\) level: weighted average between \(y_t\) and one-step ahead forecast for time \(t\), \((\ell_{t-1} + b_{t-1} = \hat{y}_{t|t-1})\)
• \(b_t\) slope: weighted average of \((\ell_{t} - \ell_{t-1})\) and \(b_{t-1}\), current and previous estimate of slope.
• Choose \(\alpha, \beta^*, \ell_0, b_0\) to minimize SSE.

ETS(A,A,N)

Holt’s linear method with additive errors.

• Assume \(\varepsilon_t = y_t - \ell_{t-1} - b_{t-1} \sim \text{NID}(0, \sigma^2)\).
• Substituting into the error correction equations for Holt’s linear method:

ETS(M,A,N)

Holt’s linear method with multiplicative errors.

• Assume \(\varepsilon_t=\frac{y_t-(\ell_{t-1}+b_{t-1})}{(\ell_{t-1}+b_{t-1})}\)
• Following a similar approach as above, the innovations state space model underlying Holt’s linear method with multiplicative errors is specified as

Example: Australian population (1/4)

aus_economy <- global_economy |> filter(Code == "AUS") |>
  mutate(Pop = Population / 1e6)
fit <- aus_economy |>
  model(AAN = ETS(Pop ~ error("A") + trend("A") + season("N")))
report(fit)

Series: Pop 
Model: ETS(A,A,N) 
  Smoothing parameters:
    alpha = 0.9999 
    beta  = 0.3266366 

  Initial states:
     l[0]      b[0]
 10.05414 0.2224818

  sigma^2:  0.0041

      AIC      AICc       BIC 
-76.98569 -75.83184 -66.68347 

Example: Australian population (2/4)

components(fit) |> autoplot()

Example: Australian population (3/4)

components(fit) |>
  left_join(fitted(fit), by = c("Country", ".model", "Year"))

# A dable: 59 x 8 [1Y]
# Key:     Country, .model [1]
# :        Pop = lag(level, 1) + lag(slope, 1) + remainder
   Country   .model  Year   Pop level slope remainder .fitted
   <fct>     <chr>  <dbl> <dbl> <dbl> <dbl>     <dbl>   <dbl>
 1 Australia AAN     1959  NA    10.1 0.222 NA           NA  
 2 Australia AAN     1960  10.3  10.3 0.222 -0.000145    10.3
 3 Australia AAN     1961  10.5  10.5 0.217 -0.0159      10.5
 4 Australia AAN     1962  10.7  10.7 0.231  0.0418      10.7
 5 Australia AAN     1963  11.0  11.0 0.223 -0.0229      11.0
 6 Australia AAN     1964  11.2  11.2 0.221 -0.00641     11.2
 7 Australia AAN     1965  11.4  11.4 0.221 -0.000314    11.4
 8 Australia AAN     1966  11.7  11.7 0.235  0.0418      11.6
 9 Australia AAN     1967  11.8  11.8 0.206 -0.0869      11.9
10 Australia AAN     1968  12.0  12.0 0.208  0.00350     12.0
# ℹ 49 more rows

Example: Australian population (4/4)

fit |>
  forecast(h = 10) |>
  autoplot(aus_economy) +
  labs(y = "Millions", title = "Population: Australia")

Damped trend method

Component form

\[\begin{align*} \hat{y}_{t+h|t} &= \ell_{t} + (\phi + \phi^2 + \dots + \phi^{h})b_{t} \\ \ell_{t} &= \alpha y_{t} + (1 - \alpha)(\ell_{t-1} + \phi b_{t-1}) \\ b_{t} &= \beta^*(\ell_{t} - \ell_{t-1}) + (1 - \beta^*)\phi b_{t-1} \end{align*}\]

• Damping parameter \(0<\phi<1\).
• If \(\phi=1\), identical to Holt’s linear trend.
• As \(h\rightarrow\infty\), \(\hat{y}{T+h}{T}\rightarrow \ell_T+\phi b_T/(1-\phi)\).
• Short-run forecasts trended, long-run forecasts constant.

Example: Australian population ct. (1/3)

aus_economy |>
  model(holt = ETS(Pop ~ error("A") + trend("Ad") + season("N"))) |>
  forecast(h = 20) |>
  autoplot(aus_economy)

Example: Australian population ct. (2/3)

fit <- aus_economy |>
  filter(Year <= 2010) |>
  model(
    ses = ETS(Pop ~ error("A") + trend("N") + season("N")),
    holt = ETS(Pop ~ error("A") + trend("A") + season("N")),
    damped = ETS(Pop ~ error("A") + trend("Ad") + season("N"))
  )

tidy(fit)
accuracy(fit)

Example: Australian population ct. (3/3)

Your turn

prices contains the price of a dozen eggs in the United States from 1900–1993

1. Use SES and Holt’s method (with and without damping) to forecast “future” data.

2. Which method gives the best training RMSE?
3. Are these RMSE values comparable?
4. Do the residuals from the best fitting method look like white noise?

Holt-Winters additive method

Holt and Winters extended Holt’s method to capture seasonality.

Component form

\[\begin{align*} \hat{y}_{t+h|t} &= \ell_{t} + hb_{t} + s_{t+h-m(k+1)} \\ \ell_{t} &= \alpha(y_{t} - s_{t-m}) + (1 - \alpha)(\ell_{t-1} + b_{t-1}) \\ b_{t} &= \beta^*(\ell_{t} - \ell_{t-1}) + (1 - \beta^*)b_{t-1} \\ s_{t} &= \gamma (y_{t} - \ell_{t-1} - b_{t-1}) + (1 - \gamma)s_{t-m} \end{align*}\]

• \(k=\) integer part of \((h-1)/m\). Ensures estimates from the final year are used for forecasting.
• Parameters:  \(0\le \alpha\le 1\),  \(0\le \beta^*\le 1\),  \(0\le \gamma\le 1-\alpha\)  and \(m=\) period of seasonality (e.g. \(m=4\) for quarterly data).

Holt-Winters additive method

• Seasonal component is usually expressed as \(s_{t} = \gamma^* (y_{t}-\ell_{t})+ (1-\gamma^*)s_{t-m}.\)
• Substitute in for \(\ell_t\): \(s_{t} = \gamma^*(1-\alpha) (y_{t}-\ell_{t-1}-b_{t-1})+ [1-\gamma^*(1-\alpha)]s_{t-m}\)
• We set \(\gamma=\gamma^*(1-\alpha)\).
• The usual parameter restriction is \(0\le\gamma^*\le1\), which translates to \(0\le\gamma\le(1-\alpha)\).

ETS(A,A,A)

Holt-Winters additive method with additive errors.

• Forecast errors: \(\varepsilon_{t} = y_t - \hat{y}_{t|t-1}\)
• \(k\) is integer part of \((h-1)/m\).

Holt-Winters multiplicative method

Seasonal variations change in proportion to the level of the series.

Component form

\[\begin{align*} \hat{y}_{t+h|t} &= (\ell_{t} + hb_{t})s_{t+h-m(k+1)} \\ \ell_{t} &= \alpha \frac{y_{t}}{s_{t-m}} + (1 - \alpha)(\ell_{t-1} + b_{t-1}) \\ b_{t} &= \beta^*(\ell_{t} - \ell_{t-1}) + (1 - \beta^*)b_{t-1} \\ s_{t} &= \gamma \frac{y_{t}}{(\ell_{t-1} + b_{t-1})} + (1 - \gamma)s_{t-m} \end{align*}\]

• \(k\) is integer part of \((h-1)/m\).
• Additive method: \(s_t\) in absolute terms — within each year \(\sum_i s_i \approx 0\).
• Multiplicative method: \(s_t\) in relative terms — within each year \(\sum_i s_i \approx m\).

ETS(M,A,M)

Holt-Winters multiplicative method with multiplicative errors.

\[\begin{align*} \text{Forecast equation} && \hat{y}_{t+h|t} &= (\ell_{t} + hb_{t}) s_{t+h-m(k+1)} \\ \text{Observation equation} && y_t &= (\ell_{t-1} + b_{t-1})s_{t-m}(1 + \varepsilon_t) \\ \text{State equations} && \ell_t &= (\ell_{t-1} + b_{t-1})(1 + \alpha \varepsilon_t) \\ && b_t &= b_{t-1} + \beta(\ell_{t-1} + b_{t-1}) \varepsilon_t \\ && s_t &= s_{t-m}(1 + \gamma \varepsilon_t) \end{align*}\]

• Forecast errors: \(\varepsilon_{t} = (y_t - \hat{y}_{t|t-1})/\hat{y}_{t|t-1}\)
• \(k\) is integer part of \((h-1)/m\).

Example: Australian holiday tourism (1/2)

aus_holidays <- tourism |>
  filter(Purpose == "Holiday") |>
  summarise(Trips = sum(Trips))
fit <- aus_holidays |>
  model(
    additive = ETS(Trips ~ error("A") + trend("A") + season("A")),
    multiplicative = ETS(Trips ~ error("M") + trend("A") + season("M"))
  )
fc <- fit |> forecast()

Example: Australian holiday tourism (2/2)

fc |>
  autoplot(aus_holidays, level = NULL) +
  labs(y = "Thousands", title = "Overnight trips")

Holt-Winters damped method

Often the single most accurate forecasting method for seasonal data:

\[\begin{align*} \hat{y}_{t+h|t} &= [\ell_{t} + (\phi + \phi^2 + \dots + \phi^{h})b_{t}]s_{t+h-m(k+1)} \\ \ell_{t} &= \alpha \left(\frac{y_{t}}{s_{t-m}}\right) + (1 - \alpha)(\ell_{t-1} + \phi b_{t-1}) \\ b_{t} &= \beta^*(\ell_{t} - \ell_{t-1}) + (1 - \beta^*)\phi b_{t-1} \\ s_{t} &= \gamma \left(\frac{y_{t}}{(\ell_{t-1} + \phi b_{t-1})}\right) + (1 - \gamma)s_{t-m} \end{align*}\]

Your turn

Apply Holt-Winters’ multiplicative method to the Gas data from aus_production.

1. Why is multiplicative seasonality necessary here?
2. Experiment with making the trend damped.
3. Check that the residuals from the best method look like white noise.

Exponential smoothing methods

Exponential smoothing methods

ETS models

Additive Error

Estimating ETS models

• Smoothing parameters \(\alpha\), \(\beta\), \(\gamma\) and \(\phi\), and the initial states \(\ell_0\), \(b_0\), \(s_0,s_{-1},\dots,s_{-m+1}\) are estimated by maximizing the “likelihood” = the probability of the data arising from the specified model.
• For models with additive errors equivalent to minimizing SSE.
• For models with multiplicative errors, not equivalent to minimizing SSE.

Innovations state space models

Estimation

\[L^*(\theta,\textbf{x}_0) = T\log\!\left(\sum_{t=1}^T \varepsilon^2_t\!\right) + 2\sum_{t=1}^T \log|k(\textbf{x}_{t-1})|\]

\[= -2\log(\text{Likelihood}) + \text{constant}\]

• Estimate parameters \(\theta = (\alpha,\beta,\gamma,\phi)\) and initial states \(\textbf{x}_0 = (\ell_0,b_0,s_0,s_{-1},\dots,s_{-m+1})\) by minimizing \(L^*\).

Parameter restrictions

Usual region

• Traditional restrictions in the methods \(0< \alpha,\beta^*,\gamma^*,\phi<1\)(equations interpreted as weighted averages).
• In models we set \(\beta=\alpha\beta^*\) and \(\gamma=(1-\alpha)\gamma^*\).
• Therefore \(0< \alpha <1\),    \(0 < \beta < \alpha\)    and \(0< \gamma < 1-\alpha\).
• \(0.8<\phi<0.98\) — to prevent numerical difficulties.

Model selection

Akaike’s Information Criterion

\[AIC = -2\log(\text{L}) + 2k\]

where \(L\) is the likelihood and \(k\) is the number of parameters initial states estimated in the model.

Automatic forecasting

From Hyndman et al. (IJF, 2002):

• Apply each model that is appropriate to the data. Optimize parameters and initial values using MLE (or some other criterion).
• Select best method using AICc:
• Produce forecasts using best method.
• Obtain forecast intervals using underlying state space model.

Residuals

Response residuals \[\hat{e}_t = y_t - \hat{y}_{t|t-1}\]

Example: Australian holiday tourism (7/9)

aus_holidays <- tourism |>
  filter(Purpose == "Holiday") |>
  summarise(Trips = sum(Trips))
fit <- aus_holidays |>
  model(ets = ETS(Trips)) |>
  report()

Series: Trips 
Model: ETS(M,N,M) 
  Smoothing parameters:
    alpha = 0.3578226 
    gamma = 0.0009685565 

  Initial states:
     l[0]      s[0]     s[-1]    s[-2]    s[-3]
 9666.501 0.9430367 0.9268433 0.968352 1.161768

  sigma^2:  0.0022

     AIC     AICc      BIC 
1331.372 1332.928 1348.046 

Example: Australian holiday tourism (9/9)

fit |>
  augment()

# A tsibble: 80 x 6 [1Q]
# Key:       .model [1]
   .model Quarter  Trips .fitted .resid   .innov
   <chr>    <qtr>  <dbl>   <dbl>  <dbl>    <dbl>
 1 ets    1998 Q1 11806.  11230.  576.   0.0513 
 2 ets    1998 Q2  9276.   9532. -257.  -0.0269 
 3 ets    1998 Q3  8642.   9036. -393.  -0.0435 
 4 ets    1998 Q4  9300.   9050.  249.   0.0275 
 5 ets    1999 Q1 11172.  11260.  -88.0 -0.00781
 6 ets    1999 Q2  9608.   9358.  249.   0.0266 
 7 ets    1999 Q3  8914.   9042. -129.  -0.0142 
 8 ets    1999 Q4  9026.   9154. -129.  -0.0140 
 9 ets    2000 Q1 11071.  11221. -150.  -0.0134 
10 ets    2000 Q2  9196.   9308. -111.  -0.0120 
# ℹ 70 more rows

Some unstable models

• Some of the combinations of (Error, Trend, Seasonal) can lead to numerical difficulties; see equations with division by a state.
• These are: ETS(A,N,M), ETS(A,A,M), ETS(A,\(A_d\),M).
• Models with multiplicative errors are useful for strictly positive data, but are not numerically stable with data containing zeros or negative values. In that case only the six fully additive models will be applied.

Exponential smoothing models

Additive Error

Multiplicative Error

Forecasting with ETS models

Traditional point forecasts: iterate the equations for \(t=T+1,T+2,\dots,T+h\) and set all \(\varepsilon_t=0\) for \(t>T\).

• Not the same as \(\text{E}(y_{t+h} | \textbf{x}_t)\) unless seasonality is additive.
• fable uses \(\text{E}(y_{t+h} | \textbf{x}_t)\).
• Point forecasts for ETS(A,*,*) are identical to ETS(M,*,*) if the parameters are the same.

Forecasting with ETS models

Prediction intervals: can only be generated using the models.

• The prediction intervals will differ between models with additive and multiplicative errors.
• Exact formulae for some models.
• More general to simulate future sample paths, conditional on the last estimate of the states, and to obtain prediction intervals from the percentiles of these simulated future paths.

Prediction intervals

PI for most ETS models: \(\hat{y}_{T+h|T} \pm c \sigma_h\), where \(c\) depends on coverage probability and \(\sigma_h\) is forecast standard deviation.

We can see them in the book HERE.