fable::ISF2020

class: inverse
background-image: linear-gradient(to right, rgba(150, 150, 150, .1), rgba(150, 150, 150, .4)), url("resources/hourglass.jpg")
background-size: cover

## Probabilistic cross-temporal hierarchies

.bottom[
### Mitchell O'Hara-Wild (<svg style="height:0.8em;top:.04em;position:relative;fill:#1da1f2;" viewBox="0 0 512 512"><path d="M459.37 151.716c.325 4.548.325 9.097.325 13.645 0 138.72-105.583 298.558-298.558 298.558-59.452 0-114.68-17.219-161.137-47.106 8.447.974 16.568 1.299 25.34 1.299 49.055 0 94.213-16.568 130.274-44.832-46.132-.975-84.792-31.188-98.112-72.772 6.498.974 12.995 1.624 19.818 1.624 9.421 0 18.843-1.3 27.614-3.573-48.081-9.747-84.143-51.98-84.143-102.985v-1.299c13.969 7.797 30.214 12.67 47.431 13.319-28.264-18.843-46.781-51.005-46.781-87.391 0-19.492 5.197-37.36 14.294-52.954 51.655 63.675 129.3 105.258 216.365 109.807-1.624-7.797-2.599-15.918-2.599-24.04 0-57.828 46.782-104.934 104.934-104.934 30.213 0 57.502 12.67 76.67 33.137 23.715-4.548 46.456-13.32 66.599-25.34-7.798 24.366-24.366 44.833-46.132 57.827 21.117-2.273 41.584-8.122 60.426-16.243-14.292 20.791-32.161 39.308-52.628 54.253z"/></svg>[@mitchoharawild](https://twitter.com/mitchoharawild))
### 27 October 2020
### Slides @ [slides.mitchelloharawild.com/isf2020](https://slides.mitchelloharawild.com/isf2020)
]

---
class: inverse
background-image: linear-gradient(to right, rgba(150, 150, 150, .1), rgba(150, 150, 150, .4)), url("resources/timezone.jpg")
background-size: cover

---
class: inverse
background-image: linear-gradient(to right, rgba(50, 50, 50, .1), rgba(200, 200, 200, .4)), url("resources/sparks.jpg")
background-size: cover

# Presentation driven development

A 'proof of concept' for what is to come.

**Everything** is subject to change, suggestions encouraged.

---
class: center

## What is fable?

.animated.fadeIn[
.sticker[![tsibble](resources/tsibble.svg)]
.sticker[![tsibbledata](resources/tsibbledata.svg)]
.sticker[![feasts](resources/feasts.svg)]

## [tidyverts.org](http://www.tidyverts.org)
]

---

# Tidy temporal data structure

### Domestic tourism in Australia

```r
library(tsibble)
tourism
```

```
#> # A tsibble: 24,320 x 5 [1Q]
#> # Key:       Region, State, Purpose [304]
#>    Quarter Region   State           Purpose  Trips
#>      <qtr> <chr>    <chr>           <chr>    <dbl>
#>  1 1998 Q1 Adelaide South Australia Business  135.
#>  2 1998 Q2 Adelaide South Australia Business  110.
#>  3 1998 Q3 Adelaide South Australia Business  166.
#>  4 1998 Q4 Adelaide South Australia Business  127.
#>  5 1999 Q1 Adelaide South Australia Business  137.
#>  6 1999 Q2 Adelaide South Australia Business  200.
#>  7 1999 Q3 Adelaide South Australia Business  169.
#>  8 1999 Q4 Adelaide South Australia Business  134.
#>  9 2000 Q1 Adelaide South Australia Business  154.
#> 10 2000 Q2 Adelaide South Australia Business  169.
#> # … with 24,310 more rows
```

---

# A peak at the features

### Compute features relating to STL decompositions

```r
library(feasts)
tourism_features <- tourism %>% 
  features(Trips, feature_set(tags = "stl"))
```

