Chapter-05-Notes.Rmd

---
title: "Chapter 05 Notes"
output: html_document
---

```{r setup, include=FALSE}

library(knitr)
library(fpp3)

knitr::opts_chunk$set(echo = TRUE)
```

## 5.1 A Tidy Forecasting Workflow

![](workflow-1.png)

#### __Data Prep (Tidy)__

* Review the data to ensure correct format
* Pre-process using `tsibble` or `tidyverse` packages

__Example__
Modeling GDP per capita over time requires transformation
```{r}
gdppc <- global_economy %>%
  mutate(GDP_per_capita = GDP / Population)
```

#### __Plot the Data (Visualize)__

* Visualization increases understanding of the data
* Identify patterns and specify model

```{r}
gdppc %>%
  filter(Country == 'Sweden') %>%
  autoplot(GDP_per_capita) +
  labs(y="$US", title = "GDP per capita for Sweden")
```

#### __Define a Model (Specify)__

* This is likely the most important step as the majority of the book is written
on it, but no real information is given except for the `fable` library.

__fable__

* Most models in `fable` are specified using model functions with the follwing syntax:
`(y ~  x)`
  + where `y` is the response variable and `x` is the model structure. Example:
  `TSLM(GDP_per_capita ~ trend())`
  + In this case, the model `TSLM()` (Time Series Linear Model), the response 
  variable is `GDP_per_capita`, and it is being modelled using `trend()`. 
* The special functions like `trend()` vary between models, check the docs under
"Specials" section for each model
* Left side of formula also supports transformations

#### __Train the Model (Estimate)__

* Next we train the model on the data
* One or more models specifications can be estimated in the `model()` function

__Example__
```{r}
fit <- gdppc %>%
  model(trend_model = TSLM(GDP_per_capita ~ trend()))
```
```{r}
fit
```

`trend_model` column contains info about the fitted model for each country. This
is expanded upon later.

#### __Check Model Performance (Evaluate)__

Once fit, a model must be evaluated on the data.

#### __Produce Forecasts (Forecast)__

Using `forecast()` you can specify the forecast interval using the `h` argument 
(which stands for horizon). You can use integers or language such as "2 years".

```{r}
fit %>% 
  forecast(h= "3 years") %>%
  filter(Country == "Sweden") %>%
  autoplot(gdppc) +
  labs(y = '$US', title = 'GDP Per Capita for Sweden')
```

***

## 5.2 Simple Forecasting Methods

For this chapter we will use Australian brick production using the `filter_index()`
method for selecting a subsection of a time series:
```{r}
bricks <- aus_production %>%
  filter_index("1970 Q1" ~ "2004 Q4")
```

#### __Average Method__

The forecast is equal to the mean of the historical data

```{r}
bricks_fit <- bricks %>% model(MEAN(Bricks))

bricks_fc <- bricks_fit %>%
  forecast(h=20)

bricks_fc %>%
  autoplot(bricks)
```

#### __Naive Method__

For this method, we set all forecasts to be the value of the last observation.
Works well for many economic and financial time series.

```{r}
bricks_fit <- bricks %>% model(NAIVE(Bricks))

bricks_fc <- bricks_fit %>%
  forecast(h=20)

bricks_fc %>%
  autoplot(bricks)
```

#### __Seasonal Naive Method__

* Useful for highly seasonal data
* Set each forecast point to be equal to the last observed value from the same 
season of the previous year

```{r}
bricks_fit <- bricks %>% model(SNAIVE(Bricks))

bricks_fc <- bricks_fit %>%
  forecast(h=20)

bricks_fc %>%
  autoplot(bricks)
```

#### __Drift Method__

* Variation on Naive method
* Allows forecasts to increase/decrease over time
* Drift: amount of change over time
* Equivalent to draing a line between first and last observation and exrapolating
into the future

```{r}
bricks_fit <- bricks %>% model(RW(Bricks ~ drift()))

bricks_fc <- bricks_fit %>%
  forecast(h=20)

bricks_fc %>%
  autoplot(bricks)
```

