%% Forecast LSTAR With Exogenous Threshold Variable
% Consider a logistic threshold-switching model (LSTAR) 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 logistic transition function, and an unknown, estimable mid-level and rate. Label the states "Contaction" and "Expansion".

tt = threshold(NaN,Type="logistic",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.

p = 4;
mdl = [arima(ARLags=p); arima(ARLags=p)];
Mdl = tsVAR(tt,mdl);

% Mdl is a partially specified tsVAR object specifying the sturcture of the model and which parameters are estimable.
% Load the quarterly US macroeconomic data set Data_USEconModel. Remove leading missing values in the series. 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);            % Response data
gUNRATE = diff(DataTable.UNRATE);   % Exogenous thresold data
dts = DataTable.Time(2:end);
T = numel(dts);

% Hold out the following sets from estimation:
% The first four observations as a presample for estimation.
% The final 25% as a forecast horizon to compare with the forecasts.

fh = ceil(0.25*T);
idxPre = 1:p;
idxEst = (idxPre(end)+1):(T-fh);
idxF = (T-fh+1):T;

% To initialize the estimation procedure, fully specify a threshold transition that has the same structure as tt, but set the mid-level to 0 and a rate of 1 (the default).

tt0 = threshold(0,Type=tt.Type);

% Fit the LSTAR 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.

EstMdl = estimate(Mdl,tt0,rRGDP(idxEst),Y0=rRGDP(idxPre),...
    Type="exogenous",Z=gUNRATE(idxEst));

% Forecast the real GDP growth rate series into the forecast horizon. Initialize the forecasts by specifying all responses in the estimation period (forecast uses only the latest, required observations). Characterize the threshold variable and provide its data in the forecast horizon.

yF = forecast(EstMdl,rRGDP(idxEst),fh,Type="exogenous",Z=gUNRATE(idxF));

% Plot the real GDP growth rate series with the forecasts.

figure
hold on
h1 = plot(dts(idxEst),rRGDP(idxEst),'b');
h2 = plot(dts(idxF),rRGDP(idxF),'b--');
h3 = plot(dts(idxF),yF,'r');
yfill = [ylim,fliplr(ylim)];
xfill = dts([idxEst(end) idxEst(end) idxF(end) idxF(end)]);
fill(xfill,yfill,'k',FaceAlpha=0.05)
legend([h1 h2 h3],["Data (estimation)" "Data (hold out)" "Forecasts"],...
    Location="NorthWest")
title("Real GDP Growth Rate Point Forecasts")
hold off

%
% Copyright 2021 The MathWorks, Inc.