```
#> # A tibble: 304 x 12
#>    Region State Purpose trend_strength seasonal_streng… seasonal_peak_y… seasonal_trough… spikiness linearity curvature
#>    <chr>  <chr> <chr>            <dbl>            <dbl>            <dbl>            <dbl>     <dbl>     <dbl>     <dbl>
#>  1 Adela… Sout… Busine…          0.451            0.380                3                1   1.62e+2    -5.51     71.4  
#>  2 Adela… Sout… Holiday          0.541            0.601                1                3   9.95e+0    48.9      78.2  
#>  3 Adela… Sout… Other            0.743            0.189                2                1   2.28e+0    95.0      43.4  
#>  4 Adela… Sout… Visiti…          0.433            0.446                1                3   6.01e+1    34.9      71.1  
#>  5 Adela… Sout… Busine…          0.453            0.140                3                1   1.18e-1     0.944    -3.30 
#>  6 Adela… Sout… Holiday          0.512            0.244                2                1   1.92e-1    10.4      23.8  
#>  7 Adela… Sout… Other            0.584            0.374                2                2   5.02e-4     4.29      3.21 
#>  8 Adela… Sout… Visiti…          0.481            0.228                0                3   7.56e+0    34.3      -0.587
#>  9 Alice… Nort… Busine…          0.526            0.224                0                1   1.79e-1    23.8      19.6  
#> 10 Alice… Nort… Holiday          0.377            0.827                3                1   8.13e-1   -19.6      10.3  
#> # … with 294 more rows, and 2 more variables: stl_e_acf1 <dbl>, stl_e_acf10 <dbl>
```

---

# A peak at the features

```r
library(ggplot2)
tourism_features %>% 
  ggplot(aes(x = trend_strength, y = seasonal_strength_year, colour = Purpose)) + 
  geom_point() + stat_density_2d(aes(fill = Purpose, alpha = ..level..), bins = 5, geom = "polygon") + 
  facet_wrap(vars(Purpose), nrow = 1) + coord_equal() + xlim(c(0,1)) + ylim(c(0,1))
```

<img src="figure/tourism-features-plot-1.svg" style="display: block; margin: auto;" />
]

---

# Total holidays across Australia

### Starting from the top

```r
library(dplyr)
aus_holiday <- tourism %>% 
  filter(Purpose == "Holiday") %>% 
  summarise(Trips = sum(Trips))

aus_holiday %>% 
  autoplot(Trips)
```

```
#> # A tsibble: 80 x 2 [1Q]
#>    Quarter  Trips
#>      <qtr>  <dbl>
#>  1 1998 Q1 11806.
#>  2 1998 Q2  9276.
#>  3 1998 Q3  8642.
#>  4 1998 Q4  9300.
#>  5 1999 Q1 11172.
#>  6 1999 Q2  9608.
#>  7 1999 Q3  8914.
#>  8 1999 Q4  9026.
#>  9 2000 Q1 11071.
#> 10 2000 Q2  9196.
#> # … with 70 more rows
```
.full-width[
<img src="figure/tourism-total-plot-1.svg" style="display: block; margin: auto;" />
]
]

---

# Forecasting with fable

### Look at the data

<br>

```r
*aus_holiday
```
]
.right-plot[

---

# Forecasting with fable

### Specify and estimate a model

<br>

```r
aus_holiday %>% 
* model(ETS(Trips))
```
]
.right-plot[

```
#> # A mable: 1 x 1
#>   `ETS(Trips)`
#>        <model>
#> 1 <ETS(M,N,M)>
```
]

---

# Forecasting with fable

### Make some forecasts

<br>

```r
aus_holiday %>% 
  model(ETS(Trips)) %>% 
* forecast(h = "3 years")
```
]
.right-plot[

```
#> # A fable: 12 x 4 [1Q]
#> # Key:     .model [1]
#>    .model     Quarter            Trips  .mean
#>    <chr>        <qtr>           <dist>  <dbl>
#>  1 ETS(Trips) 2018 Q1 N(13088, 368663) 13088.
#>  2 ETS(Trips) 2018 Q2 N(10909, 288993) 10909.
#>  3 ETS(Trips) 2018 Q3 N(10442, 294873) 10442.
#>  4 ETS(Trips) 2018 Q4 N(10624, 336446) 10624.
#>  5 ETS(Trips) 2019 Q1 N(13088, 558216) 13088.
#>  6 ETS(Trips) 2019 Q2 N(10909, 420721) 10909.
#>  7 ETS(Trips) 2019 Q3 N(10442, 415586) 10442.
#>  8 ETS(Trips) 2019 Q4 N(10624, 461447) 10624.
#>  9 ETS(Trips) 2020 Q1 N(13088, 747978) 13088.
#> 10 ETS(Trips) 2020 Q2 N(10909, 552594) 10909.
#> 11 ETS(Trips) 2020 Q3 N(10442, 536433) 10442.
#> 12 ETS(Trips) 2020 Q4 N(10624, 586585) 10624.
```
]

