# Extract the distribution parameters from a forecast that produces a fable

I have used the forecast(h="1 year") function from the fabletools package to produce a forecast from an ARIMA(1,1,3)(1,1,0)[52] model where the dataset had been transformed. A snapshot of the forecast output is the following:

# A fable: 52 x 5 [1W]
# Key:     LotGroup, .model [1]
LotGroup .model     week              Rate .mean
<chr>    <chr>    <week>            <dist> <dbl>
1 123      a_n    2019 W45 t(N(-0.78, 0.23)) 32.5
2 123      a_n    2019 W46 t(N(-0.71, 0.26)) 33.8
3 123      a_n    2019 W47 t(N(-0.53, 0.26)) 37.8
4 123      a_n    2019 W48 t(N(-0.98, 0.27)) 28.6
5 123      a_n    2019 W49 t(N(-0.96, 0.28)) 28.9
6 123      a_n    2019 W50  t(N(-0.85, 0.3)) 31.2
7 123      a_n    2019 W51  t(N(-1.2, 0.31)) 25.0


What I would like to do is extract the parameter elements of the weekly distributions ie the c(-0.78, -0.71, ...) and the c(0.23, 0.26, ...). Using the hillo() function produces PIs but does not give the parameters.

I want the parameters in order to use them in a simulation and sample from the non-transformed distribution. The samples will then be inverted back and used for further calculations.

Any ideas on how to extract the future parameters?

Thank you.

• The distributions in the fable are from the distributional package: pkg.mitchelloharawild.com/distributional. You can sample (and have the samples inverted) from the distribution using the generate() function. Jul 6, 2020 at 1:46

Here is a working example for the extraction of the desired distributional components. Bear in mind that this is the distribution on a transformed scale. Maybe this link gives some more explanation on this topic: Forecasting using transformations

library(tsibbledata)
library(tsibble)
library(dplyr)
#>
#> Attache Paket: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#>     filter, lag
#> The following objects are masked from 'package:base':
#>
#>     intersect, setdiff, setequal, union
library(fable)
library(feasts)

elec_week <- vic_elec %>%
tsibble::index_by(week = tsibble::yearweek(Date)) %>%
dplyr::group_by(week) %>%
dplyr::summarize(Demand = sum(Demand))

bc_lambda <- feasts::guerrero(elec_week$Demand) mod <- elec_week %>% model(ARIMA(box_cox(Demand, bc_lambda) ~ 0 + pdq(1, 1, 3) + PDQ(1, 1, 0, period = 52))) fc <- mod %>% forecast(h = "1 year") fc #> # A fable: 52 x 4 [1W] #> # Key: .model [1] #> .model week Demand .mean #> <chr> <week> <dist> <dbl> #> 1 "ARIMA(box_cox(Demand, bc_lambda) ~ 0~ 2015 W02 t(N(9.7e+11, 4e+22)) 1.39e6 #> 2 "ARIMA(box_cox(Demand, bc_lambda) ~ 0~ 2015 W03 t(N(1.9e+12, 4.1e+22)) 1.93e6 #> 3 "ARIMA(box_cox(Demand, bc_lambda) ~ 0~ 2015 W04 t(N(1.2e+12, 4.2e+22)) 1.51e6 #> 4 "ARIMA(box_cox(Demand, bc_lambda) ~ 0~ 2015 W05 t(N(1.5e+12, 4.2e+22)) 1.73e6 #> 5 "ARIMA(box_cox(Demand, bc_lambda) ~ 0~ 2015 W06 t(N(1.6e+12, 4.2e+22)) 1.80e6 #> 6 "ARIMA(box_cox(Demand, bc_lambda) ~ 0~ 2015 W07 t(N(1.2e+12, 4.2e+22)) 1.55e6 #> 7 "ARIMA(box_cox(Demand, bc_lambda) ~ 0~ 2015 W08 t(N(1.1e+12, 4.2e+22)) 1.48e6 #> 8 "ARIMA(box_cox(Demand, bc_lambda) ~ 0~ 2015 W09 t(N(1e+12, 4.2e+22)) 1.44e6 #> 9 "ARIMA(box_cox(Demand, bc_lambda) ~ 0~ 2015 W10 t(N(1.2e+12, 4.2e+22)) 1.56e6 #> 10 "ARIMA(box_cox(Demand, bc_lambda) ~ 0~ 2015 W11 t(N(1.1e+12, 4.2e+22)) 1.44e6 #> # ... with 42 more rows # attributes: fc$$Demand[[1]] %>% attributes() #>$$names #> [1] "dist" "transform" "inverse" #> #> $$class #> [1] "dist_transformed" "dist_default" fc$$Demand[[1]]$$dist %>% attributes() #>$$names #> [1] "mu" "sigma" #> #>$class
#> [1] "dist_normal"  "dist_default"

