DASC6510/DASC4990

Unit 2c: Forecasting Tools

Erfanul Hoque, PhD

Thompson Rivers University
The note is strongly inspired by the materials shared on the Book:
Hyndman, R. J. & Athanasopoulos, G. (2021) Forecasting:
principles and practice, 3rd edition. [Link]

2
Forcaster’s toolbox

We discuss some general tools that are useful for forecasting.

• We will describe some benchmark forecasting methods,

procedures for checking whether a forecasting method has
adequately utilised the available information, and methods for
evaluating forecast accuracy.

3
A forecasting workflow

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

4
Data preparation (tidy)

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.
5
## 7 1966 Afghanistan 1399999967. 10152331 138.
Data visualisation

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

GDP per capita for Canada

50000

40000

30000
$US

20000

10000

0
1960 1980 2000
Year [1Y]

6
Model estimation

The model() function trains models to data.

fit <- gdppc %>%

model(trend_model = TSLM(GDP_per_capita ~ trend()))
fit

## # 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
Producing forecasts

fit %>% forecast(h = "3 years")

## # A fable: 789 x 5 [1Y]

## # Key: Country, .model [263]
## Country .model Year GDP_per_capita .mean
## <fct> <chr> <dbl> <dist> <dbl>
## 1 Afghanistan trend_model 2018 N(526, 9653) 526.
## 2 Afghanistan trend_model 2019 N(534, 9689) 534.
## 3 Afghanistan trend_model 2020 N(542, 9727) 542.
## 4 Albania trend_model 2018 N(4716, 476419) 4716.
## 5 Albania trend_model 2019 N(4867, 481086) 4867.
## 6 Albania trend_model 2020 N(5018, 486012) 5018.
## 7 Algeria trend_model 2018 N(4410, 643094) 4410.
## 8 Algeria trend_model 2019 N(4489, 645311) 4489.
## 9 Algeria trend_model 2020 N(4568, 647602) 4568.
8
## 10 American Samoa trend_model 2018 N(12491, 652926) 12491.
Visualising forecasts

fit %>% forecast(h = "3 years") %>%

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

GDP per capita for Sweden

60000

level
40000
$US

80
95
20000

0
1960 1980 2000 2020
Year

9
Some simple forecasting methods

10
Some simple forecasting methods

MEAN(y): Average method

• Forecast of all future values is equal to mean of

historical data {y1 , . . . , yT }.
• Forecasts: ŷT +h|T = ȳ = (y1 + · · · + yT )/T
Clay brick production in Australia
600

500
Bricks

400

300

200

1960 Q1 1980 Q1 2000 Q1

Quarter
11
Some simple forecasting methods

NAIVE(y): Naïve method

• Forecasts equal to last observed value.

• Forecasts: ŷT +h|T = yT .
• Consequence of efficient market hypothesis.
Clay brick production in Australia
500

450
Bricks

400

350

300
1995 Q1 2000 Q1 2005 Q1 2010 Q1
Quarter

12
Some simple forecasting methods

SNAIVE(y ~ lag(m)): Seasonal naïve method

• Forecasts equal to last value from same season.

• Forecasts: ŷT +h|T = yT +h−m(k+1) , where m =
seasonal period and k is the integer part of
(h − 1)/m.
Clay brick production in Australia
500

450
Bricks

400

350

300
1995 Q1 2000 Q1 2005 Q1 2010 Q1
Quarter

13
Model fitting

The model() function trains models to data.

brick_fit <- aus_production %>%

filter() %>%
model(
Seasonal_naive = SNAIVE(Bricks),
Naive = NAIVE(Bricks),
Mean = MEAN(Bricks)
)

## # A mable: 1 x 3
## Seasonal_naive Naive Mean
## <model> <model> <model>
## 1 <SNAIVE> <NAIVE> <MEAN>

14
Producing forecasts

brick_fc <- brick_fit %>%

forecast(h = "5 years")

