The process of producing forecasts can be split up into a few fundamental steps.

1. Preparing data
2. Data visualisation
3. Specifying a model
4. Model estimation
5. Accuracy & performance evaluation
6. Producing forecasts

To illustrate the process, we will fit linear trend models to national GDP data stored in global_economy.

• The first step in forecasting is to prepare data in the correct format.

• It envolves loading in data, identifying missing values, filtering the time series, and other pre-processing tasks.

• The functionality provided by tsibble and other packages in the tidyverse simplifies this step.

• Each data set is unique, so it is an essential step to understanding the data’s features and should always be done before models are estimated.

• We create GDP per capita to be able to compare GDP across countries

gdppc <- global_economy %>%
  mutate(GDP_per_capita = GDP/Population) %>%
  select(Year, Country, GDP, Population, GDP_per_capita)
gdppc

## # A tsibble: 15,150 x 5 [1Y]
## # Key:       Country [263]
##     Year Country             GDP Population GDP_per_capita
##    <dbl> <fct>             <dbl>      <dbl>          <dbl>
##  1  1960 Afghanistan  537777811.    8996351           59.8
##  2  1961 Afghanistan  548888896.    9166764           59.9
##  3  1962 Afghanistan  546666678.    9345868           58.5
##  4  1963 Afghanistan  751111191.    9533954           78.8
##  5  1964 Afghanistan  800000044.    9731361           82.2
##  6  1965 Afghanistan 1006666638.    9938414          101. 
##  7  1966 Afghanistan 1399999967.   10152331          138. 
##  8  1967 Afghanistan 1673333418.   10372630          161. 
##  9  1968 Afghanistan 1373333367.   10604346          130. 
## 10  1969 Afghanistan 1408888922.   10854428          130. 
## # … with 15,140 more rows

• Visualisation is an essential step in understanding the data (identify common patterns, and subsequently specify an appropriate model.)

gdppc %>%
  filter(Country=="Brazil") %>%
  autoplot(GDP_per_capita) +
    labs(title = "GDP per capita for Brazil", y = "$US")

• There are many different time series models that can be used for forecasting.

• Specifying an appropriate model for the data is essential for producing appropriate forecasts.

• We can specify different models using the function model in the package fable.

• This function use a formula (y ~ x) interface. The response variable(s) are specified on the left of the formula, and the structure of the model is written on the right.

gdppc %>%
  filter(Country== "Brazil") %>%
      model(TSLM(GDP_per_capita ~ trend()))

## # A mable: 1 x 2
## # Key:     Country [1]
##   Country `TSLM(GDP_per_capita ~ trend())`
##   <fct>                            <model>
## 1 Brazil                            <TSLM>

• In this case, we are using the TSLM model (time series linear model)
• The response variable is GDP_per_capita and it is being modelled using trend()
• The special functions used to define the model’s structure vary between models (as each model can support different structures).

• Once an appropriate model is specified, the next step is to estimate (train) the model given our data.

• One or more model specifications can be estimated using the model() function.

estim<-gdppc %>%
      model(trend_model = TSLM(GDP_per_capita ~ trend()))
estim

## # A mable: 263 x 2
## # Key:     Country [263]
##    Country             trend_model
##    <fct>                   <model>
##  1 Afghanistan              <TSLM>
##  2 Albania                  <TSLM>
##  3 Algeria                  <TSLM>
##  4 American Samoa           <TSLM>
##  5 Andorra                  <TSLM>
##  6 Angola                   <TSLM>
##  7 Antigua and Barbuda      <TSLM>
##  8 Arab World               <TSLM>
##  9 Argentina                <TSLM>
## 10 Armenia                  <TSLM>
## # … with 253 more rows

• Once a model has been fitted, it is important to check how well it has performed on the data.

• There are several diagnostic tools available to check model behavior, and also accuracy measures that allow one model to be compared against another.

• It is an essential step to make sure your model account for the behavior of the data !!!

• We will see four simple forecasting methods as benchmarks. To illustrate them, we will use quarterly Australian clay brick production between 1970 and 2004.

bricks <- aus_production %>%
  filter_index("1970 Q1" ~ "2004 Q4") %>%
  select(Bricks)

• The filter_index() function is a convenient shorthand for extracting a section of a time series.

• Forecast of all future values is equal to mean of historical data \(\{y_1,\dots,y_T\}\).
• Forecasts: \(\hat{y}_{T+h|T} = \bar{y} = (y_1+\dots+y_T)/T\)

