Introduction

This module discusses the basics of time-series based forecasting.

The module is based on a free online textbook that you can access here

Hyndman, R.J., & Athanasopoulos, G. (2018). Forecasting: principles and practice, 2^nd edition. OTexts: Melbourne, Australia. OTexts.com/fpp2.

1. Getting started

Read chapter 1!

Points to take away:

The predictability of an event or a quantity depends on several factors including:

how well we understand the factors that contribute to it
how much data is available
whether the forecasts can affect the thing we are trying to forecast.

Forecasting is about predicting the future as accurately as possible, given all of the information available, including historical data and knowledge of any future events that might impact the forecasts.

Modern organisations require short-term, medium-term and long-term forecasts, depending on the specific application.

If there are no data available, or if the data available are not relevant to the forecasts, then qualitative forecasting methods must be used. These methods are not purely guesswork—there are well-developed structured approaches to obtaining good forecasts without using historical data. These methods are discussed in Chapter 4.

Quantitative forecasting can be applied when two conditions are satisfied:

numerical information about the past is available
it is reasonable to assume that the past patterns will continue into the future.

If there are no relevant data available then qualitative forecasting methods must be used.

Anything that is observed sequentially over time, is a time series. The aim of time series based forecasting, is to estimate how the sequence of observations will continue into the future.

The simplest time series forecasting methods use only information on the variable to be forecast. They will extrapolate trend and seasonal patterns, but ignore all other information. This module will focus on this type of forecasting methods:

Decomposition models
Exponential smoothing models, and
ARIMA models.

These models are discussed in Chapters 6-8, after introducing the forecaster's tools.

Explanatory models with predictor variables, and models that combine time series based forecaating models and explanatory models are discussed in a separate module on advanced techniques (Hyndman et al, chapter 9-11). These include:

Dynamic regression models
Panel data models and longitudinal models.

A forecasting task involves five steps.

Step 1: Problem definition
Step 2: Gathering information
Step 3: Preliminary (exploratory) analysis
Step 4: Choosing and fitting models
Step 5: Using and evaluating a forecasting model.

Do exercises 1 and 2, in section 1.8 of the textbook!

2. Time series graphics

The textbook comes with the fpp2 package, which can be installed in R. The package includes functions and data sets. Invoke the package using library().

You have to install it only once, on your computer. After installing a new version of R, you have to install all User Library packages, like fpp2, again. Check the installed packages in RStudio!

# install.packages("fpp2") # Install once; check "Packages in RStudio"
library(fpp2)

2.1 ts() function, and ts objects

The ts() function in R, creates time series. Following the example in the textbook:

y <- ts(c(123,39,78,52,110), start=2012)
plot(y)

autoplot(y)

Here is what you can do to create ts objects for monthly data. The data starts with 1, for the first month; 2, for the second month, up to 72 for the last months. The time series is just a straight line (e.g. a linear trend).

# save a numeric vector containing 72 monthly observations
myvector <- 1:72
# from Jan 2009 to Dec 2014 as a time series object
myts <- ts(myvector, start=c(2009, 1), end=c(2014, 12), frequency=12)
head(myts,24) # print the first 24 months

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2009   1   2   3   4   5   6   7   8   9  10  11  12
## 2010  13  14  15  16  17  18  19  20  21  22  23  24

autoplot(myts)

You can take a part of the time series using the windows() function. You can choose a start and an end period. By default, the starting period is the first observation, and the ending period is the last observation. In the example below, myts3 then is equivalent to myts2.

# subset the time series (June 2014 to December 2014)
myts2 <- window(myts, start=c(2013, 6), end=c(2014, 12))
myts3 <- window(myts, start=c(2013, 6))
myts2; myts3

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2013                      54  55  56  57  58  59  60
## 2014  61  62  63  64  65  66  67  68  69  70  71  72

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2013                      54  55  56  57  58  59  60
## 2014  61  62  63  64  65  66  67  68  69  70  71  72

# plot series
autoplot(myts2)

2.2 Time Plots

As above, you can use the autoplot() function for displaying time series.

The str() functions shows the structure of the data set, including the names of the three series.

The plot shows the three time series that are part of the melsyd data set. The plot will already make you wonder what happened in some periods!

str(melsyd)

##  Time-Series [1:283, 1:3] from 1987 to 1993: 1.91 1.85 1.86 2.14 2.12 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:3] "First.Class" "Business.Class" "Economy.Class"

melsyd[1:5,] # prints the first five observations

##      First.Class Business.Class Economy.Class
## [1,]       1.912             NA        20.167
## [2,]       1.848             NA        20.161
## [3,]       1.856             NA        19.993
## [4,]       2.142             NA        20.986
## [5,]       2.118             NA        20.497

autoplot(melsyd)

You can single out one time series, using indexing. make sure that you use the exact name of the time series. The names can be revealed by, among others, the str() function. You can also use colnames(), as below. You can add geom's from ggplot, to edit the graph.

colnames(melsyd)

## [1] "First.Class"    "Business.Class" "Economy.Class"

autoplot(melsyd[,"Economy.Class"]) +
  ggtitle("Economy class passengers: Melbourne-Sydney") +
  xlab("Year") +
  ylab("Thousands")

Another example:

autoplot(a10) +
  ggtitle("Antidiabetic drug sales") +
  ylab("$ million") +
  xlab("Year")

Apart from presenting the raw data to others, as a time series analyst and forecaster, you will need to plot the time series in order to select a proper model. We encourage you to always do so, even if the textbook and the module will introduce some automatic functions that select the best model for you.

2.3 Time Series Patterns

Make sure that you are familiar with the main components of time series, and their definitions:

Trend
Cycles
Seasonality
Error (random fluctuations; white noise)

2.4 Seasonal Plots

As one of the first steps in the analytic process, it really helps to display the time series using a variety of tools.

For the data in a10:

autoplot(a10)

Seasonal Plot

ggseasonplot(a10, year.labels=TRUE, year.labels.left=TRUE) +
  ylab("$ million") +
  ggtitle("Seasonal plot: antidiabetic drug sales")

Polar Seasonal Plot

ggseasonplot(a10, polar=TRUE) +
  ylab("$ million") +
  ggtitle("Polar seasonal plot: antidiabetic drug sales")

Both graphs show:

