library(fabletools)
library(fable)
library(tsibble)
library(dplyr)

Models built with fabletools (such as those in the fable package) are specified with a model formula. Much like lm(), the formula has two sides:

MODEL(response ~ specials + regressors)

The left hand side (LHS) says what is being modelled: the response variable, any transformations of it, and (for multivariate models) several responses.
The right hand side (RHS) says how it is modelled: model-specific specials (like pdq() for ARIMA()), deterministic terms (like trend() and season()), and exogenous regressors from the data.

Both sides are optional. If the RHS is omitted, the model uses its default (often automatically selected) specification. If the response is omitted and the data has exactly one measured variable, it is used as the response.

lung_deaths <- as_tsibble(cbind(mdeaths, fdeaths), pivot_longer = TRUE)
lung_deaths |>
  model(
    auto = ETS(value),
    manual = ETS(value ~ error("A") + trend("N") + season("A"))
  )
#> # A mable: 2 x 3
#> # Key:     key [2]
#>   key             auto       manual
#>   <chr>        <model>      <model>
#> 1 fdeaths <ETS(M,N,M)> <ETS(A,N,A)>
#> 2 mdeaths <ETS(M,A,A)> <ETS(A,N,A)>

The left hand side

Response variables

The simplest LHS is a variable from the data, which is the variable being forecast. Forecasts, fitted values and residuals are all given on the scale of this variable.

Transformations

Transformations can be written directly in the formula, using the response variable in an expression.

fit <- lung_deaths |>
  model(ETS(log(value)))

fabletools works out the inverse of the transformation, so the model is estimated on the transformed scale and the results are automatically back-transformed to the original scale of the response. Forecasts are distributions of the back-transformed response, shown as t(<distribution>). The point forecast .mean is the mean of this distribution, which is bias-adjusted.

fit |> forecast(h = 3)
#> # A fable: 6 x 5 [1M]
#> # Key:     key, .model [2]
#>   key     .model             index             value .mean
#>   <chr>   <chr>              <mth>            <dist> <dbl>
#> 1 fdeaths ETS(log(value)) 1980 Jan  t(N(6.7, 0.013))  809.
#> 2 fdeaths ETS(log(value)) 1980 Feb  t(N(6.7, 0.013))  792.
#> 3 fdeaths ETS(log(value)) 1980 Mar  t(N(6.6, 0.012))  747.
#> 4 mdeaths ETS(log(value)) 1980 Jan t(N(7.5, 0.0095)) 1832.
#> 5 mdeaths ETS(log(value)) 1980 Feb t(N(7.5, 0.0095)) 1854.
#> 6 mdeaths ETS(log(value)) 1980 Mar t(N(7.5, 0.0094)) 1732.

Transformations can be combined, for example log(value + 1) or box_cox(value, 0.3) + 10. Supported functions include arithmetic (+, -, *, /, ^), log() and its variants, exp(), sqrt() and box_cox(). Other transformations can be defined with new_transformation(), which pairs a function with its inverse.

Identifying the response with `resp()`

A transformation can involve other values besides the response, such as the 0.3 in box_cox(value, 0.3). fabletools needs to know which of these is the response variable, so it can be back-transformed.

By default, the response is identified as the only argument that is as long as the data. Other arguments must be length 1 (constant), such as 0.3 above. Any constant arguments are saved with the model and reused for back-transforming fitted values and forecasts.

When a transformation uses more than one variable from the data, it is ambiguous which is the response. For example, to forecast GDP per capita, GDP / Population could be a model of either variable. When this happens, the whole expression is treated as the response variable (in this case named GDP/Population), and it isn’t back-transformed.

Wrap the response variable in resp() to remove this ambiguity. The other data variables are then treated as time-varying parameters of the transformation.

aus_economy <- tsibbledata::global_economy |>
  filter(Code == "AUS")
fit <- aus_economy |>
  model(ETS(resp(GDP) / Population))
fit
#> # A mable: 1 x 2
#> # Key:     Country [1]
#>   Country   `ETS(resp(GDP)/Population)`
#>   <fct>                         <model>
#> 1 Australia                <ETS(M,A,N)>

The model is estimated on GDP per capita, while fitted values and forecasts are back-transformed to GDP. Since Population changes over time, its future values are needed to forecast GDP. These are given with the new_data argument of forecast():

