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) |

*y*is the univariate response series._{t}*X*is row_{t}*t*of*X*, which is the matrix of concatenated predictor series. That is,*X*is observation_{t}*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 Multiplicative ARIMA Model.$${\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.*ε*is a white noise innovation process._{t}

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) |

*u*is the unconditional disturbances process._{t}$$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.

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}$$

Substitute

*u*into the regression equation._{t}$$\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}$$

Solve for

*y*._{t}

In Equation 3,$$\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)***A*(*L*)*c*= (1 –*a*_{1}–*a*_{2}–...–*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*, the lag operator polynomial_{t}β*A*(*L*) filters the*T*-by-1 vector*X*, which is the linear combination of the predictors weighted by the regression coefficients. This filtering process requires_{t}β*P*presample observations of the predictor series.`arima`

constructs the matrix*Z*as follows:_{t}Each column of

*Z*corresponds to each term in_{t}*A*(*L*).The first column of

*Z*is the vector_{t}*X*._{t}βThe second column of

*Z*is a sequence of_{t}*d*_{2}`NaN`

s (*d*_{2}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*_{2}`NaN`

s at the beginning of the*T*-by-1 column, attaches*X*after the_{t}β`NaN`

s, but truncates the end of that product by*d*_{2}observations.The

*j*th column of*Z*is a sequence of_{t}*d*_{j}`NaN`

s (*d*is the degree of the_{j}*j*th term in*A*(*L*)), followed by the product $${L}^{{d}_{j}}{X}_{t}\beta $$. That is, the software attaches*d*_{j}`NaN`

s at the beginning of the*T*-by-1 column, attaches*X*after the_{t}β`NaN`

s, but truncates the end of that product by*d*observations._{j}

.

*Γ*= [1*–a*_{1}*–a*_{2}...*–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*, nor does it include corresponding, zero-valued coefficients in_{t}*Γ*.

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.5{X}_{t}+{u}_{t}\\ \left(1-0.8L+0.4{L}^{2}\right){u}_{t}=\left(1+0.3L\right){\epsilon}_{t}.\end{array}$$ | (4) |

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

or

$$(1-0.8L+0.4{L}^{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}\\ \vdots & \vdots & \vdots \\ {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.

You can illustrate this further through simulation and estimation.

Specify the regression model with ARIMA errors in Equation 4.

Mdl1 = 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(Mdl1.P,Mdl1.Q); numObs = T + maxPQ;... % Adjust number of observations to account for presample X1 = 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

`Mdl1`

.rng(100) % For reproducibility [y1,e1,u1] = simulate(Mdl1,T,'U0',u0,... 'E0',e0,'X',X1);

Convert

`Mdl1`

to an ARIMAX model.`[Mdl2,X2] = arima(Mdl1,'X',X1); Mdl2`

Mdl2 = arima with properties: Description: "ARIMAX(2,0,1) Model (Gaussian Distribution)" 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

`Mdl1`

. Simulate data from`Mdl2`

.y0 = Mdl1.Intercept + X1(1:maxPQ,:)*Mdl1.Beta' + u0; rng(100) y2 = simulate(Mdl2,T,'Y0',y0,'E0',e0,'X',X2); figure plot(y1,'LineWidth',3) hold on plot(y2,'r:','LineWidth',2.5) hold off title('{\bf Simulated Paths for Both Models}') legend('regARIMA Model','ARIMAX Model','Location','Best')

The simulated paths are equivalent 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.

ToEstMdl1 = regARIMA('ARLags',[1 2],'MALags',1); EstMdl1 = estimate(ToEstMdl1,y1,'E0',e0,'U0',u0,'X',X1);

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.5928e-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.

ToEstMdl2 = arima('ARLags',[1 2],'MALags',1); EstMdl2 = estimate(ToEstMdl2,y2,'E0',e0,'Y0',... y0,'X',X2);

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.3634e-10

Convert

`EstMdl1`

to an ARIMAX model.`ConvertedMdl2 = arima(EstMdl1,'X',X1)`

ConvertedMdl2 = arima with properties: Description: "ARIMAX(2,0,1) Model (Gaussian Distribution)" 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 equivalent to the ARIMAX model constant converted from the regression model with ARIMA errors. In other words,

`EstMdl2.Constant = 0.0849961`

and`ConvertedMdl2.Constant = 0.087737`

. This is because`estimate`

does not enforce the nonlinear constraint that the`arima`

converter enforces. As a result, the other estimates are not equivalent either, albeit close.

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

.