%% 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 - 1}}}\\
% {{y_{2,t - 1}}}\\
% {{y_{3,t - 1}}}
% \end{array}} \right] + \left[\begin{array}{*{20}{c}}\varepsilon_{1t}\\
% \varepsilon_{2t}\\ \varepsilon_{3t}\end{array}\right]
%% Preprocess the 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.
%
% Flag all periods corresponding to missing values in the model variables.

idx = all(~ismissing(DataTable(:,{'UNRATE' 'M1SL'})),2);

% For the rest of the example, consider only those values that of the series indicated by a true in idx.
%
% 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(idx)./(DataTable.GDPDEF(idx)/100);
rM1SL = DataTable.M1SL(idx)./(DataTable.GDPDEF(idx)/100);

dLRGDP = diff(log(rGDP));              % rGDP growth rate
dLRM1SL = diff(log(rM1SL));            % rM1SL growth rate
d3MTB = diff(DataTable.TB3MS(idx));    % Change in short-term interest rate (3MTB)
dUNRATE = diff(DataTable.UNRATE(idx)); % 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 = dates((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 VARX(4) model using varm.

Mdl = varm(n,4);

% Mdl is a varm model object serving as a template for estimation. Currently, Mdl does know have the structure in place for the regression component. However, MATLAB® creates the required structure during estimation.
%
%% Estimate the VAR(4) Model
% Estimate the parameters of the VARX(4) model using estimate. Display the parameter estimates.

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

% EstMdl is a varm model object containing the estimated parameters.
%
%% Simulate Out-Of-Sample Response Paths
% 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 pages 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);
    h1 = plot(dates(oosIdx),squeeze(YSim(:,j,:)),'Color',0.75*ones(3,1));
    hold on;   
    h2 = plot(dates(oosIdx),YSimBar(:,j),'.-k','LineWidth',2);
    h3 = plot(dates(oosIdx),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}));
    datetick;
    axis tight;
    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 1000, 3 year paths 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),...
    'NumPaths',numPaths);

% XSim is a 12-by-1000 numeric matrix of simulated exogenous paths. The rows correspond to periods and the columns correspond to paths.
%
% Simulate 3 years of 1000 future response series paths from the estimated model using the simulated exogenous data.  simulate does not accept multiple paths of predictor data, so you must simulate responses in a loop. Since the model contains 4 lags per endogenous variable, specify the last 4 observations in the estimation sample as presample data.

YSimRX = zeros(oosT,n,numPaths); % Preallocate

for j = 1:numPaths
    YSimRX(:,:,j) = simulate(EstMdl,oosT,'X',XSim(:,j),'Y0',Y0);
end

% 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);
    h1 = plot(dates(oosIdx),squeeze(YSimRX(:,j,:)),'Color',0.75*ones(3,1));
    hold on;   
    h2 = plot(dates(oosIdx),YSimBarRX(:,j),'.-k','LineWidth',2);
    h3 = plot(dates(oosIdx),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 with Simulated Unemployment Rate',respNames{j}));
    datetick;
    axis tight;
    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.