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. Note that the following example uses some functions found in the Financial Toolbox™, and therefore to run the example the Financial Toolbox must be installed. In addition, although the Datafeed Toolbox™ is not required to run this example, access to it allows the inclusion of the most recent economic data available from the Federal Reserve Economic Database. Additional Requirements: ContentsDescription of the ModelThe SmetsWouters 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" empiricallymotivated regression models, the DSGE approach focuses on "small" theoreticallyderived 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 SmetsWouters framework), unemployment tracks a widelyperceived 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 ReserveIf you have the Datafeed Toolbox, 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.matwhich contains time series for the period from 31Mar1947 to 31Mar2009. 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 Tenyear treasury bond yield
% HOANBS Nonfarm 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 Threemonth treasury bill yield
% UNRATE Unemployment rate
Since the data from FRED have different periodicities and date ranges, we use the time series features found in timetable 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 startofperiod 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 "trycatch" and load Data_USEconModel.mat if problems occur
try
Universe = timetable;
% Open a Datafeed Toolbox connection to FRED
c = fred;
for i = 1:numel(series)
fprintf('Started loading %s ... ',series{i});
% Fetch data from FRED
FredData = fetch(c, series{i});
% Set up dates
dates = FredData.Data(:,1);
dates = datetime(dates, 'ConvertFrom', 'datenum');
% Set up data
Data = FredData.Data(:,2);
% Create financial time series
fts = timetable(dates, Data, 'VariableNames', series(i));
if strcmpi(strtrim(FredData.Frequency), 'Quarterly')
fprintf('Quarterly ... ');
elseif strncmp('Mar', strtrim(FredData.Frequency), 3)
fprintf('Quarterly ... ');
else
fprintf('Monthly ... ');
end
% Combine time series
Universe = synchronize(Universe, fts, 'union');
fprintf('Finished loading %s ...\n',series{i});
end
close(c);
Universe.Properties.Description = 'U.S. Quarterly Macroeconomic Data';
% Convert to quarterly periodicity
DataTable = retime(Universe, 'quarterly', 'lastvalue');
dates = DataTable.Time;
% Dates are startofperiod dates so move to endofperiod date
dates = dates + calquarters(1)  caldays(1);
dates = lbusdate(dates.Year, dates.Month,[], [], 'datetime'); % Last business date of month
DataTable.Time = dates;
DataTable.Time.Format = 'QQQyy';
% Trim date range to period from 1947 to present
StartDate = datetime('31Mar1947', 'InputFormat', 'ddMMyyyy');
EndDate = Universe.Time(end);
TR = timerange(StartDate, EndDate, 'closed');
DataTable = DataTable(TR, :);
catch E
% Case with no internet connection
fprintf('Unable to connect to FRED. Will use local data.\n');
fprintf('Loading data from Data_USEconModel.mat ...\n');
load Data_USEconModel
end
else
% Case with no Datafeed Toolbox
fprintf('Loading data from Data_USEconModel.mat ...\n');
load Data_USEconModel
end
Loading time series data from St. Louis Federal Reserve (FRED) ...
Unable to connect to FRED. Will use local data.
Loading data from Data_USEconModel.mat ...
Business Cycle Dates from National Bureau of Economic ResearchTo 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 as follows: Start Date End Date
 
