I get this question a lot, so I thought it might help to explain some issues with AIC calculation.

First, the equation for the AIC is given by where is the likelihood of the model and is the number of parameters that are estimated (including the error variance). For a linear regression model with iid errors, fitted to observations, the log-likelihood can be written as where is the residual for the th observation. The AIC is then Since we don’t know , we estimate it using the mean squared error (the maximum likelihood estimator), giving where is a constant that depends only on the sample size and not on the model. This constant is often ignored. Thus, different software implementations can lead to different AIC values for the same model, since they may include or exclude the constant .

Now, let’s look at what R returns in a simple case using the lm() function.

set.seed(2023)
library(fpp3)
df <- tibble(
  time = seq(100),
  x = rnorm(100),
  y = x + rnorm(100)
)
fit <- lm(y ~ x, data = df)
AIC(fit)

[1] 275.6267

We can check how this is calculated by computing it ourselves.

mse <- mean(residuals(fit)^2)
n <- length(residuals(fit))
k <- length(fit$coefficients) + 1
# With constant
2*k + n*log(mse) + n*log(2*pi) + n

[1] 275.6267

# Without constant
2*k + n*log(mse)

[1] -8.161047

Clearly, AIC() applied to the output from lm() is using the version with the constant.

Now compare that with what we obtain using the TSLM() function from the fable package.

df |>
  as_tsibble(index = time) |>
  model(TSLM(y ~ x)) |>
  glance() |>
  pull(AIC)

[1] -8.161047

This is the AIC without the constant.

The situation is even more confusing with ARIMA models, and some other model classes, because some functions use approximations to the likelihood, rather than the exact likelihood.

Thus, AIC values can be compared across models fitted using the same functions, but not necessarily when models have been fitted using different functions.

I received this email today:

My team recently used some techniques found in your writings to perform forecasts … Our work has been well received by reviewers, but one commenter asked two questions that I was hoping you may be able to provide insight on.

First, they wanted to know if we could provide P-values for our prediction intervals. In our work, we said, “Observed rates were deemed significantly different from expected rates when they did not fall within the 95% PI.” This same language has been used by others published in the same journal. I am curious to hear your thoughts on giving P-values for these PIs and what the appropriate method for doing so would be (if any).

Second, they asked about making a correction for multiple comparisons. … I believe we could apply a Bonferroni correction to the PIs, but that feels too liberal. Moreover, I am curious if this is even called for given our statement of what is deemed significant and the fact that our prediction interval construction relies on a non-parametric method.

Here is my reply:

1. I don’t think this makes any sense. A p-value is the probability of obtaining observations at least as extreme as those observed given a null hypothesis. What’s the hypothesis here? In forecasting, we don’t usually have a hypothesis. Instead, we fit a model to the data, and make predictions based on the model. I guess you could make the null hypothesis “The future observations come from the forecast distributions”, and then the p-value for each future time period would be the probability of the tails beyond the observations. But it is well-known that the estimated prediction intervals are almost always too narrow due to them not taking into account all sources of variance. So the size of this test would not be well-calibrated. I think you’re better off pushing back rather than trying to meet the request.

2. A Bonferroni correction assumes independence between the intervals, and that is not true for PIs from a forecasting model. The future forecast errors are all correlated (with the strength of the correlation depending on the model and the DGP). Usually we just say that these are pointwise PI, and so we expect 5% of observations to fall outside the 95% prediction intervals. It is possible to generate uniform PI, which contain 95% of all future sample paths, but this is a little tricky due to the correlations between horizons. It could be done via simulation – simulate a 1000 future sample paths and compute the envelope that contains 950 of them.

It sounds like the reviewers are only familiar with inferential statistics, and not with predictive modelling. You could point them to Shmueli’s excellent 2010 paper “To explain or predict”, highlighting the differences between the two paradigms.

Forecast combinations have flourished remarkably in the forecasting community and, in recent years, have become part of the mainstream of forecasting research and activities. Combining multiple forecasts produced from the single (target) series is now widely used to improve accuracy through the integration of information gleaned from different sources, thereby mitigating the risk of identifying a single “best” forecast. Combination schemes have evolved from simple combination methods without estimation, to sophisticated methods involving time-varying weights, nonlinear combinations, correlations among components, and cross-learning. They include combining point forecasts, and combining probabilistic forecasts. This paper provides an up-to-date review of the extensive literature on forecast combinations, together with reference to available open-source software implementations. We discuss the potential and limitations of various methods and highlight how these ideas have developed over time. Some important issues concerning the utility of forecast combinations are also surveyed. Finally, we conclude with current research gaps and potential insights for future research.

This question arises most time I teach a forecasting workshop, and it was raised again in the following email I received today:

I have a time series that I have split into training and test datasets with an 80%-20% ratio. I fit a series of different models (ETS, BATS, ARIMA, NN etc) to the training data and generate my forecasts from each model. When evaluating the forecasts against the test set I find the model that gives the best outcome is an ARIMA(1,1,1) that was selected using the auto.arima function. My question is this, should I proceed to fit an ARIMA(1,1,1) to the whole data set, or should I use the auto.arima function again which may give me a slightly different (p,d,q) order as there is now an extra 20% of unseen data available to the forecast model? Any guidance would be greatly received.

This is a common issue. If you only have one class of model to consider (e.g., only ETS or only ARIMA), then it is easy enough to select the model on all available data using the AIC, and use the selected model to forecast the future. But if you are selecting between model classes, then you need to use either a training/test split, or (preferably) a time-series cross-validation procedure.

If you use time-series cross-validation, then there would usually be different models selected for each training set, and the cross-validated error is a measure of how well the model class works for your data. In that case, there is no single model for the training data, and you are selecting the model class rather than a specific model. This makes it clear that you should then apply the selected model class to all the data, when forecasting beyond the end of the available data. In other words, if you choose ARIMA over ETS, then you would then fit an ARIMA model to all the data, and use that model to forecast the future.

You can think of a simple training/test split as a special case of time-series cross-validation, where there is a single fold. So the same argument applies. That is, you are selecting the model class that works best for your data, and so you should apply that model class to all the data, when forecasting beyond the end of the available data.

This example also illustrates why it is important to use a time-series cross-validation procedure, rather than a simple training/test split. In this case, the ARIMA model was selected because it happened to work best for the particular training/test split that was used. But if a different split had been used, then a different model might have been selected. So the model selection is not stable. By averaging over multiple folds using a time-series cross-validation procedure, you can get a more stable estimate of the model class that works best for your data.

One of the most challenging aspects for managers when building a forecasting system is choosing how to aggregate the data at different levels. This is frequently done without the manager knowing how these choices can compromise the system’s accuracy. This paper illustrates these compromises by comparing different structures and aggregation criteria. Our paper proposes and empirically tests a framework on how to build a coherent and more accurate forecasting system. The framework’s first phase compares different time series forecasting methods, including statistical, “standard” machine learning, and deep learning. Results show that one of the statistical methods (autoregressive integrated moving average, or, for short, ARIMA) outperforms machine and deep learning methods. The second phase compares different combinations of aggregation criteria, structures of the forecasting system, and coherent forecast methods (i.e., adjustments to the forecasts at different levels of aggregation). The results show that using different criteria and structures indeed impacts predictions’ accuracy. When it is necessary to disaggregate the forecast, our results show that it is best to add more information in a grouped structure, adjusted by a bottom-up method. This combination provides the best performance, i.e., the lowest mean absolute scaled error (MASE) in most nodes, compared to the other structures and coherent forecast methods used. The results also suggest that aggregating the time series further by geographical regions is essential to improve accuracy when forecasting products’ and channels’ sales.

This paper discusses the use of forecast reconciliation with stock price time series and the corresponding stock index. The individual stock price series may be grouped using known meta-data or other clustering methods. We propose a novel forecasting framework that combines forecast reconciliation and clustering, to lead to better forecasts of both the index and the individual stock price series. The proposed approach is applied to the Dow Jones Industrial Average Index and its component stocks. The results demonstrate empirically that reconciliation improves forecasts of the stock market index and its constituents.

We have taught from the book many times, but this year we decided to pre-record short videos for each section. Our students often prefer a video explanation than reading the textbook, and we thought other readers might appreciate hearing from us as well.

These videos are embedded in most sections of the book. So far, we’ve covered the sections that we include in our own courses, but we hope to eventually have videos for all sections. Most of these were done in a single take, so they are sometimes a little rough, but hopefully still useful.

You can view the entire playlist on YouTube.

(Postponed until further notice)

Distinguished Lecture Series for the International Institute of Forecasters

Global Forecasting Models (GFM) that are trained across a set of multiple time series have shown superior results in many forecasting competitions and real-world applications compared with univariate forecasting approaches. One aspect of the popularity of statistical forecasting models such as ETS and ARIMA is their relative simplicity and interpretability (in terms of relevant lags, trend, seasonality, and others), while GFMs typically lack interpretability, especially towards particular time series. This reduces the trust and confidence of the stakeholders when making decisions based on the forecasts without being able to understand the predictions. To mitigate this problem, in this work, we propose a novel local model-agnostic interpretability approach to explain the forecasts from GFMs. We train simpler univariate surrogate models that are considered interpretable (e.g., ETS) on the predictions of the GFM on samples within a neighbourhood that we obtain through bootstrapping or straightforwardly as the one-step-ahead global black-box model forecasts of the time series which needs to be explained. After, we evaluate the explanations for the forecasts of the global models in both qualitative and quantitative aspects such as accuracy, fidelity, stability and comprehensibility, and are able to show the benefits of our approach.