## # A fable: 60 x 4 [1Q]
## # Key: .model [3]
## .model Quarter Bricks .mean
## <chr> <qtr> <dist> <dbl>
## 1 Seasonal_naive 2005 Q3 N(428, 2336) 428
## 2 Seasonal_naive 2005 Q4 N(397, 2336) 397
## 3 Seasonal_naive 2006 Q1 N(355, 2336) 355
## 4 Seasonal_naive 2006 Q2 N(435, 2336) 435
## # ... with 56 more rows

15
Visualising forecasts

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

Clay brick production in Australia

600

500
Millions of bricks

Forecast
400 Mean
Naive
Seasonal_naive
300

200

1960 Q1 1980 Q1 2000 Q1

Quarter

16
Facebook closing stock price

# Extract training data

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

# Specify, estimate and forecast

fb_stock %>%
model(
Mean = MEAN(Close),
Naive = NAIVE(Close)
) %>%
forecast(h=42) %>%
autoplot(fb_stock, level = NULL) +
labs(title = "Facebook closing stock price", y="$US") +
guides(colour=guide_legend(title="Forecast"))

17
Facebook closing stock price

Facebook closing stock price

200

150 Forecast
$US

Mean
Naive

100

50
0 500 1000
trading_day

18
Exercise to do!

• Produce forecasts using an appropriate benchmark method for

household wealth (hh_budget). Plot the results using
autoplot().
• Produce forecasts using an appropriate benchmark method for
Australian takeaway food turnover (aus_retail). Plot the
results using autoplot().

19
Residual diagnostics

20
Fitted values

• ŷt|t−1 is the forecast of yt based on observations y1 , . . . , yt−1

(we call these “fitted values”).
• Sometimes we write it as: ŷt ≡ ŷt|t−1 .
• Often not true forecasts since parameters are estimated on all
data.

For example:
• ŷt = ȳ for average method.

21
Forecasting residuals

Residuals in forecasting: difference between observed value and

its fitted value: et = yt − ŷt|t−1 .

Assumptions:
1. et ’s are uncorrelated. If they aren’t, then information
left in residuals that should be used in computing
forecasts.
2. et ’s have mean zero. If they don’t, then forecasts are
biased.
Properties (for distributions & prediction intervals):
3. et ’s have constant variance.
4. et ’s are normally distributed.
22
Facebook closing stock price

fb_stock <- gafa_stock %>%

filter(Symbol == "FB") %>%
mutate(trading_day = row_number()) %>%
update_tsibble(index = trading_day, regular = TRUE)
fb_stock %>% autoplot(Close)

200

150
Close

100

50
0 400 800 1200
trading_day [1]
23
Facebook closing stock price
fit <- fb_stock %>% model(NAIVE(Close))
augment(fit)

## # A tsibble: 1,258 x 7 [1]

## # Key: Symbol, .model [1]
## Symbol .model trading~1 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, and abbreviated variable name
## # 1: trading_day

Naïve forecasts:
ŷt|t−1 = yt−1
et = yt − ŷt|t−1 = yt − yt−1 24
Facebook closing stock price

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

200

150 colour
Close

Data
Fitted

100

50
0 400 800 1200
trading_day
25
Facebook closing stock price

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"))

200

colour
Close

175
Data
Fitted

150

125

1100 1150 1200 1250

trading_day
26
Facebook closing stock price

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

Residuals from naïve method

0
$US

−20

−40

0 400 800 1200

trading_day [1]

27
Facebook closing stock price

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

Histogram of residuals
150

100
count

50

0
−40 −20 0
.resid
28
Facebook closing stock price

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

ACF of residuals
0.06

0.03
acf

0.00

−0.03

−0.06
10 20 30
lag [1]

29
gg_tsresiduals() function

gg_tsresiduals(fit)
Innovation residuals

−20

−40
0 400 800 1200
trading_day
0.06 150

0.03
100
count
acf

0.00
50
−0.03

−0.06 0
10 20 30 −40 −20 0
lag [1] .resid

30
ACF of residuals

• 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.

