tidyverts / fabletools

General fable features useful for extension packages
http://fabletools.tidyverts.org/
89 stars 31 forks source link

How to find the mu and sigma associated with an snaive fit #383

Closed alexhallam closed 12 months ago

alexhallam commented 12 months ago

I would like to know the mu and sigma that was estimated and also the mu and sigma in the forecast horizon given a snaive fit.

library(fable)
library(tsibble)
library(tsibbledata)
library(lubridate)
library(dplyr)
aus_retail %>%
  filter(
    State %in% c("New South Wales", "Victoria"),
    Industry == "Department stores"
  ) %>% 
  model(
    snaive = SNAIVE(Turnover)
  ) |>
  glance()
# A tibble: 2 × 4
  State           Industry          .model sigma2
  <chr>           <chr>             <chr>   <dbl>
1 New South Wales Department stores snaive   492.
2 Victoria        Department stores snaive   237.

I currently see sigma2, but I am not sure what that is referring to. If I calculated backward looking var I get

aus_retail %>%
  as_tibble() |>
  filter(
    State %in% c("New South Wales", "Victoria"),
    Industry == "Department stores"
  ) %>% 
  group_by(State) |>
  summarise(sigma2 = sd(Turnover))
  State           sigma2
  <chr>            <dbl>
1 New South Wales 23111.
2 Victoria        14336.

The sigma2 and var are not matching so I probably do not understand what is going on.

mitchelloharawild commented 12 months ago

sigma2 in the glance() output is the variance of the innovation residuals.

library(fable)
library(tsibble)
library(tsibbledata)
library(lubridate)
library(dplyr)
aus_retail %>%
  filter(
    State %in% c("New South Wales", "Victoria"),
    Industry == "Department stores"
  ) %>% 
  model(
    snaive = SNAIVE(Turnover)
  ) |>
  augment() |> 
  as_tibble() |> 
  group_by(State) |> 
  summarise(sigma2 = var(.innov, na.rm = TRUE))
#> # A tibble: 2 x 2
#>   State           sigma2
#>   <chr>            <dbl>
#> 1 New South Wales   492.
#> 2 Victoria          237.

The forecast distribution parameters can be obtained using distributional::parameters():

aus_retail %>%
  filter(
    State %in% c("New South Wales", "Victoria"),
    Industry == "Department stores"
  ) %>% 
  model(
    snaive = SNAIVE(Turnover)
  ) |>
  forecast() |> 
  mutate(distributional::parameters(Turnover))
#> # A fable: 48 x 8 [1M]
#> # Key:     State, Industry, .model [2]
#>    State           Industry        .model    Month    Turnover .mean    mu sigma
#>    <chr>           <chr>           <chr>     <mth>      <dist> <dbl> <dbl> <dbl>
#>  1 New South Wales Department sto~ snaive 2019 Jan N(452, 570)  452.  452.  23.9
#>  2 New South Wales Department sto~ snaive 2019 Feb N(350, 570)  350.  350.  23.9
#>  3 New South Wales Department sto~ snaive 2019 Mar N(457, 570)  457   457   23.9
#>  4 New South Wales Department sto~ snaive 2019 Apr N(454, 570)  454.  454.  23.9
#>  5 New South Wales Department sto~ snaive 2019 May N(502, 570)  502.  502.  23.9
#>  6 New South Wales Department sto~ snaive 2019 Jun N(531, 570)  531.  531.  23.9
#>  7 New South Wales Department sto~ snaive 2019 Jul N(446, 570)  446.  446.  23.9
#>  8 New South Wales Department sto~ snaive 2019 Aug N(429, 570)  429.  429.  23.9
#>  9 New South Wales Department sto~ snaive 2019 Sep N(462, 570)  462.  462.  23.9
#> 10 New South Wales Department sto~ snaive 2019 Oct N(503, 570)  503.  503.  23.9
#> # i 38 more rows

Similar can be done for the fitted values, obtained with the fitted() function.

Created on 2023-07-07 with reprex v2.0.2