%% Fit Model to Simulated Data
% Assess estimation accuracy using simulated data from a known data-generating process (DGP). This example uses arbitrary parameter values.
%
%% Create Model for DGP
% Create a fully specified, two-state discrete-time Markov chain model for the switching mechanism.

P = [0.7 0.3; 0.1 0.9];
mc = dtmc(P);

% For each state, create a fully specified AR(1) model for the response process.

% Constants
C1 = 5;
C2 = -2;

% Autoregression coefficients
AR1 = 0.4;
AR2 = 0.2;

% Variances
V1 = 4;
V2 = 2;

% AR Submodels
dgp1 = arima('Constant',C1,'AR',AR1,'Variance',V1);
dgp2 = arima('Constant',C2,'AR',AR2,'Variance',V2);

% Create a fully specified Markov-switching dynamic regression model for the DGP.

DGP = msVAR(mc,[dgp1,dgp2]);

%
%% Simulate Response Paths from DGP
% Generate 10 random response paths of length 1000 from the DGP.

rng(1); % For reproducibility
N = 10;
n = 1000;
Data = simulate(DGP,n,'Numpaths',N);

% Data is a 1000-by-10 matrix of simulated responses.
%
%% Create Model for Estimation
% Create a partially specified Markov-switching dynamic regression model that has the same structure as the data-generating process, but specify an  unknown transition matrix and unknown submodel coefficients.

PEst = NaN(2);
mcEst = dtmc(PEst);
mdl = arima(1,0,0);
Mdl = msVAR(mcEst,[mdl; mdl]);

%
%% Create Model Containing Initial Values
% Create a fully specified Markov-switching dynamic regression model that has the same structure as Mdl, but set all estimable parameters to initial values.

P0 = 0.5*ones(2);
mc0 = dtmc(P0);
mdl01 = arima('Constant',1,'AR',0.5,'Variance',2);
mdl02 = arima('Constant',-1,'AR',0.5,'Variance',1);
Mdl0 = msVAR(mc0,[mdl01,mdl02]);

%
%% Estimate Models
% Fit the model to each simulated path.  For each path, plot the loglikelihood at each iteration of the EM algorithm.

c1 = zeros(N,1);
c2 = zeros(N,1);
v1 = zeros(N,1);
v2 = zeros(N,1);
ar1 = zeros(N,1);
ar2 = zeros(N,1);
PStack = zeros(2,2,N);

figure
hold on

for i = 1:N

    EstModel = estimate(Mdl,Mdl0,Data(:,i),'IterationPlot',true);
    
    c1(i) = EstModel.Submodels(1).Constant;
    c2(i) = EstModel.Submodels(2).Constant;
    v1(i) = EstModel.Submodels(1).Covariance;
    v2(i) = EstModel.Submodels(2).Covariance;
    ar1(i) = EstModel.Submodels(1).AR{1};
    ar2(i) = EstModel.Submodels(2).AR{1};
    PStack(:,:,i) = EstModel.Switch.P;
    
end

hold off

%
%% Assess Accuracy
% Compute the Monte Carlo mean of each estimated parameter.

c1Mean = mean(c1);
c2Mean = mean(c2);
v1Mean = mean(v1);
v2Mean = mean(v2);
ar1Mean = mean(ar1);
ar2Mean = mean(ar2);
PMean = mean(PStack,3);

% Compare population parameters to the corresponding Monte Carlo estimates.

DGPvsEstimate = [...
    C1  c1Mean
    C2  c2Mean
    V1  v1Mean
    V2  v2Mean
    AR1 ar1Mean
    AR2 ar2Mean]
P
PEstimate = PMean

% Copyright 2019 The MathWorks, Inc.