__Summary__
In most cases, there will be better forecasting methods available, but these
are good to use for benchmarks. If a different method doesn't perform better
than these simple methods, they should be abandoned.

***


## 5.3 Fitted Values and Residuals

#### __Fitted Values__

* Definition: all previous observations in a time series that we're using to forecast
* Typically denoted as $\hat{y}$

#### __Residuals__

* Definition: what's lefter over after fitting a model. 
* Innovation residuals: looking at the residuals on a transformed scale if a 
transformation has been used in the model
* Fitted values and residuals can be obtained by using `augment()`

__Example__
```{r}
augment(bricks_fit)
```

Three new columns have been added to the data. .resid and .innov are the same
as there was no transformation used in model.

If patterns are observable in the innovation residuals, then the model could be 
improved.

***


## 5.4 Residual Diagnostics

A good forecast will yield innovation residuals with these properties:

1. Residuals are uncorrelated. 
    a. If correlation exists, then there is information in the residuals that
  should be used to compute the forecast.
2. Residuals must have a mean of zero.
    a. If the mean is other than zero, the forecast is biased.
3. Innovation residuals are homoscedastic, meaning they have a constant variance.
4. Innovation residuals are normally distributed.
 
Numbers 1 and 2 are essential while numbers 3 and 4 are useful but not necessary.
Fixing issues with number 1 is hard and is not addressed until Ch 10. Fixing issues
with number 2 is easy: for mean $m$, subtract $m$ from all forecasts to remove the bias.

__Example__

For stock market prices and indexes, the best forecasting method is often the Naive
method.
```{r}
# Re-index based on trading days
google_stock <- gafa_stock %>%
  filter(Symbol == "GOOG", year(Date) >= 2015) %>%
  mutate(day = row_number()) %>%
  update_tsibble(index = day, regular = TRUE)

# Filtering to 2015 only
google_2015 <- google_stock %>% 
  filter(year(Date) == 2015)

autoplot(google_2015, Close) +
  labs(y= '$US',
       title = 'Google Daily Closing Stock Prices in 2015')
  
```

To get the residuals from the above series:
```{r}
aug <- google_2015 %>%
  model(NAIVE(Close)) %>%
  augment() 

autoplot(aug, .innov) +
  labs(y = '$US',
       title = 'Residuals from the Naive Method')
```

```{r}
aug %>%
  ggplot(aes(x = .innov)) +
  geom_histogram() +
  labs(title = 'Histogram of Residuals')
```

```{r}
aug %>%
  ACF(.innov) %>%
  autoplot() +
  labs(title = "Residuals from the Naive Method")
``` 

The above graphs show the Naive method produces a forecast that accounts for all available information. 

* Mean of residuals is near zero
* Residuals exhibit homoskedasticity 
* Histogram suggests the residuals may not be normal (would need a statistical
check to verify)
* Forecasts from this method would likely be quite good, but prediction intervals
that are computed assuming a normal distribution may be innaccurate.

A shortcut to the above three graphs is `gg_tsresiduals()`
```{r}
google_2015 %>%
  model(NAIVE(Close)) %>%
  gg_tsresiduals()
```

#### __Portmanteau Tests for Autocorrelation__

More formal test for autocorrelation by considering a whole set of $r_k$ values 
as a group, recalling $r_k$ is the autocorrelation for lag $k$.

When viewing an ACF plot, it is essentially multiple hypothesis tests, of which
any one could provide a false positive, and in all likelihood, one may do so due
to the sheer number. This would lead to the incorrect assumption that residuals
have remaining autocorrelation.

In order to overcome this, we use portmanteau tests (a French word describing
a suitcase or coat rack carrying several items of clothing.)

2 portmanteau tests recommended:

1. Box-Pierce Test
2. Ljung-Box Test
    a. This is a more accurate test than the Box-Pierce
    
```{r}
aug %>% features(.innov, ljung_box, lag=10, dof=0)
```

