%% Filter Multiple Disturbance 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 quarter and the last four quarters of government consumption expenditures and investment (GCE). Pass multiple multivariate Gaussian disturbance paths through 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.

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

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

rgcelag4 = lagmatrix(Data.GCE,0:4);
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);

%
% Generate 1000 paths of numobs observations from a 3-D Gaussian distribution.  numobs is the number of observations in the data without any missing values.

numpaths = 1000;
numseries = Mdl.NumSeries;
idx = all(~ismissing([Data array2table(rgcelag4)]),2);
numobs = sum(idx);
rng(1);
Z = mvnrnd(zeros(Mdl.NumSeries,1),eye(Mdl.NumSeries),numobs*numpaths);
Z = reshape(Z,[numobs,3,numpaths]);

%
% Filter the disturbances through the estimated model. Supply the predictor data. Return the innovations (scaled disturbances).

[Y,E] = filter(EstMdl,Z,'X',rgcelag4);

% Y and E are 244-by-3-by-1000 matrices of filtered responses and scaled disturbances, respectively. The columns correspond to the CPI growth rate, unemployment rate, and GDP growth rate, respectively.  filter applies the same predictor data to all paths.
%
% For each time point, compute the mean vector of the filtered responses among all paths.

MeanFilt = mean(Y,3);

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

Data = Data(idx,:);

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,MeanFilt(:,j));
    hold off
end

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

% Copyright 2012 The MathWorks, Inc.