%% Simulate Responses of Estimated VARX Model
% This example shows how to estimate a multivariate time series model that contains lagged endogenous and exogenous variables, and how to simulate responses. The response series are the quarterly:
% Changes in real gross domestic product (rGDP) rates (y_{1t})
% Real money supply (rM1SL) rates (y_{2t})
% Short-term interest rates (i.e., three-month treasury bill yield, y_{3t})
% from March, 1959 through March, 2009 . The exogenous series is the quarterly changes in the unemployment rates (x_t).
% Suppose that a model for the responses is this VARX(4,3) model
% \left[ {\begin{array}{*{20}{c}}
% {{y_{1t}}}\\
% {{y_{2t}}}\\
% {{y_{3t}}}
% \end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
% {{a_1}}\\
% {{a_2}}\\
% {{a_3}}
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% {{x_t}}&0&0\\
% 0&{{x_t}}&0\\
% 0&0&{{x_t}}
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{b_1}}\\
% {{b_2}}\\
% {{b_3}}
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% {\phi _{11}^{(1)}}&{\phi _{12}^{(1)}}&{\phi _{13}^{(1)}}\\
% {\phi _{21}^{(1)}}&{\phi _{22}^{(1)}}&{\phi _{23}^{(1)}}\\
% {\phi _{31}^{(1)}}&{\phi _{32}^{(1)}}&{\phi _{33}^{(1)}}
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{y_{1,t - 1}}}\\
% {{y_{2,t - 1}}}\\
% {{y_{3,t - 1}}}
% \end{array}} \right]+...+  \left[ {\begin{array}{*{20}{c}}
% {\phi _{11}^{(4)}}&{\phi _{12}^{(4)}}&{\phi _{13}^{(4)}}\\
% {\phi _{21}^{(4)}}&{\phi _{22}^{(4)}}&{\phi _{23}^{(4)}}\\
% {\phi _{31}^{(4)}}&{\phi _{32}^{(4)}}&{\phi _{33}^{(4)}}
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{y_{1,t - 4}}}\\
% {{y_{2,t - 4}}}\\
% {{y_{3,t - 4}}}
% \end{array}} \right] + \left[\begin{array}{*{20}{c}}\varepsilon_{1t}\\
% \varepsilon_{2t}\\ \varepsilon_{3t}\end{array}\right]
%% Preprocess Data
% Load the U.S. macroeconomic data set. Flag the series and their periods that contain missing values (indicated by NaN values).

load Data_USEconModel
varNaN = any(ismissing(DataTable),1); % Variables containing NaN values
seriesWithNaNs = series(varNaN)

% In this data set, the variables that contain missing values entered the sample later than the other variables. There are no missing values after sampling started for a particular variable.
% estimate accommodates missing values for responses and predictors by removing corresponding rows from the data.
%
% Compute rGDP and rM1SL, and the growth rates of rGDP, rM1SL, short-term interest rates, and the unemployment rate. Description contains a description of the data and the variable names. Reserve the last three years of data to investigate the out-of-sample performance of the estimated model.

rGDP = DataTable.GDP./(DataTable.GDPDEF/100);
rM1SL = DataTable.M1SL./(DataTable.GDPDEF/100);

dLRGDP = diff(log(rGDP));         % rGDP growth rate
dLRM1SL = diff(log(rM1SL));       % rM1SL growth rate
d3MTB = diff(DataTable.TB3MS);    % Change in short-term interest rate (3MTB)
dUNRATE = diff(DataTable.UNRATE); % Change in unemployment rate

T = numel(d3MTB);    % Total sample size
oosT = 12;           % Out-of-sample size
estT = T - oosT;     % Estimation sample size
estIdx = 1:estT;     % Estimation sample indices
oosIdx = (T - 11):T; % Out-of-sample indices
dates = DataTable.Time((end - T + 1):end);

EstY = [dLRGDP(estIdx) dLRM1SL(estIdx) d3MTB(estIdx)]; % In-sample responses
estX = dUNRATE(estIdx);                                % In-sample exogenous data
n = size(EstY,2);

OOSY = [dLRGDP(oosIdx) dLRM1SL(oosIdx) d3MTB(oosIdx)]; % Out-of-sample responses
oosX = dUNRATE(oosIdx);                                % Out-of-sample exogenous data

%
%% Create VARX Model
% Create a VAR model that characterizes the VARX(4) model using varm.

Mdl = varm(n,4);

% Mdl is a varm model object that represents a template specifying the structure of the VAR(4) model.
%
%% Estimate VAR Model
% Estimate the parameters of the VARX(4) model using estimate. Display a summary of the estimation.

