Alternative ARIMA Model Representations

Mathematical Development of regARIMA to ARIMAX Model Conversion

ARIMAX models and regression models with ARIMA errors are closely related, and the choice of which to use is generally dictated by your goals for the analysis. If your objective is to fit a parsimonious model to data and forecast responses, then there is very little difference between the two models.

If you are more interested in preserving the usual interpretation of a regression coefficient as a measure of sensitivity, i.e., the effect of a unit change in a predictor variable on the response, then use a regression model with ARIMA errors. Regression coefficients in ARIMAX models do not possess that interpretation because of the dynamic dependence on the response [1].

Suppose that you have the parameter estimates from a regression model with ARIMA errors, and you want to see how the model structure compares to ARIMAX model. Or, suppose you want some insight as to the underlying relationship between the two models.

The ARIMAX model is (t = 1,...,T):

$${\rm H}(L)y_t=c+X_t\beta +{\rm N}(L){\epsilon }_t,$$

(1)

where

• y_t is the univariate response series.

• X_t is row t of X, which is the matrix of concatenated predictor series. That is, X_t is observation t of each predictor series.

• β is the regression coefficient.

• c is the regression model intercept.

• $${\rm H}(L)=\varphi (L){(1-L)}^D\Phi (L)(1-L^s)=1-{\eta }_{1}L-{\eta }_2{L}^2-\mathrm{...}-{\eta }_P{L}^P,$$ which is the degree P lag operator polynomial that captures the combined effect of the seasonal and nonseasonal autoregressive polynomials, and the seasonal and nonseasonal integration polynomials. For more details on notation, see What Are Multiplicative ARIMA Models?.

• $${\rm N}(L)=\theta (L)\Theta (L)=1+{\nu }_{1}L+{\nu }_2{L}^2+\mathrm{...}+{\nu }_Q{L}^Q,$$ which is the degree Q lag operator polynomial that captures the combined effect of the seasonal and nonseasonal moving average polynomials.

• ε_t is a white noise innovation process.

The regression model with ARIMA errors is (t = 1,...,T)

$$\begin{array}{l}{y}_t=c+X_t\beta +u_t\\ A(L)u_t=B(L){\epsilon }_t,\end{array}$$

(2)

where

• u_t is the unconditional disturbances process.

• $$A(L)=\varphi (L){(1-L)}^D\Phi (L)(1-L^s)=1-a_{1}L-a_2{L}^2-\mathrm{...}-a_P{L}^P,$$ which is the degree P lag operator polynomial that captures the combined effect of the seasonal and nonseasonal autoregressive polynomials, and the seasonal and nonseasonal integration polynomials.

• $$B(L)=\theta (L)\Theta (L)=1+b_{1}L+b_2{L}^2+\mathrm{...}+b_Q{L}^Q,$$ which is the degree Q lag operator polynomial that captures the combined effect of the seasonal and nonseasonal moving average polynomials.

The values of the variables defined in Equation 2 are not necessarily equivalent to the values of the variables in Equation 1, even though the notation might be similar.

Consider Equation 2, the regression model with ARIMA errors. Use the following operations to convert the regression model with ARIMA errors to its corresponding ARIMAX model.

1. Solve for u_t.

   $$\begin{array}{l}{y}_t=c+X_t\beta +u_t\\ u_t=\frac{B(L)}{A(L)}{\epsilon }_t.\end{array}$$

2. Substitute u_t into the regression equation.

   $$\begin{array}{c}{y}_t=c+X_t\beta +\frac{B(L)}{A(L)}{\epsilon }_t\\ A(L)y_t=A(L)c+A(L)X_t\beta +B(L){\epsilon }_t.\end{array}$$