---

# Forecasting with fable

### See the results!

<br>

```r
aus_holiday %>% 
  model(ETS(Trips)) %>% 
  forecast(h = "3 years") %>% 
* autoplot(aus_holiday)
```
]
.right-plot[
<img src="figure/tourism-total-ets-plot-output-1.svg" style="display: block; margin: auto;" />
]

---
class: inverse, center
background-image: linear-gradient(to right, rgba(50, 50, 50, .5), rgba(50, 50, 50, .5)), url("resources/hierarchy.jpg")
background-size: cover

# Forecasting hierarchical series

`aggregate_key()`

---

# Forecasting hierarchical series

```r
data %>% 
  aggregate_key(grp_a / grp_b, y = sum(y))
```

![](resources/graph_hierarchy.png)

---

# Consider tourism by state

# Consider tourism by state

```r
tourism_state <- tourism %>% 
* aggregate_key(State, Trips = sum(Trips))
```

```
#> # A tsibble: 720 x 3 [1Q]
#> # Key:       State [9]
#>    Quarter State         Trips
#>      <qtr> <chr*>        <dbl>
#>  1 1998 Q1 <aggregated> 23182.
#>  2 1998 Q2 <aggregated> 20323.
#>  3 1998 Q3 <aggregated> 19827.
#>  4 1998 Q4 <aggregated> 20830.
#>  5 1999 Q1 <aggregated> 22087.
#>  6 1999 Q2 <aggregated> 21458.
#>  7 1999 Q3 <aggregated> 19914.
#>  8 1999 Q4 <aggregated> 20028.
#>  9 2000 Q1 <aggregated> 22339.
#> 10 2000 Q2 <aggregated> 19941.
#> # … with 710 more rows
```
]

---

# Forecasting each series separately

```r
fc <- tourism_state %>% 
  model(ets = ETS(Trips)) %>% 
  forecast(h = "3 years")
```

```
#> # A fable: 108 x 5 [1Q]
#> # Key:     State, .model [9]
#>    State  .model Quarter         Trips .mean
#>    <chr*> <chr>    <qtr>        <dist> <dbl>
#>  1 ACT    ets    2018 Q1  N(701, 7651)  701.
#>  2 ACT    ets    2018 Q2  N(717, 8032)  717.
#>  3 ACT    ets    2018 Q3  N(734, 8440)  734.
#>  4 ACT    ets    2018 Q4  N(750, 8882)  750.
#>  5 ACT    ets    2019 Q1  N(767, 9368)  767.
#>  6 ACT    ets    2019 Q2  N(784, 9905)  784.
#>  7 ACT    ets    2019 Q3 N(800, 10503)  800.
#>  8 ACT    ets    2019 Q4 N(817, 11171)  817.
#>  9 ACT    ets    2020 Q1 N(833, 11919)  833.
#> 10 ACT    ets    2020 Q2 N(850, 12756)  850.
#> # … with 98 more rows
```

]

.pull-right[
<img src="figure/tourism-agg-base-output-1.svg" style="display: block; margin: auto;" />

]

---

# Reconciling forecasts by state

```r
fc_coherent <- tourism_state %>% 
  model(ets = ETS(Trips)) %>% 
* reconcile(mint_ets = min_trace(ets)) %>%
  forecast(h = "3 years")
```

```
#> # A fable: 216 x 5 [1Q]
#> # Key:     State, .model [18]
#>    State  .model Quarter         Trips .mean
#>    <chr*> <chr>    <qtr>        <dist> <dbl>
#>  1 ACT    ets    2018 Q1  N(701, 7651)  701.
#>  2 ACT    ets    2018 Q2  N(717, 8032)  717.
#>  3 ACT    ets    2018 Q3  N(734, 8440)  734.
#>  4 ACT    ets    2018 Q4  N(750, 8882)  750.
#>  5 ACT    ets    2019 Q1  N(767, 9368)  767.
#>  6 ACT    ets    2019 Q2  N(784, 9905)  784.
#>  7 ACT    ets    2019 Q3 N(800, 10503)  800.
#>  8 ACT    ets    2019 Q4 N(817, 11171)  817.
#>  9 ACT    ets    2020 Q1 N(833, 11919)  833.
#> 10 ACT    ets    2020 Q2 N(850, 12756)  850.
#> # … with 206 more rows
```

]

