Forecasting with Exogenous Regressors

This notebook provides examples of the accepted data structures for passing the expected value of exogenous variables when these are included in the mean. For example, consider an AR(1) with 2 exogenous variables. The mean dynamics are

\[Y_t = \phi_0 + \phi_1 Y_{t-1} + \beta_0 X_{0,t} + \beta_1 X_{1,t} + \epsilon_t.\]

The \(h\)-step forecast, \(E_{T}[Y_{t+h}]\), depends on the conditional expectation of \(X_{0,T+h}\) and \(X_{1,T+h}\),

\[E_{T}[Y_{T+h}] = \phi_0 + \phi_1 E_{T}[Y_{T+h-1}] + \beta_0 E_{T}[X_{0,T+h}] +\beta_1 E_{T}[X_{1,T+h}]\]

where \(E_{T}[Y_{T+h-1}]\) has been recursively computed.

In order to construct forecasts up to some horizon \(h\), it is necessary to pass \(2\times h\) values (\(h\) for each series). If using the features of forecast that allow many forecast to be specified, it necessary to supply \(n \times 2 \times h\) values.

There are two general purpose data structures that can be used for any number of exogenous variables and any number steps ahead:

  • dict - The values can be pass using a dict where the keys are the variable names and the values are 2-dimensional arrays. This is the most natural generalization of a pandas DataFrame to 3-dimensions.

  • array - The vales can alternatively be passed as a 3-d NumPy array where dimension 0 tracks the regressor index, dimension 1 is the time period and dimension 2 is the horizon.

When a model contains a single exogenous regressor it is possible to use a 2-d array or DataFrame where dim0 tracks the time period where the forecast is generated and dimension 1 tracks the horizon.

In the special case where a model contains a single regressor and the horizon is 1, then a 1-d array of pandas Series can be used.

[1]:
# initial setup
%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn
from arch.__future__ import reindexing

seaborn.set_style("darkgrid")
plt.rc("figure", figsize=(16, 6))
plt.rc("savefig", dpi=90)
plt.rc("font", family="sans-serif")
plt.rc("font", size=14)

Simulating data

Two \(X\) variables are simulated and are assumed to follow independent AR(1) processes. The data is then assumed to follow an ARX(1) with 2 exogenous regressors and GARCH(1,1) errors.

[2]:
from arch.univariate import ARX, GARCH, ZeroMean, arch_model

burn = 250

x_mod = ARX(None, lags=1)
x0 = x_mod.simulate([1, 0.8, 1], nobs=1000 + burn).data
x1 = x_mod.simulate([2.5, 0.5, 1], nobs=1000 + burn).data

resid_mod = ZeroMean(volatility=GARCH())
resids = resid_mod.simulate([0.1, 0.1, 0.8], nobs=1000 + burn).data

phi1 = 0.7
phi0 = 3
y = 10 + resids.copy()
for i in range(1, y.shape[0]):
    y[i] = phi0 + phi1 * y[i - 1] + 2 * x0[i] - 2 * x1[i] + resids[i]

x0 = x0.iloc[-1000:]
x1 = x1.iloc[-1000:]
y = y.iloc[-1000:]
y.index = x0.index = x1.index = np.arange(1000)

Plotting the data

[3]:
ax = pd.DataFrame({"ARX": y}).plot(legend=False)
ax.legend(frameon=False)
_ = ax.set_xlim(0, 999)
../_images/univariate_univariate_forecasting_with_exogenous_variables_5_0.png

Forecasting the X values

The forecasts of \(Y\) depend on forecasts of \(X_0\) and \(X_1\). Both of these follow simple AR(1), and so we can construct the forecasts for all time horizons. Note that the value in position [i,j] is the time-i forecast for horizon j+1.

[4]:
x0_oos = np.empty((1000, 10))
x1_oos = np.empty((1000, 10))
for i in range(10):
    if i == 0:
        last = x0
    else:
        last = x0_oos[:, i - 1]
    x0_oos[:, i] = 1 + 0.8 * last
    if i == 0:
        last = x1
    else:
        last = x1_oos[:, i - 1]
    x1_oos[:, i] = 2.5 + 0.5 * last

