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 or pandas Series can be used.

[1]:
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn

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([3.86931675, 4.43465838, 4.71732919, 4.85866459, 4.9293323 ,
       4.96466615, 4.98233307, 4.99116654, 4.99558327, 4.99779163])

Fitting the model

Next, the model is fit. The parameters are precisely 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.989
Mean Model:                      AR-X   Adj. R-squared:                  0.989
Vol Model:                      GARCH   Log-Likelihood:               -1476.10
Distribution:                  Normal   AIC:                           2966.20
Method:            Maximum Likelihood   BIC:                           3000.55
                                        No. Observations:                  999
Date:                Tue, Oct 07 2025   Df Residuals:                      995
Time:                        12:45:45   Df Model:                            4
                               Mean Model
========================================================================
                 coef    std err          t      P>|t|  95.0% Conf. Int.
------------------------------------------------------------------------
Const          2.7503      0.169     16.311  8.274e-60 [  2.420,  3.081]
data[1]        0.6991  3.628e-03    192.716      0.000 [  0.692,  0.706]
x0             2.0189  2.057e-02     98.146      0.000 [  1.979,  2.059]
x1            -1.9620  2.802e-02    -70.026      0.000 [ -2.017, -1.907]
                             Volatility Model
==========================================================================
                 coef    std err          t      P>|t|    95.0% Conf. Int.
--------------------------------------------------------------------------
omega          0.1065  2.923e-02      3.643  2.693e-04 [4.920e-02,  0.164]
alpha[1]       0.1835  3.261e-02      5.628  1.821e-08   [  0.120,  0.247]
beta[1]        0.7390  4.148e-02     17.815  5.387e-71   [  0.658,  0.820]
==========================================================================

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]:
<Axes: >
../_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 2.811408 5.634993 7.603617 8.909734 9.734279 10.223992 10.489269 10.609375 10.639334 10.616290
991 15.895393 18.088062 18.654597 18.361222 17.647005 16.762042 15.845959 14.973362 14.180623 13.481866
992 14.062677 14.469561 14.450813 14.181869 13.782869 13.332214 12.878483 12.449795 12.060714 11.717157
993 13.781067 13.441605 12.844247 12.238383 11.713973 11.291455 10.963674 10.714716 10.527849 10.388444
994 14.938496 14.307436 13.466112 12.672153 12.009441 11.488372 11.091967 10.796118 10.577759 10.417571
995 10.111196 9.386117 9.057210 8.949350 8.961330 9.034393 9.134552 9.242466 9.347507 9.444244
996 11.209626 10.952390 10.981166 11.051404 11.082498 11.060972 10.997019 10.906132 10.802192 10.695520
997 12.066999 12.767649 13.086712 13.124136 12.977669 12.725378 12.423376 12.109021 11.805405 11.525493
998 16.291368 16.700134 16.337664 15.658149 14.888627 14.138408 13.455887 12.858239 12.347035 11.916476
999 18.996030 18.969988 18.069960 16.907205 15.753456 14.717169 13.830444 13.091483 12.485015 11.991839

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]:
np.float64(0.0)

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.

WARNING

Other sizes are not 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]:
np.float64(0.0)

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.949
Mean Model:                      AR-X   Adj. R-squared:                  0.949
Vol Model:                      GARCH   Log-Likelihood:               -2279.80
Distribution:                  Normal   AIC:                           4571.61
Method:            Maximum Likelihood   BIC:                           4601.05
                                        No. Observations:                  999
Date:                Tue, Oct 07 2025   Df Residuals:                      996
Time:                        12:45:46   Df Model:                            3
                               Mean Model
========================================================================
                 coef    std err          t      P>|t|  95.0% Conf. Int.
------------------------------------------------------------------------
Const         -6.5634      0.256    -25.629 7.266e-145 [ -7.065, -6.061]
data[1]        0.7693  9.555e-03     80.514      0.000 [  0.751,  0.788]
x0             1.7509  5.490e-02     31.895 3.118e-223 [  1.643,  1.859]
                             Volatility Model
==========================================================================
                 coef    std err          t      P>|t|    95.0% Conf. Int.
--------------------------------------------------------------------------
omega          2.9611      0.590      5.023  5.088e-07   [  1.806,  4.117]
alpha[1]       0.1289  4.477e-02      2.879  3.994e-03 [4.113e-02,  0.217]
beta[1]        0.3538      0.105      3.369  7.531e-04   [  0.148,  0.560]
==========================================================================

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([[2.04083966, 2.35116082, 2.71759643, 2.51314901, 2.57913983,
        2.27996525, 2.04258321, 1.94847879, 1.79868043, 1.93939794]])

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.41738901, 5.35826138, 6.62153769, 7.83981849, 8.59893642,
        8.93089675, 9.32331466, 9.08742099, 9.45813794, 9.72138886]])

Conditional Mean Alignment vs. Forecast Alignment

When fitting a model with exogenous variables, the data are aligned so that the values in x[j] are used to compute the conditional mean of y[j]. For example, in the case of an AR(1)-X, the model is

\[Y_t = \phi_0 + \phi_1 Y_{t-1} + \beta X_{t} + \epsilon_t.\]

We can recover the conditional mean by subtracting the residuals from the original data. When we do this we see that the conditional mean of observation 0 is missing since we need one lag of \(Y\) to fit the model.

[17]:
mod = arch_model(y, x=exog, mean="ARX", lags=1)
res = mod.fit(disp="off")
y - res.resid
[17]:
0            NaN
1       2.709181
2       8.108807
3      11.568810
4      14.927899
         ...