Notes on the above:

* For non-seasonal data, lag value of 10 is recommended. 
* For seasonal data, a lag of double the seasonal period is recommended. 
* `dof` stands for "degrees of freedom" and is based on the number of parameters 
used, of which the naive method uses none.
* The p-value here is > 0.05, which is not significant, and this is not 
distinguishable from white noise

Alternative approach would be the drift method
```{r}
fit <- google_2015 %>% model(RW(Close ~ drift()))
tidy(fit)
```

```{r}
augment(fit) %>% features(.innov, ljung_box, lag=10, dof=1)
```

Using the Ljung-Box test with 1 degree of freedom for the estimated parameter.
Same with the Naive method, the Drift method is also indistinguishable from white
noise.

***

## 5.5 Distributional Forecasts and Prediction Intervals

#### __Forecast Distributions__

Due to the inherent uncertainty of the future, we must account for that uncertainty
in our forecasts using a probability distribution.

#### __Prediction Intervals__

Prediction intervals is a probability range within which we expect $y_t$ to lie.
Usually this is an 80% or 95% confidence interval, but you can use another 
multiplier via the below formula and chart where $c$ is the multiplier and $\hat{\sigma}_h$
is an estimate of the standard deviation of the $h$-step forecast distribution.

$$\hat{y}_{T+h} \pm c\hat{\sigma}_h$$

::: l-body-outset
|Pct Confidence Interval |	Multiplier ($c$)|
|-----------|-----------|
|50 |	0.67|
|55 |	0.76|
|60 |	0.84|
|65 |	0.93|
|70 |	1.04|
|75 |	1.15|
|80 |	1.28|
|85 |	1.44|
|90 |	1.64|
|95 |	1.96|
|96 | 2.05|
|97 |	2.17|
|98 |	2.33|
|99 |	2.58|
:::

#### __One-step Prediction Intervals__

When forecasting one step ahead, the standard deviation of the forecast distribution
can be estimated using the standard deviation of the residuals.

__Example__
For `google_2015` stock price data, the last value observed is 758.88. So for a Naive
forecast, the forecasted value is 758.88. The standard deviation of the residuals 
is 11.19. A 95% prediction interval (1.96 critical value) for the next value is 
$$758.88 \pm 1.96(11.19) = [736.9, 780.8]$$
If a different prediction interval is needed, use the above conversion chart.

#### __Multi-step Prediction Intervals__

The further into the future we forecast, the more inherent uncertainty, and the 
wider the prediction interval. So $\sigma_h$ increases with $h$.

#### __Benchmark Methods__

When producing multi-step forecasts for the four benchmark methods, it is possible 
to mathematically derive the forecast standard deviation and thus the prediction 
intervals, but prediction intervals can easily be computed for you using the 
`fable` package.

__Example__
```{r}
google_2015 %>%
  model(NAIVE(Close)) %>%
  forecast(h=10) %>%
  hilo()
```

the `hilo()` function converts the forecast distributions into 80% and 95% 
confidence intervals by default. Other options possible via the `level` argument.

Prediction intervals shown as shaded regions when plotted
```{r}
google_2015 %>%
  model(NAIVE(Close)) %>%
  forecast(h=10) %>%
  autoplot(google_2015) +
  labs(title="Google Daily Closing Stock Price", y='$US')
```

#### __Prediction Intervals from Bootstrapped Residuals__

When a normal distribution of residuals is an unreasonable assumption, an alternative
is to us bootstrapping which only assumes residuals are uncorrelated and have 
constant variance. 

Assuming future errors will be similar to past errors, we can sample from the 
collection of errors already observed (residuals) and use that sample to simulate
the next observation in the time series and repeat the process. Doing this multiple
times will simulate an entire set of future values for the time series. 

Repeating this process will give us many possible futures which we can view using
the `generate()` function
```{r}
fit <- google_2015 %>%
  model(NAIVE(Close))
sim <- fit %>% generate(h = 30, times = 5, bootstrap = TRUE)

sim
```

