%% Perform Monte Carlo Simulation Using Markov-Switching Dynamic Regression Model
% This example shows how to characterize the distribution of a multivariate response series, modeled by a Markov-switching dynamic regression model, by summarizing the draws of a Monte Carlo simulation.
% 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] + {\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}} + {\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}} + {\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.

% 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 Multiple Paths
% Simulate 1000 separate, independent paths of responses from the model. Specify a 50-period simulation horizon. Start all simulations in the first state (that is, the state of the system at time 0 is state 1), by specifying a distribution such that state 1 has all mass.

rng('default') % For reproducibility
S0 = [1 0 0];
Y = simulate(Mdl,50,'NumPaths',1000,'S0',S0);

% Y is a 50-by-2-by-1000 array of simulated responses. Each page is a simulated path from Mdl.
%
%% Summarize Monte Carlo Draws
% The processes are stationary. Therefore, for each path, obtain the typical value of each variable by computing the sample mean.

mus = mean(Y,1);

% mus is a 1-by-2-by-1000 array.  Each page is the sample mean vector of the simulated path.
% Plot the Monte Carlo probability distribution of the response means.

figure
h1 = histogram(mus(1,1,:),'Normalization',"probability",...
    'BinWidth',0.1);
hold on
h2 = histogram(mus(1,2,:),'Normalization',"probability",...
    'BinWidth',0.1);
legend("y_1","y_2")
title('Process Means')
hold off

% For each path, obtain the variability of each variable by computing the sample standard deviation. Plot the Monte Carlo probability distribution of the response standard deviations.

sigmas = std(Y,0,1);

figure
h1 = histogram(sigmas(1,1,:),'Normalization',"probability",...
    'BinWidth',0.05);
hold on
h2 = histogram(sigmas(1,2,:),'Normalization',"probability",...
    'BinWidth',0.05);
legend("y_1","y_2")
title('Process Standard Deviations')
hold off

% Copyright 2019 The MathWorks, Inc.