bricks %>% 
  model(MEAN(Bricks))%>%
    forecast(h = 20) %>%
      autoplot(bricks, level = NULL)

# Level= Null remove the Predict Interval

• Forecasts equal to last observed value.

• Forecasts: \(\hat{y}_{T+h|T} =y_T\).

• Consequence of efficient market hypothesis.

• Naïve forecast is optimal when data follow a random walk (see Section 9.1), these are also called random walk forecasts and the RW() function can be used instead of NAIVE.

bricks %>%
  model(NAIVE(Bricks)) %>%
  forecast(h = "5 years") %>%
  autoplot(bricks, level = NULL)

• A similar method is useful for highly seasonal data.
• We set each forecast to be equal to the last observed value from the same season (e.g., the same month of the previous year).
• Forecasts: \(\hat{y}_{T+h|T} =y_{T+h-m(k+1)}\), where \(m=\) seasonal period and \(k\) is the integer part of \((h-1)/m\). (i.e., the number of complete years in the forecast period prior to time \(T+h\))
• For example, with monthly data, the forecast for all future February values is equal to the last observed February value, etc

bricks %>%
  model(SNAIVE(Bricks ~ lag("year"))) %>%
  forecast(h = "5 years") %>%
  autoplot(bricks, level = NULL) +
  labs(title = "Clay brick production in Australia")

• A variation on the naïve method is to allow the forecasts to increase or decrease over tim.

• The amount of change over time (called the drift) is set to be the average change seen in the historical data.

• Forecasts equal to last value plus average change.

\(\hat{y}_{T+h|T} = y_{T} + \frac{h}{T-1}\sum_{t=2}^T (y_t-y_{t-1}) = y_T + \frac{h}{T-1}(y_T -y_1)\)

• Equivalent to extrapolating a line drawn between first and last observations.

bricks %>%
  model(RW(Bricks ~ drift())) %>%
  forecast(h = "5 years") %>%
  autoplot(bricks, level = NULL) +
  geom_line(data = slice(bricks, range(cumsum(!is.na(Bricks)))),
            aes(y=Bricks), linetype = "dashed", colour = "#0072B2") +
  labs(title = "Clay brick production in Australia")

brick_fit <- aus_production %>%
  filter_index("1970 Q1" ~ "2004 Q4") %>%
  model(
    Seasonal_naive = SNAIVE(Bricks),
    Naive = NAIVE(Bricks),
    Drift = RW(Bricks ~ drift()),
    Mean = MEAN(Bricks)
  )

## # A mable: 1 x 4
##   Seasonal_naive   Naive         Drift    Mean
##          <model> <model>       <model> <model>
## 1       <SNAIVE> <NAIVE> <RW w/ drift>  <MEAN>

brick_fc <- brick_fit %>%
  forecast(h = "5 years")

## # A fable: 80 x 4 [1Q]
## # Key:     .model [4]
##   .model         Quarter       Bricks .mean
##   <chr>            <qtr>       <dist> <dbl>
## 1 Seasonal_naive 2005 Q1 N(409, 3026)   409
## 2 Seasonal_naive 2005 Q2 N(423, 3026)   423
## 3 Seasonal_naive 2005 Q3 N(428, 3026)   428
## 4 Seasonal_naive 2005 Q4 N(397, 3026)   397
## # … with 76 more rows

brick_fc %>%
  autoplot(aus_production, level = NULL) +
  labs(title = "Clay brick production in Australia",
       y = "Millions of bricks") +
  guides(colour = guide_legend(title = "Forecast"))

# Extract training data
apple_stock <- gafa_stock %>%
  filter(Symbol == "AAPL") %>%
  mutate(trading_day = row_number()) %>%
  update_tsibble(index=trading_day, regular=TRUE)

• Because stock prices are not observed every day, we first set up a new time index based on the trading days rather than calendar days.

• “Residuals” in a time series model are what is left over after fitting a model.

• The residuals are equal to the difference between the observations and the corresponding fitted values: \(e_t = y_t-\hat{y}_{t|t-1}\).

• A good forecasting method will yield innovation residuals with the following properties:

1. The innovation residuals are uncorrelated. If there are correlations between innovation residuals, then there is information left in the residuals which should be used in computing forecasts.

2. The innovation residuals have zero mean. If they have a mean other than zero, then the forecasts are biased.

