Time series beer with Holt-Winters algorithms

The time series beer is in the fma package:

# make sure to type in Rstudio console: install.packages("fma")
require("fma")
beer

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1991 164 148 152 144 155 125 153 146 138 190 192 192
## 1992 147 133 163 150 129 131 145 137 138 168 176 188
## 1993 139 143 150 154 137 129 128 140 143 151 177 184
## 1994 151 134 164 126 131 125 127 143 143 160 190 182
## 1995 138 136 152 127 151 130 119 153

Description of beer: Monthly Australian beer production: Jan 1991 to Aug 1995.

Visualisation of time series

Time plot

You can also embed plots, for example the time plot:

plot(beer)

Time plot with ACF and PACF

The ACF and PACF can be displayed with the time plot:

tsdisplay(beer,lag.max=70)

Other packages exist for visualisation. For instance you may prefer using packages ggplot2

# make sure to type in Rstudio console: install.packages("ggplot2")
require(ggplot2)
ggtsdisplay(beer,lag.max = 70)

The timeplot shows seasonality, no real trend and some noise.
For the time series beer, the ACF displays maxima at lags 12, 24, 36 (…) characteristic of seasonality with a period of 12 months (as it is monthly data).

Exercise: Show how to use the R functions acf() and pacf() for the time series beer. Explain the the differences with tsdisplay() (hint: look at the x-axis)

Exercise: Comment on the information provided by the PACF for the time series beer.

seasonplot()

The season plot confirms there is no trend.

seasonplot(beer)

Function ts()

Note that beer time series is stored in fma package explicitly with this information of periodicity 12. For the mink time series which has no frequency information stored, then you will need to use the function ts to indicate a period (if you see one!) for this time series e.g.:

tsdisplay(ts(mink,freq=10))

seasonplot(ts(mink,freq=10))

Function decompose()

Note that the function decompose can be visually misleading:

plot(decompose(beer))

It may look like you have trend!

Exercise: Read chapter 6 at
https://otexts.com/fpp2/decomposition.html . What other R function can be used for decomposing time series ?

Fitting Holt Winters algorithms

Only seasonal algorithms are appropriate for fitting the time series beer.

HoltWinters(beer, seasonal="additive")$SSE

## [1] 4164.097

HoltWinters(beer, seasonal="multiplicative")$SSE

## [1] 3851.269

Using SSE (Sum of Square Errors), the seasonal Holt Winters algorithm (SHWx) algorithm fits best the time series beer (lowest SSE).

Exercise: How to compute MAPE with R ? Would the same Holt-Winters algorithm selected with MAPE ?

Computing forecasts

HoltWinters(beer, seasonal="multiplicative")

## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
## 
## Call:
## HoltWinters(x = beer, seasonal = "multiplicative")
## 
## Smoothing parameters:
##  alpha: 0.02225387
##  beta : 0.6221156
##  gamma: 0.6728938
## 
## Coefficients:
##            [,1]
## a   146.0859896
## b     0.2459431
## s1    0.9794233
## s2    1.0969580
## s3    1.2794748
## s4    1.2586982
## s5    0.9702512
## s6    0.9337021
## s7    1.0637137
## s8    0.8905903
## s9    0.9948355
## s10   0.8842001
## s11   0.8398353
## s12   1.0240620

predict(HoltWinters(beer, seasonal="multiplicative"))

##           Sep
## 1995 143.3209

Using the formula \[ F_{n+k} = (L_n+k \times b_n)\times S_{n+k-s} \] with \(k=1\), \(s=12\) (frequency for beer data), and with the notation used in the lecturenotes https://www.scss.tcd.ie/Rozenn.Dahyot/ST3010/RzDTimeSeriesForecasting.pdf (chapter 8), the prediction can be computed (instead of using the R command predict()) as follow

k=1  # prediction computed 1 month ahead
s=12 # frequency 12 month 
Ln=HoltWinters(beer, seasonal="multiplicative")$coefficients[1]
bn=HoltWinters(beer, seasonal="multiplicative")$coefficients[2]
snk= HoltWinters(beer, seasonal="multiplicative")$coefficients[3]
(Ln+k*bn)*snk

##        a 
## 143.3209

Likewise prediction can be computed for \(k=5\) months ahead from the last observation of the time series beer :

k=5  # prediction computed 1 month ahead
s=12 # frequency 12 month 
Ln=HoltWinters(beer, seasonal="multiplicative")$coefficients[1]
bn=HoltWinters(beer, seasonal="multiplicative")$coefficients[2]
snk= HoltWinters(beer, seasonal="multiplicative")$coefficients[7]
(Ln+k*bn)*snk

##        a 
## 142.9332

predict(HoltWinters(beer, seasonal="multiplicative"),n.ahead=5)[5]

## [1] 142.9332

Using forecast()

forecast(HoltWinters(beer, seasonal="multiplicative"),h=1)

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Sep 1995       143.3209 135.1236 151.5182 130.7841 155.8577

plot(forecast(HoltWinters(beer, seasonal="multiplicative")))

Exercise: Explain how the R function forecast() relates to predict(). Explain how prediction intervals are computed.