31
Portmanteau tests for autocorrelation

Consider a whole set of ρk values, and develop a test to see

whether the set is significantly different from a zero set.
Box-Pierce test
X
ℓ
Q=T ρ2k
k=1
where ℓ is max lag being considered and T is number of
observations.

• If each ρk close to zero, Q will be small.

• If some ρk values large (positive or negative), Q will be large.

32
Portmanteau tests

Consider a whole set of ρk values, and develop a test to see

whether the set is significantly different from a zero set.
Ljung-Box test
X
ℓ
Q ∗ = T (T + 2) (T − k)−1 ρ2k
k=1
where ℓ is max lag being considered and T is number of
observations.

• My preferences: ℓ = 10 for non-seasonal data, ℓ = 2m for

seasonal data.
• Better performance, especially in small samples.

33
Portmanteau tests

• If data are WN, Q ∗ has χ2 distribution with (ℓ − K ) degrees of

freedom (dof) where lag = ℓ, K = no. parameters in model.
• When applied to raw data, set K = 0.

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

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

For Q ∗ , the results are not significant (i.e., the p-values are relatively large).
Thus, we can conclude that the residuals are not distinguishable from a white
noise series.
34
Exercise to do!

Compute seasonal naïve forecasts for quarterly Australian beer

production from 1992.

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

fit <- recent %>% model(SNAIVE(Beer))
fit %>% forecast() %>% autoplot(recent)

Test if the residuals are white noise.

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

gg_tsresiduals(fit)

What do you conclude?

35
Forecasts distributions and
prediction intervals

36
Forecast distributions

• A forecast ŷT +h|T is (usually) the mean of the conditional

distribution yT +h | y1 , . . . , yT .
• Most time series models produce normally distributed
forecasts.
• The forecast distribution describes the probability of observing
any future value.

37
Forecast distributions

Assuming residuals are normal, uncorrelated, sd = σ̂:

Mean: ŷT +h|T ∼ N(ȳ , (1 + 1/T )σ̂ 2 )

Naïve: ŷT +h|T ∼ N(yT , hσ̂ 2 )
Seasonal naïve: ŷT +h|T ∼ N(yT +h−m(k+1) , (k + 1)σ̂ 2 )
where k is the integer part of (h − 1)/m.
Note that when h = 1 and T is large, these all give the same
approximate forecast variance: σ̂ 2 .

38
Prediction intervals

• A prediction interval gives a region within which we expect

yT +h to lie with a specified probability.
• Assuming forecast errors are normally distributed, then a 95%
PI is

ŷT +h|T ± 1.96σ̂h

where σ̂h is the st dev of the h-step distribution.
• When h = 1, σ̂h can be estimated from the residuals.

v
u
u 1 XT
σ̂ = t e2,
T − K − M t=1 t

where K is the number of parameters, and M is the number of

missing values in the residuals (e.g. M=1 for a naive forecast,
because we can’t forecast the first observation)
39
Prediction intervals

brick_fc %>% hilo(level = 95)

## # A tsibble: 60 x 5 [1Q]
## # Key: .model [3]
## .model Quarter Bricks .mean `95%`
## <chr> <qtr> <dist> <dbl> <hilo>
## 1 Seasonal_naive 2005 Q3 N(428, 2336) 428 [333, 523]95
## 2 Seasonal_naive 2005 Q4 N(397, 2336) 397 [302, 492]95
## 3 Seasonal_naive 2006 Q1 N(355, 2336) 355 [260, 450]95
## 4 Seasonal_naive 2006 Q2 N(435, 2336) 435 [340, 530]95
## 5 Seasonal_naive 2006 Q3 N(428, 4672) 428 [294, 562]95
## 6 Seasonal_naive 2006 Q4 N(397, 4672) 397 [263, 531]95
## 7 Seasonal_naive 2007 Q1 N(355, 4672) 355 [221, 489]95
## 8 Seasonal_naive 2007 Q2 N(435, 4672) 435 [301, 569]95
## 9 Seasonal_naive 2007 Q3 N(428, 7008) 428 [264, 592]95
## 10 Seasonal_naive 2007 Q4 N(397, 7008) 397 [233, 561]95
40
## # ... with 50 more rows
Prediction intervals

