%% Forecast Model Containing Regression Component
% Consider including regression components for exogenous variables in each submodel of the threshold-switching dynamic regression model in docid:econ_ug#mw_8e7ca71d-cbea-4a0e-88f0-e67742f1ba55.
%
%% Fully Specify LSETAR Model
% Create logistic threshold transitions at mid-levels -1 and 1, with rates 3.5 and 1.5 respectively. Label the states.

t = [-1 1];
r = [3.5 1.5];
stateNames = ["Low" "Med" "High"];
tt = threshold(t,Type="logistic",Rates=[3.5 1.5],StateNames=stateNames)

% Assume the following VARX models describe the response processes of the system:
% State 1: {y_t} = \left[ {\begin{array}{*{20}{c}}
% 1\\
% { - 1}
% \end{array}} \right] +
% \left[ {\begin{array}{*{20}{c}}
% 1\\
% { - 1}
% \end{array}} \right]x_{1,t} +
% {\varepsilon _{1,t}}, where {\varepsilon _{1,t}}\sim{\rm N}\left( {\left[ {\begin{array}{*{20}{c}}
% 0\\
% 0
% \end{array}} \right],\left[ {\begin{array}{*{20}{c}}
% 1&{ - 0.1}\\
% { - 0.1}&1
% \end{array}} \right]} \right).
% State 2: {y_t} = \left[ {\begin{array}{*{20}{c}}
% 2\\
% { -2}
% \end{array}} \right]
% +
% \left[ {\begin{array}{*{20}{c}}
% 2 & 2\\
% { - 2}&-2
% \end{array}} \right]x_{2,t}
% + \left[ {\begin{array}{*{20}{c}}
% 0.5 & 0.1\\
% 0.5 & 0.5
% \end{array}} \right] y_{t-1}+ {\varepsilon _{2,t}}, where {\varepsilon _{2,t}}\sim{\rm N}\left( {\left[ {\begin{array}{*{20}{c}}
% 0\\
% 0
% \end{array}} \right],\left[ {\begin{array}{*{20}{c}}
% 2&{ - 0.2}\\
% { - 0.2}&2
% \end{array}} \right]} \right).
% State 3: {y_t} = \left[ {\begin{array}{*{20}{c}}
% 3\\
% { -3}
% \end{array}} \right]
% +
% \left[ {\begin{array}{*{20}{c}}
% 3 & 3 & 3\\
% { - 3}&-3&-3
% \end{array}} \right]x_{3,t}
% + \left[ {\begin{array}{*{20}{c}}
% 0.25 & 0\\
% 0 & 0
% \end{array}} \right] y_{t-1} + \left[ {\begin{array}{*{20}{c}}
% 0 & 0\\
% 0.25 & 0
% \end{array}} \right] y_{t-2} + {\varepsilon _{3,t}}, where {\varepsilon _{3,t}}\sim{\rm N}\left( {\left[ {\begin{array}{*{20}{c}}
% 0\\
% 0
% \end{array}} \right],\left[ {\begin{array}{*{20}{c}}
% 3&{ - 0.3}\\
% { - 0.3}&3
% \end{array}} \right]} \right).
% x_{1,t} represents a single exogenous variable, x_{2,t} represents two exogenous variables, and x_{3,t} represents three exogenous variables. Store the submodels in a vector.

% Constants (numSeries x 1 vectors)
C1 = [1; -1];
C2 = [2; -2];
C3 = [3; -3];

% Regression coefficients (numSeries x numRegressors matrices)
Beta1 = [1; -1];            % 1 regressor
Beta2 = [2 2; -2 -2];       % 2 regressors
Beta3 = [3 3 3; -3 -3 -3];  % 3 regressors

% Autoregression coefficients (numSeries x numSeries matrices)
AR1 = {};                            
AR2 = {[0.5 0.1; 0.5 0.5]};          
AR3 = {[0.25 0; 0 0] [0 0; 0.25 0]}; 

% Innovations covariances (numSeries x numSeries matrices)
Sigma1 = [1 -0.1; -0.1 1];
Sigma2 = [2 -0.2; -0.2 2];
Sigma3 = [3 -0.3; -0.3 3];