15May1937 15Jun1938
15Feb1945 15Oct1945
15Nov1948 15Oct1949
15Jul1953 15May1954
15Aug1957 15Apr1958
15Apr1960 15Feb1961
15Dec1969 15Nov1970
15Nov1973 15Mar1975
15Jan1980 15Jul1980
15Jul1981 15Nov1982
15Jul1990 15Mar1991
15Mar2001 15Nov2001
15Dec2007 15Jun2009 Transform Raw Data into Time Series for the ModelTo 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
DataTable = rmmissing(DataTable, 1);
Data = DataTable.Variables;
dates = datenum(DataTable.Time);
% Log series
CONS = log(DataTable.PCEC);
CPI = log(DataTable.CPIAUCSL);
DEF = log(DataTable.GDPDEF);
GCE = log(DataTable.GCE);
GDP = log(DataTable.GDP);
HOURS = log(DataTable.HOANBS);
INV = log(DataTable.GPDI);
WAGES = log(DataTable.COE);
% Interest rates (annual)
rFED = 0.01*(DataTable.FEDFUNDS);
rG10 = 0.01*(DataTable.GS10);
rTB3 = 0.01*(DataTable.TB3MS);
% Integrated rates
TB3 = ret2tick(0.25*rTB3);
TB3 = log(TB3(2:end));
% Unemployment rate
rUNEMP = 0.01*(DataTable.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) ];
rWAGES = [ 4*mean(diff(WAGES(1:5))); 4*diff(WAGES) ];
Display Raw DataTo 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('Investment and Output');
h = legend('GDP','INV','Location','Best');
h.FontSize = 7;
h.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('Inflation and GDP Deflator');
h = legend('CPI','DEF','Location','Best');
h.FontSize = 7;
h.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('Wages and Hours');
h = legend('WAGES','HOURS','Location','Best');
h.FontSize = 7;
h.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('Consumption');
h = legend('CONS','GCE','Location','Best');
h.FontSize = 7;
h.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('Interest Rates');
h = legend('FED','G10','TB3','Location','Best');
h.FontSize = 7;
h.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('Unemployment');
h = legend('UNEMP','Location','Best');
h.FontSize = 7;
h.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 ModelThe 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'};
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 31Mar1959 to 31Mar2009.
Main ModelThe main model is an unrestricted VAR model in the form with and Some differences between the data for our model and the SmetsWouters model should be noted. First, we use nominal series throughout since the GDP deflator is part of our model. Second, we use the 3month 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 OrderThe 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,~,logL] = estimate(varm(n,i), Y1, 'Y0', Y0);
summary = summarize(Spec);
AICtest(i) = aicbic(logL,summary.NumEstimatedParameters,summary.SampleSize);
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('Optimal Lag Order with Akaike Information Criterion');
xlabel('Lag Order');
ylabel('AIC');
hold off
AIC(1) = 8510.02
AIC(2) = 8531.85
AIC(3) = 8507.67
AIC(4) = 8496
AIC(5) = 8479.53
AIC(6) = 8459.65
AIC(7) = 8409.47
Optimal lag for model is 2.
Backtest to Assess Forecast Accuracy of the ModelThe next step is to determine the forecast accuracy of our model. To do this, we perform a MonteCarlo 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 meansquare 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
NumPaths = 500;
[T, n] = size(Y);
rng default
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
Y1 = Y(iStart:iEnd,:);
iY1 = iY(iStart:iEnd,:);
% Set up model and estimate coefficients
Spec = estimate(varm(n,nAR), Y1, 'Y0', zeros(nAR,n));
% Do forecasts
iFY = simulate(Spec, Horizon, 'Y0', Y1, 'NumPaths', 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(yysyy+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(yysyy+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
31Dec1975 0.77 0.81 2.35 2.20 0.08 0.49 16.50 0.30
31Dec1976 1.95 0.85 2.75 3.15 0.68 0.97 10.61 0.59
30Dec1977 2.05 0.81 1.38 2.62 0.49 0.74 5.48 0.29
29Dec1978 0.58 0.42 0.87 2.23 1.13 0.83 3.83 0.26
31Dec1979 1.37 0.77 1.06 1.31 0.83 1.46 6.44 1.38
31Dec1980 1.11 2.75 2.29 0.82 0.46 2.33 13.87 0.53
31Dec1981 3.27 1.78 2.50 1.72 0.65 0.49 18.31 0.89
31Dec1982 1.54 0.30 1.10 2.47 2.57 1.41 7.38 0.26
30Dec1983 2.07 0.61 1.67 1.22 1.06 0.42 7.94 0.25
31Dec1984 0.81 0.45 0.30 0.72 0.12 1.16 3.61 0.04
31Dec1985 1.12 0.48 1.35 1.38 0.16 0.63 7.64 0.30
31Dec1986 0.92 0.12 0.59 0.47 0.24 0.71 2.80 0.14
31Dec1987 0.59 0.53 0.18 0.77 0.54 1.05 3.55 0.21
30Dec1988 0.31 0.21 1.68 0.35 0.33 0.42 1.50 0.55
29Dec1989 0.91 0.67 0.71 1.36 0.39 0.61 4.13 0.35
31Dec1990 1.37 0.38 1.35 1.76 0.46 2.09 3.50 0.20
31Dec1991 0.57 0.22 0.63 1.31 0.29 0.34 2.44 0.20
31Dec1992 3.14 0.43 2.40 0.71 0.41 1.85 9.20 0.11
31Dec1993 0.87 0.38 1.12 0.41 0.52 0.80 2.15 0.26
30Dec1994 1.85 0.12 1.24 0.33 0.13 1.20 6.25 0.26
29Dec1995 0.44 0.14 0.41 0.61 0.09 0.25 4.51 0.45
31Dec1996 0.22 0.32 0.27 0.72 0.24 0.45 4.49 0.37
31Dec1997 0.32 0.81 0.68 0.42 0.10 0.18 2.41 0.30
31Dec1998 0.44 0.47 0.50 0.56 0.12 0.54 1.55 0.28
31Dec1999 0.67 0.33 0.94 0.70 0.62 0.70 2.92 0.24
29Dec2000 1.06 0.40 2.42 2.15 1.00 0.67 7.73 0.16
31Dec2001 2.30 0.06 1.64 1.98 0.21 1.85 6.51 0.17
31Dec2002 1.12 0.39 0.81 2.43 0.74 0.81 6.62 0.10
31Dec2003 0.86 0.86 1.17 1.03 0.38 0.27 2.37 0.02
31Dec2004 1.20 0.17 1.64 0.99 0.03 1.04 4.64 0.22
30Dec2005 1.23 0.67 1.28 0.42 0.19 1.20 5.53 0.49
29Dec2006 0.49 0.32 1.46 0.49 0.65 0.41 2.75 0.12
31Dec2007 3.01 0.42 2.88 3.03 1.56 3.30 11.10 0.46
Assess Forecast AccuracyThe 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 onesided. 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(['Forecast Accuracy for ' sprintf('%g',Horizon/4) 'Year Horizon']);
end
h = legend(YSeries{i},'Location','best');
h.FontSize = 7;
h.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 2008Let'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 red 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
Y1 = Y(iStart:iEnd,:); % estimation data
% Set up model and estimate coefficients
Spec = estimate(varm(n,nAR), Y1, 'Y0', zeros(nAR,n));
% Do forecasts
FT = 20;
FD = Example_QuarterlyDates(dates(iEnd), FT);
[FY, FYCov] = forecast(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(endHw+1:end),Y1(endHw+1:end,i));
hold on
scatter(dates(iEndHw+1:iEnd+Ow),Y(iEndHw+1:iEnd+Ow,i),'r.');
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(['Model Calibration to ' sprintf('%s',datestr(dates(iEnd),1))]);
end
h = legend(YSeries{i},'Location','best');
h.FontSize = 7;
h.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 DateWe 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:iStart1,:); % presample data
else
Y0 = [];
end
Y1 = Y(iStart:iEnd,:); % estimation data
% Set up model and estimate coefficients
Spec = estimate(varm(n,nAR), Y1, 'Y0', zeros(nAR,n));
% Do forecasts
FT = 20;
FD = Example_QuarterlyDates(dates(iEnd), FT);
[FY, FYCov] = forecast(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(endHw+1:end),Y1(endHw+1:end,i));
hold on
scatter(dates(iEndHw+1:iEnd+Ow),Y(iEndHw+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(['Model Calibration to ' sprintf('%s',datestr(dates(iEnd),1))]);
end
h = legend(YSeries{i},'Location','best');
h.FontSize = 7;
h.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 AnalysisFor 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;
clf;
YX = armairf(Spec.AR, {}, 'InnovCov', diag(diag(Spec.Covariance)), 'NumObs', FT)*100;
ii = 0;
for i = 1:n
for j = 1:n
ii = ii + 1;
subplot(n,n,ii,'align');
plot(YX(:,j,i));
if i == 1
title(['\bf ' Responses{j}]);
end
if j == 1
ylabel(['\bf ' Impulses{i}]);
end
ax = gca;
ax.FontSize = 7.5;
if i == n
ax.XTickMode = 'auto';
else
ax.XTick = [];
end
end
end
Real GDP ForecastThe 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 longterm 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;
rng default
iFY = simulate(Spec, FT, 'Y0', Y1, 'NumPaths', 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(['Real GDP Forecast Based on Data to End of ' sprintf('%s',datestr(dates(iEnd),1))]);
dateaxis('x',12);
legend('Forecast','LongTerm 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. ConclusionThis 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 SmetsWouters 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 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. 123162.
 FRED, St. Louis Federal Reserve, Federal Reserve Economic Database, http://research.stlouisfed.org/fred2/.
 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. 12411277.
 H. Lutkepohl (2006). New Introduction to Multiple Time Series Analysis, Springer.
 H. Lutkepohl and M. Kratzig (2004). Applied Time Series Econometrics, Cambridge University Press.
 NBER, National Bureau of Economic Research, Business Cycle Expansions and Contractions, http://www.nber.org/cycles/cyclesmain.html.
 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. 11231175.
 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. 161183.
 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. 586606.