• Models that do not satisfy these properties can be improved.

• However, that does not mean that forecasting methods that satisfy these properties cannot be improved.

• Checking these properties is important in order to see whether a method is using all of the available information, but it is not a good way to select a forecasting method.

• Adjusting for bias is easy: if the residuals have mean \(m\), then simply subtract \(m\) from all forecasts and the bias problem is solved.

• Fixing the correlation problem is harder (Ch 10).

• The residuals to also have the following two properties (but not necessary):

1. The innovation residuals have constant variance. This is known as homoscedasticity
2. The innovation residuals are normally distributed.

• These two properties make the calculation of prediction intervals easier.
• Sometimes there is usually little that you can do to ensure that your innovation residuals have constant variance and a normal distribution (alternative approach to obtaining prediction intervals is necessary).

• It is important to state the distinction between residuals and innovation residuals.

• If a transformation has been used in the model, then it is often useful to look at residuals on the transformed scale (innovation residuals)

• For example, suppose we modeled the logarithms of the data, \(w_t=log(y_t)\), the innovation residuals would be \(w_t-\hat{w_t}\), whereas the residuals are given by \(y_t-\hat{y_t}\).

• The fitted value and residuals can be checked via augment()

• Apple Stock price (Close)

apple_stock %>% autoplot(Close)

fit <- apple_stock %>% model(NAIVE(Close))
augment(fit)

## # A tsibble: 1,258 x 7 [1]
## # Key:       Symbol, .model [1]
##    Symbol .model       trading_day Close .fitted .resid .innov
##    <chr>  <chr>              <int> <dbl>   <dbl>  <dbl>  <dbl>
##  1 AAPL   NAIVE(Close)           1  79.0    NA   NA     NA    
##  2 AAPL   NAIVE(Close)           2  77.3    79.0 -1.74  -1.74 
##  3 AAPL   NAIVE(Close)           3  77.7    77.3  0.421  0.421
##  4 AAPL   NAIVE(Close)           4  77.1    77.7 -0.556 -0.556
##  5 AAPL   NAIVE(Close)           5  77.6    77.1  0.489  0.489
##  6 AAPL   NAIVE(Close)           6  76.6    77.6 -0.991 -0.991
##  7 AAPL   NAIVE(Close)           7  76.1    76.6 -0.511 -0.511
##  8 AAPL   NAIVE(Close)           8  76.5    76.1  0.399  0.399
##  9 AAPL   NAIVE(Close)           9  78.1    76.5  1.52   1.52 
## 10 AAPL   NAIVE(Close)          10  79.6    78.1  1.57   1.57 
## # … with 1,248 more rows

• We estimate a naive forecast \(\hat{y}_{t|t-1} = y_{t-1}\)

• \(e_t = y_t - \hat{y}_{t|t-1} = y_t-y_{t-1}\)

• fitted contains the fitted values;

• resid contains the residuals;

• innov contains the “innovation residuals” which, in this case, are identical to the regular residuals.

• Why is there NAs on the first row of the fit table?

augment(fit) %>%
  ggplot(aes(x = trading_day)) +
  geom_line(aes(y = Close, colour = "Data")) +
  geom_line(aes(y = .fitted, colour = "Fitted"))

• Both curves are overlapping each other, right? Let’s take a closer look

augment(fit) %>%
  filter(trading_day > 1100) %>%
  ggplot(aes(x = trading_day)) +
  geom_line(aes(y = Close, colour = "Data")) +
  geom_line(aes(y = .fitted, colour = "Fitted"))

augment(fit) %>%
  autoplot(.resid) +
  labs(y = "$US",
       title = "Residuals from naïve method")

augment(fit) %>%
  ggplot(aes(x = .resid)) +
  geom_histogram(bins = 150) +
  labs(title = "Histogram of residuals")

augment(fit) %>%
  ACF(.resid) %>%
  autoplot() + labs(title = "ACF of residuals")

• The naïve method produces forecasts might account for all available information.

• The mean of the residuals is close to zero (histogram) but may not be normal — the right tail seems a little too long, even when we ignore the outlier.

• There are some significant correlation in the residual series (ACF)

• The variation of the residuals does not stay much the same across the historical data (end of the sample)

• This forecast method might not be good, and the prediction intervals that are computed assuming a normal distribution may be inaccurate.

• gg_tsresiduals() function

gg_tsresiduals(fit)

