In this chapter, we briefly discuss four more advanced forecasting methods that build on the models discussed in earlier chapters.
12.1 Complex seasonality
So far, we have mostly considered relatively simple seasonal patterns such as quarterly and monthly data. However, higher frequency time series often exhibit more complicated seasonal patterns. For example, daily data may have a weekly pattern as well as an annual pattern. Hourly data usually has three types of seasonality: a daily pattern, a weekly pattern, and an annual pattern. Even weekly data can be challenging to forecast as there are not a whole number of weeks in a year, so the annual pattern has a seasonal period of 365.25/7\approx 52.179 on average. Most of the methods we have considered so far are unable to deal with these seasonal complexities.
We don’t necessarily want to include all the possible seasonal periods in our models — just the ones that are likely to be present in the data. For example, if we have only 180 days of data, we may ignore the annual seasonality. If the data are measurements of a natural phenomenon (e.g., temperature), we can probably safely ignore any weekly seasonality.
Figure 12.1 shows the number of calls to a North American commercial bank per 5-minute interval between 7:00am and 9:05pm each weekday over a 33 week period. Figure 12.2 shows the first four weeks of the same time series. There is a strong daily seasonal pattern with period 169 (there are 169 5-minute intervals per day), and a weak weekly seasonal pattern with period 169 \times 5=845. (Call volumes on Mondays tend to be higher than the rest of the week.) If a longer series of data were available, we may also have observed an annual seasonal pattern.
Figure 12.1: Five-minute call volume handled on weekdays between 7:00am and 9:05pm in a large North American commercial bank. Data from 3 March – 24 October 2003.
Code
df = bank_calls_5min.loc[lambda x: x["ds"] <="2003-03-29"]fig, ax = plt.subplots(figsize=(8, 3))plot_series(df, xlabel="DateTime [5m]", ylabel="Calls", title="Five-minute call volume to bank over 4 weeks", ax=ax,)
Figure 12.2: First four weeks of data from the North American commercial bank.
Apart from the multiple seasonal periods, this series has the additional complexity of missing values between the working periods.
MSTL with multiple seasonal periods
The MSTL() function is designed to deal with multiple seasonalities, and it can be used to decompose a series into its trend, multiple seasonal, and remainder components.
Figure 12.3: STL decomposition with multiple seasonality for the call volume data.
There are two seasonal patterns shown in Figure 12.3, one for the time of day (the third panel), and one for the time of week (the fourth panel). Note that in this case, there is little trend seen in the data, and the weekly seasonality is weak.
The decomposition can also be used in forecasting, with each of the seasonal components forecast using a Seasonal Naïve method, and the seasonally adjusted data forecast using ETS.
The MSTL method can also be used for forecasting, as shown below and in Figure 12.4. Since StatsForecast always produces continuous dates, we need to account for the bank being closed before 7:00am, after 9:05pm, and on weekends. We will do this by generating the correct dates directly, using the knowledge that the dataset we used to generate the forecast ends on Friday the 24th at 9:05 PM.
Figure 12.4: Forecasts of the call volume data using an STL decomposition with the seasonal components forecast using a Seasonal Naïve method, and the seasonally adjusted data forecast using ETS.
Dynamic harmonic regression with multiple seasonal periods
With multiple seasonalities, we can use Fourier terms as we did in earlier chapters (see Sections 7.4 and 10.5). Because there are multiple seasonalities, we need to add Fourier terms for each seasonal period. In this case, the seasonal periods are 169 and 845, so the Fourier terms are of the form
\sin\left(\frac{2\pi kt}{169}\right), \quad
\cos\left(\frac{2\pi kt}{169}\right), \quad
\sin\left(\frac{2\pi kt}{845}\right), \quad \text{and} \quad
\cos\left(\frac{2\pi kt}{845}\right),
for k=1,2,\dots. As usual, the fourier() function can generate these for you.
We will fit a dynamic harmonic regression model with an ARIMA error structure. The total number of Fourier terms for each seasonal period could be selected to minimise the AICc. However, for high seasonal periods, this tends to over-estimate the number of terms required, so we will use a more subjective choice with 10 terms for the daily seasonality and 5 for the weekly seasonality. We set D=d=0 in order to handle the non-stationarity through the regression terms, and P=Q=0 in order to handle the seasonality through the regression terms (this is achieved by using seasonal=False while instantiating AutoARIMA() as follows).
Figure 12.5: Forecasts from a dynamic harmonic regression applied to the call volume data.
This is a large model, containing 33 parameters: 4 ARMA coefficients, 20 Fourier coefficients for period 169, 10 Fourier coefficients for period 845, and an intercept. Two of the Fourier terms for period 845 must be dropped because they overlap with the terms of period 169 (since 845=5\times169).
Example: Electricity demand
One common application of such models is electricity demand modelling. Figure 12.6 shows half-hourly electricity demand (MWh) in Victoria, Australia, during 2012–2014, along with temperatures (degrees Celsius) for the same period for Melbourne (the largest city in Victoria).
Figure 12.6: Half-hourly electricity demand and corresponding temperatures in 2012–2014, Victoria, Australia.
Plotting electricity demand against temperature (Figure 12.7) shows that there is a nonlinear relationship between the two, with demand increasing for low temperatures (due to heating) and increasing for high temperatures (due to cooling).
Figure 12.7: Half-hourly electricity demand for Victoria, plotted against temperatures for the same times in Melbourne, the largest city in Victoria.
We will fit a regression model with a piecewise linear function of temperature (containing a knot at 18 degrees), and harmonic regression terms to allow for the daily seasonal pattern. Again, we set the orders of the Fourier terms subjectively, while using the AICc to select the order of the ARIMA errors.
Forecasting with such models is difficult because we require future values of the predictor variables. Future values of the Fourier terms are easy to compute, but future temperatures are, of course, unknown. If we are only interested in forecasting up to a week ahead, we could use temperature forecasts obtained from a meteorological model. Alternatively, we could use scenario forecasting (Section 6.5) and plug in possible temperature patterns. In the following example, we have used a repeat of the last two days of temperatures to generate future possible demand values.
Figure 12.8: Forecasts from a dynamic harmonic regression model applied to half-hourly electricity demand data.
Although the short-term forecasts look reasonable, this is a crude model for a complicated process. The residuals, plotted in Figure 12.9, demonstrate that there is a lot of information that has not been captured with this model.
Figure 12.9: Residual diagnostics for the dynamic harmonic regression model.
More sophisticated versions of this model which provide much better forecasts are described in Hyndman & Fan (2010) and Fan & Hyndman (2012).
12.2 Prophet model
A recent proposal is the Prophet model, available via PyPI. This model was introduced by Facebook (Taylor & Letham, 2018), originally for forecasting daily data with weekly and yearly seasonality, plus holiday effects. It was later extended to cover more types of seasonal data. It works best with time series that have strong seasonality and several seasons of historical data.
Prophet can be considered a nonlinear regression model (Chapter 7), of the form
y_t = g(t) + s(t) + h(t) + \varepsilon_t,
where g(t) describes a piecewise-linear trend (or “growth term”), s(t) describes the various seasonal patterns, h(t) captures the holiday effects, and \varepsilon_t is a white noise error term.
The knots (or changepoints) for the piecewise-linear trend are automatically selected if not explicitly specified. Optionally, a logistic function can be used to set an upper bound on the trend.
The seasonal component consists of Fourier terms of the relevant periods. By default, order 10 is used for annual seasonality and order 3 is used for weekly seasonality.
Holiday effects are added as simple dummy variables.
The model is estimated using a Bayesian approach to allow for automatic selection of the changepoints and other model characteristics.
We illustrate the approach using two data sets: a simple quarterly example, and then the electricity demand data described in the previous section.
Example: Quarterly cement production
For this simple quarterly example, we will repeat the analysis from Section 9.10 in which we compared an ARIMA and ETS model, but we will add in a Prophet model for comparison. By default, Prophet fits weekly and yearly seasonalities. Since we have quarterly data, we need to specify this directly using add_seasonality.
Figure 12.11: Components of a Prophet model fitted to the Victorian electricity demand
Figure 12.11 shows the trend and seasonal components of the fitted model.
The model specification is very similar to the DHR model in the previous section, although the result is different in several important ways. The Prophet model adds a piecewise linear time trend which is not really appropriate here as we don’t expect the long term forecasts to continue to follow the downward linear trend at the end of the series.
There is also substantial remaining autocorrelation in the residuals.
Prophet rarely gives better forecast accuracy than the alternative approaches, as these two examples have illustrated. However, it is a popular model for forecasting, with thousands of stars on GitHub and used as a dependency in hundreds of repositories. Its popularity can partly be explained by its open-source nature and ease of use. However, several studies, such as Garza et al. (2022) and Menculini et al. (2021), have shown that it does not perform well in terms of time and accuracy.
StatsForecast offers an adapter to Prophet with improved performance, as shown by Garza et al. (2022). It returns the best ARIMA model using Prophet’s feature preprocessing. To use it, simply replace Prophet with AutoARIMAProphet, as shown below using the previous example.
There is a more advanced model called NeuralProphet that combines neural networks with traditional time series algorithms and was inspired by Prophet. This model is available through PyPI and can be set up in a way similar to Prophet.
12.3 Vector autoregressions
One limitation of the models that we have considered so far is that they impose a unidirectional relationship — the forecast variable is influenced by the predictor variables, but not vice versa. However, there are many cases where the reverse should also be allowed for — where all variables affect each other. In previous chapters, the changes in personal consumption expenditure (C_t) were forecast based on the changes in personal disposable income (I_t). However, in this case a bi-directional relationship may be more suitable: an increase in I_t will lead to an increase in C_t and vice versa.
An example of such a situation occurred in Australia during the Global Financial Crisis of 2008–2009. The Australian government issued stimulus packages that included cash payments in December 2008, just in time for Christmas spending. As a result, retailers reported strong sales and the economy was stimulated. Consequently, incomes increased.
Such feedback relationships are allowed for in the vector autoregressive (VAR) framework. In this framework, all variables are treated symmetrically. They are all modelled as if they all influence each other equally. In more formal terminology, all variables are now treated as “endogenous”. To signify this, we now change the notation and write all variables as ys: y_{1,t} denotes the tth observation of variable y_1, y_{2,t} denotes the tth observation of variable y_2, and so on.
A VAR model is a generalisation of the univariate autoregressive model for forecasting a vector of time series.1 It comprises one equation per variable in the system. The right hand side of each equation includes a constant and lags of all the variables in the system. To keep it simple, we will consider a two variable VAR with one lag. We write a 2-dimensional VAR(1) model as
\begin{aligned}
y_{1,t} &= c_1+\phi _{11,1}y_{1,t-1}+\phi _{12,1}y_{2,t-1}+\varepsilon_{1,t} \\
y_{2,t} &= c_2+\phi _{21,1}y_{1,t-1}+\phi _{22,1}y_{2,t-1}+\varepsilon_{2,t},
\end{aligned}
\tag{12.1} where \varepsilon_{1,t} and \varepsilon_{2,t} are white noise processes that may be contemporaneously correlated. The coefficient \phi_{ii,\ell} captures the influence of the \ellth lag of variable y_i on itself, while the coefficient \phi_{ij,\ell} captures the influence of the \ellth lag of variable y_j on y_i.
If the series are stationary, we forecast them by fitting a VAR to the data directly (known as a “VAR in levels”). If the series are non-stationary, we take differences of the data in order to make them stationary, then fit a VAR model (known as a “VAR in differences”). In both cases, the models are estimated equation by equation using the principle of least squares. For each equation, the parameters are estimated by minimising the sum of squared \varepsilon_{i,t} values.
The other possibility, which is beyond the scope of this book and therefore we do not explore here, is that the series may be non-stationary but cointegrated, which means that there exists a linear combination of them that is stationary. In this case, a VAR specification that includes an error correction mechanism (usually referred to as a vector error correction model) should be included, and alternative estimation methods to least squares estimation should be used.2
Forecasts are generated from a VAR in a recursive manner. The VAR generates forecasts for each variable included in the system. To illustrate the process, assume that we have fitted the 2-dimensional VAR(1) model described in Equations 12.1, for all observations up to time T. Then the one-step-ahead forecasts are generated by \begin{align*}
\hat y_{1,T+1|T} &=\hat{c}_1+\hat\phi_{11,1}y_{1,T}+\hat\phi_{12,1}y_{2,T} \\
\hat y_{2,T+1|T} &=\hat{c}_2+\hat\phi _{21,1}y_{1,T}+\hat\phi_{22,1}y_{2,T}.
\end{align*} This is the same form as Equations 12.1, except that the errors have been set to zero and parameters have been replaced with their estimates. For h=2, the forecasts are given by \begin{align*}
\hat y_{1,T+2|T} &=\hat{c}_1+\hat\phi_{11,1}\hat y_{1,T+1|T}+\hat\phi_{12,1}\hat y_{2,T+1|T}\\
\hat y_{2,T+2|T}&=\hat{c}_2+\hat\phi_{21,1}\hat y_{1,T+1|T}+\hat\phi_{22,1}\hat y_{2,T+1|T}.
\end{align*} Again, this is the same form as Equations 12.1, except that the errors have been set to zero, the parameters have been replaced with their estimates, and the unknown values of y_1 and y_2 have been replaced with their forecasts. The process can be iterated in this manner for all future time periods.
There are two decisions one has to make when using a VAR to forecast, namely how many variables (denoted by K) and how many lags (denoted by p) should be included in the system. The number of coefficients to be estimated in a VAR is equal to K+pK^2 (or 1+pK per equation). For example, for a VAR with K=5 variables and p=3 lags, there are 16 coefficients per equation, giving a total of 80 coefficients to be estimated. The more coefficients that need to be estimated, the larger the estimation error entering the forecast.
In practice, it is usual to keep K small and include only variables that are correlated with each other, and therefore useful in forecasting each other. Information criteria are commonly used to select the number of lags to be included. Care should be taken when using the AICc as it tends to choose large numbers of lags; instead, for VAR models, we often use the BIC instead. A more sophisticated version of the model is a “sparse VAR” (where many coefficients are set to zero); another approach is to use “shrinkage estimation” (where coefficients are smaller).
A criticism that VARs face is that they are atheoretical; that is, they are not built on some economic theory that imposes a theoretical structure on the equations. Every variable is assumed to influence every other variable in the system, which makes a direct interpretation of the estimated coefficients difficult. Despite this, VARs are useful in several contexts:
forecasting a collection of related variables where no explicit interpretation is required;
testing whether one variable is useful in forecasting another (the basis of Granger causality tests);
impulse response analysis, where the response of one variable to a sudden but temporary change in another variable is analysed;
forecast error variance decomposition, where the proportion of the forecast variance of each variable is attributed to the effects of the other variables.
Example: A VAR model for forecasting US consumption
You can use a VAR model in Python via the statsmodels.tsa.api module. Note that you need to set the dates as the index of the DataFrame before fitting the model. We will fit two models: one using the AIC and another using the BIC.
A VAR(5) model is selected using the AICc (the default), while a VAR(1) model is selected using the BIC. This is not unusual — the BIC will always select a model that has fewer parameters than the AICc model as it imposes a stronger penalty for the number of parameters.
Figure 12.14: ACF of the residuals from the two VAR models. A VAR(5) model is selected by the AICc, while a VAR(1) model is selected using the BIC.
From Figure 12.14, we see that the residuals from the VAR(1) model (BIC) have significant autocorrelation for Consumption, while the VAR(5) model has effectively captured all the information in the data.
The forecasts generated by the VAR(5) model are plotted in Figure 12.15. Note that forecast_interval generates both the mean forecast and the (1-\alpha)\% prediction interval.
Figure 12.15: Forecasts for US consumption and income generated from a VAR(5) model.
12.4 Bootstrapping and bagging
Bootstrapping time series
In Section 5.5, we bootstrap the residuals of a time series in order to simulate future values of a series using a model. More generally, we can generate new time series that are similar to our observed series, using another type of bootstrap.
First, the time series is transformed if necessary, and then decomposed into trend, seasonal and remainder components using STL. Then we obtain shuffled versions of the remainder component to get bootstrapped remainder series. Because there may be autocorrelation present in an STL remainder series, we cannot simply use the re-draw procedure that was described in Section 5.5. Instead, we use a “blocked bootstrap”, where contiguous sections of the time series are selected at random and joined together. These bootstrapped remainder series are added to the trend and seasonal components, and the transformation is reversed to give variations on the original time series.
Consider the quarterly cement production in Australia from 1988 Q1 to 2010 Q2. First we check (see Figure 12.16) that the decomposition has adequately captured the trend and seasonality, and that there is no obvious remaining signal in the remainder series. (In StatsForecast, an STL decomposition is defined by specifying a single seasonality in MSTL()).
Figure 12.16: STL decomposition of quarterly Australian cement production.
Now we can generate several bootstrapped versions of the data. We will use a block size of 8 to cover two years of data and the blocked_bootstrap() function to generate the simulated data from the STL decomposition we have just obtained.
Code
def blocked_bootstrap(dcmp, block_size, num_sims): sim = {} num_blocks = np.ceil(len(dcmp) / block_size).astype(int)for i inrange(num_sims): blocks = []for _ inrange(num_blocks): start = np.random.randint(0, len(dcmp) - block_size +1) end = start + block_size blocks.append(dcmp["remainder"].iloc[start:end]) sim[str(i)] = ( dcmp["trend"]+ dcmp["seasonal"]+ pd.concat(blocks, ignore_index=True)[0 : len(dcmp)] ) sim_df = pd.DataFrame(sim)return sim_dfsim_df = blocked_bootstrap(dcmp, block_size=8, num_sims=10).assign( unique_id=cement.unique_id.iloc[0], ds=cement["ds"].to_numpy(),)fig, ax = plt.subplots()plot_series(cement, sim_df, xlabel="Quarter [1Q]", ylabel="Tonnes ('000)", palette="tab10", ax=ax)# Re-draw original `cement` time series on top and in blackax.plot("ds", "y", "k-", data=cement)fig
Figure 12.17: Ten bootstrapped versions of quarterly Australian cement production (coloured), along with the original data (black).
Bagged forecasts
One use for these bootstrapped time series is to improve forecast accuracy. If we produce forecasts from each of the additional time series, and average the resulting forecasts, we get better forecasts than if we simply forecast the original time series directly. This is called “bagging” which stands for “bootstrap aggregating”.
We demonstrate the idea using the cement data. First, we simulate many time series that are similar to the original data using bootstrapping.
For each of these series, we fit an ETS model. A different ETS model is selected for each case, although all models will be similar because the series themselves are similar. However, since the estimated parameters will be slightly different, the forecasts will also be different. This can be a time-consuming process due to the large number of series involved.
Figure 12.18: Forecasts of 100 bootstrapped series obtained using AutoETS models.
Finally, we average these forecasts for each time period to obtain the “bagged forecasts” for the original data. We also generate an ETS model applied directly to the data for comparison.
Figure 12.19: Comparing bagged AutoETS forecasts (the average of 100 bootstrapped forecasts) and AutoETS applied directly to the data.
Bergmeir et al. (2016) show that, on average, bagging gives better forecasts than just applying AutoETS() directly. Of course, it is slower because a lot more computation is required.
12.5 Exercises
Compare STL and Dynamic Harmonic Regression forecasts for one of the series in the pedestrian data set.
Try modifying the order of the Fourier terms to minimise the AICc value.
Check the residuals for each model. Do they capture the available information in the data?
Which of the two sets of forecasts are best? Explain.
Consider the weekly data on US finished motor gasoline products supplied (millions of barrels per day) (series us_gasoline).
Use the Prophet model to forecast the next 13 weeks of data.
Repeat the previous exercise using the AutoARIMA() and the AutoETS() classes from StatsForecast.
Compare the accuracy of the models you generated.
12.6 Further reading
The Prophet model is described in Taylor & Letham (2018).
Pfaff (2008) provides a book-length overview of VAR modelling and other multivariate time series models.
A current survey of the use of recurrent neural networks for forecasting is provided by Hewamalage et al. (2021).
Bootstrapping for time series is discussed in Lahiri (2003).
Bagging for time series forecasting is relatively new. Bergmeir et al. (2016) is one of the few papers which addresses this topic.
Athanasopoulos, G, Poskitt, DS, & Vahid, F. (2012). Two canonical VARMA forms: Scalar component models vis-à-vis the echelon form. Econometric Reviews, 31(1), 60–83.
Bergmeir, C, Hyndman, RJ, & Benítez, JM. (2016). Bagging exponential smoothing methods using STL decomposition and Box-Cox transformation. International Journal of Forecasting, 32(2), 303–312.
Fan, S, & Hyndman, RJ. (2012). Short-term load forecasting based on a semi-parametric additive model. IEEE Transactions on Power Systems, 27(1), 134–141.
Garza, A, Mergenthaler Canseco, M, Challú, C, & Olivares, KG. (2022). StatsForecast: Lightning fast forecasting with statistical and econometric models. PyCon Salt Lake City, Utah, US 2022.
Hamilton, JD. (1994). Time series analysis. Princeton University Press, Princeton.
Hewamalage, H, Bergmeir, C, & Bandara, K. (2021). Recurrent neural networks for time series forecasting: Current status and future directions. International Journal of Forecasting, 37(1), 388–427.
Hyndman, RJ, & Fan, S. (2010). Density forecasting for long-term peak electricity demand. IEEE Transactions on Power Systems, 25(2), 1142–1153.
Lahiri, SN. (2003). Resampling methods for dependent data. Springer Science & Business Media.
Lütkepohl, H. (2007). General-to-specific or specific-to-general modelling? An opinion on current econometric terminology. Journal of Econometrics, 136(1), 234–319.
Menculini, L, Marini, A, Proietti, M, Garinei, A, Bozza, A, Moretti, C, & Marconi, M. (2021). Comparing prophet and deep learning to ARIMA in forecasting wholesale food prices. Forecasting, 3(3), 644–662.
Pfaff, B. (2008). Analysis of integrated and cointegrated time series with R. Springer Science & Business Media.
Taylor, SJ, & Letham, B. (2018). Forecasting at scale. The American Statistician, 72(1), 37–45.
A more flexible generalisation would be a Vector ARMA process. However, the relative simplicity of VARs has led to their dominance in forecasting. Interested readers may refer to Athanasopoulos et al. (2012).↩︎
Interested readers should refer to Hamilton (1994) and Lütkepohl (2007).↩︎
