%% Estimate ARIMA Model Containing Exogenous Predictors (ARIMAX)
% Fit an ARIMAX model to simulated time series data.
%
%% Simulate Predictor and Response Data
% Create the ARIMAX(2,1,0) model for the DGP, represented by y_t in the equation
% (1 - 0.5L + 0.3{L^2}){(1 - L)^1}{y_t} = 2 + 1.5{x_{1,t}} + 2.6{x_{2,t}} - 0.3{x_{3,t}} + {\varepsilon _t},
% where \varepsilon_{t} is a series of iid Gaussian random variables with mean 0 and variance 0.1.

DGP = arima('AR',{0.5,-0.3},'D',1,'Constant',2,...
    'Variance',0.1,'Beta',[1.5 2.6 -0.3]);

% Assume that the exogenous variables x_{1,t}, x_{2,t}, and x_{3,t} are represented by the AR(1) processes
% \begin{array}{*{20}{c}}
% {{x_{1,t}} = 0.1{x_{1,t - 1}} + {\eta _{1,t}}}\\
% {{x_{2,t}} = 0.2{x_{2,t - 1}} + {\eta _{2,t}}}\\
% {{x_{3,t}} = 0.3{x_{3,t - 1}} + {\eta _{3,t}},}
% \end{array}
% where \eta_{i,t} follows a Gaussian distribution with mean 0 and variance 0.01 for i\in \left\lbrace 1,2,3\right\rbrace. Create ARIMA models that represent the exogenous variables.

MdlX1 = arima('AR',0.1,'Constant',0,'Variance',0.01);
MdlX2 = arima('AR',0.2,'Constant',0,'Variance',0.01);
MdlX3 = arima('AR',0.3,'Constant',0,'Variance',0.01);

% Simulate length 1000 exogenous series from the AR models. Store the simulated data in a matrix.

T = 1000;
rng(10); % For reproducibility
x1 = simulate(MdlX1,T);
x2 = simulate(MdlX2,T);
x3 = simulate(MdlX3,T);
X = [x1 x2 x3];

% X is a 1000-by-3 matrix of simulated time series data.  Each row corresponds to an observation in the time series, and each column corresponds to an exogenous variable.
% Simulate a length 1000 series from the DGP. Specify the simulated exogenous data.

y = simulate(DGP,T,'X',X);

% y is a 1000-by-1 vector of response data.
%
%% Estimate Model
% Create an ARIMA(2,1,0) model template for estimation.

Mdl = arima(2,1,0)

% The model description (Description property) and value of Beta suggest that the partially specified arima model object Mdl is agnostic of the exogenous predictors.
%
% Estimate the ARIMAX(2,1,0) model; specify the exogenous predictor data. Because estimate backcasts for presample responses (a process that requires presample predictor data for ARIMAX models), fit the model to the latest T – Mdl.P responses. (Alternatively, you can specify presample responses by using the 'Y0' name-value pair argument.)

EstMdl = estimate(Mdl,y((Mdl.P + 1):T),'X',X);

% EstMdl is a fully specified, estimated arima model object.
% When you estimate the model by using estimate and supply the exogenous data by specifying the 'X' name-value pair argument, MATLAB® recognizes the model as an ARIMAX(2,1,0) model and includes a linear regression component for the exogenous variables.
% The estimated model is
% \left(1-0\ldotp 56L+0\ldotp 27L^2 \right){\left(1-L\right)}^1 y_t =1\ldotp 75+{1\ldotp 48x}_{1,t} +2\ldotp 56x_{2,t} -0\ldotp 34x_{3,t} +\varepsilon_t ,
% which resembles the DGP represented by Mdl0. Because MATLAB returns the AR coefficients of the model expressed in difference-equation notation, their signs are opposite in the equation.
% Copyright 2015 The MathWorks, Inc.