Forecasting: Principles and Practice 3

library(fpp3)

-- Attaching packages --------------------- fpp3 0.3 --
v tibble      3.0.6       v tsibble     0.9.3  
v dplyr       1.0.3       v tsibbledata 0.2.0  
v tidyr       1.1.2       v feasts      0.1.6  
v lubridate   1.7.9.2     v fable       0.2.1  
v ggplot2     3.3.3       
-- Conflicts ------------------------ fpp3_conflicts --
x lubridate::date()   masks base::date()
x dplyr::filter()     masks stats::filter()
x tsibble::interval() masks lubridate::interval()
x dplyr::lag()        masks stats::lag()

https://robjhyndman.com/seminars/uwa2017/

https://robjhyndman.com/seminars/isi2019workshop/

https://github.com/robjhyndman/ETC3550Slides/tree/fable

https://github.com/rstudio-conf-2020/time-series-forecasting/blob/master/materials/labs.R

http://www.gradaanwr.net/

1 Getting started

1.1 What can be forecast?

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

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

Good forecasts capture the genuine patterns and relationships which exist in the historical data, but do not replicate past events that will not occur again.

Forecasts rarely assume that the environment is unchanging. What is normally assumed is that the way in which the environment is changing will continue into the future.

1.2 Forecasting, planning and goals

They are three different things.

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

• Goals: are what you would like to have happen. Goals should be linked to forecasts and plans, but this does not always occur. Too often, goals are set without any plan for how to achieve them, and no forecasts for whether they are realistic.

• Planning: is a response to forecasts and goals. Planning involves determining the appropriate actions that are required to make your forecasts match your goals.

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

1.3 Determining what to forecast

In the early stages of a forecasting project, decisions need to be made about what should be forecast. It is also necessary to consider the forecasting horizon.

Once it has been determined what forecasts are required, it is then necessary to find or collect the data on which the forecasts will be based.

1.4 Forecasting data and methods

If there are no data available, or if the data available are not relevant to the forecasts, then qualitative forecasting methods must be used.

Quantitative forecasting can be applied when two conditions are satisfied:

• numerical information about the past is available;
• it is reasonable to assume that some aspects of the past patterns will continue into the future.

1.4.1 Time series forecasting

The simplest time series forecasting methods use only information on the variable to be forecast, and make no attempt to discover the factors that affect its behaviour. Therefore they will extrapolate trend and seasonal patterns, but they ignore all other information such as marketing initiatives, competitor activity, changes in economic conditions, and so on.

1.4.2 Predictor variables and time series forecasting

Predictor variables are often useful in time series forecasting. For example, suppose we wish to forecast the hourly electricity demand (ED) of a hot region during the summer period. A model with predictor variables might be of the form

\[\begin{align*} \text{ED} = & f(\text{current temperature, strength of economy, population,}\\ & \qquad\text{time of day, day of week, error}). \end{align*}\]

The “error” term on the right allows for random variation and the effects of relevant variables that are not included in the model. We call this an explanatory model because it helps explain what causes the variation in electricity demand.

Because the electricity demand data form a time series, we could also use a time series model for forecasting. In this case, a suitable time series forecasting equation is of the form