# parameter extraction:
pars <- lapply(fc\$Demand, "[[", "dist")
mus <- sapply(pars, "[[", "mu")
sigmas <- sapply(pars, "[[", "sigma")

cbind(mus, sigmas)
#>                mus       sigmas
#>  [1,] 9.707294e+11 200309806003
#>  [2,] 1.872558e+12 203537898280
#>  [3,] 1.154797e+12 204131493254
#>  [4,] 1.504034e+12 204362403497
#>  [5,] 1.619694e+12 204823657928
#>  [6,] 1.209110e+12 204979142526
#>  [7,] 1.099486e+12 205339457957
#>  [8,] 1.043497e+12 205441028956
#>  [9,] 1.230285e+12 205723983361
#> [10,] 1.052427e+12 205787298940
#> [11,] 1.067510e+12 206010684952
#> [12,] 9.921512e+11 206047129343
#> [13,] 1.174856e+12 206224434053
#> [14,] 1.045946e+12 206242263253
#> [15,] 9.852037e+11 206383752697
#> [16,] 9.347329e+11 206388930518
#> [17,] 1.198340e+12 206502450303
#> [18,] 1.214986e+12 206499261092
#> [19,] 1.098114e+12 206590830908
#> [20,] 1.081210e+12 206582329842
#> [21,] 1.190340e+12 206656587029
#> [22,] 1.184583e+12 206644931572
#> [23,] 1.235165e+12 206705463367
#> [24,] 1.307883e+12 206692157392
#> [25,] 1.391260e+12 206741750660
#> [26,] 1.366485e+12 206727824195
#> [27,] 1.414352e+12 206768653706
#> [28,] 1.443564e+12 206754795155
#> [29,] 1.502585e+12 206788565941
#> [30,] 1.363670e+12 206775219170
#> [31,] 1.382080e+12 206803274467
#> [32,] 1.400818e+12 206790709906
#> [33,] 1.294506e+12 206814113223
#> [34,] 1.167685e+12 206802479778
#> [35,] 1.202383e+12 206822077472
#> [36,] 1.090755e+12 206811440257
#> [37,] 1.176906e+12 206827909398
#> [38,] 1.003286e+12 206818277011
#> [39,] 1.004477e+12 206832162095
#> [40,] 1.024545e+12 206823506205
#> [41,] 1.079534e+12 206835247423
#> [42,] 1.086830e+12 206827516715
#> [43,] 1.053291e+12 206837471789
#> [44,] 9.669133e+11 206830601778
#> [45,] 1.073288e+12 206839062898
#> [46,] 1.037378e+12 206832982734
#> [47,] 1.075157e+12 206840189746
#> [48,] 1.150829e+12 206834826842
#> [49,] 1.079385e+12 206840977566
#> [50,] 1.083486e+12 206836260661
#> [51,] 8.613049e+11 206841518994
#> [52,] 1.664814e+11 206837380086


Created on 2020-07-05 by the reprex package (v0.3.0)

• Hello. Thank you for your reply. From the link that you kindly sent I assume that the extracted parameters are not yet bias-adjusted. Is there a way to get the bias-adjusted parameters? Again, thank you for your time. Jul 6, 2020 at 13:21
• Do you mean the forecast? As it is stated in the online book: "When we use the mean, rather than the median, we say the point forecasts have been bias-adjusted.". So the mean is the default in the function forecast. But you can also get the median. Try forecast(h = "1 year", point_forecast = list(.mean = mean, .median = median)) Jul 7, 2020 at 4:39
• Thank you again. My main challenge is the standard deviation so that I can get samples from the forecasted distribution. Jul 7, 2020 at 18:27
• Forecasts are made by t(N(mu, sigma)) but the point forecasts are backtransformed. @Mitchell O'Hara-Wild: Is there a handy function in the tidyverts framework or distributional package to extract t(N(mu, sigma)) to draw samples? Jul 8, 2020 at 9:46