Accelerating the pace of engineering and science

# Econometrics Toolbox

## Modeling the United States Economy

Economists and policymakers are concerned with understanding the dynamics of economies, especially during periods of significant macroeconomic shocks. Although numerous approaches are possible, we will develop a small macroeconomic model in the style of Smets and Wouters.

Our descriptive macroeconomic model offers a nice starting point to examine the impact of various shocks on the United States economy, particularly around the period of the 2008 fiscal crisis. We will use the multiple time series tools from the Econometrics Toolbox™ to gain some insight.

Description of the Model

The Smets-Wouters model (2002, 2004, 2007) is a nonlinear system of equations in the form of a Dynamic Stochastic General Equilibrium (DSGE) model that seeks to characterize an economy derived from economic first principles. The basic model works with 7 time series: output, prices, wages, hours worked, interest rates, consumption, and investment.

Whereas a common approach in macroeconomics has been to create "large" empirically-motivated regression models, the DSGE approach focuses on "small" theoretically-derived models. It is this combination of normative rigor and parsimony that is one of the main attractions of the DSGE approach.

Armed with a small model of this sort, the linearized form can be cast as a VAR model that can be handled with standard methods of multiple time series analysis. It is an unrestricted form of this VAR model that we will examine subsequently.

For illustrative purposes, we will add an eighth time series: unemployment. Although this is not a part of the basic model (and is actually superfluous within the Smets-Wouters framework), unemployment tracks a widely-perceived measure of the "health" of an economy.

Whereas the original model was developed to model both country and aggregate European economies, we will apply the model to the United States economy as in Smets and Wouters (2007).

Obtain Economic Data from St. Louis Federal Reserve

If you have the Datafeed Toolbox™ and are online, this example will download data from the St. Louis Federal Reserve Economic Database (see link to FRED in the references) so that the analysis will incorporate the most recent available data. If not, this example will load data from the file Data_USEconModel.mat which contains time series for the period from 31-Mar-1947 to 31-Mar-2009. The series are listed in the following table.

% FRED Series     Description
% -----------     ----------------------------------------------------------------
% COE             Paid compensation of employees in $billions % CPIAUCSL Consumer price index % FEDFUNDS Effective federal funds rate % GCE Government consumption expenditures and investment in$ billions
% GDP             Gross domestic product in $billions % GDPDEF Gross domestic product price deflator % GPDI Gross private domestic investment in$ billions
% GS10            Ten-year treasury bond yield
% HOANBS          Non-farm business sector index of hours worked
% M1SL            M1 money supply (narrow money)
% M2SL            M2 money supply (broad money)
% PCEC            Personal consumption expenditures in \$ billions
% TB3MS           Three-month treasury bill yield
% UNRATE          Unemployment rate


Since the data from FRED have different periodicities and date ranges, we use the financial time series tools to build a universe of our time series at a quarterly periodicity.

if license('test', 'datafeed_toolbox')

% Load data from FRED and convert to quarterly periodicity
% Note that dates are start-of-period dates in the FRED database

fprintf('Loading time series data from St. Louis Federal Reserve (FRED) ...\n');

% FRED time series to be used for our analysis

series = { 'COE', 'CPIAUCSL', 'FEDFUNDS', 'GCE', 'GDP', 'GDPDEF', 'GPDI', ...
'GS10', 'HOANBS', 'M1SL', 'M2SL', 'PCEC', 'TB3MS', 'UNRATE' };

% Obtain data with "try-catch" and load Data_USEconModel.mat if problems occur

try
Universe = [];

% Open a Datafeed Toolbox connection to FRED

c = fred;

for i = 1:numel(series)

% Fetch data from FRED

FredData = fetch(c, series{i});

% Dates are start-of-period dates so move to end-of-period date

offset = 1;
if strcmpi(strtrim(FredData.Frequency),'Quarterly')
offset = 2;
elseif strncmp('Mar', strtrim(FredData.Frequency), 3)
offset = 2;
else
offset = 0;
end

% Set up dates

dates = FredData.Data(:,1);

mm = month(dates) + offset;
yy = year(dates);
for t = 1:numel(dates)
if mm(t) > 12
mm(t) = mm(t) - 12;
yy(t) = yy(t) + 1;
end
end
dates = lbusdate(yy, mm);

% Set up data

Data = FredData.Data(:,2);

% Create financial time series

fts = fints(dates, Data, series{i});

% Convert to quarterly periodicity