\[\begin{align*} {ED}_{t+1} = f(\text{ED}_{t}, \text{ED}_{t-1}, \text{ED}_{t-2}, \text{ED}_{t-3},\dots, \text{error} \end{align*}\]

Here, prediction of the future is based on past values of a variable, but not on external variables that may affect the system. Again, the “error” term on the right allows for random variation and the effects of relevant variables that are not included in the model.

There is also a third type of model which combines the features of the above two models.

\[\begin{align*} \text{ED}_{t+1} = f(\text{ED}_{t}, \text{current temperature, time of day, day of week, error}). \end{align*}\]

They are known as dynamic regression models, panel data models, longitudinal models, transfer function models, and linear system models (assuming that \(f\) is linear).

The model to be used in forecasting depends on the resources and data available, the accuracy of the competing models, and the way in which the forecasting model is to be used.

1.5 The basic steps in a forecasting task

1.5.1 Step 1: Problem definition.

Defining the problem carefully requires an understanding of the way the forecasts will be used, who requires the forecasts, and how the forecasting function fits within the organisation requiring the forecasts. A forecaster needs to spend time talking to everyone who will be involved in collecting data, maintaining databases, and using the forecasts for future planning.

1.5.2 Step 2: Gathering information.

There are always at least two kinds of information required: (a) statistical data, and (b) the accumulated expertise of the people who collect the data and use the forecasts. Occasionally, old data will be less useful due to structural changes in the system being forecast; then we may choose to use only the most recent data. However, remember that good statistical models will handle evolutionary changes in the system; don’t throw away good data unnecessarily.

1.5.3 Step 3: Preliminary (exploratory) analysis.

Always start by graphing the data. Are there consistent patterns? Is there a significant trend? Is seasonality important? Is there evidence of the presence of business cycles? Are there any outliers in the data that need to be explained by those with expert knowledge? How strong are the relationships among the variables available for analysis?

1.5.4 Step 4: Choosing and fitting models

The best model to use depends on the availability of historical data, the strength of relationships between the forecast variable and any explanatory variables, and the way in which the forecasts are to be used. It is common to compare two or three potential models. Each model is itself an artificial construct that is based on a set of assumptions (explicit and implicit) and usually involves one or more parameters which must be estimated using the known historical data.

1.5.5 Step 5: Using and evaluating a forecasting model.

Once a model has been selected and its parameters estimated, the model is used to make forecasts. The performance of the model can only be properly evaluated after the data for the forecast period have become available.

1.6 The statistical forecasting perspective

The further ahead we forecast, the more uncertain we are. When we obtain a forecast, we are estimating the middle of the range of possible values the random variable could take. Often, a forecast is accompanied by a prediction interval giving a range of values the random variable could take with relatively high probability.

\(y_t\) will denote the observation at time t

Suppose we denote all the information we have observed as \(\mathcal{I}\) and we want to forecast \(y_t\). We then write \(y_{t} |\mathcal{I}\) meaning “the random variable \(y_t\) given what we know in \(\mathcal{I}\).” The set of values that this random variable could take, along with their relative probabilities, is known as the “probability distribution” of \(y_{t} |\mathcal{I}\) . In forecasting, we call this the forecast distribution.

When we talk about the “forecast,” we usually mean the average value of the forecast distribution, and we put a “hat” over \(y\) to show this (\(\hat{y}_t\)).

2 Time series graphics

2.1 tsibble objects

https://otexts.com/fpp3/tsibbles.html

tsibble objects extend tidy data frames (tibble objects) by introducing temporal structure.

y <- tsibble(Year = 2015:2019, Observation = c(123, 39, 78, 52, 110), index = Year)

We have set the time series index to be the Year column, which associates the measurements (Observation) with the time of recording (Year).

For observations that are more frequent than once per year, we need to use a time class function on the index.

2.2 Time plots

melsyd_economy <- ansett %>%
  filter(Airports == "MEL-SYD", Class == "Economy")
autoplot(melsyd_economy, Passengers) +
  labs(title = "Ansett economy class passengers",
       subtitle = "Melbourne-Sydney")

The time plot immediately reveals some interesting features.

• There was a period in 1989 when no passengers were carried — this was due to an industrial dispute.
• There was a period of reduced load in 1992. This was due to a trial in which some economy class seats were replaced by business class seats.
• A large increase in passenger load occurred in the second half of 1991. There are some large dips in load around the start of each year. These are due to holiday effects.
• There is a long-term fluctuation in the level of the series which increases during 1987, decreases in 1989, and increases again through 1990 and 1991.
• There are some periods of missing observations.

PBS %>%
  filter(ATC2 == "A10") %>%
  select(Month, Concession, Type, Cost) %>%
  summarise(TotalC = sum(Cost)) %>%
  mutate(Cost = TotalC / 1e6) -> a10

autoplot(a10, Cost) +
  labs(y = "$ million", title = "Antidiabetic drug sales")

Here, there is a clear and increasing trend. There is also a strong seasonal pattern that increases in size as the level of the series increases. The sudden drop at the start of each year is caused by a government subsidisation scheme that makes it cost-effective for patients to stockpile drugs at the end of the calendar year. Any forecasts of this series would need to capture the seasonal pattern, and the fact that the trend is changing slowly.

2.3 Time series patterns

• Trend: A trend exists when there is a long-term increase or decrease in the data. Sometimes we will refer to a trend as “changing direction,” when it might go from an increasing trend to a decreasing trend.

• Seasonal: A seasonal pattern occurs when a time series is affected by seasonal factors such as the time of the year or the day of the week. Seasonality is always of a fixed and known period.

• Cyclic: A cycle occurs when the data exhibit rises and falls that are not of a fixed frequency. These fluctuations are usually due to economic conditions, and are often related to the “business cycle.”

2.4 Seasonal plots

a10 %>%
  gg_season(Cost, labels = "both") +
  labs(y = "$ million",
       title = "Seasonal plot: antidiabetic drug sales")

It is clear that there is a large jump in sales in January each year. Actually, these are probably sales in late December as customers stockpile before the end of the calendar year, but the sales are not registered with the government until a week or two later. The graph also shows that there was an unusually small number of sales in March 2008 (most other years show an increase between February and March). The small number of sales in June 2008 is probably due to incomplete counting of sales at the time the data were collected.

2.4.1 Multiple seasonal periods

Where the data has more than one seasonal pattern, the period argument can be used to select which seasonal plot is required.

vic_elec %>% gg_season(Demand, period = "day") + theme(legend.position = "none")

vic_elec %>% gg_season(Demand, period = "week") + theme(legend.position = "none")

vic_elec %>% gg_season(Demand, period = "year")

2.5 Seasonal subseries plots

a10 %>%
  gg_subseries(Cost) +
  labs(y = "$ million",
       title = "Seasonal subseries plot: antidiabetic drug sales")

The blue horizontal lines indicate the means for each month. This form of plot enables the underlying seasonal pattern to be seen clearly, and also shows the changes in seasonality over time. It is especially useful in identifying changes within particular seasons. In this example, the plot is not particularly revealing; but in some cases, this is the most useful way of viewing seasonal changes over time.

2.5.1 Example: Australian holiday tourism

total visitor nights spent on Holiday by State for each quarter

holidays <- tourism %>%
  filter(Purpose == "Holiday") %>%
  group_by(State) %>%
  summarise(Trips = sum(Trips))

head(holidays)

Strong seasonality for most states but seasonal peaks do not coincide.

autoplot(holidays, Trips) +
  labs(y = "thousands of trips",
       title = "Australian domestic holiday nights")

Southern states of Australia (Tasmania, Victoria and South Australia) have strongest tourism in Q1 (their summer), while the northern states (Queensland and the Northern Territory) have the strongest tourism in Q3 (their dry season).

gg_season(holidays, Trips) +
  labs(y = "thousands of trips",
       title = "Australian domestic holiday nights")

Corresponding subseries plots:

holidays %>%
  gg_subseries(Trips) +
  labs(y = "thousands of trips",
       title = "Australian domestic holiday nights")

2.6 Scatterplots

Explore relationships between time series.

vic_elec %>%
  filter(year(Time) == 2014) %>%
  autoplot(Demand) +
  labs(y = "Demand (GW)",
       title = "Half-hourly electricity demand: Victoria, Australia")

vic_elec %>%
  filter(year(Time) == 2014) %>%
  autoplot(Temperature) +
  labs(y = "Temperature (degrees Celsius)",
       title = "Half-hourly temperatures: Melbourne, Australia")

We can study the relationship between demand and temperature by plotting one series against the other.

vic_elec %>%
  filter(year(Time) == 2014) %>%
  ggplot(aes(x = Temperature, y = Demand)) +
  geom_point() +
  labs(y = "Demand (GW)", x = "Temperature (degrees Celsius)")

It is clear that high demand occurs when temperatures are high due to the effect of air-conditioning. But there is also a heating effect, where demand increases for very low temperatures.

2.6.1 Correlation

Correlation coefficients to measure the strength of the relationship between two variables.

\[\begin{align*} r = \frac{\sum (x_{t} - \bar{x})(y_{t}-\bar{y})}{\sqrt{\sum(x_{t}-\bar{x})^2}\sqrt{\sum(y_{t}-\bar{y})^2}}. \end{align*}\]

The value of \(r\) always lies between \(−1\) and \(1\) with negative values indicating a negative relationship and positive values indicating a positive relationship.

2.6.2 Scatterplot matrices

When there are several potential predictor variables, it is useful to plot each variable against each other variable.

visitors <- tourism %>%
  group_by(State) %>%
  summarise(Trips = sum(Trips))
visitors %>%
  ggplot(aes(x = Quarter, y = Trips)) +
  geom_line() +
  facet_grid(vars(State), scales = "free_y") +
  labs(y = "Number of visitor nights each quarter (millions)")

visitors %>%
  pivot_wider(values_from=Trips, names_from=State) %>%
  GGally::ggpairs(columns = 2:9)

In this example, mostly positive relationships are revealed, with the strongest relationships being between the neighboring states located in the south and south east coast of Australia, namely, New South Wales, Victoria and South Australia. Some negative relationships are also revealed between the Northern Territory and other regions. The Northern Territory is located in the north of Australia famous for its outback desert landscapes visited mostly in winter. Hence, the peak visitation in the Northern Territory is in the July (winter) quarter in contrast to January (summer) quarter for the rest of the regions.

2.7 Lag plots

horizontal axis shows lagged values of the time series. Each graph shows \(y_t\) plotted against \(y_{t−k}\) for different values of \(k\).

recent_production <- aus_production %>%
  filter(year(Quarter) >= 2000)
recent_production %>% gg_lag(Beer, geom = "point")

The relationship is strongly positive at lags 4 and 8, reflecting the strong seasonality in the data.

2.8 Autocorrelation

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

\[\begin{align*} r_{k} = \frac{\sum\limits_{t=k+1}^T (y_{t}-\bar{y})(y_{t-k}-\bar{y})} {\sum\limits_{t=1}^T (y_{t}-\bar{y})^2}, \end{align*}\]

The autocorrelation coefficients make up the autocorrelation function or ACF.

recent_production %>% ACF(Beer, lag_max = 9)

We usually plot the ACF to see how the correlations change with the lag \(k\). The plot is sometimes known as a correlogram.

recent_production %>%
  ACF(Beer) %>%
  autoplot()

• \(r_4\) is higher than for the other lags. This is due to the seasonal pattern in the data: the peaks tend to be four quarters apart and the troughs tend to be four quarters apart.
• \(r_2\) is more negative than for the other lags because troughs tend to be two quarters behind peaks.
• The dashed blue lines indicate whether the correlations are significantly different from zero.

2.8.1 Trend and seasonality in ACF plots

When data have a trend, the autocorrelations for small lags tend to be large and positive because observations nearby in time are also nearby in size. So the ACF of trended time series tend to have positive values that slowly decrease as the lags increase.

a10 %>%
  ACF(Cost, lag_max = 48) %>%
  autoplot()

The slow decrease in the ACF as the lags increase is due to the trend, while the “scalloped” shape is due to the seasonality.

2.9 White noise

Time series that show no autocorrelation are called white noise.

set.seed(30)
y <- tsibble(sample = 1:50, wn = rnorm(50), index = sample)
y %>% autoplot(wn) + labs(title = "White noise")

y %>%
  ACF(wn) %>%
  autoplot()

We expect each autocorrelation to be close to zero. . If one or more large spikes are outside these bounds, or if substantially more than 5% of spikes are outside these bounds, then the series is probably not white noise.