x1_oos[-1]
[4]:
array([5.30640496, 5.15320248, 5.07660124, 5.03830062, 5.01915031,
       5.00957515, 5.00478758, 5.00239379, 5.00119689, 5.00059845])

Fitting the model

Next, the most is fit. The parameters are accurately estimated.

[5]:
exog = pd.DataFrame({"x0": x0, "x1": x1})
mod = arch_model(y, x=exog, mean="ARX", lags=1)
res = mod.fit(disp="off")
print(res.summary())
                          AR-X - GARCH Model Results
==============================================================================
Dep. Variable:                   data   R-squared:                       0.991
Mean Model:                      AR-X   Adj. R-squared:                  0.991
Vol Model:                      GARCH   Log-Likelihood:               -1381.17
Distribution:                  Normal   AIC:                           2776.33
Method:            Maximum Likelihood   BIC:                           2810.68
                                        No. Observations:                  999
Date:                Mon, Mar 15 2021   Df Residuals:                      995
Time:                        17:46:38   Df Model:                            4
                               Mean Model
========================================================================
                 coef    std err          t      P>|t|  95.0% Conf. Int.
------------------------------------------------------------------------
Const          3.1333      0.148     21.155  2.508e-99 [  2.843,  3.424]
data[1]        0.7049  3.843e-03    183.404      0.000 [  0.697,  0.712]
x0             2.0062  2.360e-02     85.002      0.000 [  1.960,  2.052]
x1            -2.0337  2.557e-02    -79.541      0.000 [ -2.084, -1.984]
                             Volatility Model
==========================================================================
                 coef    std err          t      P>|t|    95.0% Conf. Int.
--------------------------------------------------------------------------
omega          0.0809  3.650e-02      2.215  2.673e-02 [9.326e-03,  0.152]
alpha[1]       0.0790  2.100e-02      3.761  1.693e-04 [3.782e-02,  0.120]
beta[1]        0.8373  5.260e-02     15.919  4.647e-57   [  0.734,  0.940]
==========================================================================

Covariance estimator: robust

Using a dict

The first approach uses a dict to pass the two variables. The key consideration here is the the keys of the dictionary must exactly match the variable names (x0 and x1 here). The dictionary here contains only the final row of the forecast values since forecast will only make forecasts beginning from the final in-sample observation by default.

Using DataFrame

While these examples make use of NumPy arrays, these can be DataFrames. This allows the index to be used to track the forecast origination point, which can be a helpful device.

[6]:
exog_fcast = {"x0": x0_oos[-1:], "x1": x1_oos[-1:]}
forecasts = res.forecast(horizon=10, x=exog_fcast)
forecasts.mean.T.plot()
[6]:
<AxesSubplot:>
../_images/univariate_univariate_forecasting_with_exogenous_variables_11_1.png

Using an array

An array can alternatively be used. This frees the restriction on matching the variable names although the order must match instead. The forecast values are 2 (variables) by 1 (forecast) by 10 (horizon).

[7]:
exog_fcast = np.array([x0_oos[-1:], x1_oos[-1:]])
print(f"The shape is {exog_fcast.shape}")
array_forecasts = res.forecast(horizon=10, x=exog_fcast)
print(array_forecasts.mean - forecasts.mean)
The shape is (2, 1, 10)
     h.01  h.02  h.03  h.04  h.05  h.06  h.07  h.08  h.09  h.10
999   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0

Producing multiple forecasts

forecast can produce multiple forecasts using the same fit model. Here the model is fit to the first 500 observations and then forecasting for the remaining values are produced. It must be the case that the x values passed for forecast have the same number of rows as the table of forecasts produced.