%VARX submodels
mdl1 = varm(Constant=C1,AR=AR1,Beta=Beta1,Covariance=Sigma1);
mdl2 = varm(Constant=C2,AR=AR2,Beta=Beta2,Covariance=Sigma2);
mdl3 = varm(Constant=C3,AR=AR3,Beta=Beta3,Covariance=Sigma3);

mdl = [mdl1; mdl2; mdl3];

% Create an LSETAR model from the switching mechanism tt and the state-specific submodels mdl. Label the series Y1 and Y2.

Mdl = tsVAR(tt,mdl,SeriesNames=["Y1" "Y2"]);

%
%% Forecast Responses Ignoring Regression Component
% If you do not supply exogenous data, forecast ignores the regression components in the submodels.
% Obtain a presample, from which to forecast, by simulating one 2-D path of 100 observations from the model. Specify the endogenous threshold variable and its delay, y_{2,t-4}.

rng(1)  % For reproducibility
numObs = 100;
Y0F = simulate(Mdl,numObs,Delay=4,Index=2);

% Forecast the LSETAR model into a 50-period horizon. Specify the endogenous threshold variable and its delay.

fh = 50;
YF = forecast(Mdl,Y0F,fh,Index=2,Delay=4);

% Plot the simulated path and forecasts for each variable on separate plots.

figure;
plot(1:numObs,Y0F(:,1),"b",(numObs+1):(numObs+fh),YF(:,1),"r--")
title("Y1")
legend(["Presample" "Forecasts"])
nexttile
plot(1:numObs,Y0F(:,2),"b",(numObs+1):(numObs+fh),YF(:,2),"r--")
legend(["Presample" "Forecasts"])
title("Y2")

%
%% Simulate Data Including Regression Component
% forecast requires exogenous data in order to generate random paths from the model.
% Assume the following AR models for the regressors:
% x_{1,t}=1+0.1x_{1,t-1}+\nu_{1,t}, where \nu_{1,t} \sim N\left(0,1\right).
% x_{2,t}=2+0.2x_{2,t-2}+\nu_{2,t}, where \nu_{2,t} \sim N\left(0,2\right).
% x_{3,t}=3+0.3x_{3,t-3}+\nu_{3,t}, where \nu_{1,t} \sim N\left(0,3\right).

MdlX1 = arima(Constant=1,AR={0.1},Variance=1);
MdlX2 = arima(Constant=2,AR={0 0.2},Variance=2);
MdlX3 = arima(Constant=3,AR={0 0 0.3},Variance=3);

% Simulate 100 observations from each model to represent in-sample data for the exogenous predictors.

X10F = simulate(MdlX1,numObs);
X20F = simulate(MdlX2,numObs);
X30F = simulate(MdlX3,numObs);
X0F = [X10F X20F X30F];

% Forecast 50 observations from each model to represent exogenous data in the forecast horizon. Specify the simulated in-sample data as a presample.

X1F = forecast(MdlX1,numObs,Y0=X10F);
X2F = forecast(MdlX2,numObs,Y0=X20F);
X3F = forecast(MdlX3,numObs,Y0=X30F);
XF = [X1F X2F X3F];

% Obtain presample response data, from which to forecast, by simulating one 2-D path of 100 observations from the model. Specify the endogenous threshold variable and its delay, and specify the simulated exogenous data.

Y0FX = simulate(Mdl,numObs,Delay=4,Index=2,X=X0F);

% Forecast the LSETAR model into a 50-period horizon. Specify the endogenous threshold variable and its delay, ans specify the exogenous data in the forecast horizon.

YFX = forecast(Mdl,Y0FX,fh,Index=2,Delay=4,X=XF);

figure;
plot(1:numObs,Y0FX(:,1),"b",(numObs+1):(numObs+fh),YFX(:,1),"r--")
title("Y1")
legend(["Presample" "Forecasts"])
nexttile
plot(1:numObs,Y0FX(:,2),"b",(numObs+1):(numObs+fh),YFX(:,2),"r--")
legend(["Presample" "Forecasts"])
title("Y2")

% Copyright 2021 The MathWorks, Inc.