3. Solve for y_t.

   $$\begin{array}{c}{y}_t=A(L)c+A(L)X_t\beta +{\displaystyle \sum _{k=1}^P{a}_k}{y}_{t-k}+B(L){\epsilon }_t\\ =A(L)c+Z_t\Gamma +{\displaystyle \sum _{k=1}^P{a}_k}{y}_{t-k}+B(L){\epsilon }_t.\end{array}$$

   (3)

   In Equation 3,

   • A(L)c = (1 – a₁ – a₂ –...– a_P)c. That is, the constant in the ARIMAX model is the intercept in the regression model with ARIMA errors with a nonlinear constraint. Though applications, such as simulate, handle this constraint, estimate cannot incorporate such a constraint. In the latter case, the models are equivalent when you fix the intercept and constant to 0.

   • In the term A(L)X_tβ, the lag operator polynomial A(L) filters the T-by-1 vector X_tβ, which is the linear combination of the predictors weighted by the regression coefficients. This filtering process requires P presample observations of the predictor series.

   • arima constructs the matrix Z_t as follows:

     • Each column of Z_t corresponds to each term in A(L).

     • The first column of Z_t is the vector X_tβ.

     • The second column of Z_t is a sequence of d₂ NaNs (d₂ is the degree of the second term in A(L)), followed by the product $$L^{d_j}{X}_t\beta$$. That is, the software attaches d₂ NaNs at the beginning of the T-by-1 column, attaches X_tβ after the NaNs, but truncates the end of that product by d₂ observations.

     • The jth column of Z_t is a sequence of d_j NaNs (d_j is the degree of the jth term in A(L)), followed by the product $$L^{d_j}{X}_t\beta$$. That is, the software attaches d_j NaNs at the beginning of the T-by-1 column, attaches X_tβ after the NaNs, but truncates the end of that product by d_j observations.

     .

   • Γ = [1 –a₁ –a₂ ... –a_P]'.

     The arima converter removes all zero-valued autoregressive coefficients of the difference equation. Subsequently, the arima converter does not associate zero-valued autoregressive coefficients with columns in Z_t, nor does it include corresponding, zero-valued coefficients in Γ.

4. Rewrite Equation 3,

   $$y_t=(1-{\displaystyle \sum _{k=1}^P{a}_k})c+X_t\beta -{\displaystyle \sum _{k=1}^P{a}_k}{X}_{t-k}\beta +{\displaystyle \sum _{k=1}^P{a}_k{y}_{t-k}}+{\epsilon }_t+{\displaystyle \sum _{k=1}^Q{\epsilon }_{t-k}}.$$

For example, consider the following regression model whose errors are ARMA(2,1):

$$\begin{array}{c}{y}_t=0.2+0.5X_t+u_t\\ \left(1-0.8L+0.4L^2\right)u_t=\left(1+0.3L\right){\epsilon }_t.\end{array}$$

(4)

The equivalent ARMAX model is:

$$\begin{array}{c}{y}_t=0.12+\left(0.5-0.4L+0.2L^2\right)X_t+0.8y_{t-1}-0.4y_{t-2}+(1+0.3L){\epsilon }_t\\ =0.12+Z_t\Gamma +0.8y_{t-1}-0.4y_{t-2}+(1+0.3L){\epsilon }_t,\end{array}$$

or

$$(1-0.8L+0.4L^2)y_t=0.12+Z_t\Gamma +(1+0.3L){\epsilon }_t,$$

where Γ = [1 –0.8 0.4]' and

$$Z_t=0.5\left[\begin{array}{ccc}{x}_1& NaN& NaN\\ x_2& x_1& NaN\\ x_3& x_2& x_1\\ ⋮& ⋮& ⋮\\ x_T& x_{T-1}& x_{T-2}\end{array}\right].$$

This model is not integrated because all of the eigenvalues associated with the AR polynomial are within the unit circle, but the predictors might affect the otherwise stable process. Also, you need presample predictor data going back at least 2 periods to, for example, fit the model to data.

Show Conversion in MATLAB®

Illustrate the conversion in MATLAB® by model simulation and estimation.

Specify the regression model with ARIMA errors in Equation 4.

MdlregARIMA0 = regARIMA('Intercept',0.2,'AR',{0.8 -0.4}, ...
               'MA',0.3,'Beta',[0.3 -0.2],'Variance',0.2);

Generate presample observations and predictor data.

rng(1); % For reproducibility
T = 100;
maxPQ = max(MdlregARIMA0.P,MdlregARIMA0.Q);
numObs  = T + maxPQ;            % Adjust number of observations to account for presample
XregARIMA = randn(numObs,2);    % Simulate predictor data
u0 = randn(maxPQ,1);            % Presample unconditional disturbances u(t)
e0 = randn(maxPQ,1);            % Presample innovations e(t)

Simulate data from the regression model with ARIMA errors MdlregARIMA0.

rng(100) % For consistent seed with later call
[y1,e1,u1] = simulate(MdlregARIMA0,T,'U0',u0, ...
    'E0',e0,'X',XregARIMA);

Convert the regression model with ARIMA errors to an ARIMAX model.

[MdlARIMAX0,XARIMAX] = arima(MdlregARIMA0,'X',XregARIMA);
MdlARIMAX0