future_population <- new_data(aus_economy, 3) |>
  mutate(Population = last(aus_economy$Population) * 1.015^(1:3))
fit |> forecast(new_data = future_population)
#> # A fable: 3 x 6 [1Y]
#> # Key:     Country, .model [1]
#>   Country   .model                  Year                  GDP   .mean Population
#>   <fct>     <chr>                  <dbl>               <dist>   <dbl>      <dbl>
#> 1 Australia ETS(resp(GDP)/Populat…  2018 t(N(53331, 3.1e+07)) 1.33e12  24967917.
#> 2 Australia ETS(resp(GDP)/Populat…  2019 t(N(52863, 7.9e+07)) 1.34e12  25342436.
#> 3 Australia ETS(resp(GDP)/Populat…  2020 t(N(52395, 1.5e+08)) 1.35e12  25722572.

Transformation parameters for each series

A common case of this is a transformation parameter that differs between series. The Guerrero method from the feasts package gives a Box-Cox \(\lambda\) for each series. After joining these to the data, the lambda_guerrero column is as long as Count, so resp() or a length-1 parameter is needed.

livestock <- tsibbledata::aus_livestock |>
  filter(State == "Victoria", Animal %in% c("Cows and heifers", "Lambs", "Pigs"))
lambda <- livestock |>
  features(Count, features = feasts::guerrero)
lambda
#> # A tibble: 3 × 3
#>   Animal           State    lambda_guerrero
#>   <fct>            <fct>              <dbl>
#> 1 Cows and heifers Victoria           0.551
#> 2 Lambs            Victoria           0.917
#> 3 Pigs             Victoria           0.247
livestock <- livestock |>
  left_join(lambda, by = c("Animal", "State"))

There are two ways to specify this model, depending on whether the parameter should vary over time.

A time-varying parameter: identify the response with resp(). The parameter’s future values must be given in new_data when forecasting.

fit <- livestock |>
  model(ETS(box_cox(resp(Count), lambda_guerrero)))
future <- new_data(livestock, 12) |>
  left_join(lambda, by = c("Animal", "State"))
fit |> forecast(new_data = future)
#> # A fable: 36 x 7 [1M]
#> # Key:     Animal, State, .model [3]
#>    Animal          State .model    Month            Count  .mean lambda_guerrero
#>    <fct>           <fct> <chr>     <mth>           <dist>  <dbl>           <dbl>
#>  1 Cows and heife… Vict… ETS(b… 2019 Jan  t(N(971, 3587)) 90506.           0.551
#>  2 Cows and heife… Vict… ETS(b… 2019 Feb  t(N(993, 4681)) 94333.           0.551
#>  3 Cows and heife… Vict… ETS(b… 2019 Mar  t(N(989, 5774)) 93704.           0.551
#>  4 Cows and heife… Vict… ETS(b… 2019 Apr  t(N(969, 6868)) 90406.           0.551
#>  5 Cows and heife… Vict… ETS(b… 2019 May t(N(1023, 7962)) 99661.           0.551
#>  6 Cows and heife… Vict… ETS(b… 2019 Jun  t(N(966, 9055)) 89935.           0.551
#>  7 Cows and heife… Vict… ETS(b… 2019 Jul t(N(959, 10149)) 88968.           0.551
#>  8 Cows and heife… Vict… ETS(b… 2019 Aug t(N(935, 11243)) 85004.           0.551
#>  9 Cows and heife… Vict… ETS(b… 2019 Sep t(N(934, 12336)) 84882.           0.551
#> 10 Cows and heife… Vict… ETS(b… 2019 Oct t(N(994, 13430)) 95081.           0.551
#> # ℹ 26 more rows

A time-invariant parameter: reduce the parameter to length 1, for example with first(). Models are estimated for each series separately, so each series uses its own \(\lambda\). The value is saved with the model, so new_data isn’t needed.

fit <- livestock |>
  model(ETS(box_cox(Count, first(lambda_guerrero))))
