1. Overview
   
2. Simple methods
   
3. Trend methods
   
4. Seasonal methods
5. Taxonomy of Exponential Smoothing Methods
6. Innovations State Space Models
7. Estimation and Model Selection
8. Forecasting with ETS models

• When making predictions, weight recent observations more than older observations
  
• Specifically, decrease weights geometrically / exponentially when moving backward
  
• Window function (a.k.a. apodization function or tapering function)
  • zero-valued outside of some chosen interval
    
  • normally symmetric around middle of interval
    
  • typically near maximum in the middle and tapering away from there
    
• Easily learned and applied (i.e. signal processing community in 1940’s)

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

Random walk forecast

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

Average forecast

\[ \begin{aligned} \hat{y}_{T+h|T} = \frac1T\sum_{t=1}^T y_t \end{aligned} \]

Consider a happy median that privileges recent information…

Forecast equation \[ \begin{aligned} \hat{y}_{T+1|T} = \alpha y_T + \alpha(1-\alpha) y_{T-1} + \alpha(1-\alpha)^2 y_{T-2}+ \cdots \end{aligned} \] where \(0 \le \alpha \le1\)

Weighted moving average with weights that decrease exponentially

• for \(y_{T}\) and \(\alpha = .2\), the weighting is \(.2\)
• for \(y_{T-1}\), the weighting is \(.2 * (1 - .2) = .16\)
• for \(y_{T-2}\), \(.2 * (1 - .2)^2 = .128\)
• etc.

Component form: forecast equation \[ \begin{align*} \hat{y}_{t+h|t} = \ell_{t} \end{align*} \]

Component form: smoothing equation \[ \begin{align*} \ell_{t} = \alpha y_{t} + (1 - \alpha)\ell_{t-1} \end{align*} \]

\(\ell_t\) is level / smoothed value of series at time \(t\)

\(\hat{y}_{t+1|t} = \alpha y_t + (1-\alpha) \hat{y}_{t|t-1}\)
Iterate for exponentially weighted moving average form

Weighted average form \[ \begin{align*} \hat{y}_{T+1|T}=\sum_{j=0}^{T-1} \alpha(1-\alpha)^j y_{T-j}+(1-\alpha)^T \ell_{0} \end{align*} \]

To choose value for \(\alpha\) and \(\ell_0\), minimize SSE: \[ \begin{aligned} \text{SSE}=\sum_{t=1}^T(y_t - \hat{y}_{t|t-1})^2. \end{aligned} \]

Unlike regression, no closed-form solution so use numerical optimization

oildata <- window(oil, start=1996)
fc <- ses(oildata, h=5)
autoplot(fc) +
  autolayer(fitted(fc), series="Fitted") +
  ylab("Oil (millions of tonnes)") + xlab("Year")

Component form: forecast \[ \begin{align*} \hat{y}_{t+h|t} &= \ell_{t+h}b_{t} \end{align*} \]

Component form: level \[ \begin{align*} \ell_{t} &= \alpha y_{t} + (1 - \alpha)(\ell_{t-1} + b_{t-1}) \end{align*} \]

Component form: trend \[ \begin{align*} 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^*\le1)\)
  
• \(\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
  
• Select \(\alpha, \beta^*, \ell_0, b_0\) to minimise SSE

window(ausair, start=1990, end=2004) %>%
  holt(h=5, PI=FALSE) %>%
  autoplot()

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

Holt’s method extended 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\) is integer part of \((h-1)/m\), ensuring final year estimates in forecast
• Parameters: \(0\le \alpha\le 1\), \(0\le \beta^*\le 1\), \(0\le \gamma\le 1-\alpha\) and \(m\) for period of seasonality (e.g. \(m=4\) for quarterly data)

• Seasonal component usually expressed as \(s_{t} = \gamma^* (y_{t}-\ell_{t}) + (1-\gamma^*)s_{t-m}\)
• Substitute for \(\ell_t\): \(s_{t} = \gamma^* (1-\alpha) (y_{t}-\ell_{t-1}-b_{t-1}) + [1-\gamma^*(1-\alpha)]s_{t-m}\)
• Set \(\gamma=\gamma^*(1-\alpha)\)
• Usual parameter restriction is \(0\le\gamma^*\le1\) or \(0\le\gamma\le(1-\alpha)\)

When seasonal variations change in proportion to level of 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*} \]

• For additive method, \(s_t\) is in absolute terms: within each year \(\sum_i s_i \approx 0\)
• For multiplicative method, \(s_t\) is in relative terms: within each year \(\sum_i s_i \approx m\)

aust <- window(austourists,start=2005)
fit1 <- hw(aust,seasonal="additive")
fit2 <- hw(aust,seasonal="multiplicative")

tmp <- cbind(Data=aust,
  "HW additive forecasts" = fit1[["mean"]],
  "HW multiplicative forecasts" = fit2[["mean"]])

autoplot(tmp) + xlab("Year") +
  ylab("International visitor night in Australia (millions)") +
  scale_color_manual(name="",
    values=c('#000000','#1b9e77','#d95f02'),
    breaks=c("Data","HW additive forecasts","HW multiplicative forecasts"))

addstates <- fit1$model$states[,1:3]
multstates <- fit2$model$states[,1:3]
colnames(addstates) <- colnames(multstates) <-
  c("level","slope","season")
p1 <- autoplot(addstates, facets=TRUE) + xlab("Year") +
  ylab("") + ggtitle("Additive states")
