%% Specify Presample Data For STAR Model Estimation
% Consider a smooth threshold-switching (STAR) model for the real US GDP growth rate, where each submodel is AR(4) and the threshold variable is the unemployment growth rate.
% Create a partially specified threshold transition for the unemployment growth rate. Specify the normal cdf transition function, and an unknown, estimable mid-level and rate. Label the states "Contaction" and "Expansion".

tt = threshold(NaN,Type="normal",Rates=NaN,...
    StateNames=["Contraction" "Expansion"]);

% tt is a partially specified threshold object, and it is agnostic of the variable and data it represents.
% Create a threshold-switching model for the real US GDP growth rate. Label the series "rRGDP".

mdl = arima(4,0,0);
submdls = [mdl; mdl];
Mdl = tsVAR(tt,submdls,SeriesNames="rRGDP");

% Load the quarterly US macroeconomic data set Data_USEconModel.mat. Compute the real GDP percent growth and unemployment growth.

load Data_USEconModel

DataTable = rmmissing(DataTable,DataVariables=["GDP" "GDPDEF" "UNRATE"]);
RGDP = DataTable.GDP./DataTable.GDPDEF;
rRGDP = price2ret(RGDP)*100;        % Response data
gUNRATE = diff(DataTable.UNRATE);   % Exogenous thresold data
dates = DataTable.Time(2:end);

% Suppose the estimation period includes growth rates from the first quarter of 1950. Identify the estimation period in the date base.

startEst = datetime(1950,1,1);
idxEst = dates >= startEst;

% To initialize the estimation procesdure, fully specify a threshold transition that has the same structure as tt, but set the mid-level to -0.5 and a rate of 100.

tt0 = threshold(-0.5,Type=tt.Type,Rates=100,StateNames=tt.StateNames);

% Fit the STAR model to the estimation period of the US GDP growth rate series. Specify the following parameters:
% Set Y0 to the responses before the estimation period to initialize the AR submodel components.
% Set Type to "exogenous" to characterize the threshold variable.
% Set Z to the threshold variable data gUNRATE in the estimation period.
% Set MaxRate to 150 to expand the search for the optimal transition function rate.
% Plot the evolution of the loglikelihood.

EstMdl = estimate(Mdl,tt0,rRGDP(idxEst),Y0=rRGDP(~idxEst),...
    Z=gUNRATE(idxEst),Type="exogenous",MaxRate=150,IterationPlot=true);

% Display an estimation summary and plot the estimated switching mechanism with the threshold data.

summarize(EstMdl)
figure
ttplot(EstMdl.Switch,Data=gUNRATE(idxEst))
dEst = dates(idxEst);
xticklabels(string(dEst(xticks)))

% The large rate suggests little mixing occurs between the models of each regime. When the quarterly unempoyment growth is less than 0.206%, the dominant model for the US GDP growth rate is in EstMdl.Submodels(1). Otherwise, the dominant submodel is in EstMdl.Submodels(2).
% Many of the coefficients are insignificant, which can suggest a search for a simpler model.
% Copyright 2021 The MathWorks, Inc.