fit |> forecast(h = 12)
#> # A fable: 36 x 6 [1M]
#> # Key:     Animal, State, .model [3]
#>    Animal           State    .model                Month            Count  .mean
#>    <fct>            <fct>    <chr>                 <mth>           <dist>  <dbl>
#>  1 Cows and heifers Victoria ETS(box_cox(Count… 2019 Jan  t(N(971, 3587)) 90506.
#>  2 Cows and heifers Victoria ETS(box_cox(Count… 2019 Feb  t(N(993, 4681)) 94333.
#>  3 Cows and heifers Victoria ETS(box_cox(Count… 2019 Mar  t(N(989, 5774)) 93704.
#>  4 Cows and heifers Victoria ETS(box_cox(Count… 2019 Apr  t(N(969, 6868)) 90406.
#>  5 Cows and heifers Victoria ETS(box_cox(Count… 2019 May t(N(1023, 7962)) 99661.
#>  6 Cows and heifers Victoria ETS(box_cox(Count… 2019 Jun  t(N(966, 9055)) 89935.
#>  7 Cows and heifers Victoria ETS(box_cox(Count… 2019 Jul t(N(959, 10149)) 88968.
#>  8 Cows and heifers Victoria ETS(box_cox(Count… 2019 Aug t(N(935, 11243)) 85004.
#>  9 Cows and heifers Victoria ETS(box_cox(Count… 2019 Sep t(N(934, 12336)) 84882.
#> 10 Cows and heifers Victoria ETS(box_cox(Count… 2019 Oct t(N(994, 13430)) 95081.
#> # ℹ 26 more rows

The second approach is usually simpler when the parameter is constant within each series.

Multivariate responses with `vars()`

Multivariate models, such as VAR(), model several response variables jointly. These are given on the LHS with vars(), and each response can have its own transformation.

lung_deaths_wide <- as_tsibble(cbind(mdeaths, fdeaths), pivot_longer = FALSE)
fit <- lung_deaths_wide |>
  model(var = VAR(vars(log(mdeaths), fdeaths)))
fit
#> # A mable: 1 x 1
#>        var
#>    <model>
#> 1 <VAR(5)>

Each response is back-transformed separately. Multivariate forecast distributions can’t yet be back-transformed directly, so forecasts from a transformed multivariate model are produced from simulated sample paths using simulate_iid().

fit |>
  mutate(var = simulate_iid(var)) |>
  forecast(h = 3)
#> # A fable: 3 x 4 [1M]
#> # Key:     .model [1]
#>   .model    index .distribution .mean[,"mdeaths"] [,"fdeaths"]
#>   <chr>     <mth>        <dist>             <dbl>        <dbl>
#> 1 var    1980 Jan  sample[5000]             1619.         696.
#> 2 var    1980 Feb  sample[5000]             1778.         751.
#> 3 var    1980 Mar  sample[5000]             1668.         696.

The right hand side

Specials

Specials are functions used in the formula to set the structure of a model. Each model has its own specials, which are documented in the “Specials” section of the model’s help file (for example, ?ETS and ?ARIMA).

lung_deaths |>
  model(
    ets = ETS(value ~ error("M") + trend("Ad") + season("M")),
    arima = ARIMA(value ~ 1 + pdq(1, 0, 0) + PDQ(0, 1, 1))
  )
#> # A mable: 2 x 3
#> # Key:     key [2]
#>   key               ets                              arima
#>   <chr>         <model>                            <model>
#> 1 fdeaths <ETS(M,Ad,M)> <ARIMA(1,0,0)(0,1,1)[12] w/ drift>
#> 2 mdeaths <ETS(M,Ad,M)> <ARIMA(1,0,0)(0,1,1)[12] w/ drift>

A special that isn’t specified uses its default. Many specials search over a set of options by default, so a partially specified formula like ARIMA(value ~ PDQ(0, 1, 1)) fixes the seasonal part and automatically selects the rest. Some models also use the intercept (1 or 0) as part of their specification, such as whether ARIMA() includes a constant.

Deterministic terms

Models that support regressors (such as TSLM() and ARIMA()) can use deterministic terms of time, which are computed from the time index:

trend(knots = NULL, origin = NULL) gives a linear trend, which is piecewise linear if knots are given.
season(period = NULL) gives seasonal dummy variables.
fourier(period = NULL, K) gives Fourier terms for smooth seasonal patterns.

lung_deaths |>
  model(
    dummy = TSLM(value ~ trend() + season()),
    harmonic = TSLM(value ~ trend() + fourier(K = 2))
  ) |>
  forecast(h = 12)
