Tidy time series analysis and forecasting

Modelling and forecasting

10th June 2024 @ UseR! 2024

Mitchell O’Hara-Wild, Nectric

Useful links

social.mitchelloharawild.com

workshop.nectric.com.au/user2024/

mitchelloharawild/workshop-fable-user2024

Forecasting can be difficult!

Many forecasts are catastrophically wrong, and better described as ‘hopecasting’.

Which is easier to forecast?

  1. daily electricity demand in 3 days time
  2. timing of next Halley’s comet appearance
  3. time of sunrise this day next year
  4. Google stock price tomorrow
  5. Google stock price in 6 months time
  6. maximum temperature tomorrow
  7. exchange rate of $US/AUS next week
  8. total sales of drugs in Australian pharmacies next month

Factors affecting forecastability

Something is easier to forecast if…

  • we have a good understanding of the factors that contribute to it,
  • there is lots of data available,
  • the forecasts cannot affect the thing we are trying to forecast,
  • there is relatively low natural/unexplainable random variation,
  • the future is somewhat similar to the past.

How are good forecasts made?

Statistically! (or computationally!)

Consider possible futures for this series:

What makes forecasts good?

Good forecasts…

  • Are uncertain (probabilistic)
  • Useful for decision making
  • Capture patterns in the data

Good forecasts aren’t necessarily accurate!

Usually accuracy is important, but not always.

Forecasts (usually) assume that interventions won’t occur, e.g. COVID-19 measures.

Some simple forecasting models

Your turn

How would you forecast visitors to Australia?

Hint: Consider the patterns seen in the data, and how statistics/code can mimic the trend, seasonality, etc.

Some simple forecasting models

The naive model

Simply use the most recent observation.

Some simple forecasting models

The seasonal naive model

Use the most recent ‘seasonal’ observation.

Some simple forecasting models

The naive with drift model

Extrapolate a line between the first and last observation.

Some simple forecasting models

The seasonal naive with drift model

Extrapolate between first and last with seasonality.

Estimating a model

Models are estimated using model():

pbs_scripts |> 
  model(
    NAIVE(Scripts),
    SNAIVE(Scripts),
    NAIVE(Scripts ~ drift()),
    SNAIVE(Scripts ~ drift()),
  )
# A mable: 1 x 4
  `NAIVE(Scripts)` `SNAIVE(Scripts)` `NAIVE(Scripts ~ drift())`
           <model>           <model>                    <model>
1          <NAIVE>          <SNAIVE>              <RW w/ drift>
# ℹ 1 more variable: `SNAIVE(Scripts ~ drift())` <model>

The mable

A mable (model table) contains 1 row per time series, and 1 column for each specified model.

Producing forecasts

Forecasts are made using forecast():

pbs_scripts <- PBS |>
  summarise(Scripts = sum(Scripts))
pbs_scripts |> 
  model(
    NAIVE(Scripts),
    SNAIVE(Scripts),
    NAIVE(Scripts ~ drift()),
    SNAIVE(Scripts ~ drift()),
  ) |> 
  forecast(h = "2 years")
# A fable: 96 x 4 [1M]
# Key:     .model [4]
   .model            Month             Scripts    .mean
   <chr>             <mth>              <dist>    <dbl>
 1 NAIVE(Scripts) 2008 Jul N(1.2e+07, 2.2e+12) 12123769
 2 NAIVE(Scripts) 2008 Aug N(1.2e+07, 4.4e+12) 12123769
 3 NAIVE(Scripts) 2008 Sep N(1.2e+07, 6.6e+12) 12123769
 4 NAIVE(Scripts) 2008 Oct N(1.2e+07, 8.8e+12) 12123769
 5 NAIVE(Scripts) 2008 Nov N(1.2e+07, 1.1e+13) 12123769
 6 NAIVE(Scripts) 2008 Dec N(1.2e+07, 1.3e+13) 12123769
 7 NAIVE(Scripts) 2009 Jan N(1.2e+07, 1.5e+13) 12123769
 8 NAIVE(Scripts) 2009 Feb N(1.2e+07, 1.8e+13) 12123769
 9 NAIVE(Scripts) 2009 Mar   N(1.2e+07, 2e+13) 12123769
10 NAIVE(Scripts) 2009 Apr N(1.2e+07, 2.2e+13) 12123769
# ℹ 86 more rows

Producing forecasts

# A fable: 96 x 4 [1M]
# Key:     .model [4]
   .model            Month             Scripts    .mean
   <chr>             <mth>              <dist>    <dbl>
 1 NAIVE(Scripts) 2008 Jul N(1.2e+07, 2.2e+12) 12123769
 2 NAIVE(Scripts) 2008 Aug N(1.2e+07, 4.4e+12) 12123769
 3 NAIVE(Scripts) 2008 Sep N(1.2e+07, 6.6e+12) 12123769
 4 NAIVE(Scripts) 2008 Oct N(1.2e+07, 8.8e+12) 12123769
 5 NAIVE(Scripts) 2008 Nov N(1.2e+07, 1.1e+13) 12123769
 6 NAIVE(Scripts) 2008 Dec N(1.2e+07, 1.3e+13) 12123769
 7 NAIVE(Scripts) 2009 Jan N(1.2e+07, 1.5e+13) 12123769
 8 NAIVE(Scripts) 2009 Feb N(1.2e+07, 1.8e+13) 12123769
 9 NAIVE(Scripts) 2009 Mar   N(1.2e+07, 2e+13) 12123769