The `.rep` variable provides the key for each of the 5 futures created with the
`times` argument. To view each predicted future:
```{r}
google_2015 %>%
  ggplot(aes(x = day)) +
  geom_line(aes(y = Close)) +
  geom_line(aes(y = .sim, colour = as.factor(.rep)),
            data = sim) +
  labs(title="Google Daily Closing Stock Price with Predictions",
       y="$US") +
  guides(col=FALSE)
```

Prediction intervals are obtained by calculating percentiles of the future
samples paths for each forecast horizon using the `forecast()` function. 
The number of samples can be controlled by the `times` argument. 
```{r}
fc <- fit %>% forecast(h = 30, bootstrap=TRUE, times=100)
fc
```
```{r}
autoplot(fc, google_2015) +
  labs(title="Google daily closing stock price", y="$US" )
```

***

## 5.6 Forecasting Using Transformations

#### Prediction Intervals with Transformations

Essentially the `fable` package handles back transformations internally, but
the process is:

* Transform data
* Make prediction
* Back-transform data to original scale

This approach preserves the probability coverage, but it will no longer be symmetric
around the point forecast.

#### Bias Adjustments

One issue with using mathematical transformations such as the Box-Cox, is that 
the back-transformed point forecast will not be the man of the forecast, but will
usually be the median. The mean is preferable, but median is acceptable.

This fails in a situation where you want to add up sales forecasts from various
regions to form a forecast for the whole country. Medians do not add up, but means do.

The difference between the simple back-transformed forecast and the mean is called
the __bias__. When we use the mean rather than the median, we say the point forecasts 
have been __bias adjusted__. 

__Example__
```{r}
prices %>%
  filter(!is.na(eggs)) %>%
  model(RW(log(eggs) ~ drift())) %>%
  forecast(h = 50) %>%
  autoplot(prices %>% filter(!is.na(eggs)),
    level = 80, point_forecast = lst(mean, median)
  ) +
  labs(title = "Annual egg prices",
       y = "$US (in cents adjusted for inflation) ")
```

Notice how the skewed forecast distribution pulls up the forecast distribution's
mean. This is a result of the added term from the bias adjustment.

You can obtain the point forecast prior to bias adjustment by using the `median()`
function on the distribution column.

***

## 5.7 Forecasting with Decomposition

This chapter uses decomposition from Ch 3 in producing forecasts.

__Example__
```{r}
us_retail_employment <- us_employment %>%
  filter(year(Month) >= 1990, Title == "Retail Trade")

dcmp <- us_retail_employment %>%
  model(STL(Employed ~ trend(window = 7), robust = TRUE)) %>%
  components() %>%
  select(-.model)

dcmp %>%
  model(NAIVE(season_adjust)) %>%
  forecast() %>%
  autoplot(dcmp) +
  labs(y="Number of People",
       title = "US Retail Employment")
```

The above shows the Naive forecast for seasonally adjusted data. This is then
"reseasonalized" by adding in the Seasonal Naive forecasts of the seasonal 
componenet:

```{r}
fit_dcmp <- us_retail_employment %>%
  model(stlf = decomposition_model(
    STL(Employed ~ trend(window=7), robust = TRUE),
    NAIVE(season_adjust)
  ))
fit_dcmp %>%
  forecast %>%
  autoplot(us_retail_employment) +
  labs(y = "Number of People",
       title = "Monthly US Retail Employment")
```

In the above code: 

* `decomposition_model()` allows you to compute forecasts via
any additive decomp whil using other model functions to forecast each of the 
decomp's components.
    * Seasonal components will automatically be forecast using `SNAIVE()`
    * The function also will reseasonalize automatically

```{r}
fit_dcmp %>% gg_tsresiduals()
```

The ACF of the residuals show significant autocorrelation because the Naive method 
does not capture the changing trend in the seasonally adjusted series. Better
methods in future chapters.

***

## 5.8 Evaluating Point Forecast Accuracy