[8]:
res = mod.fit(disp="off", last_obs=500)
exog_fcast = {"x0": x0_oos[-500:], "x1": x1_oos[-500:]}
multi_forecasts = res.forecast(start=500, horizon=10, x=exog_fcast)
multi_forecasts.mean.tail(10)
[8]:
h.01 h.02 h.03 h.04 h.05 h.06 h.07 h.08 h.09 h.10
990 4.057899 5.479729 6.526760 7.317176 7.924497 8.396850 8.767275 9.059388 9.290611 9.474107
991 8.416610 8.947636 9.169862 9.279946 9.357938 9.432114 9.509247 9.588269 9.666102 9.739891
992 3.074113 3.336329 3.928940 4.652387 5.396820 6.104208 6.746561 7.313185 7.803127 8.220671
993 3.540834 4.224427 5.038955 5.840762 6.568873 7.202884 7.741647 8.192534 8.566082 8.873392
994 6.667926 8.739973 10.187609 11.102274 11.614331 11.844837 11.890638 11.823283 11.692855 11.532874
995 12.839440 15.123608 16.253661 16.601616 16.455163 16.021744 15.443408 14.812952 14.187915 13.601524
996 21.780113 23.260720 23.202158 22.341980 21.111465 19.759851 18.427024 17.186929 16.073929 15.099005
997 22.173220 21.882728 20.887431 19.625010 18.322030 17.090708 15.981171 15.010068 14.176426 13.470503
998 22.506450 21.507597 20.210612 18.855677 17.564663 16.394081 15.364505 14.477150 13.723291 13.089638
999 20.446967 19.172556 18.018315 16.953479 15.978468 15.099800 14.321296 13.642141 13.057505 12.559846

The plot of the final 5 forecast paths shows the the mean reversion of the process.

[9]:
_ = multi_forecasts.mean.tail().T.plot()
../_images/univariate_univariate_forecasting_with_exogenous_variables_17_0.png

The previous example made use of dictionaries where each of the values was a 500 (number of forecasts) by 10 (horizon) array. The alternative format can be used where x is a 3-d array with shape 2 (variables) by 500 (forecasts) by 10 (horizon).

[10]:
exog_fcast = np.array([x0_oos[-500:], x1_oos[-500:]])
print(exog_fcast.shape)
array_multi_forecasts = res.forecast(start=500, horizon=10, x=exog_fcast)
np.max(np.abs(array_multi_forecasts.mean - multi_forecasts.mean))
(2, 500, 10)
[10]:
h.01    0.0
h.02    0.0
h.03    0.0
h.04    0.0
h.05    0.0
h.06    0.0
h.07    0.0
h.08    0.0
h.09    0.0
h.10    0.0
dtype: float64

x input array sizes

While the natural shape of the x data is the number of forecasts, it is also possible to pass an x that has the same shape as the y used to construct the model. The may simplify tracking the origin points of the forecast. Values are are not needed are ignored. In this example, the out-of-sample values are 2 by 1000 (original number of observations) by 10. Only the final 500 are used.

<h3><b>WARNING</b></h3>
Other sizes are <b>not</b> allowed. The size of the out-of-sample data must either match the original data size or the number of forecasts.
[11]:
exog_fcast = np.array([x0_oos, x1_oos])
print(exog_fcast.shape)
array_multi_forecasts = res.forecast(start=500, horizon=10, x=exog_fcast)
np.max(np.abs(array_multi_forecasts.mean - multi_forecasts.mean))
(2, 1000, 10)
[11]:
h.01    0.0
h.02    0.0
h.03    0.0
h.04    0.0
h.05    0.0
h.06    0.0
h.07    0.0
h.08    0.0
h.09    0.0
h.10    0.0
dtype: float64

Special Cases with a single x variable

When a model consists of a single exogenous regressor, then x can be a 1-d or 2-d array (or Series or DataFrame).

[12]:
mod = arch_model(y, x=exog.iloc[:, :1], mean="ARX", lags=1)
res = mod.fit(disp="off")
print(res.summary())
                          AR-X - GARCH Model Results
==============================================================================
Dep. Variable:                   data   R-squared:                       0.939
Mean Model:                      AR-X   Adj. R-squared:                  0.939
Vol Model:                      GARCH   Log-Likelihood:               -2329.43
Distribution:                  Normal   AIC:                           4670.87
Method:            Maximum Likelihood   BIC:                           4700.31
                                        No. Observations:                  999