#> # A fable: 48 x 5 [1M]
#> # Key:     key, .model [4]
#>    key     .model    index
#>    <chr>   <chr>     <mth>
#>  1 fdeaths dummy  1980 Jan
#>  2 fdeaths dummy  1980 Feb
#>  3 fdeaths dummy  1980 Mar
#>  4 fdeaths dummy  1980 Apr
#>  5 fdeaths dummy  1980 May
#>  6 fdeaths dummy  1980 Jun
#>  7 fdeaths dummy  1980 Jul
#>  8 fdeaths dummy  1980 Aug
#>  9 fdeaths dummy  1980 Sep
#> 10 fdeaths dummy  1980 Oct
#> # ℹ 38 more rows
#> # ℹ 2 more variables: value <dist>, .mean <dbl>

These terms are computed from the time index, so they can be forecast without extra data.

Exogenous regressors

Other variables from the data can be used as regressors by adding them to the RHS, in the same way as in lm(). This includes transformations of regressors (log(x)), polynomials (I(x^2)) and interactions (x:z). The intercept is controlled with 1 and 0 (or -1), and . uses all other measured variables in the data.

elec <- tsibbledata::vic_elec |>
  index_by(Date = as.Date(Time)) |>
  summarise(
    Demand = sum(Demand) / 1e3,
    Temperature = max(Temperature),
    Holiday = any(Holiday)
  ) |>
  mutate(Workday = !Holiday & !(lubridate::wday(Date) %in% c(1, 7)))

fit <- elec |>
  model(
    ARIMA(Demand ~ Temperature + I(Temperature^2) + Workday)
  )
fit
#> # A mable: 1 x 1
#>   `ARIMA(Demand ~ Temperature + I(Temperature^2) + Workday)`
#>                                                      <model>
#> 1                      <LM w/ ARIMA(0,1,4)(2,0,0)[7] errors>

Regressors given on the RHS are passed to the model’s xreg() special, which can also be used directly. Some models give xreg() extra options; for example, ARIMA(y ~ xreg(x, fixed = list(x = 1))) fixes a coefficient’s value.

Unlike deterministic terms, the future values of exogenous regressors must be given to forecast() with new_data. This also allows forecasts under different scenarios.

future_elec <- new_data(elec, 7) |>
  mutate(
    Temperature = c(26, 30, 34, 27, 25, 24, 26),
    Workday = !(lubridate::wday(Date) %in% c(1, 7))
  )
fit |>
  forecast(new_data = future_elec)
#> # A fable: 7 x 6 [1D]
#> # Key:     .model [1]
#>   .model                                                   Date      
#>   <chr>                                                    <date>    
#> 1 ARIMA(Demand ~ Temperature + I(Temperature^2) + Workday) 2015-01-01
#> 2 ARIMA(Demand ~ Temperature + I(Temperature^2) + Workday) 2015-01-02
#> 3 ARIMA(Demand ~ Temperature + I(Temperature^2) + Workday) 2015-01-03
#> 4 ARIMA(Demand ~ Temperature + I(Temperature^2) + Workday) 2015-01-04
#> 5 ARIMA(Demand ~ Temperature + I(Temperature^2) + Workday) 2015-01-05
#> 6 ARIMA(Demand ~ Temperature + I(Temperature^2) + Workday) 2015-01-06
#> 7 ARIMA(Demand ~ Temperature + I(Temperature^2) + Workday) 2015-01-07
#> # ℹ 4 more variables: Demand <dist>, .mean <dbl>, Temperature <dbl>,
#> #   Workday <lgl>

Summary

Formula	Meaning
`y`	Model `y` with the default specification
`log(y)`	Model `log(y)`, back-transformed to `y`
`box_cox(y, 0.3)`	A transformation with a constant parameter
`box_cox(y, first(lambda))`	A constant parameter from the data (one per series)
`box_cox(resp(y), lambda)`	A time-varying parameter from the data, `lambda` needed in `new_data`
`resp(y) / x`	The response is `y`, transformed using another variable `x`
`vars(y, log(z))`	A multivariate model of `y` and `z`
`y ~ pdq(1, 0, 0)`	Model-specific specials
`y ~ trend() + season()`	Deterministic terms, no `new_data` needed
`y ~ x + I(x^2)`	Exogenous regressors, `x` needed in `new_data`

The model formula