Months in which sales are relatively high (January) and low (February; March). Note that this is less clear from the plot of the series.
An increasing trend, as the colored lines in the seasonal plots are at a higher level, for more recent years, and the circles are expanding in the polar plot.

2.5 Seasonal Subseries Plots

A clearer way to detect seasonal patterns (here: the months of the year with low or high sales), is the seasonal subseries plot. For each month of the year, the mean is indicated with a horizontal line. Since for each and every month the time series is increasing (the January sales are higher for in 2007 than in 2006; and so on), there's a clear overall upward trend.

ggsubseriesplot(a10) +
  ylab("$ million") +
  ggtitle("Seasonal subseries plot: antidiabetic drug sales")

2.6 Scatterplots

Two series can be plotted separately in one graph using the facet option. Since the scales for the two series are totally different in this example (volumes demanded versus temperature), it doesn't make sense to combine them.

autoplot(elecdemand[,c("Demand","Temperature")], facets=TRUE) +
  xlab("Year: 2014") + ylab("") +
  ggtitle("Half-hourly electricity demand: Victoria, Australia")

Since demand for electricity will go up if the temperatures are especially low (heaters) or high (airco's), the two series are related in a U-shaped manner.

qplot(Temperature, Demand, data=as.data.frame(elecdemand)) +
  ylab("Demand (GW)") + xlab("Temperature (Celsius)")

Intermezzo - Anscombe data

The correlation coefficient measures the strength of the linear relationship, and can be misleading:

if the relationship is not linear. For example, the correlation for the electricity demand and temperature above is U-shaped.
in case of influential cases, like outliers.

Anscombe's quartet comprises four data sets that have identical correlation coefficients, but appear very different when graphed.

Each dataset consists of eleven $(x,y)$ points. They were constructed in 1973 by the statistician Francis Anscombe to demonstrate both the importance of graphing data before analyzing it and the effect of outliers and other influential observations on statistical properties.

Below, the Anscombe data are read from a csv-file on GitHub. We have graphed the data, and calculate the correlation coefficients for all pairs (x_i,y_i), (i=1..4).

The data set looks as follows.

library (readr)

## Warning: package 'readr' was built under R version 4.0.3

urlfile="https://raw.githubusercontent.com/ssmresearch/datasets/main/anscombecsv.csv"
anscombe<-as.data.frame(read_csv(url(urlfile)))

## 
## -- Column specification --------------------------------------------------------
## cols(
##   Observation = col_double(),
##   x1 = col_double(),
##   y1 = col_double(),
##   x2 = col_double(),
##   y2 = col_double(),
##   x3 = col_double(),
##   y3 = col_double(),
##   x4 = col_double(),
##   y4 = col_double()
## )

anscombe

The graphs are scatterplots, with fitted regression lines, using ggplot(). The correlation coefficient for (x₁,y₁) is 0.82.

anscombe[,2:3]                # x1 and y1 are in the columns 2 and 3

anscombe[,c("x1","y1")]       # is equivalent to the above

round(cor(anscombe[,2:3]),2)  # (Pearson) correlation, rounded in two decimals

##      x1   y1
## x1 1.00 0.82
## y1 0.82 1.00

library(ggplot2)
ggplot(anscombe, aes(x=x1, y=y1)) +
  geom_point(shape=1) +    # Use hollow circles
  geom_smooth(method=lm,   # Add linear regression line
              se=FALSE)    # Don't show confidence intervals

## `geom_smooth()` using formula 'y ~ x'

The regression lines and correlation coefficients are the same for all other pairs $(x,y)$!

The pair (x₂,y₂) is a non-linear relationship.

round(cor(anscombe[,4:5]),2)

##      x2   y2
## x2 1.00 0.82
## y2 0.82 1.00

ggplot(anscombe, aes(x=x2, y=y2)) +
  geom_point(shape=1) +    
  geom_smooth(method=lm, se=FALSE)

## `geom_smooth()` using formula 'y ~ x'

The pair (x₃,y₃) has an outlier: a very high value for $y$, for one observation.

round(cor(anscombe[,6:7]),2)

##      x3   y3
## x3 1.00 0.82
## y3 0.82 1.00

ggplot(anscombe, aes(x=x3, y=y3)) +
  geom_point(shape=1) +    
  geom_smooth(method=lm, se=FALSE)

## `geom_smooth()` using formula 'y ~ x'

The pair (x₄,y₄) too, has an outlier. Essentially $x$ has no variation for most observations except one. There is one observation with high values for both $x$ and $y$.

Such an example is not uncommon in economic research on small enterprises (SEs). While all SEs in the sample have a small number of employees (the $x$ variable), say in the range from 1 to 10, but diverging profits ($y$), one large company with a lot of employees and substantial profits may have entered our sample by mistake.

round(cor(anscombe[,8:9]),2)

##      x4   y4
## x4 1.00 0.82
## y4 0.82 1.00

ggplot(anscombe, aes(x=x4, y=y4)) +
  geom_point(shape=1) +    
  geom_smooth(method=lm, se=FALSE)

## `geom_smooth()` using formula 'y ~ x'

The important message is: graph your data before crunching numbers!

In case you have several series to analyze simultaneously, like in the data set visnights a matrix of scaterplots helps.

# install.packages("GGally")
library(GGally)
autoplot(visnights[,1:5], facets=TRUE) +
  ylab("Number of visitor nights each quarter (millions)")

GGally::ggpairs(as.data.frame(visnights[,1:5]))

2.7 Lag Plots

Lag plots can also help in identifying seasonal patterns.

For each quarter, in these quartely series, the values are plotted against lagged values. The positive correlations at lags 4 and 8 confirm the seasonal pattern, with peaks (purple) in Q4 and troughs in Q2.

The graph is somewhat confusing. You can use it for analysis (detection of seasonal patterns), but your readers (or managers) may not find it very insightful.

str(ausbeer)

##  Time-Series [1:218] from 1956 to 2010: 284 213 227 308 262 228 236 320 272 233 ...

beer2 <- window(ausbeer, start=1992)
gglagplot(beer2)

2.8 Autocorrelation

Autocorrelation is extremely important, especially in discussing and understanding the techniques in chapters 6 and following.

Autocorrelation measures the linear relationship between lagged values of a time series.

In the graph below, the quarterly seasonal pattern becomes very clear, as the correlations between values and lagged values, at lag 4, are significantly positive. The dashed blue line represents the boundary between significant correlations and random correlations. The negative correlations at lags 2, 6 and so on, are caused by peaks and troughs which are 2 quarters apart (peaks in Q4, and troughs in Q2).

ggAcf(beer2)

For data with a trend, the autocorrelations for small lags are large and positive because observations nearby in time are also nearby in size. The ACF of trended time series has positive values that slowly decrease as the lags increase.

When data are seasonal, the autocorrelations will be larger for the seasonal lags than for other lags.

The combined effect occurs in the

aelec <- window(elec, start=1980)
autoplot(aelec) + xlab("Year") + ylab("GWh")

ggAcf(aelec, lag=48)

2.9 White Noise

Time series that show no autocorrelation are called white noise.

Below, $y$ is a set of values randomly drawn from a normal distribution. The set.seed(#) sees to it that you can replicate the same graph. The # (here 30) is just the starting point in a series of random numbers, which are not that random after all!

Any other value then 30 will generate a somewhat different graph - but still white noise.

The ACF will show no or only a few significant correlations. Since significance is tested at 95% confidence, random numbers may produce a significant correlation in 1 out of 20 cases!

Try some other seed numbers (≠30), and longer series ($>50$) to see the effect.

set.seed(30)
y <- ts(rnorm(50))
autoplot(y) + ggtitle("White noise")

ggAcf(y)

2.10 Exercises

Complete exercises 2.3; 2.7 and 2.8 from the textbook.

3. The Forecaster's Toolbox

3.1 Simple Forecasting Methods

In the average method forecasts of all future values are equal to the average (or mean) of the historical data.

For naïve forecasts, we simply set all forecasts to be the value of the last observation.

A similar method is useful for seasonal data. Here, we set each forecast to be equal to the last observed value from the same season of the year. The function snaive() is not applied to the small series of annual data, below.

A variation on the naïve method is to allow the forecasts to increase or decrease over time, where the amount of change over time (called the drift) is set to be the average change observed in the historical data.

y <- ts(c(123,39,78,52,110), start=2012)
autoplot(y)

meanf(y, h=2)         # Average method

##      Point Forecast    Lo 80    Hi 80     Lo 95   Hi 95
## 2017           80.4 19.74314 141.0569 -29.44201 190.242
## 2018           80.4 19.74314 141.0569 -29.44201 190.242

naive(y, h=2)         # Naive method

##      Point Forecast    Lo 80    Hi 80        Lo 95    Hi 95
## 2017            110 38.02459 181.9754  -0.07688897 220.0769
## 2018            110  8.21140 211.7886 -45.67222928 265.6722

rwf(y, 2)             # Naive method, alternative experession

##      Point Forecast    Lo 80    Hi 80        Lo 95    Hi 95
## 2017            110 38.02459 181.9754  -0.07688897 220.0769
## 2018            110  8.21140 211.7886 -45.67222928 265.6722

rwf(y, 2, drift=TRUE) # Naive method with drift

##      Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## 2017         106.75  23.77923 189.7208 -20.14285 233.6428
## 2018         103.50 -27.68831 234.6883 -97.13521 304.1352

In a more elaborate example, we use the ausbeer data set in fpp2.

The quarterly data start in 1956 (Q1) and end in 2010 (Q2). An easy way to check this, is by using the head() and tail() functions to display the first and last (here, 4) observations.

We use windows() to create a subset of fairly recent data starting in 1992 and ending in 2007 to train our models. Data after 2007 will be used to test the data.

The autolayer() adds the various forecasts to the training data.

head(ausbeer,4);tail(ausbeer,4)

##      Qtr1 Qtr2 Qtr3 Qtr4
## 1956  284  213  227  308

##      Qtr1 Qtr2 Qtr3 Qtr4
## 2009            419  488
## 2010  414  374

# Set training data from 1992 to 2007
beer2 <- window(ausbeer,start=1992,end=c(2007,4))
# Plot some forecasts
autoplot(beer2) +
  autolayer(meanf(beer2, h=11),
            series="Mean", PI=FALSE) +
  autolayer(naive(beer2, h=11),
            series="Na?ve", PI=FALSE) +
  autolayer(snaive(beer2, h=11),
            series="Seasonal na?ve", PI=FALSE) +
  ggtitle("Forecasts for quarterly beer production") +
  xlab("Year") + ylab("Megalitres") +
  guides(colour=guide_legend(title="Forecast"))

The data from 2008 to 2010 (Q2) are known. They can be used as a test set. Displaying training and test data, along with the forecasts, follows the following code.

Obviously, given the seasonal pattern, the forecasts of the seasonally naive method outperform those from the other methods. Still, the forecasts for the low seasons seem to be too low. It may be that using information from previous years, rather than from the last year only, might have resulted in more accurate forecasts.

These simple methods are often used to benchmark the forecasts from more complex methods. The use of complex methods is only justified if they clearly and systematically outperform simple methods.

beer3 <- window(ausbeer,start=1992) # Training plus test data
# Plot some forecasts
autoplot(beer3) +                   # Plot training plus test data
  autolayer(meanf(beer2, h=11),
            series="Mean", PI=FALSE) +
  autolayer(naive(beer2, h=11),
            series="Na?ve", PI=FALSE) +
  autolayer(snaive(beer2, h=11),
            series="Seasonal na?ve", PI=FALSE) +
  ggtitle("Forecasts for quarterly beer production") +
  xlab("Year") + ylab("Megalitres") +
  guides(colour=guide_legend(title="Forecast"))

Using data on Google closing stock prices:

Which of the methods you think is best?

autoplot(goog200) +
  autolayer(meanf(goog200, h=40),
            series="Mean", PI=FALSE) +
  autolayer(rwf(goog200, h=40),
            series="Naïve", PI=FALSE) +
  autolayer(rwf(goog200, drift=TRUE, h=40),
            series="Drift", PI=FALSE) +
  ggtitle("Google stock (daily ending 6 Dec 2013)") +
  xlab("Day") + ylab("Closing Price (US$)") +
  guides(colour=guide_legend(title="Forecast"))

3.2 Transformations

Calendar Adjustments

Some of the variation in seasonal data may be due to simple calendar effects.

In such cases, it is usually much easier to remove the variation before fitting a forecasting model. The monthdays() function will compute the number of days in each month or quarter.

milk2 <- cbind(Monthly = milk, Days = monthdays(milk), 
         DailyAverage = milk/monthdays(milk))
autoplot(milk2, facet=TRUE) +
    xlab("Years") + ylab("Pounds") +
    ggtitle("Milk production per cow")

Population Adjustments

Any data affected by population changes can be adjusted to give per-capita data.

Inflation Adjustments

Data which are affected by the value of money are best adjusted before modelling.

To make these adjustments, a price index is used.

Mathematical Transformations

If the data show variation that increases or decreases with the level of the series, then a transformation can be useful. Logarithmic transformations are quite common, especially in economic data.

A specific family of transformations that includes both logarithms and power transformations, are Box-Cox transformations which depend on a parameter $λ$.

The logarithm in a Box-Cox transformation is always a natural logarithm (base $e$, $e=2.71828..$).

If $λ=0$, natural logarithms are used. If $λ≠0$, a power transformation is used.

The $BoxCox()$ computes an optimal value for $λ$. The outcome will be a very precise number (e.g. $λ=0.2654076$. Since the impact of small changes in $λ$ on the results are negligible, most researchers use a rounded number (e.g. $λ=0.25$ or $λ=0.30$).

autoplot(elec)

(lambda <- BoxCox.lambda(elec))

## [1] 0.2654076

An issue with using mathematical transformations such as Box-Cox transformations is that the back-transformed point forecast will not be the mean of the forecast distribution. This is called bias.

Forecasts based on transformed data have to be adjusted for this bias. Bias adjustment is not done by default in the forecast package. Whenever applicable, use the argument biasadj=TRUE. We will come back to this later.

To see the impact of bias adjustments:

fc <- rwf(eggs, drift=TRUE, lambda=0, h=50, level=80)
fc2 <- rwf(eggs, drift=TRUE, lambda=0, h=50, level=80,
           biasadj=TRUE)
autoplot(eggs) +
  autolayer(fc, series="Simple back transformation") +
  autolayer(fc2, series="Bias adjusted", PI=FALSE) +
  guides(colour=guide_legend(title="Forecast"))

3.3 Residual Diagnostics

Residuals (the difference between the observations and the model-based predictions) can be obtained using the res() function.

Ideally, the residuals:

Are normally distributed, with a mean of zero (to be checked with a histogram)
Are not autocorrelated, to be checked with an ACF of the residual series.

autoplot(goog200) +
  xlab("Day") + ylab("Closing Price (US$)") +
  ggtitle("Google Stock (daily ending 6 December 2013)")

res <- residuals(naive(goog200))
autoplot(res) + xlab("Day") + ylab("") +
  ggtitle("Residuals from naive method")

gghistogram(res) + ggtitle("Histogram of residuals")

## Warning: Removed 1 rows containing non-finite values (stat_bin).

ggAcf(res)+ggtitle("ACF of residuals")

Portmanteau Tests

Rather than inspecting the ACF to detect significant correlations at specific lags, it is better to test a set of autocorrelations at various lags simultaneously.

There are various such portmanteau tests, of which the Ljung-Box Test seems to be the most accurate. We are testing whether there are significant correlations in the group. Small $p$-values would indicate that this is the case. A value $p>0.05$ indicates that there are no significant autocorrelations, and that the series of residuals can be considered white noise.

The default test in Box.test() is Box-Pierce. For the more accurate Box-Ljung test, use the option type="Lj". The Ljung–Box test is a statistical test of whether a group of autocorrelations of a time series are different from zero. Instead of testing randomness at each distinct lag, it tests the overall randomness based on a number of lags, and is therefore a portmanteau test. Portmanteau is a French word for a traveling bag (symbolic of something containing many items).

We can use:

# lag=h and fitdf=K
Box.test(res, lag=10, fitdf=0)

## 
##  Box-Pierce test
## 
## data:  res
## X-squared = 10.611, df = 10, p-value = 0.3886

Box.test(res,lag=10, fitdf=0, type="Lj")

## 
##  Box-Ljung test
## 
## data:  res
## X-squared = 11.031, df = 10, p-value = 0.3551

The test-statistic is not that unlikely to happen ($p=0.36$) under the assumption that all autocorrelations are zero. Therefore, we accept the null hypohesis that all autocorrelations in the group are zero. If $p<0.05$ we start doubting the white noise hypothesis. Another model may be considered.

Earlier on, we considered naive models, even for time series with a clear seasonal pattern. We would expect the residuals to be large, and autocorrelated. The Box-Ljung would inform us of the same.

For example, for the time series on beer consumption:

autoplot(ausbeer)

rb <- residuals(naive(ausbeer))
ggAcf(rb)

Box.test(rb,lag=10, fitdf=0, type="Lj")

## 
##  Box-Ljung test
## 
## data:  rb
## X-squared = 612.47, df = 10, p-value < 2.2e-16

3.4 Evaluating Forecast Accuracy

Forecast Errors

The proof of the pudding is in the eating.

The residuals tell us how well the model fits the data. However there is no guarantee that a model with a good fit, produces superior forecasts.

It is important to evaluate forecast accuracy using real forecasts. Since we don't want to wait for future events, we split our data into a training set, and a test set containing recent observations. We can then see how well the model predicts the data in the test set. This procedure is sometimes referred to as backcasting, for obvious reasons.

Commonly used measures of accuracy are:

Mean absolute error (MAE)
Root mean squared error (RMSE)
Mean absolute percentage error (MAPE)
Symmetric MAPE (sMAPE).

beer2 <- window(ausbeer,start=1992,end=c(2007,4))
beerfit1 <- meanf(beer2,h=10)
beerfit2 <- rwf(beer2,h=10)
beerfit3 <- snaive(beer2,h=10)
autoplot(window(ausbeer, start=1992)) +
  autolayer(beerfit1, series="Mean", PI=FALSE) +
  autolayer(beerfit2, series="Na?ve", PI=FALSE) +
  autolayer(beerfit3, series="Seasonal na?ve", PI=FALSE) +
  xlab("Year") + ylab("Megalitres") +
  ggtitle("Forecasts for quarterly beer production") +
  guides(colour=guide_legend(title="Forecast"))

beer3 <- window(ausbeer, start=2008)
accuracy(beerfit1, beer3)

##                   ME     RMSE      MAE        MPE     MAPE     MASE        ACF1
## Training set   0.000 43.62858 35.23438 -0.9365102 7.886776 2.463942 -0.10915105
## Test set     -13.775 38.44724 34.82500 -3.9698659 8.283390 2.435315 -0.06905715
##              Theil's U
## Training set        NA
## Test set      0.801254

accuracy(beerfit2, beer3)

##                       ME     RMSE      MAE         MPE     MAPE     MASE
## Training set   0.4761905 65.31511 54.73016  -0.9162496 12.16415 3.827284
## Test set     -51.4000000 62.69290 57.40000 -12.9549160 14.18442 4.013986
##                     ACF1 Theil's U
## Training set -0.24098292        NA
## Test set     -0.06905715  1.254009

accuracy(beerfit3, beer3)

##                     ME     RMSE  MAE        MPE     MAPE      MASE       ACF1
## Training set -2.133333 16.78193 14.3 -0.5537713 3.313685 1.0000000 -0.2876333
## Test set      5.200000 14.31084 13.4  1.1475536 3.168503 0.9370629  0.1318407
##              Theil's U
## Training set        NA
## Test set      0.298728

You can limit the output to the measures of your choice.

accuracy(beerfit3, beer3)[,c("RMSE","MAE","MPE","MAPE")]

##                  RMSE  MAE        MPE     MAPE
## Training set 16.78193 14.3 -0.5537713 3.313685
## Test set     14.31084 13.4  1.1475536 3.168503

acc3 <- accuracy(beerfit3, beer3)
acc3[,c("RMSE","MAE","MPE","MAPE")]

##                  RMSE  MAE        MPE     MAPE
## Training set 16.78193 14.3 -0.5537713 3.313685
## Test set     14.31084 13.4  1.1475536 3.168503

An elegant way to compare the results of various models, is to rearrange the output using cbind().

rmse1 <- accuracy(beerfit1, beer3)[,c("RMSE")]
rmse2 <- accuracy(beerfit2, beer3)[,c("RMSE")]
rmse3 <- accuracy(beerfit3, beer3)[,c("RMSE")]
cbind(rmse1,rmse2,rmse3)

##                 rmse1    rmse2    rmse3
## Training set 43.62858 65.31511 16.78193
## Test set     38.44724 62.69290 14.31084

Using the goog200 data:

googfc1 <- meanf(goog200, h=40)
googfc2 <- rwf(goog200, h=40)
googfc3 <- rwf(goog200, drift=TRUE, h=40)
autoplot(subset(goog, end = 240)) +
  autolayer(googfc1, PI=FALSE, series="Mean") +
  autolayer(googfc2, PI=FALSE, series="Na?ve") +
  autolayer(googfc3, PI=FALSE, series="Drift") +
  xlab("Day") + ylab("Closing Price (US$)") +
  ggtitle("Google stock price (daily ending 6 Dec 13)") +
  guides(colour=guide_legend(title="Forecast"))

googtest <- window(goog, start=201, end=240)
accuracy(googfc1, googtest)

##                         ME      RMSE       MAE        MPE     MAPE      MASE
## Training set -4.296286e-15  36.91961  26.86941 -0.6596884  5.95376  7.182995
## Test set      1.132697e+02 114.21375 113.26971 20.3222979 20.32230 30.280376
##                   ACF1 Theil's U
## Training set 0.9668981        NA
## Test set     0.8104340  13.92142

accuracy(googfc2, googtest)

##                      ME      RMSE       MAE       MPE      MAPE     MASE
## Training set  0.6967249  6.208148  3.740697 0.1426616 0.8437137 1.000000
## Test set     24.3677328 28.434837 24.593517 4.3171356 4.3599811 6.574582
##                     ACF1 Theil's U
## Training set -0.06038617        NA
## Test set      0.81043397  3.451903

accuracy(googfc3, googtest)

##                         ME      RMSE       MAE         MPE      MAPE     MASE
## Training set -5.998536e-15  6.168928  3.824406 -0.01570676 0.8630093 1.022378
## Test set      1.008487e+01 14.077291 11.667241  1.77566103 2.0700918 3.119002
##                     ACF1 Theil's U
## Training set -0.06038617        NA
## Test set      0.64732736  1.709275

Below we have replicated the model comparison as it apperars in the textbook:

We only use the results for the test set (our main interest!)
We select four measures
We combinee them using rbind()

MeanMethod <- accuracy(googfc1, googtest)[2,c("RMSE","MAE","MAPE","MASE")]
NaiveMethod <- accuracy(googfc2, googtest)[2,c("RMSE","MAE","MAPE","MASE")]
DriftMethod <- accuracy(googfc3, googtest)[2,c("RMSE","MAE","MAPE","MASE")]
round(rbind(MeanMethod,NaiveMethod,DriftMethod),2)

##               RMSE    MAE  MAPE  MASE
## MeanMethod  114.21 113.27 20.32 30.28
## NaiveMethod  28.43  24.59  4.36  6.57
## DriftMethod  14.08  11.67  2.07  3.12

Cross Validation

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

In this procedure, there are a series of test sets, each consisting of a single observation. The training set consists of all observations that occurred prior to the observation in the test set.

The forecast accuracy is computed by averaging over the test sets.

One-step forecasts may not be as relevant as multi-step forecasts. The cross-validation procedure can be modified to allow multi-step errors to be used.

In the example below, we use one-step and four-step forecasts. As can be expected, the errors (or rather, the mean of the squared errors) increases, as it's harder to predict further into the future. Also, the forecast errors are larger than the residuals (which are based on observed data).

sqrt(mean(residuals(rwf(goog200, drift=TRUE))^2, na.rm=TRUE))

## [1] 6.168928

e1 <- tsCV(goog200, rwf, drift=TRUE, h=1)
sqrt(mean(e1^2, na.rm=TRUE))

## [1] 6.233245

e2 <- tsCV(goog200, rwf, drift=TRUE, h=2)
sqrt(mean(e2^2, na.rm=TRUE))

## [1] 7.527005

3.5 Prediction Intervals

Forecasts are normally point forecasts. However, we also want to make statements like with p% certainty the future value will lie in the interval low to high.

Prediction intervals for 80% and 90% are automatically generated, by default for the next 10 periods. Intuitively, the interval widens as the forecast is further into the future.

Graphically, this is displayed by the shaded regions around the point forecasts.

x <- naive(goog200); x

##     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 201       531.4783 523.5222 539.4343 519.3105 543.6460
## 202       531.4783 520.2267 542.7298 514.2705 548.6861
## 203       531.4783 517.6980 545.2586 510.4031 552.5534
## 204       531.4783 515.5661 547.3904 507.1428 555.8138
## 205       531.4783 513.6880 549.2686 504.2704 558.6862
## 206       531.4783 511.9900 550.9666 501.6735 561.2830
## 207       531.4783 510.4285 552.5280 499.2854 563.6711
## 208       531.4783 508.9751 553.9814 497.0627 565.8939
## 209       531.4783 507.6101 555.3465 494.9750 567.9815
## 210       531.4783 506.3190 556.6375 493.0005 569.9561

autoplot(x)

3.6 The forecast() Function

The forecast() function produces forecasts based on either:

Time series, or
Time series models.

If the input is a time series, then the forecast is based on a model automatically selected using the ETS algorithm (discussed in chapter 7).

Preferably, we build our own model, and base our forecasts on that model.

The automatic forecast of ausbeer is illustrated below.

forecast(ausbeer, h=4)

##         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2010 Q3       404.6025 385.8695 423.3355 375.9529 433.2521
## 2010 Q4       480.3982 457.5283 503.2682 445.4216 515.3748
## 2011 Q1       417.0367 396.5112 437.5622 385.6456 448.4277
## 2011 Q2       383.0996 363.5063 402.6930 353.1341 413.0651

3.7 Exercises

Complete exercises 3.5 and 3.6, in the textbook.

5. The Linear Model

5.1 Simple Linear Regression

5.2 Least squares estimation

The regression model describes a linear relationship between the forecast variable $y$ and a single predictor variable $x$.

Below we show the time series of quarterly growth rates of real personal consumption expenditure $y$, and real personal disposable income $x$,

head(uschange,6)

##         Consumption     Income Production   Savings Unemployment
## 1970 Q1   0.6159862  0.9722610 -2.4527003 4.8103115          0.9
## 1970 Q2   0.4603757  1.1690847 -0.5515251 7.2879923          0.5
## 1970 Q3   0.8767914  1.5532705 -0.3587079 7.2890131          0.5
## 1970 Q4  -0.2742451 -0.2552724 -2.1854549 0.9852296          0.7
## 1971 Q1   1.8973708  1.9871536  1.9097341 3.6577706         -0.1
## 1971 Q2   0.9119929  1.4473342  0.9015358 6.0513418         -0.1

autoplot(uschange[,c("Consumption","Income")]) +
  ylab("% change") + xlab("Year")

uschange %>%
  as.data.frame() %>%
  ggplot(aes(x=Income, y=Consumption)) +
  ylab("Consumption (quarterly % change)") +
  xlab("Income (quarterly % change)") +
  geom_point() +
  geom_smooth(method="lm", se=FALSE)

## `geom_smooth()` using formula 'y ~ x'

The equation is estimated in R using the tslm() function:

model <- tslm(Consumption ~ Income, data=uschange)
summary(model)

## 
## Call:
## tslm(formula = Consumption ~ Income, data = uschange)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.40845 -0.31816  0.02558  0.29978  1.45157 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.54510    0.05569   9.789  < 2e-16 ***
## Income       0.28060    0.04744   5.915 1.58e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6026 on 185 degrees of freedom
## Multiple R-squared:  0.159,  Adjusted R-squared:  0.1545 
## F-statistic: 34.98 on 1 and 185 DF,  p-value: 1.577e-08

The fitted line has a positive slope, reflecting the positive relationship between income and consumption.

The slope coefficient shows that a one unit increase in $x% (a one percentage point increase in income) results in an increase of 0.28 percentage points in consumption expenditure).

Multiple Regression

If we have more than one predictor, then we apply multiple regression.

We can first inspect visually if there are correlations in the data.

autoplot(uschange, facet=TRUE)

uschange %>%
  as.data.frame() %>%
  GGally::ggpairs()

Since there are correlations between consuption and all four other variables, we can fit the complete model.

model2 <- tslm(
  Consumption ~ Income + Production + Unemployment + Savings,
  data=uschange)
summary(model2)

## 
## Call:
## tslm(formula = Consumption ~ Income + Production + Unemployment + 
##     Savings, data = uschange)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.88296 -0.17638 -0.03679  0.15251  1.20553 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.26729    0.03721   7.184 1.68e-11 ***
## Income        0.71449    0.04219  16.934  < 2e-16 ***
## Production    0.04589    0.02588   1.773   0.0778 .  
## Unemployment -0.20477    0.10550  -1.941   0.0538 .  
## Savings      -0.04527    0.00278 -16.287  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3286 on 182 degrees of freedom
## Multiple R-squared:  0.754,  Adjusted R-squared:  0.7486 
## F-statistic: 139.5 on 4 and 182 DF,  p-value: < 2.2e-16

The results suggest that only income and savings have a significant effect on consumption. The three asteriks (***) indicate that the probability that these effects are zero is very small ($p<.001$). There's a negative effect of unemployment on consumption (as expected: if unemployment goes up, consuption goes down). The effect, however, is not significant ($0.05<p<0.10$). A challenge with these data is that the predictors are intercorrelated. The regression coefficients control for the effect of all other predictors, which may result in small and/or insignificant effects (for production and unemployment), even though the correlations between consumption on the one hand and production and unemployment on the other hand, are strong.

A measure of how well a linear regression model fits the data, is R². It is called the coefficient of determination, and reflects the proportion of variation in the forecast variable explained by the model.

For forecasting, insignificant coefficients do not necessarily have to be dropped. If they contain some information that leads to better forecasts, then it pays off to keep them!

The fitted data (in green) by and large follow the pattern of the data (in orange). In some cases, the fitted data are substantiaaly higher than the observed data.

autoplot(uschange[,'Consumption'], series="Data") +
  autolayer(model2$fitted.values, series="Fitted") +
  xlab("Year") + ylab("") +
  ggtitle("Percent change in US consumption expenditure") +
  guides(colour=guide_legend(title="Legend"))

cbind(Data = uschange[,"Consumption"],
      Fitted = fitted(model2)) %>%
  as.data.frame() %>%
  ggplot(aes(x=Data, y=Fitted)) +
  geom_point() +
  ylab("Fitted (predicted values)") +
  xlab("Data (actual values)") +
  ggtitle("Percent change in US consumption expenditure") +
  geom_abline(intercept=0, slope=1)

5.3 Evaluating the Regression Model

The main information for evaluating the model is obtained with checkresiduals().

In addition to a plot of the residuals, the ACF and a histogram, the Breusch-Godfrey test for serial correlation is calculated. The p-value is just above the critical value of 0.05, so we can accept the null hypothesis of no serial correlation - even though there are some high positive residuals especially in the first half of the series. This is also reflected in the histogram which resembles a normal distribution but is slight positive skew and outliers on the right of the distribution.

checkresiduals(model2)

## 
##  Breusch-Godfrey test for serial correlation of order up to 8
## 
## data:  Residuals from Linear regression model
## LM test = 14.874, df = 8, p-value = 0.06163

Scatterplots of residuals against fitted values, and against each of the predictors can help in detecting any non-random patterns.

Note that, in the code below:

we have used the function as.data.frame() to convert uschange from a time series (ts object) to a data frame (df)
we have added a column to df, containing the residuals

# Residuals against predictors
df <- as.data.frame(uschange)
df[,"Residuals"] <- as.numeric(residuals(model2))
head(df)

p1 <- ggplot(df, aes(x=Income, y=Residuals)) +
  geom_point()
p2 <- ggplot(df, aes(x=Production, y=Residuals)) +
  geom_point()
p3 <- ggplot(df, aes(x=Savings, y=Residuals)) +
  geom_point()
p4 <- ggplot(df, aes(x=Unemployment, y=Residuals)) +
  geom_point()
gridExtra::grid.arrange(p1, p2, p3, p4, nrow=2)

# Residuals against fitted values
cbind(Fitted = model2$fitted.values, Residuals=residuals(model2)) %>% as.data.frame() %>%
  ggplot(aes(x=Fitted, y=Residuals)) + geom_point()

SPURIOUS CORRELATION

Some models fit very well, even though the variable to be forecasted and the predictors have no theoretic relation to one another.

From the scatterplot and the fitted regression line, it can be seen that residuals are negative if $aussies<15$, then positive if $15<aussies<60$ (with an exception for aussies ∼ 40), and negative for high values of aussies. Note that given the upward trend in both series, high values correspond to recent observations. In short, the residuals do not randomly fluctuate around zero.

Such patterns may indicate non-linear relationships between the variables, or - like here, most likely - spurious correlation. The series just move in the same, upward direction, but in slightly different patterns; fluctuations in the short and medium run in one variable, cannot be explained by the other variable.

aussies <- window(ausair, end=2011)
autoplot(cbind(aussies, guinearice), facet=TRUE)

new <- as.data.frame(cbind(aussies, guinearice))
ggplot(new, aes(x=aussies, y=guinearice)) +
  geom_point() + geom_smooth(method=lm, se=FALSE)

## `geom_smooth()` using formula 'y ~ x'

Test

5.4 Some Useful Predictors

Dummy Variables

Categorical variable taking only two values (e.g., yes or no) can be handled within the framework of multiple regression models by creating a dummy variable which take of 0 or 1. For example, in time series, a dummy for public holiday can be added to the equation.

Dummy variables can also be used for outliers, or improbable observations, in the data. Rather than omitting the outlier, dummy variables can "correct" for the effects of outliers.

For interventions (competitor activity; advertising expenditure; ...) that may have affected the variable to be forecast may have lasted only for one or a few periods, we can use a spike variable that takes 1 in the period of the intervention and 0 elsewhere. A spike variable is equivalent to a dummy variable for outliers.

If there are more than two categories, then the variable can be coded using several dummy variables (one fewer than the total number of categories). For example, for public holidays there are two values (yes or no), and therefore 2-1=1 dummy. If we want to distinguish normal days (including Saturdays) versus public holidays and Sundays, we would add two dummies, one for public holidays and one for Sundays. Note that both variables can be coded 1, in case a public holiday happens on a Sunday (like Easter).

Seasonal Dummy Variables

We can use dummy variables for seasonal data.

The tslm() function automatically handles dummy variables. For quarterly, we would need three dummies (for quarters 2, 3 and 4, if we want to estimate the effects of these quarters versus quarter 1 as the base category).

For example:

beer2 <- window(ausbeer, start=1992)
autoplot(beer2) + xlab("Year") + ylab("Megalitres")

fit.beer <- tslm(beer2 ~ trend + season)
summary(fit.beer)

## 
## Call:
## tslm(formula = beer2 ~ trend + season)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -42.903  -7.599  -0.459   7.991  21.789 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 441.80044    3.73353 118.333  < 2e-16 ***
## trend        -0.34027    0.06657  -5.111 2.73e-06 ***
## season2     -34.65973    3.96832  -8.734 9.10e-13 ***
## season3     -17.82164    4.02249  -4.430 3.45e-05 ***
## season4      72.79641    4.02305  18.095  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.23 on 69 degrees of freedom
## Multiple R-squared:  0.9243, Adjusted R-squared:  0.9199 
## F-statistic: 210.7 on 4 and 69 DF,  p-value: < 2.2e-16

There is an average downward trend of -0.34 megalitres per quarter.

On average:

Production in the second quarter is 34.7 megalitres lower than the first quarter.
The third quarter is 17.8 megalitres lower than the first quarter.
Production is highest in the fourth quarter: 72.8 megalitres higher than the first quarter (and $72.8 + 34.7 = 107.5$ higher than in the second quarter).

autoplot(beer2, series="Data") +
  autolayer(fitted(fit.beer), series="Fitted") +
  xlab("Year") + ylab("Megalitres") +
  ggtitle("Quarterly Beer Production")

cbind(Data=beer2, Fitted=fitted(fit.beer)) %>%
  as.data.frame() %>%
  ggplot(aes(x = Data, y = Fitted, colour = as.factor(cycle(beer2)))) +
  geom_point() +
  ylab("Fitted") + xlab("Actual values") +
  ggtitle("Quarterly beer production") +
  scale_colour_brewer(palette="Dark2", name="Quarter") +
  geom_abline(intercept=0, slope=1)

Fourier Series

Fourier terms are an alternative to using seasonal dummy variables, especially for long seasonal periods.

With Fourier terms, we need fewer predictors than with dummy variables, especially when m is large. This makes them useful for weekly data, for example, where m ≈ 52. For short seasonal periods (e.g., quarterly data), there is little advantage in using Fourier terms over seasonal dummy variables.

These Fourier terms are produced using the fourier() function. For example:

fourier.beer <- tslm(beer2 ~ trend + fourier(beer2, K=2))
summary(fourier.beer)

## 
## Call:
## tslm(formula = beer2 ~ trend + fourier(beer2, K = 2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -42.903  -7.599  -0.459   7.991  21.789 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               446.87920    2.87321 155.533  < 2e-16 ***
## trend                      -0.34027    0.06657  -5.111 2.73e-06 ***
## fourier(beer2, K = 2)S1-4   8.91082    2.01125   4.430 3.45e-05 ***
## fourier(beer2, K = 2)C1-4  53.72807    2.01125  26.714  < 2e-16 ***
## fourier(beer2, K = 2)C2-4  13.98958    1.42256   9.834 9.26e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.23 on 69 degrees of freedom
## Multiple R-squared:  0.9243, Adjusted R-squared:  0.9199 
## F-statistic: 210.7 on 4 and 69 DF,  p-value: < 2.2e-16

5.5 Selecting Predictors

For selecting the best predictors to use in a regression model, we can make use of five measures.

fit.consMR <- tslm(
  Consumption ~ Income + Production + Unemployment + Savings,
  data=uschange)
CV(fit.consMR)

##           CV          AIC         AICc          BIC        AdjR2 
##    0.1163477 -409.2980298 -408.8313631 -389.9113781    0.7485856

Which measure should we use?

While R² is widely used, its tendency to select too many predictor variables makes it less suitable for forecasting.

We recommend that one of the AICc, AIC, or CV statistics be used, each of which has forecasting as their objective. If the value of T is large, they will all lead to the same model.

The best strategy is to gradually enter new variables, and see whether the newly added variable leads to a better model fit.

Suppose we want to test whether savings lead to a better model fit, starting from a model with only income as a predictor.

model1 <- tslm(
  Consumption ~ Income,
  data=uschange)
model2 <- tslm(
  Consumption ~ Income + Savings,
  data=uschange)
CV(model1); CV(model2)

##           CV          AIC         AICc          BIC        AdjR2 
##    0.3721815 -185.4377314 -185.3065838 -175.7444055    0.1544792

##           CV          AIC         AICc          BIC        AdjR2 
##    0.1286616 -388.7271972 -388.5074170 -375.8027628    0.7163990

Since the CV measure and the information criteria (AIC; AICc; and BIC) are smaller (or more negative) and the R² higher for model2 that includes savings as a predictor, we conclude that adding savings to the model leads to a better fit (and, hence, a better level of understanding and ability predict).

To show you what happens if we include a variable uncorrelated to consumption, we add the production figures of some other country (Prod2) as a random variable in the same range as those of the country at hand.

usc2 <- as.data.frame(uschange)
summary(usc2$Production);nrow(usc2)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -6.85104  0.05568  0.65793  0.50806  1.30572  4.14957

## [1] 187

set.seed(376)
usc2$Prod2 <- runif(nrow(usc2), min=-5, max=+6)
head(usc2)

usc2t <- ts(usc2,start=1970, frequency=4)
str(usc2t)

##  Time-Series [1:187, 1:6] from 1970 to 2016: 0.616 0.46 0.877 -0.274 1.897 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:6] "Consumption" "Income" "Production" "Savings" ...

model3 <- tslm(Consumption ~ Income + Savings, data=usc2t)
model4 <- tslm(Consumption ~ Income + Savings + Prod2, data=usc2t)
round(CV(model3),4)

##        CV       AIC      AICc       BIC     AdjR2 
##    0.1287 -388.7272 -388.5074 -375.8028    0.7164

round(CV(model4),4)

##        CV       AIC      AICc       BIC     AdjR2 
##    0.1287 -388.4452 -388.1137 -372.2896    0.7175

The (adjusted) R² is slightly higher for model 4. Adding variables to the model always leads to a better R²! This makes it a poor measure for comparing models.

The BIC seems to be the measure most sensitive to adding poor predictors. The value increases (becomes less negative), suggesting that model 3 is a better model by dropping marginal or irrelevant predictors.

modelfinal <- tslm(Consumption ~ Income + Production + Savings + Unemployment, data=uschange)
autoplot(uschange[,'Consumption'], series="Data") +
  autolayer(modelfinal$fitted.values, series="Fitted") +
  xlab("Year") + ylab("") +
  ggtitle("Percent change in US consumption expenditure") +
  guides(colour=guide_legend(title="Legend"))

5.6 Regression Forecasts

fit.consBest <- tslm(Consumption ~ Income + Savings + Unemployment,data = uschange)
h <- 4
newdata <- data.frame(
  Income = c(1, 1, 1, 1),
  Savings = c(0.5, 0.5, 0.5, 0.5),
  Unemployment = c(0, 0, 0, 0))
fcast.up <- forecast(fit.consBest, newdata = newdata)
newdata <- data.frame(
  Income = rep(-1, h),
  Savings = rep(-0.5, h),
  Unemployment = rep(0, h))
fcast.down <- forecast(fit.consBest, newdata = newdata)

autoplot(uschange[, 1]) +
  ylab("% change in US consumption") +
  autolayer(fcast.up, PI = TRUE, series = "increase") +
  autolayer(fcast.down, PI = TRUE, series = "decrease") +
  guides(colour = guide_legend(title = "Scenario"))

fit.cons <- tslm(Consumption ~ Income, data = uschange)
h <- 4
fcast.ave <- forecast(fit.cons, newdata = data.frame(Income = rep(mean(uschange[,"Income"]), h)))
fcast.up <- forecast(fit.cons, newdata = data.frame(Income = rep(5, h)))
autoplot(uschange[, "Consumption"]) +
  ylab("% change in US consumption") +
  autolayer(fcast.ave, series = "Average increase",PI = TRUE) +
  autolayer(fcast.up, series = "Extreme increase", PI = TRUE) +
  guides(colour = guide_legend(title = "Scenario"))

Forecasting: Principles and Practice