%% Initialize Multivariate Model Simulation from Multiple Starting Conditions
% This example shows how to initialize simulated paths from presample responses and initial states. The example uses arbitrary parameter values.
%
%% Fully Specify LSETAR Model
% Consider the following 2-D LSETAR model.
% State 1, "Low": {y_t} = \left[ {\begin{array}{*{20}{c}}
% 1\\
% { - 1}
% \end{array}} \right] + {\varepsilon _{1,t}}, where {\varepsilon _{1,t}}\sim{\rm N}\left( {\left[ {\begin{array}{*{20}{c}}
% 0\\
% 0
% \end{array}} \right],\left[ {\begin{array}{*{20}{c}}
% 1&{ - 0.1}\\
% { - 0.1}&1
% \end{array}} \right]} \right).
% State 2 , "Med": {y_t} = \left[ {\begin{array}{*{20}{c}}
% 2\\
% { -2}
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% 0.5 & 0.1\\
% 0.5 & 0.5
% \end{array}} \right] y_{t-1}+ {\varepsilon _{2,t}}, where {\varepsilon _{2,t}}\sim{\rm N}\left( {\left[ {\begin{array}{*{20}{c}}
% 0\\
% 0
% \end{array}} \right],\left[ {\begin{array}{*{20}{c}}
% 2&{ - 0.2}\\
% { - 0.2}&2
% \end{array}} \right]} \right).
% State 3, "High": {y_t} = \left[ {\begin{array}{*{20}{c}}
% 3\\
% { -3}
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% 0.25 & 0\\
% 0 & 0
% \end{array}} \right] y_{t-1} + \left[ {\begin{array}{*{20}{c}}
% 0 & 0\\
% 0.25 & 0
% \end{array}} \right] y_{t-2} + {\varepsilon _{3,t}}, where {\varepsilon _{3,t}}\sim{\rm N}\left( {\left[ {\begin{array}{*{20}{c}}
% 0\\
% 0
% \end{array}} \right],\left[ {\begin{array}{*{20}{c}}
% 3&{ - 0.3}\\
% { - 0.3}&3
% \end{array}} \right]} \right).
% The system is in state 1 when y_{2,t-4} <-1, the system is in state 2 when -1\le y_{2,t-4} <1, and the system is in state 3 otherwise.
% The transition function is logistic. The transition rate from state 1 to 2 is 3.5, and the transition rate from state 1 to 3 is 1.5.
% Create logistic threshold transitions at mid-levels -1 and 1, with rates 3.5 and 1.5 respectively. Label the states.

t = [-1 1];
r = [3.5 1.5];
stateNames = ["Low" "Med" "High"];
tt = threshold(t,Type="logistic",Rates=[3.5 1.5],StateNames=stateNames);

% Create the VAR submodels by using varm. Store the submodels in a vector with order corresponding to the regimes in tt.StateNames.

% Constants (numSeries x 1 vectors)
C1 = [1; -1];
C2 = [2; -2];
C3 = [3; -3];

% Autoregression coefficients (numSeries x numSeries matrices)
AR1 = {};                            % 0 lags
AR2 = {[0.5 0.1; 0.5 0.5]};          % 1 lag
AR3 = {[0.25 0; 0 0] [0 0; 0.25 0]}; % 2 lags

% Innovations covariances (numSeries x numSeries matrices)
Sigma1 = [1 -0.1; -0.1 1];
Sigma2 = [2 -0.2; -0.2 2];
Sigma3 = [3 -0.3; -0.3 3];

% VAR Submodels
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];

% Create an LSETAR model from the switching mechanism tt and the state-specific submodels mdl. Label the series Y1 and Y2.

Mdl = tsVAR(tt,mdl,SeriesNames=["Y1" "Y2"])

% Mdl is a fully specified tsVAR object representing a multivariate three-state LSETAR model. tsVAR object functions enable you to specify threshold variable characteristics and data.
%
%% Initialize Simulation from Presample Responses
% Consider simulating 5 paths initialized from presample responses. Specify a numPreObs-by-numSeries-by-numPaths array of presample responses, where:
% numPreObs is the number of presample responses per series and path. You must specify enough presample observations to initialize all AR components in the VAR models and the endogenous threshold variable. The largest AR component order is 2 and the threshold variable delay is 4, therefore simulate requires numPreObs=4 presample observations per series and path.
% numSeries=2, the number of response series in the system.
% numPaths=5, the number of independent paths to simulate.

delay = 4;
numPaths = 5;
Y0 = zeros(delay,Mdl.NumSeries,numPaths);
for j = 2:numPaths
    Y0(:,:,j) = 10*j*ones(delay,Mdl.NumSeries);
end

% Simulate 10 paths of length 100 from the LSETAR model from the presample. Specify the endogenous threshold variable and its delay, y_{2,t-4}.

numObs = 100;
rng(1);
Y = simulate(Mdl,numObs,NumPaths=numPaths,Y0=Y0,Index=2,Delay=4);

% Y is a 100-by-2-by-5 array of simulate response paths. For example, Y(50,2,3) is the simulated response of path 3, of series Y2, at time point 50.
% Plot the simulated paths for each variable on separate plots.

tiledlayout(2,1)
nexttile
plot(squeeze(Y(:,1,:)))
title("Y1")
nexttile
plot(squeeze(Y(:,2,:)))
title("Y2")

% The system quickly settles regardless of the presample.
%
%% Initialize Simulation from States
% Simulate three paths of length 100, where each of the three states initialize a path. Specify state indices for initialization, and specify the endogenous threshold variable and its delay.

S0 = 1:Mdl.NumStates;
numPaths = numel(S0);
Y = simulate(Mdl,numObs,NumPaths=numPaths,S0=S0,Index=2,Delay=4);
tiledlayout(2,1)
nexttile
plot(squeeze(Y(:,1,:)))
title("Y1")
nexttile
plot(squeeze(Y(:,2,:)))
title("Y2")

% Copyright 2021 The MathWorks, Inc.