.pull-right[
<img src="figure/tourism-agg-coherent-output-1.svg" style="display: block; margin: auto;" />

]

---
class: inverse, center
background-image: linear-gradient(to right, rgba(50, 50, 50, .5), rgba(50, 50, 50, .5)), url("resources/many_clocks.jpg")
background-size: cover

# Forecasting temporal granularities

`aggregate_index()`

---

# Forecasting temporal granularities

```r
library(moment) # pkg.mitchelloharawild.com/moment/
data %>% 
  aggregate_index(c(tu_year(1), tu_quarter(2), tu_quarter(1), tu_month(1)), y = sum(y))
```

![](resources/graph_temporal.png)

---

# Temporal aggregations of Australian holidays

```r
library(moment)
aus_holiday_temporal <- aus_holiday %>% 
* aggregate_index(c(tu_year(1), tu_quarter(2), tu_quarter(1)), Trips = sum(Trips))
```

```
#> # A tsibble: 140 x 2 [1Y, 2Q, 1Q]
#>     Quarter  Trips
#>    <moment>  <dbl>
#>  1     1998 39024.
#>  2     1999 38719.
#>  3     2000 38599.
#>  4     2001 37990.
#>  5     2002 38159.
#>  6     2003 37277.
#>  7     2004 37216.
#>  8     2005 34833.
#>  9     2006 37218.
#> 10     2007 37601.
#> # … with 130 more rows
```
]

.pull-right[
<img src="figure/tourism-time-agg-plot-1.svg" style="display: block; margin: auto;" />
]

---

# Forecasting each series separately

```r
fc <- aus_holiday_temporal %>% 
  model(ets = ETS(Trips)) %>% 
  forecast(h = "3 years")
```

```
#> # A fable: 21 x 4 [1Y, 2Q, 1Q]
#> # Key:     .model [1]
#>    .model    Quarter             Trips  .mean
#>    <chr>    <moment>            <dist>  <dbl>
#>  1 ets          2018 N(44587, 1768165) 44587.
#>  2 ets          2019 N(44587, 3535974) 44587.
#>  3 ets          2020 N(44587, 5303783) 44587.
#>  4 ets    2018 Q3-Q4  N(24133, 490947) 24133.
#>  5 ets    2019 Q3-Q4  N(22168, 614810) 22168.
#>  6 ets    2020 Q3-Q4   N(25076, 8e+05) 25076.
#>  7 ets    2021 Q3-Q4 N(23110, 1060950) 23110.
#>  8 ets    2022 Q3-Q4 N(26019, 1408284) 26019.
#>  9 ets    2023 Q3-Q4 N(24053, 1855035) 24053.
#> 10 ets       2018 Q1  N(13088, 368663) 13088.
#> # … with 11 more rows
```

]

]

---

# Reconciling forecasts over time

```r
fc <- aus_holiday_temporal %>% 
  model(ets = ETS(Trips)) %>% 
* reconcile(mint_ets = min_trace(ets)) %>%
  forecast(h = "3 years")
```

```
#> # A fable: 42 x 4 [1Y, 2Q, 1Q]
#> # Key:     .model [2]
#>    .model      Quarter             Trips  .mean
#>    <chr>      <moment>            <dist>  <dbl>
#>  1 ets            2018 N(44587, 1768165) 44587.
#>  2 ets            2019 N(44587, 3535974) 44587.
#>  3 ets            2020 N(44587, 5303783) 44587.
#>  4 mint_ets       2018  N(45454, 461746) 45454.
#>  5 mint_ets       2019  N(46170, 743625) 46170.
#>  6 mint_ets       2020 N(46885, 1119338) 46885.
#>  7 ets      2018 Q3-Q4  N(24133, 490947) 24133.
#>  8 ets      2019 Q3-Q4  N(22168, 614810) 22168.
#>  9 ets      2020 Q3-Q4   N(25076, 8e+05) 25076.
#> 10 ets      2021 Q3-Q4 N(23110, 1060950) 23110.
#> # … with 32 more rows
```

]

.pull-right[
<img src="figure/tourism-time-agg-coherent-output-1.svg" style="display: block; margin: auto;" />

]