if strcmpi(strtrim(FredData.Frequency), 'Quarterly')
fprintf('Quarterly ... ');
elseif strncmp('Mar', strtrim(FredData.Frequency), 3)
fprintf('Quarterly ... ');
else
fprintf('Monthly ... ');
fts = toquarterly(fts);
end

% Combine time series

Universe = merge(fts, Universe);

end

close(c);

Universe.desc = 'U.S. Macroeconomic Data';
Universe.freq = 'quarterly';

% Trim date range to period from 1947 to present

StartDate = datenum('31-Mar-1947');
EndDate = datenum(Universe.dates(end));

Universe = Universe([datestr(StartDate,1) '::' datestr(EndDate,1)]);

% Convert combined time series into date and data arrays

dates = Universe.dates;
Data = fts2mat(Universe.(series));
Dataset = dataset([{Data},series],'ObsNames',cellstr(datestr(dates,'QQ-YY')));

% Uncomment next line to save data in Data_USEconModelUpdate.mat

% save Data_USEconModelUpdate series dates Data Dataset

catch E

% Case with no internet connection

fprintf('Unable to connect to FRED. Will use local data.\n');

end

else

% Case with no Datafeed Toolbox

end

Loading time series data from St. Louis Federal Reserve (FRED) ...
Unable to connect to FRED. Will use local data.


Business Cycle Dates from National Bureau of Economic Research

To examine the interplay between the business cycle and our model, we also include the dates for peaks and troughs of the economic cycle from the National Bureau of Economic Research (see link to NBER in the references). We arbitrarily set the middle of the listed month as the start or end date of a recession.

Recessions = [ datenum('15-May-1937'), datenum('15-Jun-1938');
datenum('15-Feb-1945'), datenum('15-Oct-1945');
datenum('15-Nov-1948'), datenum('15-Oct-1949');
datenum('15-Jul-1953'), datenum('15-May-1954');
datenum('15-Aug-1957'), datenum('15-Apr-1958');
datenum('15-Apr-1960'), datenum('15-Feb-1961');
datenum('15-Dec-1969'), datenum('15-Nov-1970');
datenum('15-Nov-1973'), datenum('15-Mar-1975');
datenum('15-Jan-1980'), datenum('15-Jul-1980');
datenum('15-Jul-1981'), datenum('15-Nov-1982');
datenum('15-Jul-1990'), datenum('15-Mar-1991');
datenum('15-Mar-2001'), datenum('15-Nov-2001');
datenum('15-Dec-2007'), datenum('15-Jun-2009') ];

Recessions = busdate(Recessions);


Transform Raw Data into Time Series for the Model

To work with our time series, we create two sets of time series. The first set contains differenced or rate data for each of our time series and the second set contains integrated or cumulative data. Time series that have exponential growth are also transformed into logarithmic series prior to differencing.

% Remove dates with NaN values

ii = any(isnan(Data),2);
dates(ii) = [];
Data(ii,:) = [];
Dataset(ii,:) = [];

% Log series

CONS = log(Dataset.PCEC);
CPI = log(Dataset.CPIAUCSL);
DEF = log(Dataset.GDPDEF);
GCE = log(Dataset.GCE);
GDP = log(Dataset.GDP);
HOURS = log(Dataset.HOANBS);
INV = log(Dataset.GPDI);
M1 = log(Dataset.M1SL);
M2 = log(Dataset.M2SL);
WAGES = log(Dataset.COE);

% Interest rates (annual)

rFED = 0.01*(Dataset.FEDFUNDS);
rG10 = 0.01*(Dataset.GS10);
rTB3 = 0.01*(Dataset.TB3MS);

% Integrated rates

FED = ret2tick(0.25*rFED);
FED = log(FED(2:end));
G10 = ret2tick(0.25*rG10);
G10 = log(G10(2:end));
TB3 = ret2tick(0.25*rTB3);
TB3 = log(TB3(2:end));

% Unemployment rate

rUNEMP = 0.01*(Dataset.UNRATE);

UNEMP = ret2tick(0.25*rUNEMP);
UNEMP = log(UNEMP(2:end));

% Annualized rates

rCONS = [ 4*mean(diff(CONS(1:5))); 4*diff(CONS) ];
rCPI = [ 4*mean(diff(CPI(1:5))); 4*diff(CPI) ];
rDEF = [ 4*mean(diff(DEF(1:5))); 4*diff(DEF) ];
rGCE = [ 4*mean(diff(GCE(1:5))); 4*diff(GCE) ];
rGDP = [ 4*mean(diff(GDP(1:5))); 4*diff(GDP) ];
rHOURS = [ 4*mean(diff(HOURS(1:5))); 4*diff(HOURS) ];
rINV = [ 4*mean(diff(INV(1:5))); 4*diff(INV) ];
rM1 = [ 4*mean(diff(M1(1:5))); 4*diff(M1) ];
rM2 = [ 4*mean(diff(M2(1:5))); 4*diff(M2) ];
rWAGES = [ 4*mean(diff(WAGES(1:5))); 4*diff(WAGES) ];


