MATLAB Examples

# 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.

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.

• Financial Toolbox

## 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, 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 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 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 = 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 start-of-period dates so move to end-of-period 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 = 'QQQ-yy'; % Trim date range to period from 1947 to present StartDate = datetime('31-Mar-1947', 'InputFormat', 'dd-MM-yyyy'); 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 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 as follows:

Start Date     End Date
-----------   -----------
15-May-1937   15-Jun-1938
15-Feb-1945   15-Oct-1945
15-Nov-1948   15-Oct-1949
15-Jul-1953   15-May-1954
15-Aug-1957   15-Apr-1958
15-Apr-1960   15-Feb-1961
15-Dec-1969   15-Nov-1970
15-Nov-1973   15-Mar-1975
15-Jan-1980   15-Jul-1980
15-Jul-1981   15-Nov-1982
15-Jul-1990   15-Mar-1991
15-Mar-2001   15-Nov-2001
15-Dec-2007   15-Jun-2009

## 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 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 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('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 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'}; 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,~,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 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 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(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(['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 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 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(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),'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 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 = 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(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(['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 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; 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 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; 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','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.