%% Simulate Paths from Model with VARX Submodels
% Generate random paths 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 = [10 1 1; 1 10 1; 1 1 10];
mc = dtmc(P);

% mc is a fully specified dtmc object.  dtmc normalizes the rows of P so that they sum to 1.
% For each regime, use varm to create a VAR model that describes the response process within the regime. Specify all parameter values.

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

% Autoregression coefficients
AR1 = {};                            
AR2 = {[0.5 0.1; 0.5 0.5]};          
AR3 = {[0.25 0; 0 0] [0 0; 0.25 0]}; 

% Regression coefficients
Beta1 = [1;-1];           
Beta2 = [2 2;-2 -2];      
Beta3 = [3 3 3;-3 -3 -3]; 

% Innovations covariances
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];

% 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 Ignoring Regression Component
% If you do not supply exogenous data, simulate ignores the regression components in the submodels. Simulate 3 separate, independent paths of responses, innovations, and state indices of length 5 from the model.

rng(1); % For reproducibility
[Y,E,SP] = simulate(Mdl,5,'NumPaths',3)

% Y and E are 5-by-2-by-3 arrays of simulated responses and innovations, respectively.  Rows correspond to time points, columns correspond to variables in the system, and pages correspond to paths. SP is a 5-by-3 matrix whose columns correspond to paths.
% Simulate a single path of responses, innovations, and states into a simulation horizon of length 50. Then plot each path separately.

rng0 = rng; % Store settings to reproduce state sequence.
[Y,E,SP] = simulate(Mdl,50);

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 Data Including Regression Component
% Simulate exogenous data for the three regressors by generating 50 random observations from the 3-D standard Gaussian distribution.

X = randn(50,3);

% Generate one random response, innovation, and state path of length 50. Specify the simulated exogenous data for the submodel regression components. Plot the results.

rng(rng0); % Reproduce state sequence in previous simulation.
[Y,E,SP] = simulate(Mdl,50,'X',X);

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])

% Copyright 2019 The MathWorks, Inc.