Display Raw Data

To see what our time series look like, we plot each of the differenced time series (identified with a lowercase r preceding the series mnemonic) and overlay shaded bands that identify periods of economic recession as determined by NBER.

clf;

subplot(3,2,1,'align');
plot(dates, [rGDP, rINV]);
recessionplot;
dateaxis('x');
title('\bfInvestment and Output');
h = legend('GDP','INV','Location','Best');
set(h,'FontSize',7,'Box','off');
axis([dates(1) - 600, dates(end) + 600, 0, 1]);
axis 'auto y'

subplot(3,2,2,'align');
plot(dates, [rCPI, rDEF]);
recessionplot;
dateaxis('x');
title('\bfInflation and GDP Deflator');
h = legend('CPI','DEF','Location','Best');
set(h,'FontSize',7,'Box','off');
axis([dates(1) - 600, dates(end) + 600, 0, 1]);
axis 'auto y'

subplot(3,2,3,'align');
plot(dates, [rWAGES, rHOURS]);
recessionplot;
dateaxis('x');
title('\bfWages and Hours');
h = legend('WAGES','HOURS','Location','Best');
set(h,'FontSize',7,'Box','off');
axis([dates(1) - 600, dates(end) + 600, 0, 1]);
axis 'auto y'

subplot(3,2,4,'align');
plot(dates, [rCONS, rGCE]);
recessionplot;
dateaxis('x');
title('\bfConsumption');
h = legend('CONS','GCE','Location','Best');
set(h,'FontSize',7,'Box','off');
axis([dates(1) - 600, dates(end) + 600, 0, 1]);
axis 'auto y'

subplot(3,2,5,'align');
plot(dates, [rFED, rG10, rTB3]);
recessionplot;
dateaxis('x');
title('\bfInterest Rates');
h = legend('FED','G10','TB3','Location','Best');
set(h,'FontSize',7,'Box','off');
axis([dates(1) - 600, dates(end) + 600, 0, 1]);
axis 'auto y'

subplot(3,2,6,'align');
plot(dates, rUNEMP);
recessionplot;
dateaxis('x');
title('\bfUnemployment');
h = legend('UNEMP','Location','Best');
set(h,'FontSize',7,'Box','off');
axis([dates(1) - 600, dates(end) + 600, 0, 1]);
axis 'auto y'


The shaded bands to identify recessions are plotted using the utility function recessionplot.

Set up the Main Model

The main model for our analysis uses the seven time series described in Smets and Wouters (2007) plus an appended eighth time series. These time series are listed in the following table with their relationship to raw FRED counterparts. The variable Y contains the main time series for the model and the variable iY contains integrated data from Y that will be used in forecasting analyses.

%   Model       FRED Series     Transformation from FRED Data to Model Time Series
%   --------    -----------     --------------------------------------------------
%   rGDP        GDP             rGDP = diff(log(GDP))
%   rDEF        GDPDEF          rDEF = diff(log(GDPDEF))
%   rWAGES      COE             rWAGE = diff(log(COE))
%   rHOURS      HOANBS          rWORK = diff(log(WORK))
%   rTB3        TB3MS           rTB3 = 0.01*TB3MS
%   rCONS       PCEC            rCONS = diff(log(PCEC))
%   rINV        GPDI            rINV = diff(log(GPDI))
%   rUNEMP      UNRATE          rUNEMP = 0.01*UNRATE

Y = [rGDP, rDEF, rWAGES, rHOURS, rTB3, rCONS, rINV, rUNEMP];
iY = [GDP, DEF, WAGES, HOURS, TB3, CONS, INV, UNEMP];

YSeries = {'Output (GDP)', 'Prices', 'Total Wages', 'Hours Worked', ...
'Cash Rate', 'Consumption', 'Private Investment', 'Unemployment'};
YAbbrev = {'GDP', 'DEF', 'WAGES', 'HOURS', 'TB3', 'CONS', 'INV', 'UNEMP'};

YInfo = 'U.S. Macroeconomic Model';

n = numel(YSeries);

fprintf('The date range for available data is %s to %s.\n', ...
datestr(dates(1),1),datestr(dates(end),1));

The date range for available data is 31-Mar-1959 to 31-Mar-2009.


