%% Forecast US GDP Rate
% Consider a two-state Markov-switching dynamic regression model of the postwar US real GDP growth rate. The model has the parameter estimates presented in [1].
%
%% Create Markov-Switching Dynamic Regression Model
% Create a fully specified discrete-time Markov chain model that describes the regime switching mechanism. Label the regimes.

P = [0.92 0.08; 0.26 0.74];
mc = dtmc(P,'StateNames',["Expansion" "Recession"]);

% Create separate, fully specified AR(0) models for the two regimes.

sigma = 3.34; % Homoscedastic models across states
mdl1 = arima('Constant',4.62,'Variance',sigma^2);
mdl2 = arima('Constant',-0.48,'Variance',sigma^2);
mdl = [mdl1 mdl2];

% Create the Markov-switching dynamic regression model from the switching mechanism mc and the state-specific submodels mdl.

Mdl = msVAR(mc,mdl);

% Mdl is a fully specified msVAR object.
%
%% Load and Preprocess Data
% forecast requires observations to initialize the model. Load the US GDP data set.

load Data_GDP

% Data contains quarterly measurements of the US real GDP in the period 1947:Q1–2005:Q2. The period of interest in [1] is 1947:Q2–2004:Q2. For more details on the data set, enter Description at the command line.
% Transform the data to an annualized rate series by:
% Converting the data to a quarterly rate within the estimation period
% Annualizing the quarterly rates

qrate = diff(Data(2:230))./Data(2:229); % Quarterly rate
arate = 100*((1 + qrate).^4 - 1);       % Annualized rate

% The transformation drops the first observation.
%
%% Forecast US GDP Rates
% Forecast the model over a 12-quarter forecast horizon. Initialize the model by supplying the annualized rate series.

numPeriods = 12;
yf = forecast(Mdl,arate,numPeriods);

% yf is a 12-by-1 vector of model forecasts. yf(j) is the j-step-ahead optimal point forecast.
% Plot the observed annualized GDP rate from 1980 with the model forecasts.

dates = datetime(dates(3:230),'ConvertFrom','datenum',...
    'Format','yyyy:QQQ','Locale','en_US');
dt1980Q1 = datetime("1980:Q1",'InputFormat','yyyy:QQQ',...
    'Locale','en_US'); % Specify US date format for 1980:Q1. 
idx = dates >= dt1980Q1;
figure;
plot(dates(idx),arate(idx),'k',...
    dates(end) + calquarters(1:numPeriods),yf,'r--')
xlabel("Year")
ylabel("GDP (Annualized Rate)")
recessionplot
legend("Observations","Forecasts")

% Copyright 2019 The MathWorks, Inc.
%
%