---
class: inverse, center
background-image: linear-gradient(to bottom, rgba(50, 50, 50, .8), rgba(50, 50, 50, .5)), url("resources/clock_stack.jpg")
background-size: cover

# Cross-temporal forecast reconciliation

`aggregate_key() %>% aggregate_index()`

---

# Cross-temporal forecast reconciliation

```r
data %>% 
  aggregate_key(grp_a / grp_b, y = sum(y)) %>% 
  aggregate_index(c(tu_year(1), tu_quarter(2), tu_quarter(1), tu_month(1)), y = sum(y))
```

.pull-left[
![](resources/graph_hierarchy.png)
]
.pull-right[
![](resources/graph_temporal.png)
]
---

# Purpose of travel at multiple granularities

```r
library(moment)
tourism_purpose_temporal <- tourism %>% 
* aggregate_key(Purpose, Trips = sum(Trips)) %>%
* aggregate_index(c(tu_year(1), tu_quarter(1)), Trips = sum(Trips))
```

```
#> # A tsibble: 500 x 3 [1Y, 1Q]
#> # Key:       Purpose [5]
#>     Quarter Purpose       Trips
#>    <moment> <chr*>        <dbl>
#>  1     1998 Business     15474.
#>  2     1998 Holiday      39024.
#>  3     1998 Other         2820.
#>  4     1998 Visiting     26844.
#>  5     1998 <aggregated> 84162.
#>  6     1999 Business     15940.
#>  7     1999 Holiday      38719.
#>  8     1999 Other         3320.
#>  9     1999 Visiting     25508.
#> 10     1999 <aggregated> 83488.
#> # … with 490 more rows
```
]

.pull-right[
<img src="figure/tourism-crosstemp-agg-plot-1.svg" style="display: block; margin: auto;" />
]

---

# Forecasting each series separately

```r
fc <- tourism_purpose_temporal %>% 
  model(ets = ETS(Trips)) %>% 
  forecast(h = "3 years")
```

```
#> # A fable: 75 x 5 [1Y, 1Q]
#> # Key:     Purpose, .model [5]
#>    Purpose  .model  Quarter             Trips  .mean
#>    <chr*>   <chr>  <moment>            <dist>  <dbl>
#>  1 Business ets        2018   N(22296, 3e+06) 22296.
#>  2 Business ets        2019   N(22296, 6e+06) 22296.
#>  3 Business ets        2020 N(22296, 9078586) 22296.
#>  4 Business ets     2018 Q1   N(4937, 137658)  4937.
#>  5 Business ets     2018 Q2   N(5604, 209086)  5604.
#>  6 Business ets     2018 Q3   N(5888, 268396)  5888.
#>  7 Business ets     2018 Q4   N(5575, 293305)  5575.
#>  8 Business ets     2019 Q1    N(4937, 3e+05)  4937.
#>  9 Business ets     2019 Q2   N(5604, 367634)  5604.
#> 10 Business ets     2019 Q3   N(5888, 427147)  5888.
#> # … with 65 more rows
```

]

]

---

# Reconciling forecasts over time

```r
fc <- tourism_purpose_temporal %>% 
  model(ets = ETS(Trips)) %>% 
* reconcile(mint_ets = min_trace(ets)) %>%
  forecast(h = "3 years")
```

```
#> # A fable: 150 x 5 [1Y, 1Q]
#> # Key:     Purpose, .model [10]
#>    Purpose  .model    Quarter             Trips  .mean
#>    <chr*>   <chr>    <moment>            <dist>  <dbl>
#>  1 Business ets          2018   N(22296, 3e+06) 22296.
#>  2 Business ets          2019   N(22296, 6e+06) 22296.
#>  3 Business ets          2020 N(22296, 9078586) 22296.
#>  4 Business mint_ets     2018  N(22431, 629048) 22431.
#>  5 Business mint_ets     2019 N(22803, 1115028) 22803.
#>  6 Business mint_ets     2020 N(23174, 1621665) 23174.
#>  7 Business ets       2018 Q1   N(4937, 137658)  4937.
#>  8 Business ets       2018 Q2   N(5604, 209086)  5604.
#>  9 Business ets       2018 Q3   N(5888, 268396)  5888.
#> 10 Business ets       2018 Q4   N(5575, 293305)  5575.
#> # … with 140 more rows
```

]

.pull-right[
<img src="figure/tourism-crosstemp-agg-coherent-output-1.svg" style="display: block; margin: auto;" />

]