Date:                Mon, Mar 15 2021   Df Residuals:                      996
Time:                        17:46:39   Df Model:                            3
                               Mean Model
========================================================================
                 coef    std err          t      P>|t|  95.0% Conf. Int.
------------------------------------------------------------------------
Const         -6.3855      0.275    -23.199 4.709e-119 [ -6.925, -5.846]
data[1]        0.7787  9.786e-03     79.571      0.000 [  0.759,  0.798]
x0             1.7489  6.309e-02     27.719 4.126e-169 [  1.625,  1.873]
                             Volatility Model
==========================================================================
                 coef    std err          t      P>|t|    95.0% Conf. Int.
--------------------------------------------------------------------------
omega          5.4298      1.319      4.117  3.841e-05   [  2.845,  8.015]
alpha[1]       0.1375  5.762e-02      2.387  1.699e-02 [2.460e-02,  0.250]
beta[1]        0.0000      0.248      0.000      1.000   [ -0.486,  0.486]
==========================================================================

Covariance estimator: robust

These two examples show that both formats can be used.

[13]:
forecast_1d = res.forecast(horizon=10, x=x0_oos[-1])
forecast_2d = res.forecast(horizon=10, x=x0_oos[-1:])
print(forecast_1d.mean - forecast_2d.mean)

## Simulation-forecasting

mod = arch_model(y, x=exog, mean="ARX", lags=1, power=1.0)
res = mod.fit(disp="off")
     h.01  h.02  h.03  h.04  h.05  h.06  h.07  h.08  h.09  h.10
999   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0

Simulation

forecast supports simulating paths. When forecasting a model with exogenous variables, the same value is used to in all mean paths. If you wish to also simulate the paths of the x variables, these need to generated and then passed inside a loop.

Static out-of-sample x

This first example shows that variance of the paths when the same x values are used in the forecast. There is a sense the out-of-sample x are treated as deterministic.

[14]:
x = {"x0": x0_oos[-1], "x1": x1_oos[-1]}
sim_fixedx = res.forecast(horizon=10, x=x, method="simulation", simulations=100)
sim_fixedx.simulations.values.std(1)
[14]:
array([[0.91938951, 1.17924331, 1.2191893 , 1.15076098, 1.37975046,
        1.49420146, 1.4678764 , 1.40853659, 1.53294681, 1.41081542]])

Simulating the out-of-sample x

This example simulates distinct paths for the two exogenous variables and then simulates a single path. This is then repeated 100 times. We see that variance is much higher when we account for variation in the x data.

[15]:
from numpy.random import RandomState


def sim_ar1(params: np.ndarray, initial: float, horizon: int, rng: RandomState):
    out = np.zeros(horizon)
    shocks = rng.standard_normal(horizon)
    out[0] = params[0] + params[1] * initial + shocks[0]
    for i in range(1, horizon):
        out[i] = params[0] + params[1] * out[i - 1] + shocks[i]
    return out


simulations = []
rng = RandomState(20210301)
for i in range(100):
    x0_sim = sim_ar1(np.array([1, 0.8]), x0.iloc[-1], 10, rng)
    x1_sim = sim_ar1(np.array([2.5, 0.5]), x1.iloc[-1], 10, rng)
    x = {"x0": x0_sim, "x1": x1_sim}
    fcast = res.forecast(horizon=10, x=x, method="simulation", simulations=1)
    simulations.append(fcast.simulations.values)

Finally the standard deviation is quite a bit larger. This is a most accurate value fo the long-run variance of the forecast residuals which should account for dynamics in the model and any exogenous regressors.

[16]:
joined = np.concatenate(simulations, 1)
joined.std(1)
[16]:
array([[3.043168  , 4.72751552, 6.22303597, 7.59272654, 8.43144415,
        8.69788172, 8.67148595, 8.75879958, 9.05931905, 9.27729954]])