fb_stock <- gafa_stock %>%
  filter(Symbol == "FB") %>%
  mutate(trading_day = row_number()) %>%
  update_tsibble(index=trading_day, regular=TRUE)

fit2 <- fb_stock %>% model(NAIVE(Close))
augment(fit2)

## # A tsibble: 1,258 x 7 [1]
## # Key:       Symbol, .model [1]
##    Symbol .model       trading_day Close .fitted .resid .innov
##    <chr>  <chr>              <int> <dbl>   <dbl>  <dbl>  <dbl>
##  1 FB     NAIVE(Close)           1  54.7    NA   NA     NA    
##  2 FB     NAIVE(Close)           2  54.6    54.7 -0.150 -0.150
##  3 FB     NAIVE(Close)           3  57.2    54.6  2.64   2.64 
##  4 FB     NAIVE(Close)           4  57.9    57.2  0.720  0.720
##  5 FB     NAIVE(Close)           5  58.2    57.9  0.310  0.310
##  6 FB     NAIVE(Close)           6  57.2    58.2 -1.01  -1.01 
##  7 FB     NAIVE(Close)           7  57.9    57.2  0.720  0.720
##  8 FB     NAIVE(Close)           8  55.9    57.9 -2.03  -2.03 
##  9 FB     NAIVE(Close)           9  57.7    55.9  1.83   1.83 
## 10 FB     NAIVE(Close)          10  57.6    57.7 -0.140 -0.140
## # … with 1,248 more rows

gg_tsresiduals(fit2)

• We assume that the residuals are white noise (uncorrelated, mean zero, constant variance). If they aren’t, then there is information left in the residuals that should be used in computing forecasts.

• So a standard residual diagnostic is to check the ACF of the residuals of a forecasting method.

• We expect these to look like white noise.

• But we can do a more formal test for autocorrelation:

\(r_k =\) autocorrelation of residual at lag \(k\)

Consider a whole set of \(r_{k}\) values, and develop a test to see whether the set is significantly different from a zero set (instead of doing it individually) - portmanteau test

Box-Pierce test

\[Q = T \sum_{k=1}^\ell r_k^2\] where \(\ell\) is max lag being considered and \(T\) is number of observations.

• If each \(r_k\) close to zero, \(Q\) will be small.
• If some \(r_k\) values large (positive or negative), \(Q\) will be large.

Ljung-Box test - preferred method

\[Q^* = T(T+2) \sum_{k=1}^\ell (T-k)^{-1}r_k^2\] where \(\ell\) is max lag being considered and \(T\) is number of observations.

• My preferences: \(\ell=10\) for non-seasonal data, \(h=2m\) for seasonal data (where \(m\) is seasonal period).

• However, the test is not good when \(\ell\) is large, so if these values are larger than \(T/5\), then use \(\ell = T/5\).

• The null hypothesis for both test is the residuals are not distinguishable from a white noise series (not significant results).

• If data are WN, \(Q^*\) has \(\chi^2\) distribution with \((\ell - K)\) degrees of freedom where \(K=\) no. parameters in model.
• When applied to raw data, set \(K=0\).
• lag \(= \ell\), dof \(= K\)

#Apple Stocks 
augment(fit) %>%
  features(.resid, ljung_box, lag=10, dof=0)

## # A tibble: 1 × 4
##   Symbol .model       lb_stat lb_pvalue
##   <chr>  <chr>          <dbl>     <dbl>
## 1 AAPL   NAIVE(Close)    25.8   0.00399

• The results are significant (i.e., the p-values less than 5%). Thus, we can conclude that the residuals are distinguishable from a white noise series.

#Facebook Stocks
augment(fit2) %>%
  features(.resid, ljung_box, lag=10, dof=0)

## # A tibble: 1 × 4
##   Symbol .model       lb_stat lb_pvalue
##   <chr>  <chr>          <dbl>     <dbl>
## 1 FB     NAIVE(Close)    12.1     0.276

• The results are not significant (i.e., the p-values are large). Thus, we can conclude that the residuals are not distinguishable from a white noise series.

fit_drift <- apple_stock %>% model(Drift = RW(Close ~ drift()))
augment(fit_drift)