Main Model

The main model is an unrestricted VAR model in the form

with

and

Some differences between the data for our model and the Smets-Wouters model should be noted. First, we use nominal series throughout since the GDP deflator is part of our model. Second, we use the 3-month Treasury Bill rate in place of the Federal Funds rate due to greater coverage. Third, we use the change in hours worked rather than the integrated series. Fourth, of course, we have added unemployment. Finally, we do not detrend the data with either series trends or a common GDP trend.

Optimal Lag Order

The first step in our analysis is to determine the "optimal" number of autoregressive lags based on the Akaike Information Criterion (AIC). We set up models with up to 7 lags (7 quarters) and perform the AIC test using the toolbox function aicbic. The minimum AIC test statistic identifies the optimal number of lags. Depending upon the source, the theory would suggest that 3 lags are sufficient and practice dictates between 2 and 4 lags. We will use 2 lags subsequently.

nARmax = 7;

Y0 = Y(1:nARmax,:);
Y1 = Y(nARmax+1:end,:);

AICtest = zeros(nARmax,1);
for i = 1:nARmax
Spec = vgxset('n', n, 'Constant', true, 'nAR', i, 'Series', YSeries);
[Spec, SpecStd, LLF] = vgxvarx(Spec, Y1, [], Y0);
AICtest(i) = aicbic(LLF,Spec.NumActive,Spec.T);
fprintf('AIC(%d) = %g\n',i,AICtest(i));
end
[AICmin, nAR] = min(AICtest);

fprintf('Optimal lag for model is %d.\n',nAR);

clf;

plot(AICtest);
hold on
scatter(nAR,AICmin,'filled','b');
title('\bfOptimal Lag Order with Akaike Information Criterion');
xlabel('Lag Order');
ylabel('AIC');
hold off

AIC(1) = -8438.02
AIC(2) = -8459.85
AIC(3) = -8435.67
AIC(4) = -8424
AIC(5) = -8407.53
AIC(6) = -8387.65
AIC(7) = -8337.47
Optimal lag for model is 2.


Backtest to Assess Forecast Accuracy of the Model

The next step is to determine the forecast accuracy of our model. To do this, we perform a Monte-Carlo simulation with 500 sample paths for each year from 1975 to the most recent prior year. Given 500 sample paths for each year, we estimate the root mean-square error (RMSE) between subsequent realizations and forecasts over the time horizon. For this analysis, the forecast horizon is 1 year.

The RMSE forecasts work with integrated simulated forecast data to compute forecast accuracy because integrated forecasts provides a better measure of where the model is going than to work with differenced data.

The results appear in the following table. Each row contains results for the end date of the estimation period which is also the start date for the forecast period. Following the date, each row contains the forecast RMSE for each series over the forecast horizon.

nAR = 2;

syy = 1975;                 % start year for backtest
eyy = 2008;                 % end year for backtest

Horizon = 4;                % number of quarters for forecast horizon

[T, n] = size(Y);

FError = NaN(eyy - syy + 1, n);
FDates = zeros(eyy - syy + 1, 1);

fprintf('RMSE of Actual vs Model Forecast (x 100) with Horizon of %d Quarters\n',Horizon);
fprintf('%12s','ForecastDate');
for i = 1:n
fprintf('  %7s',YAbbrev{i});
end
fprintf('\n');

for yy = syy:eyy

StartDate = lbusdate(1959,3);
EndDate = lbusdate(yy,12);

if StartDate < dates(1)
error(message('econ:Demo_USEconModel:EarlyStartDate', datestr( dates( 1 ), 1 )));
end
if EndDate > dates(end)
error(message('econ:Demo_USEconModel:LateStartDate', datestr( dates( end ), 1 )));
end

% Locate indexes in data for specified date range

iStart = find(StartDate <= dates,1);
if iStart < 1
iStart = 1;
end
iEnd = find(EndDate <= dates,1);
if iEnd > numel(dates)
iEnd = numel(dates);
end

if iStart > 1
Y0 = Y(1:iStart-1,:);
else
Y0 = [];
end
Y1 = Y(iStart:iEnd,:);
iY1 = iY(iStart:iEnd,:);

% Set up model and estimate coefficients

Spec = vgxset('n', n, 'Constant', true, 'nAR', nAR, 'Series', YSeries, 'Info', YInfo);
Spec = vgxvarx(Spec, Y1, [], Y0);

% Do forecasts

NumPaths = 500;
iFY = vgxsim(Spec, Horizon, [], Y1, [], NumPaths);
iFY = repmat(iY1(end,:),[Horizon,1,NumPaths]) + 0.25*cumsum(iFY);
eFY = mean(iFY,3);