EstMdl = estimate(Mdl,EstY,'X',estX);
summarize(EstMdl)

% EstMdl is an fully specified, estimated varm model object.
%
%% Simulate Out-Of-Sample Response Paths Using Same Exogenous Data per Path
% Simulate 1000, 3 year response series paths from the estimated model assuming that the exogenous unemployment rate is a fixed series. Since the model contains 4 lags per endogenous variable, specify the last 4 observations in the estimation sample as presample data.

numPaths = 1000;
Y0 = EstY((end-3):end,:);
rng(1); % For reproducibility
YSim = simulate(EstMdl,oosT,'X',oosX,'Y0',Y0,'NumPaths',numPaths);

% YSim is a 12-by-3-by-1000 numeric array of simulated responses. The rows of YSim correspond to out-of-sample periods, the columns correspond to the response series, and the leaves correspond to paths.
%
% Plot the response data and the simulated responses. Identify the 5%, 25%, 75% and 95% percentiles, and the mean and median of the simulated series at each out-of-sample period.

YSimBar = mean(YSim,3);
YSimQrtl = quantile(YSim,[0.05 0.25 0.5 0.75 0.95],3);
RepDates = repmat(dates(oosIdx),1,1000);
respNames = {'dLRGDP' 'dLRM1SL' 'd3MTB'};

figure;
for j = 1:n
    subplot(3,1,j);
    Data = array2timetable(squeeze(YSim(:,j,:)),'RowTimes',dates(oosIdx));
    h1 = plot(Data.Time,Data{:,:},'Color',0.75*ones(3,1));
    hold on;   
    h2 = plot(Data.Time,YSimBar(:,j),'.-k','LineWidth',2);
    h3 = plot(Data.Time,squeeze(YSimQrtl(:,j,:)),':r','LineWidth',1.5);
    h4 = plot(dates((end - 30):end),[EstY((end - 18):end,j); OOSY(:,j)],...
        'b','LineWidth',2);
    title(sprintf('%s',respNames{j}));
    hold off;
end
legend([h1(1) h2(1) h3(1) h4],{'Simulated Series','Simulation Mean',...
    'Simulation Quartile','Data'},'Location',[0.4 0.1 0.01 0.01],...
    'FontSize',8);

%
%% Simulate Out-Of-Sample Response Paths Using Random Exogenous Data
% Suppose that the change in the unemployment rate is an AR(4) model, and fit the model to the estimation sample data.

MdlUNRATE = arima('ARLags',1:4);
EstMdlUNRATE = estimate(MdlUNRATE,estX,'Display','off');

% EstMdlUNRATE is an arima model object containing the parameter estimates.
%
% Simulate 1, 3 year path from the estimated AR(4) model for the change in unemployment rate. Since the model contains 4 lags, specify the last 4 observations in the estimation sample as presample data.

xSim = simulate(EstMdlUNRATE,oosT,'Y0',estX(end-3:end));

% xSim is a 12-by-1 numeric vector of a simulated exogenous path.
%
% Simulate 3 years of 1000 future response series paths from the estimated model using the simulated exogenous data. Since the model contains 4 lags per endogenous variable, specify the last 4 observations in the estimation sample as presample data.

YSimRX = simulate(EstMdl,oosT,'X',xSim,'Y0',Y0,'NumPaths',numPaths);

% YSimRX is a 12-by-3-by-1000 numeric array of simulated responses.
%
% Plot the response data and the simulated responses. Identify the 5%, 25%, 75% and 95% percentiles, and the mean and median of the simulated series at each out-of-sample period.

YSimBarRX = mean(YSimRX,3);
YSimQrtlRX = quantile(YSimRX,[0.05 0.25 0.5 0.75 0.95],3);

figure;
for j = 1:n
    subplot(3,1,j);
    Data = array2timetable(squeeze(YSimRX(:,j,:)),'RowTimes',dates(oosIdx));
    h1 = plot(Data.Time,Data{:,:},'Color',0.75*ones(3,1));
    hold on;   
    h2 = plot(Data.Time,YSimBarRX(:,j),'.-k','LineWidth',2);
    h3 = plot(Data.Time,squeeze(YSimQrtlRX(:,j,:)),':r','LineWidth',1.5);
    h4 = plot(dates((end - 30):end),[EstY((end - 18):end,j); OOSY(:,j)],...
        'b','LineWidth',2);
    title(sprintf('%s',respNames{j}));
    hold off;
end
legend([h1(1) h2(1) h3(1) h4],{'Simulated Series','Simulation Mean',...
    'Simulation Quartile','Data'},'Location',[0.4 0.1 0.01 0.01],...
    'FontSize',8)

% Copyright 2012 The MathWorks, Inc.