## # A tsibble: 1,258 x 7 [1]
## # Key:       Symbol, .model [1]
##    Symbol .model trading_day Close .fitted .resid .innov
##    <chr>  <chr>        <int> <dbl>   <dbl>  <dbl>  <dbl>
##  1 AAPL   Drift            1  79.0    NA   NA     NA    
##  2 AAPL   Drift            2  77.3    79.1 -1.80  -1.80 
##  3 AAPL   Drift            3  77.7    77.3  0.359  0.359
##  4 AAPL   Drift            4  77.1    77.8 -0.618 -0.618
##  5 AAPL   Drift            5  77.6    77.2  0.426  0.426
##  6 AAPL   Drift            6  76.6    77.7 -1.05  -1.05 
##  7 AAPL   Drift            7  76.1    76.7 -0.574 -0.574
##  8 AAPL   Drift            8  76.5    76.2  0.336  0.336
##  9 AAPL   Drift            9  78.1    76.6  1.46   1.46 
## 10 AAPL   Drift           10  79.6    78.1  1.50   1.50 
## # … with 1,248 more rows

gg_tsresiduals(fit_drift)

augment(fit_drift) %>%
  features(.resid, ljung_box, lag=10, dof=1)

## # A tibble: 1 × 4
##   Symbol .model lb_stat lb_pvalue
##   <chr>  <chr>    <dbl>     <dbl>
## 1 AAPL   Drift     25.8   0.00219

• We express the uncertainty in our forecasts using a probability distribution
• A forecast \(\hat{y}_{T+h|T}\) is (usually) the mean of the conditional distribution \(y_{T+h} \mid y_1, \dots, y_{T}\).
• Most time series models produce normally distributed forecasts.
• The forecast distribution describes the probability of observing any future value.

• A prediction interval gives a region within which we expect \(y_{T+h}\) to lie with a specified probability.
• Assuming forecast errors are normally distributed, then a 95% PI is

\[\hat{y}_{T+h|T} \pm 1.96 \hat\sigma_h\]

where \(\hat\sigma_h\) is the st dev of the \(h\)-step distribution.

• When \(h=1\), \(\hat\sigma_h\) can be estimated from the residuals.

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

\[ \hat{\sigma} = \sqrt{\frac{1}{T-K-M}\sum_{t=1}^{T}e_t^2}\] - K is the number of parameters estimated in the forecasting method,

• M is the number of missing values in the residuals. (For example, M = 1 for a naive forecast, because we can’t forecast the first observation.)

• When T is large enough, T and K does not matter*

• For the four benchmark methods, it is possible to mathematically derive the forecast standard deviation under the assumption of uncorrelated residuals.

• If the \(\hat{\sigma}\) is the residual standard deviation (h=1), then \(\hat{\sigma_h}\) is:

• m seasonal period
• k is the integer part of \(\frac{h-1}{m}\)

brick_fc <- brick_fit %>%
  forecast(h = 1)

brick_fc %>% hilo(level = 95)

## # A tsibble: 1 x 5 [1Q]
## # Key:       .model [1]
##   .model Quarter       Bricks .mean        `95%`
##   <chr>    <qtr>       <dist> <dbl>       <hilo>
## 1 Naive  2005 Q1 N(397, 1960)   397 [310, 484]95

• \(\hat{\sigma} = \sqrt{1960}\)

• So \(\hat{\sigma_2} = \sqrt{1960} \times \sqrt{1} = 44.27\)

• Interval \(397 \pm 1.96\times44.27 = 397 \pm 86.7 = [310-484]\)

brick_fc <- brick_fit %>%
  forecast(h = 2)

brick_fc %>% hilo(level = 95)

## # A tsibble: 2 x 5 [1Q]
## # Key:       .model [1]
##   .model Quarter       Bricks .mean        `95%`
##   <chr>    <qtr>       <dist> <dbl>       <hilo>
## 1 Naive  2005 Q1 N(397, 1960)   397 [310, 484]95
## 2 Naive  2005 Q2 N(397, 3920)   397 [274, 520]95

• Interval \(397 \pm 1.96\times\sqrt3920 = 397 \pm 122 = [274-519]\)

• Point forecasts often useless without a measure of uncertainty (such as prediction intervals).
• Prediction intervals require a stochastic model (with random errors, etc) - Bootstrap.
• For most models, prediction intervals get wider as the forecast horizon increases.
• Use level argument to control coverage.
• Check residual assumptions before believing them.
• Usually too narrow due to unaccounted uncertainty.

• We observe some drift, and we could smooth the data by log it

fit <- eggs %>%
  model(RW(log(eggs) ~ drift()))
fit

