%% Perform Constrained Estimation
% This example shows how to specify equality constraints on estimable submodel coefficients and threshold parameters.
% Consider the model for the real US GDP growth rate in docid:econ_ug#mw_badc2965-b156-4786-893e-1cbd387d0668. Assume the submodels are AR(4), but they include only the model constant (time trend in differenced series) and fourth lag. Also, suppose the transition function rate 3 for the unemployment growth series.
% Create the partially specified threshold transition for the unemployment growth rate. Specify the known transition rate.

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

% Create a threshold-switching model for the real US GDP growth rate. Explicitly specify the AR lags to include in the submodels by using the ARLags option.

mdl = arima(ARLags=4);
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;

% Fully specify a threshold transition that has the same structure as tt, including the equality constraint on the transition function rate. Set the initial mid-level to 0.

tt0 = threshold(0,Type=tt.Type,Rates=rConstraint,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",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)))

% Copyright 2021 The MathWorks, Inc.