p2 <- autoplot(multstates, facets=TRUE) + xlab("Year") +
  ylab("") + ggtitle("Multiplicative states")
gridExtra::grid.arrange(p1,p2,ncol=2)

\[ \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(y_{t} / s_{t-m}) + (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 \frac{y_{t}}{(\ell_{t-1} + \phi b_{t-1})} + (1 - \gamma)s_{t-m} \end{align*} \]

Often single most accurate forecasting method for seasonal data

air <- window(ausair, start = 1990, end = 2004)
fit3 <- holt(air, h = 15)
fit4 <- holt(air, damped = TRUE, phi = .9, h = 15)
autoplot(air) +
  autolayer(fit3, series = 'Holt\'s method', PI = FALSE) +
  autolayer(fit4, series = 'Dampled Holt\'s method', PI = FALSE) +
  ggtitle('Forecasts from Holt\'s method') +
  xlab('Year') +
  ylab('Air passengers in Australia (millions)') +
  guides(color = guide_legend(title = 'Forecast'))

Short hand representation of methods:

• (N, N): Simple exponential smoothing
• (A, N): Holt’s linear method
• (A\(_d\), N): Additive damped trend method
• (A, A): Additive Holt-Winters’ method
• (A, M): Multiplicative Holt-Winters’ method
• (A\(_d\), M): Holt-Winters’ damped method

1. Methods: Apply algorithms to obtain point forecasts
2. Models: Apply stochastic generating process to forecast an entire distribution.
   • Generate point forecasts and forecast intervals
   • Beneficial for model selection

1. Each model contains two equations:
   • Observation equation for observed data components.
     
   • State equation for unobserved changes in level, trend, and seasonal components.

2. Each method contains two models distinguished by either additive or multiplicative errors.

3. Error, Trend, Seasonal (ETS) labeling is used to differentiate models, methods, and type.

18 variations of ETS models exist and are classified as follows:

1. Error: \(\{\)\(A\) | \(M\)\(\}\)
   • Additive Errors
   • Multipicative Errors

2. Trend: \(\{\)\(N\) | \(A\) | \(A_d\)\(\}\)

3. Seasonal: \(\{\)\(N\) | \(A\) | \(M\)\(\}\)

• \((\)\(A\), \(N\), \(N\)\()\): Simple exponential smoothing; additive errors
• \((\)\(M\), \(N\), \(N\)\()\): Simple exponential smoothing; multiplicative errors
  
• \((\)\(A\), \(A\), \(N\)\()\): Holt’s linear method; additive errors
• \((\)\(M\), \(A\), \(N\)\()\): Holt’s linear method; multiplicative errors
• \((\)\(M\), \(A\), \(M\)\()\): Multiplicative Holt-Winters’ method; multiplicative errors

• 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 maximising the “likelihood” or the probability of the data arising from the specified model.
• Constraints are placed on parameters for weighted averages and to prevent numerical difficulties in estimating the models.
• We will estimate models with the ets() function in the forecast package.

Akaike’s Information Criterion (AIC):

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

Corrected Akaike’s Information Criterion (AIC\(_c\)):

\[ \text{AIC}_{\text{c}} = \text{AIC} + \frac{2(k+1)(k+2)}{T-k} \]

Schwarz’s Bayesian Information Criterion (BIC):

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

Automatic forecasting using the ets() function:

• Automatically selects models when arguments are set to default.

• Produces forecasts using best method.

• Obtain forecast intervals using underlying state spacemodel.

Tourist example from text:

aust <- window(austourists, start=2005) 
fit <- ets(aust) 
summary(fit)

## ETS(M,A,M) 
## 
## Call:
##  ets(y = aust) 
## 
##   Smoothing parameters:
##     alpha = 0.1908 
##     beta  = 0.0392 
##     gamma = 2e-04 
## 
##   Initial states:
##     l = 32.3679 
##     b = 0.9281 
##     s = 1.0218 0.9628 0.7683 1.2471
## 
##   sigma:  0.0383
## 
##      AIC     AICc      BIC 
## 224.8628 230.1569 240.9205 
## 
## Training set error measures:
##                      ME     RMSE     MAE        MPE     MAPE     MASE
## Training set 0.04836907 1.670893 1.24954 -0.1845609 2.692849 0.409454
##                   ACF1
## Training set 0.2005962

Model evaluation:

##           AIC     AICc      BIC
## [1,] 224.8628 230.1569 240.9205

Train set error measures:

##                      ME     RMSE     MAE        MPE     MAPE     MASE
## Training set 0.04836907 1.670893 1.24954 -0.1845609 2.692849 0.409454
##                   ACF1
## Training set 0.2005962

cbind('Residuals' = residuals(fit), 
      'Forecast errors' = residuals(fit, type='response')) %>% 
  autoplot(facet=TRUE) + xlab("Year") + ylab("")

Models only! Cannot be generated using the methods.

• The prediction intervals will differ between models with additive and multiplicative errors.
  
• 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.
  
• Options are available in R using the forecast function in the forecast package.

fit %>% forecast(h=8) %>%
  autoplot() +
  ylab("International visitor night in Australia (millions)")

1. https://github.com/robjhyndman/ETC3550Slides/raw/master/7-exponentialsmoothing.pdf
2. https://en.wikipedia.org/wiki/Exponential_smoothing
3. https://machinelearningmastery.com/exponential-smoothing-for-time-series-forecasting-in-python/