%% Perform Constrained Estimation
% Fit simulated data to a Markov-switching dynamic regression model with VARX submodels.  Specify equality constraints for estimation. This example uses arbitrary parameter values.
%
%% Create Model for DGP
% Create a fully specified, three-state discrete-time Markov chain model for the switching mechanism.

P = [0.8 0.1 0.1; 0.2 0.6 0.2; 0 0.1 0.9];
mc = dtmc(P);

% For each state, create a fully specified VARX(1) model for the response process. Specify the same model constant and lag 1 AR coefficient matrix for all submodels. For each model, specify a different 2-by-1 regression coefficient for the one exogenous variable.

% Constants
C = [1;-1];

% Autoregression coefficients
AR = {[0.6 0.1; 0.4 0.2]};

% Regression coefficients
Beta1 = [0.2;-0.4];
Beta2 = [0.6;-1.0];
Beta3 = [0.9;-1.3];

% VAR Submodels
dgp = varm('Constant',C,'AR',AR,'Covariance',5*eye(2));
dgp1 = dgp;
dgp1.Beta = Beta1;
dgp2 = dgp;
dgp2.Beta = Beta2;
dgp3 = dgp;
dgp3.Beta = Beta3;

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

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

%
%% Simulate Data from DGP
% Simulate data for the exogenous series by generating 1000 observations from the Gaussian distribution with mean 0 and variance 100.

rng(1); % For reproducibility
X = 10*randn(1000,1);

% Generate a random response path of length 1000 from the DGP. Specify the simulated exogenous data for the submodel regression components.

Data = simulate(DGP,1000,'X',X);

% Data is a 1000-by-1 vector 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 regression coefficients. Specify the true values of the constants and AR coefficient matrices.

PEst = NaN(3);
mcEst = dtmc(PEst);
mdl = varm(2,1);
mdl.Constant = C;
mdl.AR = AR;
mdl.Beta = NaN(2,1);
Mdl = msVAR(mcEst,[mdl; mdl; mdl]);

% Because the values of the constants and AR coefficient matrices are specified in Mdl, estimate treats them as equality constraints for estimation.
%
%% 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 and set parameters with equality constraints to their values specified in Mdl.

P0 = (1/3)*ones(3);
mc0 = dtmc(P0);
mdl0 = varm('Constant',C,'AR',AR,'Covariance',eye(2));
mdl01 = mdl0;
mdl01.Beta = [0.1;-0.1];
mdl02 = mdl0;
mdl02.Beta = [0.5;-0.5];
mdl03 = mdl0;
mdl03.Beta = [1;-1];
Mdl0 = msVAR(mc0,[mdl01; mdl02; mdl03]);

%
%% Estimate Model
% Fit the model to the simulated data.  Specify the exogenous data for the regression component. Plot the loglikelihood at each iteration of the EM algorithm.

figure
EstMdl = estimate(Mdl,Mdl0,Data,'X',X,'IterationPlot',true);

%
%% Assess Accuracy
% Compare the estimated regression coefficient vectors and transition matrix to their true values.

Beta1
Beta1Estimate = EstMdl.Submodels(1).Beta

Beta2
Beta2Estimate = EstMdl.Submodels(2).Beta

Beta3
Beta3Estimate = EstMdl.Submodels(3).Beta

P
PEstimate = EstMdl.Switch.P

% Copyright 2019 The MathWorks, Inc.