%% Display Estimation Summary Separately for Each State
% Create the following fully specified Markov-switching model the DGP.
% State transition matrix: P = \left[ {\begin{array}{*{20}{c}}
% {0.5}&{0.2}&{0.3}\\
% {0.2}&{0.6}&{0.2}\\
% 0.2&{0.1}&{0.7}
% \end{array}} \right].
% State 1: \left[ {\begin{array}{*{20}{c}}
% {{y_{1,t}}}\\
% {{y_{2,t}}}
% \end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
% -1\\
% -1
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% {-0.5}&{0.1}\\
% {0.2}&{-0.75}
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{y_{1,t - 1}}}\\
% {{y_{2,t - 1}}}
% \end{array}} \right] + \varepsilon _{1,t}, where \varepsilon_{1,t} \sim N_2 \left(\left\lbrack \begin{array}{c}
% 0\\
% 0
% \end{array}\right\rbrack ,\left\lbrack \begin{array}{cc}
% 0\ldotp 5 & 0\\
% 0 & 1
% \end{array}\right\rbrack \right).
% State 2: \[\left[ {\begin{array}{*{20}{c}}
% {{y_{1,t}}}\\
% {{y_{2,t}}}
% \end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
% -1\\
% 2
% \end{array}} \right] +\varepsilon _{2,t}, where \varepsilon_{2,t} \sim N_2 \left(\left\lbrack \begin{array}{c}
% 0\\
% 0
% \end{array}\right\rbrack ,\left\lbrack \begin{array}{cc}
% 1 & 0\\
% 0 & 1
% \end{array}\right\rbrack \right).
% State 3:\[\left[ {\begin{array}{*{20}{c}}
% {{y_{1,t}}}\\
% {{y_{2,t}}}
% \end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
% 1\\
% 2
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% {0.5}&{0.1}\\
% {0.2}&{0.75}
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{y_{1,t - 1}}}\\
% {{y_{2,t - 1}}}
% \end{array}} \right] + \varepsilon _{3,t}, where \varepsilon_{3,t} \sim N_2 \left(\left\lbrack \begin{array}{c}
% 0\\
% 0
% \end{array}\right\rbrack ,\left\lbrack \begin{array}{cc}
% 1 & -0\ldotp 1\\
% -0\ldotp 1 & 2
% \end{array}\right\rbrack \right).

PDGP = [0.5 0.2 0.3; 0.2 0.6 0.2; 0.2 0.1 0.7];
mcDGP = dtmc(PDGP);

constant1 = [-1; -1];
constant2 = [-1; 2];
constant3 = [1; 2];
AR1 = [-0.5 0.1; 0.2 -0.75];
AR3 = [0.5 0.1; 0.2 0.75];
Sigma1 = [0.5 0; 0 1];
Sigma2 = eye(2);
Sigma3 = [1 -0.1; -0.1 2];

mdl1DGP = varm(Constant=constant1,AR={AR1},Covariance=Sigma1);
mdl2DGP = varm(Constant=constant2,Covariance=Sigma2);
mdl3DGP = varm(Constant=constant3,AR={AR3},Covariance=Sigma3);
mdlDGP = [mdl1DGP; mdl2DGP; mdl3DGP];

MdlDGP = msVAR(mcDGP,mdlDGP);

% Generate a random response path of length 1000 from the DGP.

rng(1) % For reproducibiliy
Y = simulate(MdlDGP,1000);

% Create a partially specified Markov-switching model that has the same structure as the DGP, but that the transition matrix, and all submodel coefficients and innovations covariance matrices are unknown and estimable.

mc = dtmc(nan(3));
mdlar = varm(2,1);
mdlc = varm(2,0);
Mdl = msVAR(mc,[mdlar; mdlc; mdlar]);

% Initialize the estimation procedure by fully specifying a Markov-switching model that has the same structure as Mdl, but has the following parameter values:
% A randomly drawn transition matrix
% Randomly drawn contant vectors for each model
% AR self lags of 0.1 and cross lags of 0.
% The identify matrix for the innovations covariance.

P0 = randi(10,3,3);
mc0 = dtmc(P0);
constant01 = randn(2,1);
constant02 = randn(2,1);
constant03 = randn(2,1);
AR0 = 0.1*eye(2);
Sigma0 = eye(2);
mdl01 = varm(Constant=constant01,AR={AR0},Covariance=Sigma0);
mdl02 = varm(Constant=constant02,Covariance=Sigma0);
mdl03 = varm(Constant=constant03,AR={AR0},Covariance=Sigma0);
submdl0 = [mdl01; mdl02; mdl03];
Mdl0 = msVAR(mc0,submdl0);

% Fit the Markov-switching model to the simulated series. Plot the loglikelihood after each iteration of the EM algorithm.

EstMdl = estimate(Mdl,Mdl0,Y,IterationPlot=true);

% The plot displays the evolution of the loglikelihood with increasing iterations of the EM algorithm. The procedure terminates when one of the stopping criteria is satisfied.
% Display an estimation summary of the model.

summarize(EstMdl)

% Display an estimation summary separately for each state.

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

%
%
% Copyright 2021 The MathWorks, Inc.