% Assess Forecast Quality

Ow = max(0,min(Horizon,(size(Y,1) - iEnd)));		% overlap between actual and forecast data

if Ow >= Horizon
h = Horizon;
else
h = [];
end

FDates(yy-syy+1) = lbusdate(yy,12);
if ~isempty(h)
Yerr = iY(iEnd+1:iEnd+Ow,:) - eFY(1:Ow,:);

Ym2 = Yerr(1:h,:) .^ 2;
Yrmse = sqrt(mean(Ym2,1));

fprintf('%12s',datestr(EndDate,1));
for i = 1:n
fprintf('  %7.2f',100*Yrmse(i));
end
FError(yy-syy+1,:) = 100*Yrmse';
fprintf('\n');
end
end

RMSE of Actual vs Model Forecast (x 100) with Horizon of 4 Quarters
ForecastDate      GDP      DEF    WAGES    HOURS      TB3     CONS      INV    UNEMP
31-Dec-1975     0.77     0.81     2.35     2.20     0.08     0.49    16.50     0.30
31-Dec-1976     1.95     0.85     2.75     3.15     0.68     0.97    10.61     0.59
30-Dec-1977     2.05     0.81     1.38     2.62     0.49     0.74     5.48     0.29
29-Dec-1978     0.58     0.42     0.87     2.23     1.13     0.83     3.83     0.26
31-Dec-1979     1.37     0.77     1.06     1.31     0.83     1.46     6.44     1.38
31-Dec-1980     1.11     2.75     2.29     0.82     0.46     2.33    13.87     0.53
31-Dec-1981     3.27     1.78     2.50     1.72     0.65     0.49    18.31     0.89
31-Dec-1982     1.54     0.30     1.10     2.47     2.57     1.41     7.38     0.26
30-Dec-1983     2.07     0.61     1.67     1.22     1.06     0.42     7.94     0.25
31-Dec-1984     0.81     0.45     0.30     0.72     0.12     1.16     3.61     0.04
31-Dec-1985     1.12     0.48     1.35     1.38     0.16     0.63     7.64     0.30
31-Dec-1986     0.92     0.12     0.59     0.47     0.24     0.71     2.80     0.14
31-Dec-1987     0.59     0.53     0.18     0.77     0.54     1.05     3.55     0.21
30-Dec-1988     0.31     0.21     1.68     0.35     0.33     0.42     1.50     0.55
29-Dec-1989     0.91     0.67     0.71     1.36     0.39     0.61     4.13     0.35
31-Dec-1990     1.37     0.38     1.35     1.76     0.46     2.09     3.50     0.20
31-Dec-1991     0.57     0.22     0.63     1.31     0.29     0.34     2.44     0.20
31-Dec-1992     3.14     0.43     2.40     0.71     0.41     1.85     9.20     0.11
31-Dec-1993     0.87     0.38     1.12     0.41     0.52     0.80     2.15     0.26
30-Dec-1994     1.85     0.12     1.24     0.33     0.13     1.20     6.25     0.26
29-Dec-1995     0.44     0.14     0.41     0.61     0.09     0.25     4.51     0.45
31-Dec-1996     0.22     0.32     0.27     0.72     0.24     0.45     4.49     0.37
31-Dec-1997     0.32     0.81     0.68     0.42     0.10     0.18     2.41     0.30
31-Dec-1998     0.44     0.47     0.50     0.56     0.12     0.54     1.55     0.28
31-Dec-1999     0.67     0.33     0.94     0.70     0.62     0.70     2.92     0.24
29-Dec-2000     1.06     0.40     2.42     2.15     1.00     0.67     7.73     0.16
31-Dec-2001     2.30     0.06     1.64     1.98     0.21     1.85     6.51     0.17
31-Dec-2002     1.12     0.39     0.81     2.43     0.74     0.81     6.62     0.10
31-Dec-2003     0.86     0.86     1.17     1.03     0.38     0.27     2.37     0.02
31-Dec-2004     1.20     0.17     1.64     0.99     0.03     1.04     4.64     0.22
30-Dec-2005     1.23     0.67     1.28     0.42     0.19     1.20     5.53     0.49
29-Dec-2006     0.49     0.32     1.46     0.49     0.65     0.41     2.75     0.12
31-Dec-2007     3.01     0.42     2.88     3.03     1.56     3.30    11.10     0.46


Assess Forecast Accuracy