#### __Train / Test Sets__

It is essential to evaluate forecast accuracy using genuine forecasts, and the 
accuracy of forecasts can only be determined by testing the model on unseen data.

Therefore, it is common practice to divide the available data into two portions
called the "__train/test split__" (usually 80/20). At minimum, the test set
should be as large as the desired forecast horizon (in order to accurately compare).
 
Helpful notes:

* Just because a model fits the training data well, it does not mean it will
forecast with accuracy
* If you have enough parameters, you can get a perfect fit
* Overfitting to the training data is just as bad as underfitting the data

#### __Creating the Train/Test Split__

`filter()` can easily extract a portion of a time series
```{r}
aus_production %>%
  filter(year(Quarter) >= 1995) # extracts all data from 1995 on
```

`slice()` allows the use of indices
```{r}
aus_production %>%
  slice(n()-19:0) # extracts the last 20 observations (5 years)
```

`slice()` also works with groups to subset observations from each key
```{r}
aus_retail %>%
  group_by(State, Industry) %>%
  slice(1:12)
```

#### __Forecast Errors__

An "error" is the difference between the forecast and the observed value. It 
differs from residuals in the following ways:

* Residuals are in regards to the training set while errors are in regards to the
test set
* Residuals are based on on-step forecasts while errors can involve multi-step
forecasts


#### __Scale-Dependent Measures__

Accuracy measures that are based only on the errors are scale-dependent and cannot 
be used to make comparison between series that involve different units.

Two most commonly used scale dependent measures:

1. MAE - Mean Absolute Error: $mean(\lvert{e_t}\rvert)$
2. RMSE - Root Mean Squared Error: $\sqrt{mean(e^2_t)}$


#### __Percentage Measures__

Percentage measures have the advantage of being unit-free, so can be compared 
between data sets.The most common is MAPE.

Disadvantages:

* When $y_t=0$ as they can be infinite or undefined.
* They assume the unit of measurement has a meaningful zero.
    *Example: percentage error makes no sense when measuring the accuracy of temp 
    forecasts for F$^o$ or C$^o$ as they have arbitrary zero points.
* They put a heavier penalty on negative errors than on positive errors
    * This lead to the development of sMAPE (symmetric MAPE), but this method still
    involves division by a number close to zero, and the value of sMAPE can be 
    negative, so it's not really an "absolute measure" at all
    
Percentage Measure Formulas

1. MAPE - Mean Absolute Percentage Error: $mean(\lvert{p_t}\rvert)$
2. sMAPE - symmetric MAPE: $mean(200\lvert{y_t - \hat{y}_t}\rvert/(y_t+\hat{y}_t))$

#### __Scaled Measures__

Methods scale the errors based on the training MAE from a simple forecast method.
A scaled error is < 1 if it arises from a better forecast than the average one-step
Naive forecast computed on the training data. It is > 1 if the forecast is worse.

The two scaled measures:

1. MASE - Mean Absolute Scaled Error: $mean(\lvert{q_j}\rvert)$
2. RMSSE - Root Mean Squared Scaled Error: $\sqrt{mean(q^2_j)}$

$q_j$ is a rather complicated formula you won't remember

__Example__
```{r}
recent_production <- aus_production %>%
  filter(year(Quarter) >= 1992)

beer_train <- recent_production %>%
  filter(year(Quarter) <= 2007)

beer_fit <- beer_train %>%
  model(
    Mean = MEAN(Beer),
    `Naive` = NAIVE(Beer),
    `Seasonal Naive` = SNAIVE(Beer),
    Drift = RW(Beer ~ drift())
  )

beer_fc <- beer_fit %>%
  forecast(h=10)

beer_fc %>%
  autoplot(
    recent_production,
    level = NULL
  ) +
  labs(
    y = 'Megalitres',
    title = 'Forecasts for Quarterly Beer Production'
  ) +
  guides(color = guide_legend(title = 'Forecast'))


accuracy(beer_fc, recent_production)
  
```