MdlARIMAX0 = 
  arima with properties:

     Description: "ARIMAX(2,0,1) Model (Gaussian Distribution)"
      SeriesName: "Y"
    Distribution: Name = "Gaussian"
               P: 2
               D: 0
               Q: 1
        Constant: 0.12
              AR: {0.8 -0.4} at lags [1 2]
             SAR: {}
              MA: {0.3} at lag [1]
             SMA: {}
     Seasonality: 0
            Beta: [1 -0.8 0.4]
        Variance: 0.2

Generate presample responses for the ARIMAX model to ensure consistency with the regression model with ARIMA errors. Simulate data from the ARIMAX model.

y0 = MdlregARIMA0.Intercept + XregARIMA(1:maxPQ,:)*MdlregARIMA0.Beta' + u0;
rng(100) % For consistent seed with earlier call
y2 = simulate(MdlARIMAX0,T,'Y0',y0,'E0',e0,'X',XARIMAX);

figure
plot(y1,'LineWidth',3)
hold on
plot(y2,'r:','LineWidth',2.5)
hold off
title("\bf Simulated Paths")
legend("regARIMA Model","ARIMAX Model",'Location','best')

The simulated paths are equal because the arima converter enforces the nonlinear constraint when it converts the regression model intercept to the ARIMAX model constant.

Fit a regression model with ARIMA errors to the simulated data.

MdlregARIMA0 = regARIMA('ARLags',[1 2],'MALags',1);
EstMdlregARIMA = estimate(MdlregARIMA0,y1,'E0',e0,'U0',u0,'X',XregARIMA);

 
    Regression with ARMA(2,1) Error Model (Gaussian Distribution):
 
                  Value      StandardError    TStatistic      PValue  
                 ________    _____________    __________    __________

    Intercept     0.14074        0.1014         1.3879         0.16518
    AR{1}         0.83061        0.1375         6.0407      1.5349e-09
    AR{2}        -0.45402        0.1164        -3.9007      9.5927e-05
    MA{1}         0.42803       0.15145         2.8262       0.0047109
    Beta(1)       0.29552      0.022938         12.883       5.597e-38
    Beta(2)      -0.17601      0.030607        -5.7506      8.8941e-09
    Variance      0.18231      0.027765         6.5663      5.1569e-11

Fit an ARIMAX model to the simulated data.

MdlARIMAX = arima('ARLags',[1 2],'MALags',1);
EstMdlARIMAX = estimate(MdlARIMAX,y2,'E0',e0,'Y0',...
    y0,'X',XARIMAX);

 
    ARIMAX(2,0,1) Model (Gaussian Distribution):
 
                 Value      StandardError    TStatistic      PValue  
                ________    _____________    __________    __________

    Constant    0.084996      0.064217         1.3236         0.18564
    AR{1}        0.83136       0.13634         6.0975      1.0775e-09
    AR{2}       -0.45599       0.11788        -3.8683       0.0001096
    MA{1}          0.426       0.15753         2.7043       0.0068446
    Beta(1)        1.053       0.13685         7.6949      1.4166e-14
    Beta(2)      -0.6904       0.19262        -3.5843      0.00033796
    Beta(3)      0.45399       0.15352         2.9572       0.0031047
    Variance     0.18112      0.028836          6.281      3.3635e-10

Convert the estimated regression model with ARIMA errors EstMdlregARIMA to an ARIMAX model.

ConvertedMdlARIMAX = arima(EstMdlregARIMA,'X',XregARIMA)

ConvertedMdlARIMAX = 
  arima with properties:

     Description: "ARIMAX(2,0,1) Model (Gaussian Distribution)"
      SeriesName: "Y"
    Distribution: Name = "Gaussian"
               P: 2
               D: 0
               Q: 1
        Constant: 0.087737
              AR: {0.830611 -0.454025} at lags [1 2]
             SAR: {}
              MA: {0.428031} at lag [1]
             SMA: {}
     Seasonality: 0
            Beta: [1 -0.830611 0.454025]
        Variance: 0.182313

The estimated ARIMAX model constant is not equal to the ARIMAX model constant converted from the regression model with ARIMA errors. In other words, EstMdlARIMAX.Constant is 0.084996 and ConvertedMdlARIMAX.Constant = 0.087737. The reason for the discrepancy is estimate does not enforce the nonlinear constraint that the arima converter enforces. As a result, the other estimates are close, but not equal.

References

[1] Hyndman, R. J. (2010, October). "The ARIMAX Model Muddle." Rob J. Hyndman. Retrieved May 4, 2017 from https://robjhyndman.com/hyndsight/arimax/.