library(tsibble)
library(dplyr)
library(fable)
library(fabletools)
library(vctrs)
library(purrr)
library(distributional)Calculating growth rates from fable forecasts
Intro
This is an attempt to calculate growth rates from fable forecasts of a volume variable. Data was the tourism dataset from tsibble package.
Libraries loaded
Libs needed for this:
Base fit and forecast
tsb_test <- tourism |>
filter(Region %in% c('Adelaide') ) |>
filter(Purpose %in% c('Business', 'Holiday') )
fit_test <- tsb_test |>
model( ets = ETS(log(Trips) ) )
fc_test <- fit_test |>
forecast(h=3)Base autoplot
autoplot(fc_test)
Experimental: growth rate calculation from forecast
na_distribution <- dist_transformed( dist_normal(NA,NA),
transform = function(Trips) exp(Trips),
inverse = function(Trips) log(Trips))
Trips_lagged <- fc_test |>
group_by_key() |>
mutate ( dplyr::lag(Trips, default = na_distribution), .keep = "none") |>
pull()
pull_first_element <- function(x) {x[[1]]}
y_fc <- fc_test$Trips |>
vec_data() |>
lapply (pull_first_element) |>
vec_restore(fc_test$Trips )
y_lagged <- Trips_lagged |>
vec_data() |>
lapply (pull_first_element) |>
vec_restore(fc_test$Trips )
y_diff <- map2_vec(y_fc, y_lagged, `-` )
y_diff_transformed <- dist_transformed( y_diff,
transform = function(Trips) exp(Trips),
inverse = function(Trips) log(Trips) )
fc_growth <- fc_test |>
mutate(Trips = y_diff_transformed ) |>
mutate(.mean = mean(Trips))The echo: false option disables the printing of code (only output is displayed).
Autoplot of growth rates of the forecast
autoplot(fc_growth)Warning in stats::rnorm(times, x[["mu"]], x[["sigma"]]): wyprodukowano wartości
NA
Warning in stats::rnorm(times, x[["mu"]], x[["sigma"]]): wyprodukowano wartości
NA
Warning in stats::rnorm(times, x[["mu"]], x[["sigma"]]): wyprodukowano wartości
NA
Warning in stats::rnorm(times, x[["mu"]], x[["sigma"]]): wyprodukowano wartości
NA
Warning in stats::rnorm(times, x[["mu"]], x[["sigma"]]): wyprodukowano wartości
NA
Warning in stats::rnorm(times, x[["mu"]], x[["sigma"]]): wyprodukowano wartości
NA
Warning in stats::rnorm(times, x[["mu"]], x[["sigma"]]): wyprodukowano wartości
NA
Warning in stats::rnorm(times, x[["mu"]], x[["sigma"]]): wyprodukowano wartości
NAWarning: Removed 2 rows containing missing values or values outside the scale range
(`geom_line()`).
Autoplot of growth rates of forecasts excluding NAs
fc_growth_not_na <- fc_growth |>
filter(!is.na(.mean))
autoplot(fc_growth_not_na)
Forecast data table
fc_test |> knitr::kable()| Region | State | Purpose | .model | Quarter | Trips | .mean |
|---|---|---|---|---|---|---|
| Adelaide | South Australia | Business | ets | 2018 Q1 | t(N(5, 0.05)) | 147.6086 |
| Adelaide | South Australia | Business | ets | 2018 Q2 | t(N(5.1, 0.054)) | 175.5875 |
| Adelaide | South Australia | Business | ets | 2018 Q3 | t(N(5.2, 0.057)) | 191.2471 |
| Adelaide | South Australia | Holiday | ets | 2018 Q1 | t(N(5.3, 0.019)) | 212.5453 |
| Adelaide | South Australia | Holiday | ets | 2018 Q2 | t(N(5.1, 0.02)) | 169.4951 |
| Adelaide | South Australia | Holiday | ets | 2018 Q3 | t(N(5.1, 0.021)) | 166.4275 |
fc_growth |> knitr::kable()| Region | State | Purpose | .model | Quarter | Trips | .mean |
|---|---|---|---|---|---|---|
| Adelaide | South Australia | Business | ets | 2018 Q1 | t(N(NA, NA)) | NA |
| Adelaide | South Australia | Business | ets | 2018 Q2 | t(N(0.17, 0.1)) | 1.2492463 |
| Adelaide | South Australia | Business | ets | 2018 Q3 | t(N(0.084, 0.11)) | 1.1484114 |
| Adelaide | South Australia | Holiday | ets | 2018 Q1 | t(N(NA, NA)) | NA |
| Adelaide | South Australia | Holiday | ets | 2018 Q2 | t(N(-0.23, 0.04)) | 0.8129706 |
| Adelaide | South Australia | Holiday | ets | 2018 Q3 | t(N(-0.019, 0.041)) | 1.0016028 |
Lognormal distribution
y_diff_exp <- y_diff |> exp()
fc_growth_exp <- fc_test |>
mutate(Trips = y_diff_exp ) |>
mutate(.mean = mean(Trips))
autoplot(fc_growth_exp)Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_line()`).
fc_growth_exp |> knitr::kable()| Region | State | Purpose | .model | Quarter | Trips | .mean |
|---|---|---|---|---|---|---|
| Adelaide | South Australia | Business | ets | 2018 Q1 | lN(NA, NA) | NA |
| Adelaide | South Australia | Business | ets | 2018 Q2 | lN(0.17, 0.1) | 1.2509001 |
| Adelaide | South Australia | Business | ets | 2018 Q3 | lN(0.084, 0.11) | 1.1501272 |
| Adelaide | South Australia | Holiday | ets | 2018 Q1 | lN(NA, NA) | NA |
| Adelaide | South Australia | Holiday | ets | 2018 Q2 | lN(-0.23, 0.04) | 0.8131273 |
| Adelaide | South Australia | Holiday | ets | 2018 Q3 | lN(-0.019, 0.041) | 1.0018079 |
Check growth rates by hand
Mean is calculated from formula:
\[ mean = \exp \left(\mu +{\frac {\sigma ^{2}}{2}}\right) \]
mus <- sapply( y_diff |> vec_data() , \(x) x[["mu"]] )
sigmas <- sapply( y_diff |> vec_data() , \(x) x[["sigma"]] )
growth_by_hand <- mus + (sigmas ^2)/2
growth_by_hand[1] NA 0.223863368 0.139872576 NA -0.206867547
[6] 0.001806264exp_growth_by_hand <- exp(growth_by_hand)
exp_growth_by_hand[1] NA 1.2509001 1.1501272 NA 0.8131273 1.0018079fc_augment <- fit_test |>
augment() |>
select(-any_of(c(".innov", ".fitted", ".resid"))) |>
mutate(Trips = dist_lognormal( log(Trips), 0)) |>
mutate(.mean = mean(Trips))
dimnames(fc_augment$Trips) <- "Trips"
fc_augment <- fc_augment |>
as_fable(response = "Trips", distribution = "Trips")
y_as_lnorm <- fc_test |> pull(Trips) |> vec_data() |> lapply (\(x) x[[1]] ) |> vec_restore( fc_test$Trips ) |> exp()
fc_as_lnorm <- fc_test |>
mutate( Trips = y_as_lnorm ) |>
mutate( .mean = mean(Trips) )
fc_all <- dplyr::bind_rows (fc_augment, fc_as_lnorm)autoplot(fc_all, level=90)
fc_augment |> tail() |> knitr::kable()| Region | State | Purpose | .model | Quarter | Trips | .mean |
|---|---|---|---|---|---|---|
| Adelaide | South Australia | Holiday | ets | 2016 Q3 | lN(5, 0) | 154.1169 |
| Adelaide | South Australia | Holiday | ets | 2016 Q4 | lN(5, 0) | 147.1020 |
| Adelaide | South Australia | Holiday | ets | 2017 Q1 | lN(5.4, 0) | 215.9403 |
| Adelaide | South Australia | Holiday | ets | 2017 Q2 | lN(5.3, 0) | 201.7341 |
| Adelaide | South Australia | Holiday | ets | 2017 Q3 | lN(5.1, 0) | 163.7119 |
| Adelaide | South Australia | Holiday | ets | 2017 Q4 | lN(5.4, 0) | 214.0531 |
Trips_lagged <- fc_all |>
group_by_key() |>
mutate ( dplyr::lag(Trips, default = dist_lognormal(NA,NA)), .keep = "none") |>
pull()
fc_diff <- fc_all |>
mutate(Trips = exp( log(Trips) - log(Trips_lagged) )) |>
mutate(.mean = mean(Trips))autoplot(fc_diff, level=90)Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_line()`).
fc_diff |> tsibble::group_by_key() |> slice_tail(n=6) |> knitr::kable()| Region | State | Purpose | .model | Quarter | Trips | .mean |
|---|---|---|---|---|---|---|
| Adelaide | South Australia | Business | ets | 2017 Q2 | lN(0.3, 0) | 1.3464928 |
| Adelaide | South Australia | Business | ets | 2017 Q3 | lN(0.064, 0) | 1.0665732 |
| Adelaide | South Australia | Business | ets | 2017 Q4 | lN(0.063, 0) | 1.0646371 |
| Adelaide | South Australia | Business | ets | 2018 Q1 | lN(-0.31, 0.05) | 0.7484511 |
| Adelaide | South Australia | Business | ets | 2018 Q2 | lN(0.17, 0.1) | 1.2509001 |
| Adelaide | South Australia | Business | ets | 2018 Q3 | lN(0.084, 0.11) | 1.1501272 |
| Adelaide | South Australia | Holiday | ets | 2017 Q2 | lN(-0.068, 0) | 0.9342122 |
| Adelaide | South Australia | Holiday | ets | 2017 Q3 | lN(-0.21, 0) | 0.8115231 |
| Adelaide | South Australia | Holiday | ets | 2017 Q4 | lN(0.27, 0) | 1.3074985 |
| Adelaide | South Australia | Holiday | ets | 2018 Q1 | lN(-0.017, 0.019) | 0.9930028 |
| Adelaide | South Australia | Holiday | ets | 2018 Q2 | lN(-0.23, 0.04) | 0.8131273 |
| Adelaide | South Australia | Holiday | ets | 2018 Q3 | lN(-0.019, 0.041) | 1.0018079 |