Chapter 4 Forecasting

Forecasting has always been an important part of the time series field (De Gooijer and Hyndman (2006)). Macroeconomic forecasts are done in many places: Public Administration (notably Treasuries), Central Banks, International Institutions (e.g. IMF, OECD), banks, big firms. These institutions are interested in the point estimates (\(\sim\) most likely value) of the variable of interest. They also sometimes need to measure the uncertainty (\(\sim\) dispersion of likely outcomes) associated to the point estimates.¹³

Forecasts produced by professional forecasters are available on these web pages:

• Philly Fed Survey of Professional Forecasters.
• ECB Survey of Professional Forecasters.
• IMF World Economic Outlook.
• OECD Global Economic Outlook.
• European Commission Economic Forecasts.

How to formalize the forecasting problem? Assume the current date is \(t\). We want to forecast the value that variable \(y_t\) will take on date \(t+1\) (i.e., \(y_{t+1}\)) based on the observation of a set of variables gathered in vector \(x_t\) (\(x_t\) may contain lagged values of \(y_t\)).

The forecaster aims at minimizing (a function of) the forecast error. It is usal to consider the following (quadratic) loss function: \[ \underbrace{\mathbb{E}([y_{t+1} - y^*_{t+1}]^2)}_{\mbox{Mean square error (MSE)}} \] where \(y^*_{t+1}\) is the forecast of \(y_{t+1}\) (function of \(x_t\)).

Proposition 4.1 (Smallest MSE) The smallest MSE is obtained with the expectation of \(y_{t+1}\) conditional on \(x_t\).

Proof. See Appendix 8.4.

Proposition 4.2 Among the class of linear forecasts, the smallest MSE is obtained with the linear projection of \(y_{t+1}\) on \(x_t\). This projection, denoted by \(\hat{P}(y_{t+1}|x_t):=\boldsymbol\alpha'x_t\), satisfies: \[\begin{equation} \mathbb{E}\left( [y_{t+1} - \boldsymbol\alpha'x_t]x_t \right)=\mathbf{0}.\tag{4.1} \end{equation}\]

