%% Display Estimation Summary for One State
% Create the following fully specified SETAR model for the DGP.
% State 1: \left[ {\begin{array}{*{20}{c}}
% {{y_{1,t}}}\\
% {{y_{2,t}}}
% \end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
% -1\\
% -4
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% {-0.5}&{0.1}\\
% {0.2}&{-0.75}
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{y_{1,t - 1}}}\\
% {{y_{2,t - 1}}}
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% {{\varepsilon _{1,t}}}\\
% {{\varepsilon _{2,t}}}
% \end{array}} \right].
% State 2: \[\left[ {\begin{array}{*{20}{c}}
% {{y_{1,t}}}\\
% {{y_{2,t}}}
% \end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
% 1\\
% 4
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% {{\varepsilon _{1,t}}}\\
% {{\varepsilon _{2,t}}}
% \end{array}} \right]\].
% State 3:\[\left[ {\begin{array}{*{20}{c}}
% {{y_{1,t}}}\\
% {{y_{2,t}}}
% \end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
% 1\\
% 4
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% {0.5}&{0.1}\\
% {0.2}&{0.75}
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{y_{1,t - 1}}}\\
% {{y_{2,t - 1}}}
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% {{\varepsilon _{1,t}}}\\
% {{\varepsilon _{2,t}}}
% \end{array}} \right]\].
% \left[ {\begin{array}{*{20}{c}}
% {{\varepsilon _{1,t}}}\\
% {{\varepsilon _{2,t}}}
% \end{array}} \right]\sim {\rm N}_2\left(\left[\begin{array}{*{20}{c}}
% 0\\
% 0
% \end{array}\right],\left[\begin{array}{*{20}{c}}
% 2 & -1\\
% -1 & 1
% \end{array}\right]\right).
% The system is in state 1 when y_{2,t-4} <-3, the system is in state 2 when -3\le y_{2,t-4} <3, and the system is in state 3 otherwise.

t = [-3 3];
ttDGP = threshold(t);
constant1 = [-1; -4];
constant2 = [1; 4];
constant3 = [1; 4];
AR1 = [-0.5 0.1; 0.2 -0.75];
AR3 = [0.5 0.1; 0.2 0.75];
Sigma = [2 -1; -1 1];
mdl1DGP = varm(Constant=constant1,AR={AR1});
mdl2DGP = varm(Constant=constant2);
mdl3DGP = varm(Constant=constant3,AR={AR3});
mdlDGP = [mdl1DGP; mdl2DGP; mdl3DGP];
MdlDGP = tsVAR(ttDGP,mdlDGP,Covariance=Sigma);

% DIsplay a summary of the unestimated DGP.

summarize(MdlDGP)

% summarize prints an object display.
% Generate a random response path of length 500 from the DGP. Specify that second response variable with a delay of 4 as the threshold variable.

rng(10) % For reproducibiliy
y = simulate(MdlDGP,500,Index=2,Delay=4);

% Create a partially specified threshold-switching model that has the same structure as the DGP, but specify the transition mid-level, submodel coefficients, and model-wide covariance as unknown for estimation.

tt = threshold([NaN; NaN]);
mdlar = varm(2,1);
mdlc = varm(2,0);
Mdl = tsVAR(tt,[mdlar; mdlc; mdlar],Covariance=nan(2));

% Fully specify a threshold transition that has the same structure as tt, but set the mid-levels to -1 and 1.

t0 = [-1 1];
tt0 = threshold(t0);

% Fit the threshold-switching model to the simulated series. Specify the threshold variable y_{2,t-4}. Plot the loglikelihood after each iteration of the threshold search algorithm.

EstMdl = estimate(Mdl,tt0,y,IterationPlot=true,Index=2,Delay=4);

%
% The plot displays the evolution of the loglikelihood as the estimation procedure searches for optimal levels. The procedure terminates when one of the stopping criteria is satisfied.
% Display an estimation summary for state 3 only.

summarize(EstMdl,3)

% Copyright 2020 The MathWorks, Inc.