Structural Vector Autoregression demo
This script illustrates how to work with Structural Vector Autoregression (SVAR) models with a simple example.
Contents
- Model
- Data
- Estimation
- Impulse Responses
- Historical decomposition
- Forecast Error Variance decomposition
Model
The model is a structural VAR with three variables (Inflation, Output and Federal Funds Rate) and four lags.
where variables are ranked in the following order
and with the following imposed restriction on the structural matrix A
Data
In this example, we estimate the model parameters on US data.
- Consumer Price Index for All Urban Consumers: All Items (CPIAUCSL)
- Federal Funds Rate (FEDFUNDS)
- Civilian Noninstitutional Population (CNP16OV)
- Gross Domestic Product (GDP)
- Gross Domestic Product: Implicit Price Deflator (GDPDEF)
if contains(matlabpath,['toolbox' filesep 'datafeed'])
% Fetches data from FRED and saves it in raw_data.mat
% Create connection with FRED server
c = fred('https://research.stlouisfed.org/fred2/');
start_date = '01/01/1960';
rawdata.cpi = fetch(c, 'CPIAUCSL' ,start_date,end_date);
rawdata.ffr = fetch(c, 'FEDFUNDS' ,start_date,end_date);
rawdata.pop = fetch(c, 'CNP16OV' ,start_date,end_date);
rawdata.gdp = fetch(c, 'GDP' ,start_date,end_date);
rawdata.gdp_defl = fetch(c, 'GDPDEF' ,start_date,end_date);
save raw_data rawdata start_date end_date
% Close connection with server
disp('Successfully fetched data!')
disp('Unable to download data from FRED!')
disp('Could not find the datafeed toolbox.')
end
Successfully fetched data!
If the toolbox is not available in your MATLAB license or fetching data fails, it loads previously dowloaded data (raw_data.mat)
% Load raw data from .mat file
disp('Loading previously downloaded data ...')
Data preparation
Downloaded data has to be further manipulated before it can be used in our model.
In particular, we have to:
- Convert monthly data (CPI, FEDFUNDS and CNP16OV) into quarterly data
- Create a series for Inflation from the CPI index
- Construct Real GDP per capita
- Remove the trend component to get cyclical fluctuations of GDP
Monthly data to quarterly data
quarters = rawdata.gdp.Data(:,1);
% Quarters to months legend:
% To see dates, datestr(quarters)
% Monthly to quarterly data
rawdata.cpi.Data = rawdata.cpi.Data(ismember(rawdata.cpi.Data(:,1),quarters),:);
rawdata.ffr.Data = rawdata.ffr.Data(ismember(rawdata.ffr.Data(:,1),quarters),:);
rawdata.pop.Data = rawdata.pop.Data(ismember(rawdata.pop.Data(:,1),quarters),:);
Year-on-year inflation rate
Pi(t-4) = 100*(rawdata.cpi.Data(t,2) - rawdata.cpi.Data(t-4,2))/rawdata.cpi.Data(t-4,2);
% NOTE: The first four observations are lost when constructing Pi.
% Update sample period and data
quarters = quarters(5:end);
rawdata.gdp.Data = rawdata.gdp.Data(ismember(rawdata.gdp.Data(:,1),quarters),:);
rawdata.gdp_defl.Data = rawdata.gdp_defl.Data(ismember(rawdata.gdp_defl.Data(:,1),quarters),:);
rawdata.ffr.Data = rawdata.ffr.Data(ismember(rawdata.ffr.Data(:,1),quarters),:);
rawdata.pop.Data = rawdata.pop.Data(ismember(rawdata.pop.Data(:,1),quarters),:);
Real GDP per capita
% Create population index
index_date = '01-Jan-1990';
pop_index = rawdata.pop.Data(:,2) ./ rawdata.pop.Data(rawdata.pop.Data(:,1)==datenum(index_date),2);
gdp_real = rawdata.gdp.Data(:,2) ./ rawdata.gdp_defl.Data(:,2);
% Divide by population index
gdp_real_cap = gdp_real ./ pop_index;
Remove trend in output
% Choose detrending option: 'log-linear' or 'HP filtering'
detrend_opt = 'HP filtering';
Y = detrend(100.*log(gdp_real_cap));
[~,Y] = hpfilter(100.*log(gdp_real_cap),1600);
Cleaned data
dates = datetime(datestr(rawdata.ffr.Data(:,1)));
R = rawdata.ffr.Data(:,2);
% Final data ready to be used in estimation
data = timetable(dates,Pi,Y,R,'VariableNames',{'Inflation','Output','FedFunds'});
Plot data (optional)
% Select 'true' or 'false' if you want to plot data
plot(dates,data.Variables)
legend(data.Properties.VariableNames)
Estimation
Now that data is ready, we can estimate the VAR using the VAR Toolbox.
First, we need to tell MATLAB where the toolbox is located:
% Add VAR toolbox to the MATLAB path (including all subfolders)
addpath(genpath('VAR_toolbox'))
Next, we need to specify some variables that determine the shape of our VAR, that is the number of lags and whether the VAR should have a constant, a trend or both.
Mean and trends
Before estimating the VAR, you need to choose between the following options:
- const_trend = 0 --> No constant and no trend
- const_trend = 1 --> Only constant and no trend
- const_trend = 2 --> Constant and trend
This choice depends on whether your endogenous variables (Inflation, Output and Interest rate) display a non-zero mean and a trend. Therefore, you should inspect your variables and make a choice accordingly.
In practice, it is recommended to include both constant and trend in your estimation. If it turns out that your time series does not have a non-zero mean nor a trend, the respective estimated parameters will be equal (or very close) to zero.
% Choose constant and trend
Finally, call the appropriate function to estimate the model and show parameters estimates:
[VAR, VAR_options] = VARmodel(X,nbr_lags,const_trend);
% Show estimated model parameters
VAR_options.vnames = data.Properties.VariableNames;
VARprint(VAR,VAR_options);
Reduced form VAR estimation:
Inflation Output FedFunds
c 0.3964 0.2305 0.5102
trend -0.0011 -0.0007 -0.0026
Inflation(-1) 0.9528 -0.0588 -0.0015
Output(-1) 0.2775 1.0243 0.6967
FedFunds(-1) 0.0881 -0.0527 0.6162
Inflation(-2) -0.0356 0.0757 0.0148
Output(-2) -0.1874 0.0413 -0.2558
FedFunds(-2) 0.0713 -0.0060 0.1492
Inflation(-3) 0.1342 -0.0674 0.1320
Output(-3) -0.0946 -0.2336 -0.1915
FedFunds(-3) -0.0026 0.1522 0.2315
Inflation(-4) -0.1431 0.0373 -0.0364
Output(-4) 0.0965 -0.0337 0.0312
FedFunds(-4) -0.1396 -0.1135 -0.1190
Impulse responses
In this section, we compute and plot the dynamic responses of endogenous variables to a structural shock.
Identification
But wait, what is a structural shock!? Remember in the slides we made an important distinction between innovations and structural shocks. Loosely speaking, a structural shock is a sudden movement in the error terms that we can clearly identify coming from a specific source. Conversely, innovations are sudden movements in the error terms which might be caused by several factors.
How do we tell a structural shock from an innovation? One way is to use economic theory to impose some restrictions on the parameters of the model. Given the restriction, we can now say that a sudden movement in the error term is due to a movement in a particular variable.
Remember when we introduced the model, we said that the structural matrix A was lower triangular. This is an example of one of such restrictions. This restriction implies that only some variables are able to contemporaneously influence the other variables. This identification scheme is called recursive or short-run restriction.
In particular, given the ordering of the variables and the restriction on A,
we deduce that:
- The policy rate () responds to contemporaneous changes to other variables
- Output () responds current inflation but not interest rate
- Inflation () is not contemporaneously affected by or and only responds to these variables with a lag
This ordering of the variables is reasonable since we expect the central bank to immediately respond to movements in inflation and output with changes in the policy rate.
On the other hand, since both inflation and output are assumed to not contemporaneously react to movements in the policy rate, we can say that any sudden change in the error term of is due to movements in only. That is, it is a shock to the policy rate that does not come from a reaction to movements in other variables. This is a monetary policy shock. Implementation
In practice, the restriction on the structural matrix A is imposed by doing a Cholesky decomposition of the estimated variance-covariance matrix of the residuals. % Choose the identification scheme
VAR_options.ident = 'oir'; % 'oir' selects a recursive scheme
% Choose the horizon for the impulse responses
% Apply the identification scheme and compute impulse responses
[IRF,VAR] = VARir(VAR,VAR_options);
It is usually good practice to report impulse responses along with confidence intervals.
The VAR toolbox computes a lower and an upper bound for the impulse responses via bootstrap:
- The residuals are resampled and the model parameters are re-estimated a large number of times.
- For each new estimate, impulse responses are computed.
- All impulse responses are compared to identify the lower, upper and median impulse response.
% Compute confidence intervals using bootstrap methods
[IRF_lower,IRF_upper,IRF_median] = VARirband(VAR,VAR_options);
Loop 10 / 100 draws
Loop 20 / 100 draws
Loop 30 / 100 draws
Loop 40 / 100 draws
Loop 50 / 100 draws
Loop 60 / 100 draws
Loop 70 / 100 draws
Loop 80 / 100 draws
Loop 90 / 100 draws
Loop 100 / 100 draws
-- Done!
Finally, we can plot the impulse responses:
% Figures related options
VAR_options.savefigs = false;
% Plot impulse response functions
VARirplot(IRF_median,VAR_options,IRF_lower,IRF_upper);
Impulse response to Inflation
Impulse response to Output
Impulse response to FedFunds
Historical decomposition
The historical decomposition summarises the history of each endogenous variable in the light of the VAR.
It asks the following question: given the estimated model, what is the sequence of shocks that is able to replicate the time series of Inflation, Output and Federal Funds Rate?
In other words, the historical decomposition tells us what portion of the deviation of the endogenous variables from their unconditional mean is due to the each shock.
% Compute historical decomposition
% Plot historical decomposition
VARhdplot(HistDecomp,VAR_options);
Historical decomposition of Inflation
Historical decomposition of Output
Historical decomposition of FedFunds
Forecast Error Variance decomposition
The variance decomposition of the forecast errors tells us the contribution of each shock to the total variability of each endogenous variable.
It is used to understand how much of the variability of each variable is explained by a given shock. The dynamic nature of the VAR model also implies that the contribution of each shock may change over time, for instance one shock may be important in the first few periods but less important in the long-run.
% Compute forecast error variance decomposition
[FEVD,VAR] = VARfevd(VAR,VAR_options);
% Compute confidence interval via bootstrap
[FEVDINF,FEVDSUP,FEVDMED] = VARfevdband(VAR,VAR_options);
Loop 10 / 100 draws
Loop 20 / 100 draws
Loop 30 / 100 draws
Loop 40 / 100 draws
Loop 50 / 100 draws
Loop 60 / 100 draws
Loop 70 / 100 draws
Loop 80 / 100 draws
Loop 90 / 100 draws
Loop 100 / 100 draws
-- Done!
VARfevdplot(FEVDMED,VAR_options,FEVDINF,FEVDSUP);
Variance decomposition of Inflation
Variance decomposition of Output
Variance decomposition of FedFunds