995    11.955259
996     9.434079
997    12.959983
998    15.739526
999    19.897383
Length: 1000, dtype: float64

Conditional Mean uses target alignment

When modeling the conditional mean in an AR-X, HAR-X, or LS model, the \(X\) data is target-aligned. This requires that when modeling the mean of y[t], the correct values of \(X\) must appear in x[t]. Mathematically, the \(X\) matrix used when estimating a model should have the structure (using the Python indexing convention of a T-element data set having indices 0, 1, …, T-1):

\[\begin{split}\left[\begin{array}{c} X_{0}\\ X_{1}\\ \vdots\\ X_{T-1} \end{array}\right]\end{split}\]

forecast uses origin alignment

Forecasting with \(X\) values aligns them differently. When producing a 1-step-ahead forecast for \(Y_{t+1}\) using information available at time \(t\), the \(X\) values used for this forecast must appear in for t. This is needed since when once wants to produce true out-of-sample forecasts (see below), it must be the case that the final row of x passed forecast must all occur after the final time stamp of the most recent \(Y\) value. Mathematically, the \(X\) matrix used in forecasting should have the following structure (using Python indexing convention so that a \(T\) observation dataset will have indices 0, 1, …, T-1).

\[\begin{split}\left[\begin{array}{cccc} E\left[X_{1}|\mathcal{F}_0\right] & E\left[X_{2}|\mathcal{F}_0\right] & \ldots & E\left[X_{h|\mathcal{F}_0}\right]\\ E\left[X_{2}|\mathcal{F}_1\right] & E\left[X_{3}|\mathcal{F}_0\right] & \ldots & E\left[X_{h+1}|\mathcal{F}_1\right]\\ \vdots & \vdots & \vdots & \vdots\\ E\left[X_{T}|\mathcal{F}_{T-1}\right] & E\left[X_{T+1}|\mathcal{F}_{T-1}\right] & \ldots & E\left[X_{T+h-1}|\mathcal{F}_{T-1}\right] \end{array}\right]\end{split}\]

where \(|\mathcal{F}_{s}\) is the time-\(s\) information set.

If you use the same x value in the model when forecasting, you will see different values due to this alignment difference. Naively using the same x values ie equivalent to setting

\[E\left[X_{s}|\mathcal{F}_{s-1} \right] = X_{s-1}\]

In general this would not be correct when forecasting, and will always produce forecasts that differ from the conditional mean. In order to recover the conditional mean using the forecast function, it is necessary to shift the \(X\) values by -1, so that once shifted, the x values will have the relationship

\[E\left[X_{s}|\mathcal{F}_{s-1} \right] = X_{s} .\]

Here we shift the \(X\) data by -1 so x[s] is treated as being in the information set for y[s-1]. Also, note that the final forecast is NaN. Conceptually this must be the case because the value of \(X\) at 999 should be ahead of 999 (i.e., at observation 1,000), and we do not have this value.

[18]:
exog_dict = {col: exog[[col]].shift(-1) for col in exog}
fcast = res.forecast(horizon=1, x=exog_dict, start=0)
fcast.mean
[18]:
h.1
0 2.709181
1 8.108807
2 11.568810
3 14.927899
4 18.308219
... ...
995 9.434079
996 12.959983
997 15.739526
998 19.897383
999 NaN

1000 rows × 1 columns

(Nearly) out-of-sample forecasts

These “in-sample” forecasts are not really forecasts at all but are just fitted values with a different alignment. If you want real (nearly) out-of-sample forecasts\(\dagger\), it is necessary to replace the actual values of \(X\) with their conditional expectation. This can be done by taking the fitted values from AR(1) models of the \(X\) variables.

\(\dagger\) These are not true out-of-sample since the parameters were estimated using data from that same range of indices where these forecasts target. True out-of-sample requires both using forecast \(X\) values and parameters estimated without the period being forecasted.

[19]:
res0 = ARX(exog["x0"], lags=1).fit()
res1 = ARX(exog["x1"], lags=1).fit()
forecast_x = pd.concat(
    [res0.forecast(start=0).mean, res1.forecast(start=0).mean], axis=1
)
forecast_x.columns = ["x0f", "x1f"]
in_samp_forcast_exog = {"x0": forecast_x[["x0f"]], "x1": forecast_x[["x1f"]].shift(-1)}
fcast = res.forecast(horizon=1, x=in_samp_forcast_exog, start=0)
fcast.mean
[19]:
h.1
0 1.930577
1 7.941269
2 9.155056
3 13.509107
4 18.158888
... ...
995 8.826984
996 12.112520
997 13.941321
998 17.352874
999 NaN

1000 rows × 1 columns

True out-of-sample forecasts

In order to make a true out-of-sample prediction, we need the expected values of X from the end of the data we have. These can be constructed by forecasting the two \(X\) variables and then passing these values as x to forecast.

[20]:
mod = arch_model(y, x=exog, mean="ARX", lags=1)
res = mod.fit(disp="off")

actual_x_oos = {
    "x0": res0.forecast(horizon=10).mean,
    "x1": res1.forecast(horizon=10).mean,
}
fcasts = res.forecast(horizon=10, x=actual_x_oos)
fcasts.mean
[20]:
h.01 h.02 h.03 h.04 h.05 h.06 h.07 h.08 h.09 h.10
999 19.016205 19.044782 18.209077 17.106733 16.002979 15.004062 14.142287 13.417293 12.815734 12.32041