10 NAIVE(Scripts) 2009 Apr N(1.2e+07, 2.2e+13) 12123769
# ℹ 86 more rows

The fable

A fable (forecast table) is like a tsibble, but with distributions for the forecasted variable.

Visualising forecasts

pbs_scripts |> 
  model(
    NAIVE(Scripts),
    SNAIVE(Scripts),
    NAIVE(Scripts ~ drift()),
    SNAIVE(Scripts ~ drift()),
  ) |> 
  forecast(h = "2 years") |> 
  autoplot(pbs_scripts)

The forecasting workflow

# Tidied / prepared data here, e.g.
vic_tourism |> 
  # Estimate the models using the data
  model(
    # Model specifications here, e.g.
    SNAIVE(Trips)
  ) |> 
  # Produce forecasts
  forecast(h = "2 years") |> 
  # Visualise forecasts
  autoplot(pbs_scripts)

Exercise 6

Specify appropriate model(s) for the total number of tourists arriving in Victoria.

Estimate them with model() and produce forecasts for the next 5 years with forecast().

Plot the forecasts, and visually evaluate their suitability.

Forecasting with regression

Regression can also forecast, use time as regressors for the response variable.

  • trend via the time index \(t\)
  • seasonality via seasonal dummies \(s_{it}\)

\[y_t = \beta_0 + \beta_1 t + \sum_{i=1}^m \beta_{s_i} s_{it} + \varepsilon_t\]

Specifying time series linear models

TSLM(y ~ trend() + season())

Forecasting with regression

Specifying time series linear models

TSLM(y ~ trend() + season())

Exercise 7

Produce forecasts for total Takings of Australian tourist accommodation from aus_accommodation using linear regression.

library(fpp3)
aus_accommodation_total <- fpp3::aus_accommodation |> 
  summarise(Takings = sum(Takings), Occupancy = mean(Occupancy))

Forecasting with regression

Regression can also incorporate other covariates. For example:

aus_accommodation_total |>
  model(
    TSLM(Takings ~ trend() + season() + Occupancy)
  )
# A mable: 1 x 1
  `TSLM(Takings ~ trend() + season() + Occupancy)`
                                           <model>
1                                           <TSLM>

Difficult to forecast

Forecasting with regressors requires you to forecast the regressors too! Useful regressors are easily known, predictable, or controllable.

Exponential smoothing

Time series patterns often change over time.

The fixed coefficients of regression is unsuitable for most time series.

Exponential smoothing (ETS) is similar to regression, but allows the ‘coefficients’ (states) to vary over time.

\[y_t = \ell_{t-1} + b_{t-1} + s_{t-m} + \varepsilon_t\]

Specifying ETS models

ETS(y ~ error("A") + trend("A") + season("A"))

Exponential smoothing

Exponential smoothing also works for ‘multiplicative’ patterns.

This is where the variation grows proportionately to \(\ell_t\).

For example, multiplicative seasonality:

\[y_t = (\ell_{t-1} + b_{t-1})s_{t-m} + \varepsilon_t\]

Specifying ETS models

ETS(y ~ error("A") + trend("A") + season("M"))

Exponential smoothing

The best ETS model can be chosen automatically using AICc.

Simply don’t specify the right side of the formula for automatic selection.

Specifying ETS models

ETS(y)

Exercise 8

Is the seasonality of total Australian accommodation takings from aus_accommodation_total additive or multiplicative?

Estimate an ETS model for the data, does the automatic ETS model match the patterns you see in a time plot?

ARIMA

Unlike earlier models which directly describe trends and seasonality, ARIMA forecasts patterns using autocorrelations (ACFs).

Recall the ACF plot

pbs_scripts |> 
  ACF(Scripts) |> 
  autoplot()

ARIMA

Unlike ETS, ARIMA can’t directly handle multiplicative patterns.

This isn’t a problem though, since we have a transformational trick 🪄

\[\log(x \times y) = \log(x) + \log(y)\]

Before we can use ARIMA, we need to transform our data to be additive.

ARIMA

Log transform the data to make multiplicative patterns additive.

pbs_scripts |> 
  autoplot(log(Scripts))

ARIMA

With the appropriate transformation, we can now estimate an ARIMA model:

Specifying ARIMA models

These will use automatic ARIMA selection (often described as auto.arima()).

ARIMA(y), or ARIMA(log(y))

Other transformations

Multiplicative patterns aren’t always exactly multiplicative - for this we often use power transformations via box_cox(y, lambda).

More information: https://otexts.com/fpp3/transformations.html

ARIMA

Exercise 9

Identify if a transformation is necessary for the aus_accommodation_total Takings

Then estimate an automatically selected ARIMA model for this data.

Compare ARIMA forecasts with the automatic ETS model, how do they differ?

ETS vs ARIMA

Which model is better?

It depends!

Capturing patterns

  • Both handle time varying trends and seasonality.
  • ETS directly captures multiplicative patterns.
  • ARIMA can forecast cyclical patterns.

Evaluating accuracy

We can determine which model works best on a specific dataset using accuracy evaluation…

Fortunately, that’s up next!

⏰ Time for a break!

Up next…

  • Looking for patterns in residuals
  • Evaluating forecast accuracy
  • Choosing the best forecasts

Useful links

social.mitchelloharawild.com

workshop.nectric.com.au/user2024/

mitchelloharawild/workshop-fable-user2024

Unsplash credits

Thanks to these Unsplash contributors for their photos