Proof. Consider the function \(f:\) \(\boldsymbol\alpha \rightarrow \mathbb{E}\left( [y_{t+1} - \boldsymbol\alpha'x_t]^2 \right)\). We have: \[ f(\boldsymbol\alpha) = \mathbb{E}\left( y_{t+1}^2 - 2 y_t x_t'\boldsymbol\alpha + \boldsymbol\alpha'x_t x_t'\boldsymbol\alpha] \right). \] We have \(\partial f(\boldsymbol\alpha)/\partial \boldsymbol\alpha = \mathbb{E}(-2 y_{t+1} x_t + 2 x_t x_t'\boldsymbol\alpha)\). The function is minimised for \(\partial f(\boldsymbol\alpha)/\partial \boldsymbol\alpha =0\).

Eq. (4.1) implies that \(\mathbb{E}\left( y_{t+1}x_t \right)=\mathbb{E}\left(x_tx_t' \right)\boldsymbol\alpha\). (Note that \(x_t x_t'\boldsymbol\alpha=x_t (x_t'\boldsymbol\alpha)=(\boldsymbol\alpha'x_t) x_t'\).)

Hence, if \(\mathbb{E}\left(x_tx_t' \right)\) is nonsingular, \[\begin{equation} \boldsymbol\alpha=[\mathbb{E}\left(x_tx_t' \right)]^{-1}\mathbb{E}\left( y_{t+1}x_t \right).\tag{4.2} \end{equation}\]

The MSE then is: \[ \mathbb{E}([y_{t+1} - \boldsymbol\alpha'x_t]^2) = \mathbb{E}{(y_{t+1}^2)} - \mathbb{E}\left( y_{t+1}x_t' \right)[\mathbb{E}\left(x_tx_t' \right)]^{-1}\mathbb{E}\left(x_ty_{t+1} \right). \]

Consider the regression \(y_{t+1} = \boldsymbol\beta'\mathbf{x}_t + \varepsilon_{t+1}\). The OLS estimate is: \[ \mathbf{b} = \left[ \underbrace{ \frac{1}{T} \sum_{i=1}^T \mathbf{x}_t\mathbf{x}_t'}_{\mathbf{m}_1} \right]^{-1}\left[ \underbrace{ \frac{1}{T} \sum_{i=1}^T \mathbf{x}_t'y_{t+1}}_{\mathbf{m}_2} \right]. \] If \(\{x_t,y_t\}\) is covariance-stationary and ergodic for the second moments then the sample moments (\(\mathbf{m}_1\) and \(\mathbf{m}_2\)) converges in probability to the associated population moments and \(\mathbf{b} \overset{p}{\rightarrow} \boldsymbol\alpha\) (where \(\boldsymbol\alpha\) is defined in Eq. (4.2)).

Example 4.1 (Forecasting an MA(q) process) Consider the MA(q) process: \[ y_t = \mu + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \dots + \theta_q \varepsilon_{t-q}, \] where \(\{\varepsilon_t\}\) is a white noise sequence (Def. 1.1).

We have:¹⁴ \[\begin{eqnarray*} &&\mathbb{E}(y_{t+h}|\varepsilon_{t},\varepsilon_{t-1},\dots) =\\ &&\left\{ \begin{array}{lll} \mu + \theta_h \varepsilon_{t} + \dots + \theta_q \varepsilon_{t-q+h} \quad &for& \quad h \in [1,q]\\ \mu \quad &for& \quad h > q \end{array} \right. \end{eqnarray*}\] and \[\begin{eqnarray*} &&\mathbb{V}ar(y_{t+h}|\varepsilon_{t},\varepsilon_{t-1},\dots)= \mathbb{E}\left( [y_{t+h} - \mathbb{E}(y_{t+h}|\varepsilon_{t},\varepsilon_{t-1},\dots)]^2 \right) =\\ &&\left\{ \begin{array}{lll} \sigma^2(1+\theta_1^2+\dots+\theta_{h-1}^2) \quad &for& \quad h \in [1,q]\\ \sigma^2(1+\theta_1^2+\dots+\theta_q^2) \quad &for& \quad h>q. \end{array} \right. \end{eqnarray*}\]

Example 4.2 (Forecasting an AR(p) process) (See this web interface.) Consider the AR(p) process: \[ y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \dots + \phi_p y_{t-p} + \varepsilon_t, \] where \(\{\varepsilon_t\}\) is a white noise sequence (Def. 1.1).

Using the notation of Eq. (2.4), we have: \[ \mathbf{y}_t - \boldsymbol\mu = F (\mathbf{y}_{t-1}- \boldsymbol\mu) + \boldsymbol\xi_t, \] with \(\boldsymbol\mu = [\mu,\dots,\mu]'\) (\(\mu\) is defined in Eq. (2.8)). Hence: \[ \mathbf{y}_{t+h} - \boldsymbol\mu = \boldsymbol\xi_{t+h} + F \boldsymbol\xi_{t+h-1} + \dots + F^{h-1} \boldsymbol\xi_{t+1} + F^h (\mathbf{y}_{t}- \mu). \] Therefore: \[\begin{eqnarray*} \mathbb{E}(\mathbf{y}_{t+h}|y_{t},y_{t-1},\dots) &=& \boldsymbol\mu + F^{h}(\mathbf{y}_t - \boldsymbol\mu)\\ \mathbb{V}ar\left( [\mathbf{y}_{t+h} - \mathbb{E}(\mathbf{y}_{t+h}|y_{t},y_{t-1},\dots)] \right) &=& \Sigma + F\Sigma F' + \dots + F^{h-1}\Sigma (F^{h-1})', \end{eqnarray*}\] where: \[ \Sigma = \left[ \begin{array}{ccc} \sigma^2 & 0& \dots\\ 0 & 0 & \\ \vdots & & \ddots \\ \end{array} \right]. \]

Alternative approach: Taking the (conditional) expectations of both sides of \[ y_{t+h} - \mu = \phi_1 (y_{t+h-1} - \mu) + \phi_2 (y_{t+h-2} - \mu) + \dots + \phi_p (y_{t-p} - \mu) + \varepsilon_{t+h}, \] we obtain: \[\begin{eqnarray*} \mathbb{E}(y_{t+h}|y_{t},y_{t-1},\dots) &=& \mu + \phi_1\left(\mathbb{E}[y_{t+h-1}|y_{t},y_{t-1},\dots] - \mu\right)+\\ &&\phi_2\left(\mathbb{E}[y_{t+h-2}|y_{t},y_{t-1},\dots] - \mu\right) + \dots +\\ && \phi_p\left(\mathbb{E}[y_{t+h-p}|y_{t},y_{t-1},\dots] - \mu\right), \end{eqnarray*}\] which can be exploited recursively.

The recursion begins with \(\mathbb{E}(y_{t-k}|y_{t},y_{t-1},\dots)=y_{t-k}\) (for any \(k \ge 0\)).

Example 4.3 (Forecasting an ARMA(p,q) process) Consider the process: \[\begin{equation} y_t = c + \phi_1 y_{t-1} + \dots + \phi_p y_{t-p} + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \dots + \theta_q \varepsilon_{t-q},\tag{4.3} \end{equation}\] where \(\{\varepsilon_t\}\) is a white noise sequence (Def. 1.1). We assume that the MA part of the process is invertible (see Eq. (2.24)), which implies that the information contained in \(\{y_{t},y_{t-1},y_{t-2},\dots\}\) is identical to that in \(\{\varepsilon_{t},\varepsilon_{t-1},\varepsilon_{t-2},\dots\}\).

While one could use a recursive algorithm to compute the conditional mean (as in Example 4.2), it is convenient to employ the Wold decomposition of this process (see Theorem 2.2 and Prop. 2.7 for the computation of the \(\psi_i\)’s in the context of ARMA processes): \[ y_t = \mu + \sum_{i=0}^{+\infty} \psi_i \varepsilon_{t-i}. \] This implies: \[\begin{eqnarray*} y_{t+h} &=& \mu + \sum_{i=0}^{h-1} \psi_i \varepsilon_{t+h-i} + \sum_{i=h}^{+\infty} \psi_i \varepsilon_{t+h-i}\\ &=& \mu + \sum_{i=0}^{h-1} \psi_i \varepsilon_{t+h-i} + \sum_{i=0}^{+\infty} \psi_{i+h} \varepsilon_{t-i}. \end{eqnarray*}\]

Since \(\mathbb{E}(y_{t+h}|y_t,y_{t-1},\dots)=\mu+\sum_{i=0}^{+\infty} \psi_{i+h} \varepsilon_{t-i}\), we get: \[ \mathbb{V}ar(y_{t+h}|y_t,y_{t-1},\dots) =\mathbb{V}ar\left(\sum_{i=0}^{h-1} \psi_i \varepsilon_{t+h-i}\right)= \sigma^2 \sum_{i=0}^{h-1} \psi_i^2. \]

How to use the previous formulas in practice?

One has first to select a specification and to estimate the model. Two methods to determine relevant specifications:

1. Information criteria (see Definition 2.11).
2. Box-Jenkins approach.

Box and Jenkins (1976) have proposed an approach that is now widely used.

1. Data transformation. The data should be transformed to “make them stationary”. To do so, one can e.g. take logarithms, take changes in the considered series, remove (deterministic) trends.
2. Select \(p\) and \(q\). This can be based on the PACF approach (see Section 2.4), or on selection criteria (see Definition 2.11).
3. Estimate the model parameters. See Section 2.8.
4. Check that the estimated model is consistent with the data. See below.

Assessing the performances of a forecasting model

Once one has fitted a model on a given dataset (of length \(T\), say), one compute MSE (mean square errors) to evaluate the performance of the model. But this MSE is the in-sample one. It is easy to reduce in-sample MSE. Typically, if the model is estimated by OLS, adding covariates mechanically reduces the MSE. That is, even if additional data are irrelevant, the \(R^2\) of the regression increases. Adding irrelevant variables increases the (in-sample) \(R^2\) but is bound to increase the out-of-sample MSE.

Therefore, it is important to analyse out-of-sample performances of the forecasting model:

1. Estimate a model on a sample of reduced size (\(1,\dots,T^*\), with \(T^*<T\))
2. Use the remaining available periods (\(T^*+1,\dots,T\)) to compute out-of-sample forecasting errors (and compute their MSE). In an out-of-sample exercise, it is important to make sure that the data used to produce a forecasts (as of date \(T^*\)) where indeed available on date \(T^*\).

Diebold-Mariano test

How to compare different forecasting approaches? Diebold and Mariano (1995) have proposed a simple test to address this question.

Assume that you want to compare approaches A and B. You have historical data sets and you have implemented both approaches in the past, providing you with two sets of forecasting errors: \(\{e^{A}_t\}_{t=1,\dots,T}\) and \(\{e^{B}_t\}_{t=1,\dots,T}\).

It may be the case that your forecasts serve a specific purpose and that, for instance, you dislike positive forecasting errors and you care less about negative errors. We assume you are able to formalise this by means of a loss function \(L(e)\). For instance:

• If you dislike large positive errors, you may set \(L(e)=\exp(e)\).
• If you are concerned about both positive and negative errors (indifferently), you may set \(L(e)=e^2\) (standard approach).

Let us define the sequence \(\{d_t\}_{t=1,\dots,T} \equiv \{L(e^{A}_t)-L(e^{B}_t)\}_{t=1,\dots,T}\) and assume that this sequence is covariance stationary. We consider the following null hypothesis: \(H_0:\) \(\bar{d}=0\), where \(\bar{d}\) denotes the population mean of the \(d_t\)s. Under \(H_0\) and under the assumption of covariance-stationarity of \(d_t\), we have (Theorem 1.1): \[ \sqrt{T} \bar{d}_T \overset{d}{\rightarrow} \mathcal{N}\left(0,\sum_{j=-\infty}^{+\infty} \gamma_j \right), \]
where the \(\gamma_j\)s are the autocovariances of \(d_t\).

Hence, assuming that \(\hat{\sigma}^2\) is a consistent estimate of \(\sum_{j=-\infty}^{+\infty} \gamma_j\) (for instance the one given by the Newey-West formula, see Def. (1.6)), we have, under \(H_0\): \[ DM_T := \sqrt{T}\frac{\bar{d}_T}{\sqrt{\hat{\sigma}^2}} \overset{d}{\rightarrow} \mathcal{N}(0,1). \] \(DM_T\) is the test statistics. For a test of size \(\alpha\), the critical region is:¹⁵ \[ ]-\infty,-\Phi^{-1}(1-\alpha/2)] \cup [\Phi^{-1}(1-\alpha/2),+\infty[, \] where \(\Phi\) is the c.d.f. of the standard normal distribution.

Example 4.4 (Forecasting Swiss GDP growth) We use a long historical time series of the Swiss GDP growth taken from the Jordà, Schularick, and Taylor (2017) dataset (see Figure 1.3, and Example 2.4).

We want to forecast this GDP growth. We envision two specifications : an AR(1) specification (the one advocated by the AIC criteria, see Example 2.4), and an ARMA(2,2) specification. We are interested in 2-year-ahead forecasts (i.e., \(h=2\) since the data are yearly).

Code

library(IdSS)
library(forecast)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

Code

data <- subset(JST,iso=="CHE")
T <- dim(data)[1]
y <- c(NaN,log(data$gdp[2:T]/data$gdp[1:(T-1)]))
first.date <- T-50
e1 <- NULL; e2 <- NULL;h<-2
for(T.star in first.date:(T-h)){
  estim.model.1 <- arima(y[1:T.star],order=c(1,0,0))
  estim.model.2 <- arima(y[1:T.star],order=c(2,0,2))
  e1 <- c(e1,y[T.star+h] - predict(estim.model.1,n.ahead=h)$pred[h])
  e2 <- c(e2,y[T.star+h] - predict(estim.model.2,n.ahead=h)$pred[h])
}
res.DM <- dm.test(e1,e2,h = h,alternative = "greater")
res.DM

## 
##  Diebold-Mariano Test
## 
## data:  e1e2
## DM = -0.82989, Forecast horizon = 2, Loss function power = 2, p-value =
## 0.7946
## alternative hypothesis: greater

With alternative = "greater" The alternative hypothesis is that method 2 is more accurate than method 1. Since we do not reject the null (the p-value being of 0.795), we are not led to use the more sophisticated model (ARMA(2,2)) and we keep the simple AR(1) model.

Assume now that we want to compare the AR(1) process to a VAR model (see Def. 3.1). We consider a bivariate VAR, where GDP growth is complemented with CPI-based inflation rate.

Code

library(vars)
infl <- c(NaN,log(data$cpi[2:T]/data$cpi[1:(T-1)]))
y_var <- cbind(y,infl)
e3 <- NULL
for(T.star in first.date:(T-h)){
  estim.model.3 <- VAR(y_var[2:T.star,],p=1)
  e3 <- c(e3,y[T.star+h] - predict(estim.model.3,n.ahead=h)$fcst$y[h,1])
}
res.DM <- dm.test(e1,e3,h = h,alternative = "greater")
res.DM

## 
##  Diebold-Mariano Test
## 
## data:  e1e3
## DM = -0.11555, Forecast horizon = 2, Loss function power = 2, p-value =
## 0.5458
## alternative hypothesis: greater

Again, we do not find that the alternative model (here the VAR(1) model) is better than the AR(1) model to forecast GDP growth.

4.1 Exercise

Code

render_exercises("ps3", "ex-forecast-performance")

Exercise 4.1 (Forecast comparison and the Diebold–Mariano test) PS3 | Advanced | Forecast evaluation
Review Section 4.

Three one-step forecasting models are compared on 200 observations: A, a random walk; B, an AR(1); and C, an ARMA(2,2). The out-of-sample exercise uses 100 rolling estimation windows of length 100.

	Model A	Model B	Model C
In-sample RMSE	1.03	0.96	0.94
Out-of-sample RMSE	1.11	1.02	1.09
Mean loss differential \(\bar d\) versus A	—	-0.191	-0.044
\(\hat\gamma_0\)	—	0.090	0.132
\(\hat\gamma_1\)	—	0.061	0.086
\(\hat\gamma_2\)	—	0.020	0.047
\(\hat\gamma_3\)	—	0.003	0.011

Here \(d_t=e_{j,t}^2-e_{A,t}^2\), for \(j=B,C\).

1. Discuss the in-sample and out-of-sample rankings.

2. Using a Newey–West long-run variance with maximum lag 3 and Bartlett weights, estimate \(\mathbb{V}ar(\bar d)\).

3. Test at 5% whether B and C have equal predictive accuracy to A, and interpret the sign of the test statistic.

4.1.1 Exercise from the 2024 exam

Code

render_exercises("exam-2024", "ex-exam2024-ma1-forecasting")

Exercise 4.2 (Inflation forecasts from an MA(1)) EXAM-2024 | Exam | True/false
Review Section 4.

Let annual inflation be \(\pi_t=0.02+\varepsilon_t+0.4\varepsilon_{t-1}\), where \(\varepsilon_t\overset{i.i.d.}{\sim}\mathcal N(0,\sigma^2)\) and \(\sigma=0.01\).

State whether each assertion is true or false and justify your answer.

1. For \(h>0\), \(\mathbb E(\pi_{t+h}+\cdots+\pi_{t+1}\mid\varepsilon_t,\varepsilon_{t-1},\ldots)=0.02h+0.4\varepsilon_t\).
2. \(\mathbb{V}ar(\pi_t)=0.01\).
3. Conditional on \(\varepsilon_t=0\), a 95% interval for \(\pi_{t+1}\) is \([-1.96\sigma,1.96\sigma]\).
4. \(\pi_t\) is correlated with \(\pi_{t-1}\) but uncorrelated with \(\pi_{t-2}\).
5. \(\mathbb{V}ar_t(\pi_{t+3})=\mathbb{V}ar_t(\pi_{t+2})\).

4.1.2 Exercise from the 2025 exam

Code

render_exercises("exam-2025", "ex-exam2025-ma2-gdp-forecasting")

Exercise 4.3 (GDP forecasts from an MA(2) growth model) EXAM-2025 | Exam | True/false
Review Section 4.

Let \(y_t=\log(GDP_t/GDP_{t-1})\) and suppose \[ y_t=\mu+\varepsilon_t+0.5\varepsilon_{t-1}+0.2\varepsilon_{t-2}, \qquad \varepsilon_t\overset{i.i.d.}{\sim}\mathcal N(0,\sigma^2). \] Recall that \(\log(GDP_{t+h})=\log(GDP_t)+y_{t+h}+\cdots+y_{t+1}\). Conditional moments at date \(t\) use current and past information.

State whether each assertion is true or false and justify your answer.

1. \(\mathbb E_t[\log(GDP_{t+2})]=\log(GDP_t)+2\mu+\varepsilon_{t+2}+1.5\varepsilon_{t+1}+0.7\varepsilon_t\).
2. \(\mathbb{V}ar(y_t)=0.29\sigma^2\).
3. Conditional on date-\(t\) information, a 95% interval for \(y_{t+2}\) is \([\mu-1.96\sigma,\mu+1.96\sigma]\).
4. The first-order autocorrelation of \(y_t\) is \(0.6/1.29\).
5. \(\mathbb E[\log(GDP_{t+h}/GDP_t)]=h\mu\).

References

Box, George.E. P., and Gwilym M. Jenkins. 1976. Time Series Analysis: Forecasting and Control. Holden-Day.

De Gooijer, Jan G., and Rob J. Hyndman. 2006. “25 Years of Time Series Forecasting.” International Journal of Forecasting 22 (3): 443–73. https://doi.org/https://doi.org/10.1016/j.ijforecast.2006.01.001.

Diebold, Francis, and Roberto Mariano. 1995. “Comparing Predictive Accuracy.” Journal of Business & Economic Statistics 13 (3): 253–63. https://EconPapers.repec.org/RePEc:bes:jnlbes:v:13:y:1995:i:3:p:253-63.

Jordà, Òscar, Moritz Schularick, and Alan M. Taylor. 2017. “Macrofinancial History and the New Business Cycle Facts.” NBER Macroeconomics Annual 31 (1): 213–63. https://ideas.repec.org/a/ucp/macann/doi10.1086-690241.html.

---

13. In its inflation report, the Bank of England displays charts showing the conditional distribution of future inflation, called fan charts. This fan charts show the uncertainty associated with future inflation. See this page.↩︎

14. The present reasoning relies on the assumption that the \(\varepsilon_t\)’s are observed. But this is generally not the case in practice, where the \(\varepsilon_t\)’s generally have to be estimated.↩︎

15. This ShinyApp application illustrates the notion of statistical test (illustrating the p-value and the cirtical region, in particular).↩︎