## # A mable: 1 x 1
##   `RW(log(eggs) ~ drift())`
##                     <model>
## 1             <RW w/ drift>

fc <- fit %>%
  forecast(h = 50)
fc

## # A fable: 50 x 4 [1Y]
## # Key:     .model [1]
##    .model                   year             eggs .mean
##    <chr>                   <dbl>           <dist> <dbl>
##  1 RW(log(eggs) ~ drift())  1994 t(N(4.1, 0.018))  61.8
##  2 RW(log(eggs) ~ drift())  1995 t(N(4.1, 0.036))  61.4
##  3 RW(log(eggs) ~ drift())  1996 t(N(4.1, 0.054))  61.0
##  4 RW(log(eggs) ~ drift())  1997 t(N(4.1, 0.073))  60.5
##  5 RW(log(eggs) ~ drift())  1998 t(N(4.1, 0.093))  60.1
##  6 RW(log(eggs) ~ drift())  1999    t(N(4, 0.11))  59.7
##  7 RW(log(eggs) ~ drift())  2000    t(N(4, 0.13))  59.3
##  8 RW(log(eggs) ~ drift())  2001    t(N(4, 0.15))  58.9
##  9 RW(log(eggs) ~ drift())  2002    t(N(4, 0.17))  58.6
## 10 RW(log(eggs) ~ drift())  2003    t(N(4, 0.19))  58.2
## # … with 40 more rows

fc %>% autoplot(eggs) +
  labs(title="Annual egg prices",
       y="US$ (adjusted for inflation)")

• Back-transformed point forecasts are medians.

• Back-transformed PI have the correct coverage.

• For many purposes, this is acceptable, although the mean is usually preferable ( medians do not add up, whereas means do).

• we can adjust the forecast back to the mean (bias adjusted)

fc %>%
  autoplot(eggs, level = c(80, 90), point_forecast = lst(mean, median)) +
  labs(title="Annual egg prices",
       y="US$ (adjusted for inflation)")

\[y_t = \hat{S}_t + \hat{A}_t\]

• \(\hat{A}_t\) is seasonally adjusted component \((\hat{A}_t =\hat{T}_t + \hat{R}_t)\)

• \(\hat{S}_t\) is seasonal component.

• Forecast both term separately.

• The seasonal component is (usually) unchanging (or changing extremely slowly) - forecast by simply taking the last year of the estimated component (a seasonal naïve method).

• Forecast \(\hat{A}_t\) using non-seasonal time series method.

• Combine forecasts of \(\hat{S}_t\) and \(\hat{A}_t\) to get forecasts of original data.

us_retail_employment <- us_employment %>%
  filter(year(Month) >= 1990, Title == "Retail Trade") %>%
  select(-Series_ID)

autoplot(us_retail_employment)

dcmp <- us_retail_employment %>%
  model(STL(Employed)) %>%
  components() %>% select(-.model)
dcmp

## # A tsibble: 357 x 6 [1M]
##       Month Employed  trend season_year remainder season_adjust
##       <mth>    <dbl>  <dbl>       <dbl>     <dbl>         <dbl>
##  1 1990 Jan   13256. 13288.      -33.0      0.836        13289.
##  2 1990 Feb   12966. 13269.     -258.     -44.6          13224.
##  3 1990 Mar   12938. 13250.     -290.     -22.1          13228.
##  4 1990 Apr   13012. 13231.     -220.       1.05         13232.
##  5 1990 May   13108. 13211.     -114.      11.3          13223.
##  6 1990 Jun   13183. 13192.      -24.3     15.5          13207.
##  7 1990 Jul   13170. 13172.      -23.2     21.6          13193.
##  8 1990 Aug   13160. 13151.       -9.52    17.8          13169.
##  9 1990 Sep   13113. 13131.      -39.5     22.0          13153.
## 10 1990 Oct   13185. 13110.       61.6     13.2          13124.
## # … with 347 more rows

• For the seasonally adjusted component, any non-seasonal forecasting method may be used. We will use the naive method in this example.

dcmp %>%
  model(NAIVE(season_adjust)) %>%
  forecast() %>%
  autoplot(dcmp) +
  labs(title = "Naive forecasts of seasonally adjusted data")

• The previous graphs shows naïve forecasts of the seasonally adjusted US retail employment data.

• These are then “reseasonalised” by adding in the seasonal naïve forecasts of the seasonal component.

• We can use the function decomposition_model() instead - compute forecasts via any additive decomposition, using model functions to forecast each of the decomposition’s components.

