%% Compute Smoothed State Probabilities from Model with VARX Submodels
% Compute smoothed state probabilities from a three-state Markov-switching dynamic regression model for a 2-D VARX response process. This example uses arbitrary parameter values for the DGP.
%
%% Create Fully Specified Model for DGP
% Create a three-state discrete-time Markov chain model for the switching mechanism.

P = [5 1 1; 1 5 1; 1 1 5];
mc = dtmc(P);

% mc is a fully specified dtmc object.  dtmc normalizes the rows of P so that they sum to 1.
% For each state, create a fully specified VARX(0) model (constant and regression coefficient matrix only) for the response process. Specify different constant vectors across models.  Specify the same regression coefficient for the two regressors, and specify the same covariance matrix.  Store the VARX models in a vector.

% Constants
C1 = [1;-1];
C2 = [3;-3];
C3 = [5;-5];

% Regression coefficient
Beta = [0.2 0.1;0 -0.3];

% Covariance matrix
Sigma = [1.8 -0.4; -0.4 1.8];

% VARX submodels
mdl1 = varm('Constant',C1,'Beta',Beta,...
    'Covariance',Sigma);
mdl2 = varm('Constant',C2,'Beta',Beta,...
    'Covariance',Sigma);
mdl3 = varm('Constant',C3,'Beta',Beta,...
    'Covariance',Sigma);

mdl = [mdl1; mdl2; mdl3];

% mdl contains three fully specified varm model objects.
% For the DGP, create a fully specified Markov-switching dynamic regression model from the switching mechanism mc and the submodels mdl.

Mdl = msVAR(mc,mdl);

% Mdl is a fully specified msVAR model.
%
%% Simulate Data from DGP
% Simulate data for the two exogenous series by generating 30 observations from the standard 2-D Gaussian distribution.

rng(1) % For reproducibility
X = randn(30,2);

% Generate one random response and state path, both of length 30, from the DGP. Specify the simulated exogenous data for the submodel regression components.

[Y,~,SP] = simulate(Mdl,30,'X',X);

% Y is a 30-by-2 matrix of one simulated response path.  SP is a 30-by-1 vector of one simulated state path.
%
%% Compute State Probabilities
% Compute smoothed state probabilities from the DGP given the simulated response data.

SS = smooth(Mdl,Y,'X',X);

% SS is a 30-by-2 matrix of smoothed state probabilities for each period in the simulation horizon.
% Plot the simulated state path and the smoothed state probabilities on subplots in the same figure.

figure
subplot(2,1,1)
plot(SP,'m')
yticks([1 2 3])
legend({'Simulated states'})
subplot(2,1,2)
plot(SS,'-')
legend({'Smoothed s1','Smoothed s2','Smoothed s3'})

% Copyright 2019 The MathWorks, Inc.
%
%