%% Specify Presample Data
% Consider the Markov-switching model in docid:econ_ug.mw_c02ef02d-b3af-460e-bdbc-e7e55599fe0f, but assume that the submodels are AR(1) instead. Consider fitting the model to observations in the period 1960:Q1–2004:Q2.
%
% Create the model template for estimation. Specify AR(1) submodels.

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

% Because the submodels are AR(1), each requires one presample observation to initialize its dynamic component for estimation.
% 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 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);

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

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

%
% Simulate a response path from the fitted model over the estimation period.  Specify the presample observation.

rng(1) % For reproducibility
numObs = sum(idxEst);
aratesim = simulate(EstMdl,numObs,'Y0',y0);

% Plot the observations and simulated path of annualized rates, and identify periods of recession by using recessionplot.

figure;
plot(dates(idxEst),[arate(idxEst) aratesim])
recessionplot
xlabel("Time")
ylabel("Annualized GDP Rate")
legend("Observed","Simulated");

% Copyright 2019 The MathWorks, Inc.
%
%
%