%% Specify Presample Data
% Consider the data in Compute Filtered Probabilities of Recession, but assume that the period of interest is 1960:Q1–2004:Q2. Also, consider adding an autoregressive term to each submodel.
%
%% Create Partially Specified Model for Estimation
% Create a partially specified Markov-switching dynamic regression model for estimation. Specify AR(1) submodels.

P = NaN(2);
mc = dtmc(P,'StateNames',["Expansion" "Recession"]);
mdl = arima(1,0,0);
Mdl = msVAR(mc,[mdl; mdl]);

% Because the submodels are AR(1), each requires one presample observation to initialize its dynamic component for estimation.
%
%% Create Fully Specified Model Containing Initial Values
% Create the model containing initial parameter values for the estimation procedure.

mc0 = dtmc(0.5*ones(2),'StateNames',["Expansion" "Recession"]);
submdl01 = arima('Constant',1,'Variance',1,'AR',0.001);
submdl02 = arima('Constant',-1,'Variance',1,'AR',0.001);
Mdl0 = msVAR(mc0,[submdl01; submdl02]);

%
%% Load and Preprocess Data
% Load the data. Transform the entire set to an annualized rate series.

load Data_GDP
qrate = diff(Data)./Data(1:(end - 1)); 
arate = 100*((1 + qrate).^4 - 1);

% Identify the presample and estimation sample periods using the dates associated with the annualized rate series. Because the transformation applies the first difference, you must drop the first observation date from the original sample.

dates = datetime(dates(2:end),'ConvertFrom','datenum',...
    'Format','yyyy:QQQ','Locale','en_US');
estPrd = datetime(["1960:Q2" "2004:Q2"],'InputFormat','yyyy:QQQ',...
    'Format','yyyy:QQQ','Locale','en_US');
idxEst = isbetween(dates,estPrd(1),estPrd(2));
idxPre = dates < (estPrd(1)); 

%
%% Estimate Model
% Fit the model to the estimation sample data. Specify the presample observation.

arate0 = arate(idxPre);
arateEst = arate(idxEst);
EstMdl = estimate(Mdl,Mdl0,arateEst,'Y0',arate0);

% EstMdl is a fully specified msVAR object.
%
%% Compute State Probabilities
% Compute filtered and smoothed state probabilities from the estimated model and data in the estimation period.  Specify the presample observation. Plot the estimated probabilities of recession on subplots in the same figure.

FS = filter(EstMdl,arateEst,'Y0',arate0);
SS = smooth(EstMdl,arateEst,'Y0',arate0);

figure;
subplot(2,1,1)
plot(dates(idxEst),FS(:,2),'r')
title("Current Filter Probabilities and NBER Recessions")
recessionplot
subplot(2,1,2)
plot(dates(idxEst),SS(:,2),'r')
title("Full-Sample Smoothed Probabilities and NBER Recessions")
recessionplot

% Copyright 2019 The MathWorks, Inc.
%