• Point forecasts often useless without a measure of uncertainty

(such as prediction intervals).
• Prediction intervals require a stochastic model (with random
errors, etc).
• 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.

41
Evaluating forecast accuracy

42
Training and test sets

Training data Test data

time

• A model which fits the training data well will not necessarily
forecast well.
• A perfect fit can always be obtained by using a model with
enough parameters.
• Over-fitting a model to data is just as bad as failing to
identify a systematic pattern in the data.
• 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.

43
Forecast errors

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

forecast.
eT +h = yT +h − ŷT +h|T ,

where the training data is given by {y1 , . . . , yT }

• 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 ŷT +h|T .

44
Measures of forecast accuracy

Forecasts for quarterly beer production

500

Forecast
Megalitres

Mean
450 Naive
Seasonal_naive

400

1995 Q1 2000 Q1 2005 Q1 2010 Q1

Quarter

45
Measures of forecast accuracy

yT +h = (T + h)th observation, h = 1, . . . , H
ŷT +h|T = its forecast based on data up to time T .
eT +h = yT +h − ŷT +h|T

MAE = mean(|eT +h |) q
MSE = mean(eT2 +h ) RMSE = mean(eT2 +h )
MAPE = 100mean(|eT +h |/|yT +h |)

• MAE, MSE, RMSE are all scale dependent.

• MAPE is scale independent but is only sensible if yt ≫ 0 for
all t, and y has a natural zero.

46
Measures of forecast accuracy

Mean Absolute Scaled Error

MASE = mean(|eT +h |/Q)

where Q is a stable measure of the scale of the time series {yt }.

Proposed by Hyndman and Koehler (IJF, 2006).
For non-seasonal time series,

X
T
−1
Q = (T − 1) |yt − yt−1 |
t=2

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

method.

47
Measures of forecast accuracy

Mean Absolute Scaled Error

MASE = mean(|eT +h |/Q)

where Q is a stable measure of the scale of the time series {yt }.

Proposed by Hyndman and Koehler (IJF, 2006).
For seasonal time series,

X
T
−1
Q = (T − m) |yt − yt−m |
t=m+1

works well. Then MASE is equivalent to MAE relative to a

seasonal naïve method.

48
Measures of forecast accuracy

Forecasts for quarterly beer production

500

Forecast
Megalitres

Drift

450 Mean
Naive
Seasonal_naive

400

1995 Q1 2000 Q1 2005 Q1 2010 Q1

Quarter

49
Measures of forecast accuracy

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)
)
beer_fc <- beer_fit %>%
forecast(h = 10)

50
Measures of forecast accuracy

accuracy(beer_fit)

## # A tibble: 4 x 6
## .model .type RMSE MAE MAPE MASE
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Drift Training 65.3 54.8 12.2 3.83
## 2 Mean Training 43.6 35.2 7.89 2.46
## 3 Naive Training 65.3 54.7 12.2 3.83
## 4 Seasonal_naive Training 16.8 14.3 3.31 1

accuracy(beer_fc, recent_production)

## # A tibble: 4 x 6
## .model .type RMSE MAE MAPE MASE
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Drift Test 64.9 58.9 14.6 4.12
## 2 Mean Test 38.4 34.8 8.28 2.44
## 3 Naive Test 62.7 57.4 14.2 4.01
## 4 Seasonal_naive Test 14.3 13.4 3.17 0.937

51
Time series cross-validation

52
Time series cross-validation
Traditional evaluation
Training data Test data
time

Time series cross-validation

h=1

time
53
Time series cross-validation
Traditional evaluation
Training data Test data
time

Time series cross-validation

h=2

time
54
Time series cross-validation
Traditional evaluation
Training data Test data
time

Time series cross-validation

h=3

