%% Simulate Multiple Response Paths
% Estimate a VAR(4) model of the consumer price index (CPI), the unemployment rate, and the gross domestic product (GDP). Include a linear regression component containing the current and the last 4 quarters of government consumption expenditures and investment. Simulate multiple paths from the estimated model.
% Load the Data_USEconModel data set. Compute the real GDP.

load Data_USEconModel
DataTable.RGDP = DataTable.GDP./DataTable.GDPDEF*100;

%
% Plot all variables on separate plots.

figure;
subplot(2,2,1)
plot(DataTable.Time,DataTable.CPIAUCSL);
ylabel('Index');
title('Consumer Price Index');
subplot(2,2,2)
plot(DataTable.Time,DataTable.UNRATE);
ylabel('Percent');
title('Unemployment Rate');
subplot(2,2,3)
plot(DataTable.Time,DataTable.RGDP);
ylabel('Output');
title('Real Gross Domestic Product');
subplot(2,2,4)
plot(DataTable.Time,DataTable.GCE);
ylabel('Billions of $');
title('Government Expenditures');

%
% Stabilize the CPI, GDP, and GCE by converting each to a series of growth rates. Synchronize the unemployment rate series with the others by removing its first observation. Create a new data set containing the transformed variables.

inputVariables = {'CPIAUCSL' 'RGDP' 'GCE'};
Data = varfun(@price2ret,DataTable,'InputVariables',inputVariables);
Data.Properties.VariableNames = inputVariables;
Data.UNRATE = DataTable.UNRATE(2:end);
dates = DataTable.Time(2:end);

%
% Expand the GCE rate series to a matrix that includes its current value and up through four lagged values. Remove GCE and observations containing any missing values from Data.

rgcelag4 = lagmatrix(Data.GCE,0:4);
idx = all(~ismissing([Data table(rgcelag4)]),2);
Data = Data(idx,:);
Data.GCE = [];

%
% Create a default VAR(4) model using the shorthand syntax.

Mdl = varm(3,4);
Mdl.SeriesNames = ["rcpi" "unrate" "rgdpg"];

%
% Estimate the model using the entire sample. Specify the GCE matrix as data for the regression component.

EstMdl = estimate(Mdl,Data.Variables,'X',rgcelag4);

%
% Simulate 1000 paths from the estimated model. Specify that the length of the paths is the same as the length of the data without any missing values. Supply the predictor data. Return the innovations (scaled disturbances).

numpaths = 1000;
numobs = size(Data,1);
rng(1); % For reproducibility
[Y,E] = simulate(EstMdl,numobs,'X',rgcelag4,'NumPaths',numpaths);

% Y is a 244-by-3-by-1000 matrix of simulated responses. E is a matrix whose dimensions correspond to the dimensions of Y, but represents the simulated, scaled disturbances. The columns correspond to the CPI growth rate, unemployment rate, and GDP growth rate, respectively.  simulate applies the same predictor data to all paths.
%
% For each time point, compute the mean vector of the simulated responses among all paths.

MeanSim = mean(Y,3);

% MeanSim is a 244-by-3 matrix containing the average of the simulated responses at each time point.
%
% Plot the simulated responses, their averages, and the data.

figure;
for j = 1:Mdl.NumSeries
    subplot(2,2,j)
    plot(Data.Time,squeeze(Y(:,j,:)),'Color',[0.8,0.8,0.8])
    title(Mdl.SeriesNames{j});
    hold on
    h1 = plot(Data.Time,Data{:,j});
    h2 = plot(Data.Time,MeanSim(:,j));
    hold off
end

hl = legend([h1 h2],'Data','Mean');
hl.Position = [0.6 0.25 hl.Position(3:4)];

% Copyright 2012 The MathWorks, Inc.