%% Perform Monte Carlo Estimation
% Consider the model in docid:econ_ug.mw_3756343a-4312-4918-969d-0758bf8fee0e.
%
%% Create Fully Specified Model
% Create the Markov-switching model excluding the regression component.

P = [10 1 1; 1 10 1; 1 1 10];
mc = dtmc(P);
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]}; 
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,'Covariance',Sigma1);
mdl2 = varm('Constant',C2,'AR',AR2,'Covariance',Sigma2);
mdl3 = varm('Constant',C3,'AR',AR3,'Covariance',Sigma3);
mdl = [mdl1; mdl2; mdl3];
Mdl = msVAR(mc,mdl);

%
%% Simulate Multiple Paths
% Generate 1000 random paths of responses for 50 time steps. Start all simulations at the first state.

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

% Y is a 50-by-2-by-1000 array of simulated response paths.
%
%% Compute Monte Carlo Distribution
% For each variable and path, compute the process mean.

mus = mean(Y,1);

% For each variable, plot the Monte Carlo distribution of the process mean.

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 variable and path, compute the process standard deviation.

sigmas = std(Y,0,1);

% For each variable, plot the Monte Carlo distribution of the process standard deviation.

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.