time
55
Time series cross-validation
Traditional evaluation
Training data Test data
time

Time series cross-validation

h=4

time

Also known as "evaluation on a rolling forecasting origin" 56

Time series cross-validation

Stretch with a minimum length of 3, growing by 1 each step.

fb_stretch <- fb_stock %>%
stretch_tsibble(.init = 3, .step = 1) %>%
filter(.id != max(.id))

## # A tsibble: 790,650 x 4 [1]

## # Key: .id [1,255]
## Date Close trading_day .id
## <date> <dbl> <int> <int>
## 1 2014-01-02 54.7 1 1
## 2 2014-01-03 54.6 2 1
## 3 2014-01-06 57.2 3 1
## 4 2014-01-02 54.7 1 2
## 5 2014-01-03 54.6 2 2
## 6 2014-01-06 57.2 3 2
## 7 2014-01-07 57.9 4 2
## # ... with 790,643 more rows 57
Time series cross-validation

Estimate Naïve models for each window.

fit_cv1 <- fb_stretch %>%

model(NAIVE(Close))

## # A mable: 1,255 x 3
## # Key: .id, Symbol [1,255]
## .id Symbol `NAIVE(Close)`
## <int> <chr> <model>
## 1 1 FB <NAIVE>
## 2 2 FB <NAIVE>
## 3 3 FB <NAIVE>
## 4 4 FB <NAIVE>
## # ... with 1,251 more rows

58
Time series cross-validation

Produce one step ahead forecasts from all models.

fc_cv1 <- fit_cv1 %>%

forecast(h=1)

## # A fable: 1,255 x 5 [1]

## # Key: .id, Symbol [1,255]
## .id Symbol trading_day Close .mean
## <int> <chr> <dbl> <dist> <dbl>
## 1 1 FB 4 N(57, 3.5) 57.2
## 2 2 FB 5 N(58, 2.5) 57.9
## 3 3 FB 6 N(58, 1.9) 58.2
## 4 4 FB 7 N(57, 1.7) 57.2
## # ... with 1,251 more rows

59
Time series cross-validation

# Cross-validated
fc_cv1 %>% accuracy(fb_stock)
# Training set
fb_stock %>% model(NAIVE(Close))%>% accuracy()

RMSE MAE MAPE

Cross-validation 2.418 1.469 1.266
Training 2.414 1.468 1.264

A good way to choose the best forecasting model is to find the model
with the smallest RMSE computed using time series cross-validation.

60
Exercise for Class!

Are the following statements true or false? Explain your answer.

• A model with small residuals will give good forecasts.

61
Exercise for Class!

Are the following statements true or false? Explain your answer.

• A model with small residuals will give good forecasts.

• The best measure of forecast accuracy is MAPE.

61
Exercise for Class!

Are the following statements true or false? Explain your answer.

• A model with small residuals will give good forecasts.

• The best measure of forecast accuracy is MAPE.

• Always choose the model with the best forecast accuracy as

measured on the test set.

61
Exercise to do!

tourism contains quarterly visitor nights (in thousands) from 2000

to 2017 for 76 regions of Australia.

• Extract data from the Gold Coast region using filter() and
aggregate total overnight trips (sum over Purpose) using
summarise(). Call this new dataset gc_tourism.
• Using slice() or filter(), create three training sets for this
data excluding the last 1, 2 and 3 years. For example,
gc_train_1 <- gc_tourism |> slice(1:(n()-4)).
• Compute one year of forecasts for each training set using the
seasonal naïve ( SNAIVE() ) method. Call these gc_fc_1,
gc_fc_2 and gc_fc_3, respectively.
• Use accuracy() to compare the test set forecast accuracy
using MAPE. Comment on these.

62
Next Lecture!

• In the next lecture, we learn The Regression Model and

Exponential Smoothing

Please go to Chapter 7 & 8 of text book (Hyndman, R. J. &

Athanasopoulos, G. (2021) Forecasting: principles and practice,
3rd edition. [Link] beforehand.

63