Providing forecasts for ultra-long time series plays a vital role in various activities, such as investment decisions, industrial production arrangements, and farm management. This paper develops a novel distributed forecasting framework to tackle challenges associated with forecasting ultra-long time series by using the industry-standard MapReduce framework. The proposed model combination approach facilitates distributed time series forecasting by combining the local estimators of time series models delivered from worker nodes and minimizing a global loss function. In this way, instead of unrealistically assuming the data generating process (DGP) of an ultra-long time series stays invariant, we make assumptions only on the DGP of subseries spanning shorter time periods. We investigate the performance of the proposed approach with AutoRegressive Integrated Moving Average (ARIMA) models using the real data application as well as numerical simulations. Compared to directly fitting the whole data with ARIMA models, our approach results in improved forecasting accuracy and computational efficiency both in point forecasts and prediction intervals, especially for longer forecast horizons. Moreover, we explore some potential factors that may affect the forecasting performance of our approach.

Accurate forecasts of ambulance demand are crucial inputs when planning and deploying staff and fleet. Such demand forecasts are required at national, regional and sub-regional levels, and must take account of the nature of incidents and their priorities. These forecasts are often generated independently by different teams within the organization. As a result, forecasts at different levels may be inconsistent, resulting in conflicting decisions and a lack of coherent coordination in the service. To address this issue, we exploit the hierarchical and grouped structure of the demand time series, and apply forecast reconciliation methods to generate both point and probabilistic forecasts that are coherent and use all the available data at all levels of disaggregation. The methods are applied to daily incident data from an ambulance service in Great Britain, from October 2015 to July 2019, disaggregated by nature of incident, priority, managing health board, and control area. We use an ensemble of forecasting models, and show that the resulting forecasts are better than any individual forecasting model. We validate the forecasting approach using time-series cross-validation.

Real-time monitoring using in-situ sensors is becoming a common approach for measuring water-quality within watersheds. High-frequency measurements produce big datasets that present opportunities to conduct new analyses for improved understanding of water-quality dynamics and more effective management of rivers and streams. Of primary importance is enhancing knowledge of the relationships between nitrate, one of the most reactive forms of inorganic nitrogen in the aquatic environment, and other water-quality variables. We analysed high-frequency water-quality data from in-situ sensors deployed in three sites from different watersheds and climate zones within the National Ecological Observatory Network, USA. We used generalised additive mixed models to explain the nonlinear relationships at each site between nitrate concentration and conductivity, turbidity, dissolved oxygen, water temperature, and elevation. Temporal auto-correlation was modelled with an auto-regressive–moving-average (ARIMA) model and we examined the relative importance of the explanatory variables. Total deviance explained by the models was high for all sites (99%). Although variable importance and the smooth regression parameters differed among sites, the models explaining the most variation in nitrate contained the same explanatory variables. This study demonstrates that building a model for nitrate using the same set of explanatory water-quality variables is achievable, even for sites with vastly different environmental and climatic characteristics. Applying such models will assist managers to select cost-effective water-quality variables to monitor when the goals are to gain a spatial and temporal in-depth understanding of nitrate dynamics and adapt management plans accordingly.

Commentary on Basellini, Camarda & Booth (2022) Thirty years on: A review of the Lee-Carter method for forecasting mortality, International Journal of Forecasting, in press. [preprint version]

I’ve been interviewed for several podcasts over the last few years. It’s always fun to talk about my work, and I hope there are enough differences between them to make it interesting for listeners. Here is a full list of them.

Time series often reflect variation associated with other related variables. Controlling for the effect of these variables is useful when modeling or analysing the time series. We introduce a novel approach to normalize time series data conditional on a set of covariates. We do this by modeling the conditional mean and the conditional variance of the time series with generalized additive models using a set of covariates. The conditional mean and variance are then used to normalize the time series. We illustrate the use of conditionally normalized series using two applications involving river network data. First, we show how these normalized time series can be used to impute missing values in the data. Second, we show how the normalized series can be used to estimate the conditional autocorrelation function and conditional cross-correlation functions via additive models. Finally we use the conditional cross-correlations to estimate the time it takes water to flow between two locations in a river network.

R package: conduits

Collections of time series that are formed via aggregation are prevalent in many fields. These are commonly referred to as hierarchical time series and may be constructed cross-sectionally across different variables, temporally by aggregating a single series at different frequencies, or may even be generalised beyond aggregation as time series that respect linear constraints. When forecasting such time series, a desirable condition is for forecasts to be coherent, that is to respect the constraints. The past decades have seen substantial growth in this field with the development of reconciliation methods that not only ensure coherent forecasts but can also improve forecast accuracy. This paper serves as both an encyclopaedic review of forecast reconciliation and an entry point for researchers and practitioners dealing with hierarchical time series. The scope of the article includes perspectives on forecast reconciliation from machine learning, Bayesian statistics and probabilistic forecasting as well as applications in economics, energy, tourism, retail demand and demography.