1
$\begingroup$

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.

$\endgroup$
1

1 Answer 1

1
$\begingroup$

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)
#> Lade nötiges Paket: fabletools
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)

$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – pcsksa5
    Jul 6, 2020 at 13:21
  • $\begingroup$ 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)) $\endgroup$
    – Tim-TU
    Jul 7, 2020 at 4:39
  • $\begingroup$ Thank you again. My main challenge is the standard deviation so that I can get samples from the forecasted distribution. $\endgroup$
    – pcsksa5
    Jul 7, 2020 at 18:27
  • $\begingroup$ 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? $\endgroup$
    – Tim-TU
    Jul 8, 2020 at 9:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.