The forecast errors are visualized in the following plot. On each subplot, the blue line plots the average of the RMSE forecast errors associated with each date along with the sample mean (green line) and standard deviation (dotted red lines) of these errors over all dates. A value of 1 on the plot corresponds with a 1% forecast error.

Note that the standard deviation of forecast errors is somewhat misleading since forecast errors are one-sided. Nonetheless, the standard deviation offers a qualitative guide to the variability of forecast errors.

mFError = NaN(size(FError));
sFError = NaN(size(FError));
for i = 1:n
mFError(:,i) = nanmean(FError(:,i));
sFError(:,i) = nanstd(FError(:,i));
end

for i = 1:n
subplot(ceil(n/2),2,i,'align');
plot(FDates,FError(:,i));
hold on
plot(FDates,mFError(:,i),'g');
plot(FDates,[mFError(:,i) - sFError(:,i),mFError(:,i) + sFError(:,i)],':r');
recessionplot;
dateaxis('x',12);
if i == 1
title(['\bfForecast Accuracy for ' sprintf('%g',Horizon/4) '-Year Horizon']);
end
h = legend(YSeries{i},'Location','best');
set(h,'FontSize',7,'Box','off');
hold off
end


With the exception of private investment, all forecasts tend to fall within plus or minus 2% of their realized values over a 1 year time horizon. However, around recessions, the errors tend to be larger - which implies either unmodeled or mismodeled effects that have not been captured in our linearized model.

Analysis to the End of 2008

Let's look at our forecasts in greater detail. We calibrate the model at the end of 2008 to see how the fiscal meltdown of the prior few months might play out. The model is our VAR(2) model with calibration over the period from 1959 to 2008 and with forecasts 5 years into the future.

The results are plotted below, where actual data are plotted with green dots, forecasts are identified with both a blue line for the mean forecast and red dotted lines for the standard deviations of the forecast.

nAR = 2;

% Set up date range

StartDate = lbusdate(1959,3);
EndDate = lbusdate(2008,12);

if StartDate < dates(1)
error(message('econ:Demo_USEconModel:EarlyStartDate', datestr( dates( 1 ), 1 )));
end
if EndDate > dates(end)
error(message('econ:Demo_USEconModel:LateStartDate', datestr( dates( end ), 1 )));
end

% Locate indexes in data for specified date range

iStart = find(StartDate <= dates,1);
if iStart < 1
iStart = 1;
end
iEnd = find(EndDate <= dates,1);
if iEnd > numel(dates)
iEnd = numel(dates);
end

% Set up data for estimation

D1 = dates(iStart:iEnd,:);			% dates for specified date range
if iStart > 1
Y0 = Y(1:iStart-1,:);			% presample data
else
Y0 = [];
end
Y1 = Y(iStart:iEnd,:);				% estimation data

% Set up model and estimate coefficients

Spec = vgxset('n', n, 'Constant', true, 'nAR', nAR, 'Series', YSeries, 'Info', YInfo);
Spec = vgxvarx(Spec, Y1, [], Y0);

% Do forecasts

FT = 20;
FD = Example_QuarterlyDates(dates(iEnd), FT);

[FY, FYCov] = vgxpred(Spec, FT, [], Y1);
FYSigma = zeros(size(FY));
for t = 1:FT
FYSigma(t,:) = sqrt(diag(FYCov{t}))';
end

Hw = 20;                                    % number of historical quarters to display
Fw = 20;                                    % number of forecast quarters to display
Ow = max(0,min(Fw,(size(Y,1) - iEnd)));     % overlap between historical and forecast data

clf;

for i = 1:n
subplot(ceil(n/2),2,i,'align');
plot(D1(end-Hw+1:end),Y1(end-Hw+1:end,i));
hold on
scatter(dates(iEnd-Hw+1:iEnd+Ow),Y(iEnd-Hw+1:iEnd+Ow,i),'.');
plot([D1(end); FD],[Y1(end,i); FY(:,i)],'b');
plot(FD,[FY(:,i) - FYSigma(:,i), FY(:,i) + FYSigma(:,i)],':r');
dateaxis('x',12);
if i == 1
title(['\bfModel Calibration to ' sprintf('%s',datestr(dates(iEnd),1))]);
end
h = legend(YSeries{i},'Location','best');
set(h,'FontSize',7,'Box','off');
hold off
end


This plot suggests that, from the end of 2008 onward, a recovery is likely to begin sometime in 2009 - but with two distinct results. At the level of the macro economy, the model suggests that output will increase but with an increase in both interest rates and inflation. At the household level, however, the recovery might take longer which is most evident in the gradual reduction of unemployment.

