Fit a VAR(4) model to the consumer price index (CPI) and unemployment rate data. Then, simulate responses by filtering a random series of Gaussian distributed disturbances through the estimated model.

Load the Data_USEconModel data set.

load Data_USEconModel

Plot the two series on separate plots.

figure;
plot(DataTable.Time,DataTable.CPIAUCSL);
title('Consumer Price Index');
ylabel('Index');
xlabel('Date');

figure;
plot(DataTable.Time,DataTable.UNRATE);
title('Unemployment rate');
ylabel('Percent');
xlabel('Date');

Stabilize the CPI by converting it to a series of growth rates. Synchronize the two series by removing the first observation from the unemployment rate series. Create a new data set containing the transformed variables, and do not include any rows containing at least one missing observation.

rcpi = price2ret(DataTable.CPIAUCSL);
unrate = DataTable.UNRATE(2:end);
idx = all(~ismissing([rcpi unrate]),2);
Data = array2timetable([rcpi(idx) unrate(idx)],...
    'RowTimes',DataTable.Time(idx),'VariableNames',{'rcpi','unrate'});

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

Mdl = varm(2,4);

Estimate the model using the entire data set.

EstMdl = estimate(Mdl,Data.Variables);

EstMdl is a fully specified, estimated varm model object.

Generate a numobs-by-2 series of random Gaussian distributed values, where numobs is the number of observations in the data.

numobs = size(Data,1);
rng(1) % For reproducibility
Z = mvnrnd(zeros(Mdl.NumSeries,1),eye(Mdl.NumSeries),numobs);

To simulate responses, filter the disturbances through the estimated model.

Y = filter(EstMdl,Z);

Y is a 245-by-2 matrix of simulated responses. The first and second columns contain the simulated CPI growth rate and unemployment rate, respectively.

Plot the simulated and true responses.

figure;
plot(Data.Time,Y(:,1));
hold on;
plot(Data.Time,Data.rcpi)
ylabel('Growth rate')
xlabel('Date')
title('CPI Growth Rate');
legend('Simulation','True')

figure;
plot(Data.Time,Y(:,2));
hold on;
plot(Data.Time,Data.unrate)
ylabel('Percent');
xlabel('Date')
title('Unemployment Rate');
legend('Simulation','True')

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)];