%% Estimate Model Containing Regression Component
% Consider the following data-generating process (DGP) for a 2-D time-varying TAR (TVTAR) model containing an exogenous regression component.
% State 1: {y_t} = \left[ {\begin{array}{*{20}{c}}
% 1\\
% { - 1}
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% {0.6}&{0.1}\\
% {0.4}&{0.2}
% \end{array}} \right]{y_{t - 1}} + {x_t}\left[ {\begin{array}{*{20}{c}}
% {0.2}\\
% { - 0.4}
% \end{array}} \right] + {\varepsilon _{1,t}}, where \varepsilon_{1,t} \sim N_2 \left(0,0\ldotp 5I_2 \right) and I_2 is the 2-by-2 identity matrix.
% State 2: {y_t} = \left[ {\begin{array}{*{20}{c}}
% 2\\
% { - 2}
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% {0.6}&{0.1}\\
% {0.4}&{0.2}
% \end{array}} \right]{y_{t - 1}} + {x_t}\left[ {\begin{array}{*{20}{c}}
% {0.6}\\
% { - 1.0}
% \end{array}} \right] + {\varepsilon _{2,t}}, where \varepsilon_{2,t} \sim N_2 \left(0,I_2 \right).
% State 3: {y_t} = \left[ {\begin{array}{*{20}{c}}
% 3\\
% { -3}
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% {0.6}&{0.1}\\
% {0.4}&{0.2}
% \end{array}} \right]{y_{t - 1}} + {x_t}\left[ {\begin{array}{*{20}{c}}
% {0.9}\\
% { - 1.3}
% \end{array}} \right] + {\varepsilon _{3,t}}, where \varepsilon_{3,t} \sim N_2 \left(0,1\ldotp 5I_2 \right).
% The exogenous threshold variable z_t represents a linear time trend over the sample period. In this example, z_t is the sampling period normalized by the sample size.
% The system is in state 1 when z_t <0\ldotp 3, the system is in state 2 when 0\ldotp 3\le z_t <0\ldotp 6, and the system is in state 3 otherwise.
% Create a TVTAR model that represents the DGP.

tDGP = [0.3 0.6];
ttDGP = threshold(tDGP);

% Constant vectors
C1 = [1;-1];
C2 = [2;-2];
C3 = [3;-3];

% Common Autoregression coefficient
AR = {[0.6 0.1; 0.4 0.2]};

% Regression coefficient vectors
Beta1 = [0.2;-0.4];
Beta2 = [0.6;-1.0];
Beta3 = [0.9;-1.3];

% Innovations covariance matrices
Sigma1 = 0.5*eye(2);
Sigma2 = eye(2);
Sigma3 = 1.5*eye(2);

% VAR Submodels
mdl1 = varm(Constant=C1,AR=AR,Beta=Beta1,Covariance=Sigma1);
mdl2 = varm(Constant=C2,AR=AR,Beta=Beta2,Covariance=Sigma2);
mdl3 = varm(Constant=C3,AR=AR,Beta=Beta3,Covariance=Sigma3);
mdlDGP = [mdl1; mdl2; mdl3];

DGP = tsVAR(ttDGP,mdlDGP);

% For the exogenous predictors x_t, simulate a 2-D 1000-period from the Gaussian distribution with mean 0 and standard deviation 10.

rng(1); % For reproducibility
numPeriods = 1000;
X = 10*randn(numPeriods,2);

% For the threshold variable, create a 1000-element vector of equally spaced elements from 0 to 1.

z = linspace(0,1,numPeriods)';

% Simulate a 1000-period path from the DGP.

Y = simulate(DGP,numPeriods,X=X,Type="exogenous",Z=z);

% Create a TVTAR model template for estimation. Specify all estimable parameters as unknown by using NaNs.

t = [NaN NaN];
tt = threshold(t);

mdl = varm(2,1);        % Unknown contant and AR coefficient by default.
mdl.Beta = NaN(2,1);    % Configure exogenous regression component for estimation
Mdl = tsVAR(tt,[mdl; mdl; mdl]);

% Specify initial values of 0.25 and 0.7 for the regime thresholds.

t0 = [.25 .7];
tt0 = threshold(t0);

% The largest lag among all models is 1. Therefore the the estimation procedure requires one presample observation to initialize the model.
% Identify indices for the required presample, and identify indices for the estimation sample.

p = mdl.P;
idxPre = 1:p;
idxEst = (p + 1):numPeriods;

% Fit the TVTAR model to the data. Specify presample responses Y0, characterize the threshold variable Type and provide its data Z, and specify the exogenous data X.

EstMdl = estimate(Mdl,tt0,Y(idxEst,:),Y0=Y(idxPre,:),...
    Type="exogenous",Z=z(idxEst),X=X(idxEst,:));

% Display an estimation summary separately for each state, and display the estimated threshold transitions.

summarize(EstMdl,1)
summarize(EstMdl,2)
summarize(EstMdl,3)
EstMdl.Switch

% Copyright 2021 The MathWorks, Inc.