Graph shows the three forecast methods applied and the mean. Accuracy metrics
also shown. Seasonal Naive is obviously the best performer.

Non-seasonal example:
```{r}
google_fit <- google_2015 %>%
  model(
    Mean = MEAN(Close),
    Naive = NAIVE(Close),
    Drive = RW(Close ~ drift())
  )

google_jan_2016 <- google_stock %>%
  filter(yearmonth(Date) == yearmonth("2016 Jan"))

google_fc <- google_fit %>%
  forecast(google_jan_2016)

google_fc %>%
  autoplot(bind_rows(google_2015, google_jan_2016),
           level=NULL) +
  labs(y = "$US",
       title = "Google Closing Stock Prices from Jan 2015")+
  guides(color = guide_legend(title = "Forecast"))


```

In some instances, one of the simple methods will be best forecasting method,
but will likely serve as a benchmark to compare how well other methods work.

***

## 5.9 Evaluating Distributional Forecast Accuracy

When evaluating distributional forecasts, we need to use different measures.

#### __Quantile Scores__

  * A low quantile score indicates a better estimate of the quantile
  * Text does not say what low is, but gives 4.86 as an example answer which 
  implies this is low

__Example__
```{r}
google_fc %>%
  filter(.model == 'Naive', Date == "2016-01-04") %>%
  accuracy(google_stock, list(qs=quantile_score), probs=0.10)
```

### __Scale-Free Comparisons Using Skill Scores__

  * Computing a forecast accuracy measure relative to some benchmark method.
  * CRPS is a specific skill score that is left undefined by the curriculum as 
  to its acronym or actual formula.

__Example__
```{r}
google_fc %>%
  accuracy(google_stock, list(skill = skill_score(CRPS)))
```

Above we are comparing the drift method to the Naive method. Naive is 0 because
it can't improve on itself. The other two methods have larger CRPS values, so they
are negative values. The skill value is interpreted as "The drift method is 26.6%
worse than the Naive method".

`skill_score()`
  * Will always calculate the the CRPS for the appropriate benchmark. 
    + Example: when data is seasonal, it will use seasonal naive
  * Can be used with any accuracy measure:
    + MSE
      - Ensure test set is large enough to allow reliable calculation of the error
      measure
    + MASE or RMSSE are scale-free and preferred
    
***


## __5.10 Time Series Cross-Validation__

* For a time series, the cross_validation procedure contains a series of training 
sub-sets and a series of test sets with only 1 observation.
* Forecast accuracy is computed by averaging over the test sets.

__Example__
```{r}
google_2015_tr <- google_2015 %>%
  stretch_tsibble(.init = 3, .step = 1) %>%
  relocate(Date, Symbol, .id)

google_2015_tr

```

Above:
  * `stretch_tsibble()` creates many training sets
    + `.init` determines the size of the first training set
    + `.step` increases each successive training set by that number of steps
  * for the output, the `.id` column indicates which training set
  
Using accuracy to evaluate the forecast accuracy:
```{r}
google_2015_tr %>%
  model(RW(Close ~ drift())) %>%
  forecast(h = 1) %>%
  accuracy(google_2015)

google_2015 %>%
  model(RW(Close ~ drift())) %>%
  accuracy()
```

 A good way to choose the best forecasting model is to find the model with the 
 smallest RMSE computed using time series cross-validation
 
 __Example__
 
 Below code evaluates the forecasting performance of 1 to 8-step-ahead drift
 forecasts. Forecast error increases as the forecast horizon increases.
 
```{r}
google_2015_tr <- google_2015 %>%
  stretch_tsibble(.init = 3, .step = 1)

fc <- google_2015_tr %>%
  model(RW(Close ~ drift())) %>%
  forecast(h = 8) %>%
  group_by(.id) %>%
  mutate(h = row_number()) %>%
  ungroup()

fc %>%
  accuracy(google_2015, by = c('h', '.model')) %>%
  ggplot(aes(x = h, y = RMSE)) +
  geom_point()
  
```