This section of code can be run with different start and end dates. Specifically, if you run the model to the end of 2006 (change EndDate to lbusdate(2006,12)), you can see hints of an upcoming downturn with a projected small dip in real GDP (GDP net the GDP deflator), a small drop in hours worked and a small rise in unemployment. Thus, at the end of 2006, our model predicted a slowdown or mild recession.

To set up the forecast dates for the plot, we use the following helper function which is also used in the next analysis.

type Example_QuarterlyDates.m

function qDates = Example_QuarterlyDates(startDate,numDates)
%EXAMPLE_QUARTERLYDATES Compute quarterly forecast dates
%
% Syntax:
%
%   qDates = Example_QuarterlyDates(startDate,numDates)
%
% Description:
%
%   Compute a specified number of quarterly dates after a start date.
%
% Input Arguments:
%
%   startDate - Starting date number.
%   numDates - Number of quarterly dates to compute.
%
% Output Argument:
%
%   qDates - numDates quarterly date numbers, starting with the first
%     quarter after startDate.

% Copyright 2010 The MathWorks, Inc.

qDates = zeros(numDates,1);
qq = month(startDate)/3;
yy = year(startDate);

for i = 1:numDates

qq = qq + 1;
if qq > 4
qq = 1;
yy = yy + 1;
end
qDates(i) = lbusdate(yy,3*qq);

end


Analysis to the Current Available Date

We now repeat our previous analysis with more recent data. With downloaded data from FRED, we can run our analysis on the most current available data. Otherwise, we have data to the end of March 2009.

nAR = 2;

% Set up date range

StartDate = lbusdate(1959,3);
EndDate = dates(end);

if StartDate < dates(1)
error(message('econ:Demo_USEconModel:EarlyStartDate', datestr( dates( 1 ), 1 )));
end
if EndDate > dates(end)
error(message('econ:Demo_USEconModel:LateStartDate', datestr( dates( end ), 1 )));
end

% Locate indexes in data for specified date range

iStart = find(StartDate <= dates,1);
if iStart < 1
iStart = 1;
end
iEnd = find(EndDate <= dates,1);
if iEnd > numel(dates)
iEnd = numel(dates);
end

% Set up data for estimation

D1 = dates(iStart:iEnd,:);			% dates for specified date range
if iStart > 1
Y0 = Y(1:iStart-1,:);			% presample data
else
Y0 = [];
end
Y1 = Y(iStart:iEnd,:);				% estimation data

% Set up model and estimate coefficients

Spec = vgxset('n', n, 'Constant', true, 'nAR', nAR, 'Series', YSeries, 'Info', YInfo);
Spec = vgxvarx(Spec, Y1, [], Y0);

% Do forecasts

FT = 20;
FD = Example_QuarterlyDates(dates(iEnd), FT);

[FY, FYCov] = vgxpred(Spec, FT, [], Y1);
FYSigma = zeros(size(FY));
for t = 1:FT
FYSigma(t,:) = sqrt(diag(FYCov{t}))';
end

Hw = 20;                                    % number of historical quarters to display
Fw = 20;                                    % number of forecast quarters to display
Ow = max(0,min(Fw,(size(Y,1) - iEnd)));     % overlap between historical and forecast data

clf;

for i = 1:n
subplot(ceil(n/2),2,i,'align');
plot(D1(end-Hw+1:end),Y1(end-Hw+1:end,i));
hold on
scatter(dates(iEnd-Hw+1:iEnd+Ow),Y(iEnd-Hw+1:iEnd+Ow,i),'.');
plot([D1(end); FD],[Y1(end,i); FY(:,i)],'b');
plot(FD,[FY(:,i) - FYSigma(:,i), FY(:,i) + FYSigma(:,i)],':r');
dateaxis('x',12);
if i == 1
title(['\bfModel Calibration to ' sprintf('%s',datestr(dates(iEnd),1))]);
end
h = legend(YSeries{i},'Location','best');
set(h,'FontSize',7,'Box','off');
hold off
end


How well does the model perform given this new information? Some things to consider are changes between the prior analysis to the end of 2008 and the current analysis. The backtest suggests that the model is reasonably accurate for a period of about 1 year, which would imply that the realized values ought to fall within the 1 standard deviation error bands over the next 4 quarters.

Impulse Response Analysis

For policymakers, a primary question is - what is likely to happen if a particular change occurs, especially if it is a change due to a policy decision? An impulse response analysis provides a ceteris paribus sensitivity analysis of the dynamics of a system. The following plot shows the projected dynamic responses of each time series along each column in reaction to a 1 standard deviation impulse along each row. The units for each plot are percentage deviations from the initial state for each time series.

