%% Simulate Multivariate Markov-Switching Dynamic Regression Model
% This example shows how to generate random response and state paths from a two Markov-switching dynamic regression model.
% Consider the response processes y_{1t} and y_{2t} that switch between 3 states, governed by the latent process s_t with this observed transition matrix
% \[P = \left[ {\begin{array}{*{20}{c}}
% 10&1&1\\
% 1&10&1\\
% 1&1&10
% \end{array}} \right].\]
% Suppose that \[y_t = \left[ {\begin{array}{*{20}{c}}
% y_{1t}&y_{2t}
% \end{array}} \right]^\prime\] is a VARX model in each state:
% \[{y_t} = \left\{ {\begin{array}{*{20}{l}}
% {\left[ {\begin{array}{*{20}{c}}
% 1\\
% { - 1}
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% 1\\
% { - 1}
% \end{array}} \right]{x_{1t}} + {\varepsilon _{1t}}}&{;{s_t} = 1}\\
% {\left[ {\begin{array}{*{20}{c}}
% 2\\
% { - 2}
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% {0.5}&{0.1}\\
% {0.5}&{0.5}
% \end{array}} \right]{y_{t - 1}} + \left[ {\begin{array}{*{20}{c}}
% 2&2\\
% { - 2}&{ - 2}
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{x_{1t}}}\\
% {{x_{2t}}}
% \end{array}} \right] + {\varepsilon _{2t}}}&{;{s_t} = 2}\\
% {\left[ {\begin{array}{*{20}{c}}
% 3\\
% { - 3}
% \end{array}} \right] + \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}} + \left[ {\begin{array}{*{20}{c}}
% 3&3&3\\
% { - 3}&{ - 3}&{ - 3}
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{x_{1t}}}\\
% {{x_{2t}}}\\
% {{x_{3t}}}
% \end{array}} \right] + {\varepsilon _{3t}}}&{;{s_t} = 3}.
% \end{array}} \right.\]
% where, for j = 1,2,3, \varepsilon_{\textrm{jt}} is Gaussian with mean 0 and covariance matrix
% \[\Sigma_j=j\left[ {\begin{array}{*{20}{c}}
% 1&{ - 0.1}\\
% { - 0.1}&1
% \end{array}} \right].\]
%
%% Create a Fully Specified Model
% Create the Markov-switching dynamic regression model that describes y_t and s_t (for more details, see docid:econ_ug.mw_acc27d90-3ad5-46e5-a47b-406d5df15972).

% Switching mechanism
P = [10 1 1; 1 10 1; 1 1 10];
mc = dtmc(P);

% VAR submodels 
C1 = [1;-1];
C2 = [2;-2];
C3 = [3;-3];
AR1 = {};                            
AR2 = {[0.5 0.1; 0.5 0.5]};          
AR3 = {[0.25 0; 0 0] [0 0; 0.25 0]}; 
Beta1 = [1; -1];
Beta2 = [2 2;-2 -2];
Beta3 = [3 3 3;-3 -3 -3];
Sigma1 = [1 -0.1; -0.1 1];
Sigma2 = [2 -0.2; -0.2 2];
Sigma3 = [3 -0.3; -0.3 3];
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 = msVAR(mc,[mdl1; mdl2; mdl3]);

% Mdl is a fully specified msVAR object.
%
%% Simulate One Path
% Simulate one path of responses, innovations, and state indices from the model. Specify a 50-period simulation horizon, that is, a 50-observation path for all variables.

rng('default') % For reproducibility
[Y,E,SP] = simulate(Mdl,50);

% Y and E are 50-by-2 vectors of simulated responses and innovations, respectively. SP is a 50-by-1 vector of state indices.  By default, simulate ignores the regression components.
% Plot the results in separate subplots.

figure
subplot(3,1,1); 
plot(Y);
ylabel("Response");
grid on;
legend("y_1","y_2")

subplot(3,1,2)
plot(E)
ylabel("Innovation")
grid on
legend("e_1","e_2")

subplot(3,1,3)
plot(SP,'m')
ylabel("State")
yticks([1 2 3])

%
%% Simulate Multiple Paths
% Simulate 3 separate, independent paths of responses and state indices from the model. Specify a 5-period simulation horizon.

[Y3,~,SP3] = simulate(Mdl,5,'NumPaths',3)

% Y is a 5-by-2-by-3 array of simulated responses. SP is a 5-by-3 matrix of simulated state indices.
% The simulated observation of y_{1t} \;in the second period of the first path is 1.9558.
% The simulated observation of y_{2t} \;in the second period of the first path is 3.9371.
% The corresponding simulated state in the second period of the first path is a state 3.
%
%% Simulate Model Containing Regression Component
% simulate requires exogenous data over the simulation horizon for the regression component. Generate a 50-by-3 matrix of iid standard normal observations for the required exogenous data.

X = randn(50,3);

% Simulate one path of responses, innovations, and state indices from the model. Specify a 50-period simulation horizon. Plot the results. To determine how the response and state paths compare to the model without regression components, reset the random seed to default.

rng('default');
[YReg,EReg,SPReg] = simulate(Mdl,50,'X',X);

figure
subplot(3,1,1)
plot([Y YReg])
ylabel("Response")
grid on
legend("y_1","y_2","yreg_1","yreg_2")

subplot(3,1,2)
plot([E EReg])
ylabel("Innovation")
grid on
legend("e_1","e_2","ereg_1","ereg_2")

subplot(3,1,3)
plot([SP SPReg])
ylabel("State")
legend("sp","spreg")
yticks([1 2 3]);

% Innovations and state paths are equivalent between models, but response paths differ.
% Copyright 2019 The MathWorks, Inc.