"The only function of economic forecasting is to make astrology look respectable."
A real forecaster may not like this quotation because it expresses that forecasting is some sort of dubious mystical activity. In truth, most forecasts are the result of modern, replicable and largely quantitative methods. Among others, one key field of forecasting is economics. Economic agents such as governments, businesses, central banks and financial institutions make forecasts of the major economic variables like gross domestic product (GDP), unemployment, consumption, inflation and interest rates.
The aim of this post is to estimate different models to forecast Swiss GDP growth after an appreciation of 15% of the CHF vis-à-vis the EUR as it was the case when the SNB repealed the minimum rate. The models predict considerably lower GDP growth which confirmed that the models reflect the negative relationship between GDP growth and the exchange rate. In fact an ARMAX- and a VARMAX model are elaborated. The results show that both types of models exhibit quite reasonable predictive power.
This estimation is done with the SARIMAX and the VARMAX routines of Statsmodels. First of all we import some packages and load the data, which can be found on the database of the St. Louis Fred and the OECD. In fact the VARMAX model contains besides Swiss GDP growth also inflation, unemployment rate, exports growth, a shorter and a longer yield spread. The exogenous data contains lagged GDP of the Euro Area and lagged CHF/EUR exchangerate. Conducting Granger Causality tests shows that from a statistically point of view its justified to treat the exchange rate as exogenous which is an assumption that is often used when forecasting GDP growth in a small open economy, although it does not make a lot of sense economically. Economically it can reasonably be assumed that Euro Area GDP growth is exogenous to Swiss GDP growth because Switzerland is a small open economy and thus is unlikely to affect GDP growth in the EU.
In [1]:
import numpy as np
import pandas as pd
import statsmodels.api as sm
from pandas.io.data import DataReader
import matplotlib.pyplot as plt
import warnings
warnings.simplefilter('ignore')
In [2]:
table = pd.read_excel('forecasting_data.xlsx')
idx = table.observation_date
data = table[['SwissGDP','Inflation','Unemployment','Exports','Shortspread','Longspread']]
data.index=idx
xdata = table[['EUR/CHF','EUR/CHF_L1','EUR/CHF_L2','EUR/CHF_L3','EuroGDP_L1']]
xdata.index=idx
SwissGDP = data[['SwissGDP']]
fig1 = plt.figure(figsize=(8,5))
ax1 = fig1.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(SwissGDP, lags=6, ax=ax1)
ax2 = fig1.add_subplot(212)
fig1 = sm.graphics.tsa.plot_pacf(SwissGDP, lags=6, ax=ax2)
plt.show(fig1)
The figure above shows the ACF and the PACF of Swiss GDP growth. The ACF is significantly different from zero up to two lags while the PACF is only significantly different from zero for lag one. Therefore, p = 4 and q = 4 should be according to Box and Jenkins (1976) an appropriate maximum order for the ARMAX(p,q) model. In an automated procedure several information criteria are calculated to trade-off explanatory power with parsimony.
In [3]:
aic = pd.DataFrame(np.zeros((5,5), dtype=float))
bic = pd.DataFrame(np.zeros((5,5), dtype=float))
# Iterate over all ARMA(p,q) models with p,q in [5,5]
for p in range(5):
for q in range(5):
if p == 0 and q == 0:
continue
# Estimate the model with no missing datapoints
mod = sm.tsa.statespace.SARIMAX(SwissGDP, xdata, order=(p,0,q), enforce_invertibility=True)
try:
res = mod.fit(maxiter=5000, disp = False)
aic.iloc[p,q] = res.aic
bic.iloc[p,q] = res.bic
except:
aic.iloc[p,q] = np.nan
bic.iloc[p,q] = np.nan
aic.iloc[0,0] = np.nan
bic.iloc[0,0] = np.nan
print(aic)
print(bic)
q = aic.min().idxmin()
p = aic.idxmin()[q]
As can be seen above, the AIC and the BIC are both minimized for $p = 1$ and $q = 2$, therefore an $ARMAX(1,2)$ model is used here. Before estimating the model, some notes about the general $ARMA$ model.
Due to Wold (1938), a covariance-stationary equation can be written as the so-called "Wold-Decomposition": $$\begin{equation*} Y_t = \mu + \sum_{j=0}^{\infty} \psi_j \epsilon_{t-j} \end{equation*}$$ where $\epsilon_t$ is the error from forecasting $Y_t$ solely with lagged values of $Y_t$. A Wold representation in principle requires fitting an infinite number of parameters ($\psi_1, \psi_2, ...$). With a finite number of observations this will never be possible. However, in practice a typical assumption is that $Y_t$ is a causal process and $\psi(L)$ can be approximated from its causal representation: $$\begin{equation*} \sum_{j=0}^\infty \psi_j L^j = \frac{\Theta(L)}{\Phi(L)} \equiv \frac{(1 + \theta_1L + \theta_2L^2 + ... + \theta_qL^q)}{(1 - \phi_1L - \phi_2L^2 - ... - \phi_pL^p)} \end{equation*}$$
ARMA modelling seeks a finite parameter approximation to the Wold representation. Box and Jenkins (1976) emphasise parsimony in the choice of the AR order $p$ and the MA $q$ order
The general ARMA model is: $$\begin{equation*} Y_t = \sum_{i=0}^p \phi_i Y_{t-i} + \epsilon_t + \sum_{j=0}^q \theta_j \epsilon_{t-j} \end{equation*}$$ where $\epsilon_t$ is a disturbance term(white noise) and $Y_t$ is the stationary forecast variable. The typical strategy implies:
In [4]:
# Statespace
mod = sm.tsa.statespace.SARIMAX(SwissGDP, xdata, order=(p,0,q))
res = mod.fit(maxiter=5000)
fig2 = res.plot_diagnostics()
print(res.summary())
plt.show(fig2)
The residual plot suggest that the ARMA model fits pretty well. The fitted values track the actual values and even cover the recession during the financial crisis. The residuals look reasonably independent and identically distributed (i.i.d.) according to the ACF and also resemble a normal distribution. By looking at the ACF of the residuals it’s observable that there are no lags which are autocorrelated.
Next, we generate the exogenous variables for the forecasting task. For simplicity, it is assumed that the exchange rate decreases by 15% in 2015Q1 (due to the decision of the SNB to repeal the euro minimum rate on January 15) and is zero otherwise. The lagged Euro Area GDP is set to its mean. Then the model is used to forecast 12 quaters ahead.
In [5]:
# Generate Exogenous Data for Forecasting
fx = pd.DataFrame(index=list(range(0, 12)) , columns=xdata.columns)
fx = fx.fillna(0)
fx.iloc[0,0] = -15.00
fx.iloc[1,1] = -15.00
fx.iloc[2,2] = -15.00
fx.iloc[3,3] = -15.00
fx.iloc[0:12,4] = np.mean(xdata[['EuroGDP_L1']]).values
In [6]:
# In-sample one-step-ahead predictions, and out-of-sample forecasts
nforecast = 11
predict = res.get_prediction(end=mod.nobs + nforecast, exog = fx)
""" idx = np.arange(len(predict.predicted_mean)) """
predict_ci90 = predict.conf_int(alpha=0.1)
predict_ci60 = predict.conf_int(alpha=0.4)
predict_ci30 = predict.conf_int(alpha=0.7)
fidx = predict_ci90.index
predict_ci90.iloc[138,:] = SwissGDP.iloc[138,:].values
predict_ci60.iloc[138,:] = SwissGDP.iloc[138,:].values
predict_ci30.iloc[138,:] = SwissGDP.iloc[138,:].values
# Graph
fig3, ax = plt.subplots(figsize=(12,6))
ax.grid()
ax.plot(idx,SwissGDP,'k',label='', alpha=0.7, linewidth=0.75)
ax.plot(idx,SwissGDP,'k.', label='Actual Values', alpha=0.7, linewidth=0.75)
# Plot
ax.plot( predict.predicted_mean[:-nforecast-1], 'gray',label='In Sample Prediction')
ax.plot(fidx[-nforecast-1:], predict.predicted_mean[-nforecast-1:], 'k--', linestyle='--', linewidth=1.5,label='Out of Sample Prediction')
ax.fill_between(fidx[-nforecast-2:], predict_ci90.iloc[-nforecast-2:, 0], predict_ci90.iloc[-nforecast-2:, 1], alpha=0.15, label= '')
ax.fill_between(fidx[-nforecast-2:], predict_ci60.iloc[-nforecast-2:, 0], predict_ci60.iloc[-nforecast-2:, 1], alpha=0.15, label= '')
ax.fill_between(fidx[-nforecast-2:], predict_ci30.iloc[-nforecast-2:, 0], predict_ci30.iloc[-nforecast-2:, 1], alpha=0.15, label= '')
# Retrieve and also plot the NBER recession indicators
rec = DataReader('USREC', 'fred', start=idx[0], end=idx[138])
recr = rec.iloc[0:413:3,0]
ylim = ax.get_ylim()
ax.fill_between(recr.index, ylim[0], ylim[1], recr, facecolor='k', alpha=0.1);
ax.legend()
ax.set_xlim(730000.0, 736603.0)
ax.set(title='ARMAX - Forecasting Swiss GDP');
plt.show(fig3)
We see this scenario predicts Swiss GDP to be considerably lower because of the appreciation. This reflects the positive relationship between a depreciation and GDP growth. Therefore, an appreciation of the swiss franc vis-à-vis the euro leads to lower Swiss GDP growth.
Next, we do the same, including some other variables. That is inflation, unemployment rate, exports growth, a shorter and a longer yield spread. As the automated procedure suggests to take a $VARMAX(2,0)$ model, some notes about VAR models:
A $VAR(p)$ process in it's basic form is: $$\begin{equation*} y_t = c + \Phi_1y_{t-1} + ... + \Phi_py_{t-p} + \epsilon_t \end{equation*}$$ Here $y_t$ represents a set of variables collected in a vector, c denotes a vector of constants, $\Phi_j$ a matrix of autoregressive coefficients and $\epsilon_t$ is white noise. Since the parameters of $\Phi_j$ are unknown, we have to estimate these parameters. To do so, the following is defined: $$\begin{align} Y &:= (y_1, y_2, ... , y_T) \notag \\ \notag B &:= (c, \Phi_1, \Phi_2, ... , \Phi_p) \\ \notag Z &:= \left[ \begin{array}{cccc} 1 & 1 & ... & 1\\ y_{0} & y_1 & ... & y_{T-1} \\ \vdots & \vdots & \ddots & \vdots\\ y_{-p+1} & y_{-p+2} & ... & y_{T-p} \end{array} \right] \\ \notag U &:= (\epsilon_1, \epsilon_2, ... , \epsilon_T) \\ \notag \end{align}$$ The VAR(p) can now be written compactly as $$\begin{equation*} Y = BZ + U \end{equation*}$$
Coefficients can now easily be estimated via OLS or Maximum Likelihood.
In [7]:
# Iterate over all VARMAX(p,q) models with p,q in [5,5]
for p in range(5):
for q in range(5):
if p == 0 and q == 0:
continue
# Estimate the model with no missing datapoints
mod = sm.tsa.VARMAX(data, xdata, order=(p,q), enforce_invertibility=True)
try:
res = mod.fit(maxiter=100, disp = False)
aic.iloc[p,q] = res.aic
bic.iloc[p,q] = res.bic
except:
aic.iloc[p,q] = np.nan
bic.iloc[p,q] = np.nan
aic.iloc[0,0] = np.nan
bic.iloc[0,0] = np.nan
print(aic)
print(bic)
q = aic.min().idxmin()
p = aic.idxmin()[q]
As already mentioned above, both information criteria are minimized within the $VARMAX(2,0)$ model.
In [8]:
# Statespace
mod = sm.tsa.VARMAX(data, xdata, order=(p,q), enforce_invertibility=True)
res = mod.fit(maxiter=5000, disp = False)
fig4 = res.plot_diagnostics(0)
print(res.summary())
plt.show(fig4)
The residual diagnostics plot looks quite similar to the one of the ARMAX model, therefore we conclude same things as above.
As the ARMAX model, the VARMAX model predicts the GDP to be considerably lower.
In [9]:
# In-sample one-step-ahead predictions, and out-of-sample forecasts
nforecast = 11
predict = res.get_prediction(end=mod.nobs + nforecast, exog = fx)
""" idx = np.arange(len(predict.predicted_mean)) """
predict_ci90 = predict.conf_int(alpha=0.1)
predict_ci60 = predict.conf_int(alpha=0.4)
predict_ci30 = predict.conf_int(alpha=0.7)
fidx = predict_ci90.index
predict_ci90.iloc[138,:] = SwissGDP.iloc[138,:].values
predict_ci60.iloc[138,:] = SwissGDP.iloc[138,:].values
predict_ci30.iloc[138,:] = SwissGDP.iloc[138,:].values
# Graph
fig5, ax = plt.subplots(figsize=(12,6))
ax.grid()
ax.plot(idx,SwissGDP,'k',label='', alpha=0.7, linewidth=0.75)
ax.plot(idx,SwissGDP,'k.', label='Actual Values', alpha=0.7, linewidth=0.75)
# Plot
ax.plot( predict.predicted_mean[:-nforecast-1][['SwissGDP']], 'gray',label='In Sample Prediction')
ax.plot(fidx[-nforecast-1:], predict.predicted_mean[-nforecast-1:][['SwissGDP']], 'k--', linestyle='--', linewidth=1.5,label='Out of Sample Prediction')
ax.fill_between(fidx[-nforecast-2:], predict_ci90.iloc[-nforecast-2:, 0], predict_ci90.iloc[-nforecast-2:, 6], alpha=0.15, label= '')
ax.fill_between(fidx[-nforecast-2:], predict_ci60.iloc[-nforecast-2:, 0], predict_ci60.iloc[-nforecast-2:, 6], alpha=0.15, label= '')
ax.fill_between(fidx[-nforecast-2:], predict_ci30.iloc[-nforecast-2:, 0], predict_ci30.iloc[-nforecast-2:, 6], alpha=0.15, label= '')
# Retrieve and also plot the NBER recession indicators
rec = DataReader('USREC', 'fred', start=idx[0], end=idx[138])
recr = rec.iloc[0:413:3,0]
ylim = ax.get_ylim()
ax.fill_between(recr.index, ylim[0], ylim[1], recr, facecolor='k', alpha=0.1);
ax.legend()
ax.set_xlim(730000.0, 736603.0)
ax.set(title='VARMAX - Forecasting Swiss GDP');
plt.show(fig5)