Impulses = YAbbrev;
Responses = YAbbrev;

W0 = zeros(FT, n);

clf;

ii = 0;
for i = 1:n
WX = W0;
WX(1,i) = sqrt(Spec.Q(i,i));
YX = 100*(vgxproc(Spec, WX, [], Y1) - vgxproc(Spec, W0, [], Y1));
for j = 1:n
ii = ii + 1;
subplot(n,n,ii,'align');
plot(YX(:,j));
if i == 1
title(['\bf ' Responses{j}]);
end
if j == 1
ylabel(['\bf ' Impulses{i}]);
end
if i == n
set(gca,'xtickmode','auto');
else
set(gca,'xtick',[]);
end
end
end


Real GDP Forecast

The final analysis is a forecast of real GDP based on the calibrated model to the current available date. The projected value is compared with a long-term trend value based on the past 30 years of real GDP data.

iY1 = iY(iStart:iEnd,:);

% Simulate forecasts of cumulative values of model time series

NumPaths = 1000;
iFY = vgxsim(Spec, FT, [], Y1, [], NumPaths);
iFY = repmat(iY1(end,:),[FT,1,NumPaths]) + 0.25*cumsum(iFY);
iFY = iFY(:,1,:) - iFY(:,2,:);
FGDP = mean(iFY,3);
FGDPSigma = std(iFY,0,3);
FGDP0 = GDP(iEnd) - DEF(iEnd);

w = 120;
H = [ ones(w,1) (1:w)' ];
trendParam = H \ (GDP(iEnd - w + 1:iEnd) - DEF(iEnd - w + 1:iEnd));
trendFGDP = [ ones(FT,1) w + (1:FT)' ] * trendParam;

clf;

plot(FD, [FGDP, trendFGDP]);
hold on
plot(FD, [FGDP - FGDPSigma, FGDP + FGDPSigma],':r');
title(['\bfReal GDP Forecast Based on Data to End of ' sprintf('%s',datestr(dates(iEnd),1))]);
dateaxis('x',12);
legend('Forecast','Long-Term Trend','Location','Best');
grid on;
ylabel('Log Real GDP');


If the forecast curve starts below and then approaches the trend line, this implies an expansion - and recovery - since GDP growth would have to exceed trend growth to return to the trend line. This result is the case for forecasts from late 2008 or early 2009, which would suggest a strong recovery during 2010.

Conclusion

This example shows a few things you can do with the multiple time series analysis tools in the Econometrics Toolbox. We have shown that a VAR model based loosely on the Smets-Wouters model provides reasonably accurate forecasts for most of the economic series in our model. Thus, the scripts in this example are a good point of departure to begin testing your own ideas in macroeconomic modeling and analysis.

References

1. M. Del Negro, F. Schorfheide, F. Smets, and R. Wouters (2007), "On the Fit of New Keynesian Models," Journal of Business & Economic Statistics, Vol. 25, No. 2, pp. 123-162.

2. FRED, St. Louis Federal Reserve, Federal Reserve Economic Database, http://research.stlouisfed.org/fred2/.

3. M. Kimball (1995), "The Quantitative Analytics of the Basic Neomonetarist Model", Journal of Money, Credit, and Banking, Part 2: Liquidity, Monetary Policy, and Financial Intermediation, Vol. 27, No. 4, pp. 1241-1277.

4. H. Lutkepohl (2006). New Introduction to Multiple Time Series Analysis, Springer.

5. H. Lutkepohl and M. Kratzig (2004). Applied Time Series Econometrics, Cambridge University Press.

6. NBER, National Bureau of Economic Research, Business Cycle Expansions and Contractions, http://www.nber.org/cycles/cyclesmain.html.

7. F. Smets and R. Wouters (2002), An Estimated Stochastic Dynamic General Equilibrium Model of the Euro Area, European Central Bank, Working Paper Series, No. 171. Also in Journal of the European Economic Association, Vol. 1, No. 5, 2003, pp. 1123-1175.

8. F. Smets and R. Wouters (2004), Comparing Shocks and Frictions in US and Euro Area Business Cycles: A Bayesian DSGE Approach, European Central Bank, Working Paper Series, No. 391. Also in Journal of Applied Econometrics, Vol. 20, No. 1, 2005, pp. 161-183.

9. F. Smets and R. Wouters (2007), Shocks and Frictions in US Business Cycles: A Bayesian DSGE Approach, European Central Bank, Working Paper Series, No. 722. Also in American Economic Review, Vol. 97, No. 3, 2007, pp. 586-606.