• Seasonal components will be forecast automatically using SNAIVE() (standard)

• The function will also reseasonalising for you.

• The variances from both the seasonally adjusted and seasonal forecasts are combined (predict intervals)

us_retail_employment %>%
  model(stlf = decomposition_model(
    STL(Employed ~ trend(window = 7), robust = TRUE),
    NAIVE(season_adjust)
  )) %>%
  forecast() %>%
  autoplot(us_retail_employment)

• So we divide the data into two sub-sample - training and test.

• The test set must not be used for any aspect of model development or calculation of forecasts.

• Forecast accuracy is based only on the test set.

• It should be at least as large as the maximum forecast horizon required

Forecast “error”: the difference between an observed value and its forecast.

\[ e_{T+h} = y_{T+h} - \hat{y}_{T+h|T}, \]

• training data {\(y_{1},..., y_T\)}; test data {\(y_{T+1}, y_{T+2}, ...\)}
• Residuals - Comes from training data, while forecast errors comes from test data
• Unlike residuals, forecast errors on the test set involve multi-step forecasts.
• These are true forecast errors as the test data is not used in computing \(\hat{y}_{T+h|T}\).

train <- aus_production %>%
  filter(between(year(Quarter), 1992, 2007))

beer <- aus_production %>%
  filter(year(Quarter) >= 1992)

beer_fc_plot <- train %>%
  model(
    Mean = MEAN(Beer),
    Naive = NAIVE(Beer),
    Seasonal_naive = SNAIVE(Beer),
    Drift = RW(Beer ~ drift())
  ) %>% forecast(h=10) %>% autoplot(beer, level = NULL) +
    labs(title = "Forecasts for quarterly beer production",
         y = "Megalitres") +
    guides(colour=guide_legend(title="Forecast"))

beer_fc_plot

• What is the most accurate model? How can we tell?

• For forecast (forecast errors as well) on the same scale data, we look at:

Mean Absolute error: \(MAE = mean(|e_t|)\)

Root mean squared error: \(RMSE = \sqrt{mean({e_t}^2)}\)

• MAE is easily interpretable, but it less popular
• RMSE is more sensitive to outliers (penalise large errors)

• For forecast on different scale data (unit-free), we look at:

Mean absolute percentage error: \[MAPE= 100\times mean(|e_{T+h}|/ |y_{T+h}|)\]

• The MAPE could be undefined if \(y=0\), or extreme values if \(y\) close to zero.

Mean Absolute Scaled Error \[ MASE = \text{mean}(|e_{T+h}|/Q) \] where \(Q\) is a stable measure of the scale of the time series \(\{y_t\}\).

For non-seasonal time series, \[ Q = (T-1)^{-1}\sum_{t=2}^T |y_t-y_{t-1}| \]

works well. Then MASE is equivalent to MAE relative to a naïve method.

Mean Absolute Scaled Error

\[ MASE = mean(|e_{T+h}|/Q) \]

where \(Q\) is a stable measure of the scale of the time series \(\{y_t\}\).

For seasonal time series, \[ Q = (T-m)^{-1}\sum_{t=m+1}^T |y_t-y_{t-m}| \]

Then MASE is equivalent to MAE relative to a seasonal naïve method.

If it is less than one - Better forecast than the average one-step naïve forecast. If it is more than one - Worse than the average one-step naïve forecast

recent_production <- aus_production %>%
  filter(year(Quarter) >= 1992)

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

beer_fit <- train %>%
  model(
    Mean = MEAN(Beer),
    Naive = NAIVE(Beer),
    Seasonal_naive = SNAIVE(Beer),
    Drift = RW(Beer ~ drift())
  )
beer_fc <- beer_fit %>%
  forecast(h = 10)

• We also can evaluate the distributional forecast accuracy
• We can evaluate a given quantile forecast - Quantile scores
• Evaluate a given predict interval - Winkler Score
• Evaluate the whole forecast distribution - Continuous Ranked Probability Score or CRPS

• A more sophisticated version of training/test sets is time series cross-validation.

• There are a series of test sets, each consisting of a single observation.

• The corresponding training set consists only of all observations that occurred prior to the observation that forms the test set.

• Forecast accuracy averaged over test sets.

• Also known as “evaluation on a rolling forecasting origin”

• We can also produce cross-validation for multi-step forecast (\(h=2,3,4,...\))