This page gives examples and replicates empirical results to show how to use the R-package lpirfs.
Example: Linear impulse responses with local projections
Load library:
library(lpirfs)Load data from package to estimate a simple, new-Keynesian, closed-economy model. These data are used by Jordà (2005) in chapter IV, p. 174. See the data’s help file in the package or the paper for a detailed description of the data.
endog_data <- interest_rules_var_dataEstimate linear impulse responses with function lp_lin. Note that the endogenous and exogenous data have to be a data.frame.
results_lin <- lp_lin(endog_data,
lags_endog_lin = 4, # Number of lags for endogenous data
trend = 0, # 0 = no trend, 1 = trend, 2 = trend & trend^2
shock_type = 1, # 0 = standard deviation shock, 1 = unit shock
confint = 1.96, # Width of confidence bands:
# 1 = 68%, 1.67 = 90%, 1.96 = 95%
hor = 12) # Number of cores to use. When NULL, the number of cores
# is chosen automatically Create plots for impulse responses of linear model with function plot_lin.
linear_plots <- plot_lin(results_lin)Display single impulse responses:
- The first plot shows the response of the first variable (GDP_gap) to a shock of the first variable (GDP_gap).
- The second plot shows the response of the first variable (GDP_gap) to a shock of the second variable (inflation).
- …
linear_plots[[1]]
linear_plots[[2]]
Display all plots (compare with Figure 5 in Jordà (2005), p. 176):
# Load libraries to plot all impulse responses
# The package does not depend on those packages so they have to be installed
library(ggpubr)
library(gridExtra)
# Show all plots
lin_plots_all <- sapply(linear_plots, ggplotGrob)
marrangeGrob(lin_plots_all, nrow = ncol(endog_data), ncol = ncol(endog_data), top = NULL)
Example: Nonlinear impulse responses with local projections
Load data set from package to estimate a nonlinear, new-Keynesian, closed-economy model. These data are used by Jordà (2005) in chapter IV, p.176. See the data’s help file or the paper for a detailed description of the data.
endog_data <- interest_rules_var_dataThe switching variable (\(z_t\)) can either be decomposed by the Hodrick-Prescott filter (see Auerbach and Gorodnichenko, 2013) or directly plugged into the following logistic function:
\[F({z_t}) = \ \frac{e^{(-\gamma z_t)}}{\left(1 + e^{(-\gamma z_t)}\right)},\]
where \(\gamma > 0\). To differentiate between the two regimes, the endogenous variables (\(\boldsymbol{y}_{t-p}\)) are multiplied with the values of the transition function at t − 1 where
- Regime 1 \((R_1)\): \(\boldsymbol{y}_{t-i}\cdot(1-F(z_{t-1}))\), \(\qquad i = 1, ..., p\),
and
- Regime 2 \((R_2)\): \(\boldsymbol{y}_{t-i}\cdot F(z_{t-1})\), \(\hskip 4.5em i = 1, ..., p\).
IMPORTANT: The index of \(z\) is set to t − 1 in order to avoid contemporaneous feedback (see Auerbach and Gorodnichenko 2012). You can suppress this option, however, with ‘lag_switching = FALSE’.
Nonlinear impulse responses are computed as in Ahmed and Cassou (2016). First, a reduced VAR is estimated to obtain the covariance matrix of the residuals \(\Sigma\). The Cholesky decomposition is then applied to obtain the shock matrix with columns denoted by \(d_j\). IRFs for both regimes are estimated via:
\[\hat{IR}^{R_1}(t,h,d_j) = \hat{\boldsymbol{B}}_{1, R_1}^h d_j, \ \ \ \ \ \ \ \ \ h = 1, ..., H, \] and
\[\hat{IR}^{R_2}(t,h,d_j) = \hat{\boldsymbol{B}}_{1, R_2}^h d_j, \ \ \ \ \ \ \ \ \ h = 1, ..., H, \]
with \(\hat{\boldsymbol{B}}_{1, R1}^0 = I\) and \(\hat{\boldsymbol{B}}_{1, R2}^0 = I\). The parameters are obtained by running a sequence of OLS forecasts (local projections):
\[ \boldsymbol{y}_{t+h} = \boldsymbol{\alpha}^h + \boldsymbol{B}_{1, R_1}^{h} \left(\boldsymbol{y}_{t-1}\cdot(1-F(z_{t-1})\right) \ + \ ...\ +\ \boldsymbol{B}_{p, R_1}^{h} \left(\boldsymbol{y}_{t-p}\cdot(1-F(z_{t-1})\right) + \\ \hskip 7.6em \boldsymbol{B}_{1, R_2}^{h} \left(\boldsymbol{y}_{t-1}\cdot F(z_{t-1}\right)) \ + \ ... \ + \ \boldsymbol{B}_{p, R_2}^{h} \left(\boldsymbol{y}_{t-p}\cdot F(z_{t-1}\right)) + \boldsymbol{\varepsilon}_{t+h}^h, \] with \(h = 1,..., H.\)
\[ \]
Estimate nonlinear impulse responses with package function lp_nl. Note that the endogenous and exogenous data have to be a data.frame.
results_nl <- lp_nl(endog_data,
lags_endog_lin = 4, # Number of lags for (reduced) linear VAR to obtain shock matrix.
lags_endog_nl = 4, # Number of lags for nonlinear model.
trend = 1, # 0 = no trend, 1 = trend, 2 = trend & trend^2.
shock_type = 0, # 0 = standard deviation shock, 1 = unit shock.
confint = 1.67, # Width of confidence bands:
# 1 = 68%, 1.67 = 90%, 1.96 = 95%.
hor = 12, # Length of irf horizons.
switching = endog_data$Infl, # Use inflation rate as switching variable.
use_hp = TRUE, # Use HP-filter? TRUE or FALSE.
lambda = 1600, # Ravn and Uhlig (2002):
# Anuual data = 6.25
# Quarterly data = 1600
# Monthly data = 129,600
gamma = 6) Save values from transition function and list which contains further information on the data for plot function.
fz <- results_nl$fz
specs <- results_nl$specsPlot output gap.
# Make date sequence and store data in a data.frame for ggplot.
dates <- seq(as.Date("1955/12/1"), as.Date("2003/1/1"), by = "quarter")
data_df <- data.frame(x = dates, fz = fz, switch_var = specs$switching[(specs$lags_endog_nl+1):length(endog_data$FF)])
# Plot inflation rate
ggplot(data = data_df) +
geom_line(aes(x = x, y = switch_var)) +
theme_bw() +
ylab("") +
xlab("Date") +
scale_x_date(date_breaks = "5 year", date_labels = "%Y")
# Plot transition function
ggplot(data = data_df) +
geom_line(aes(x = x, y = fz)) +
theme_bw() +
ylab("") +
xlab("Date") +
scale_x_date(date_breaks = "5 year", date_labels = "%Y")
Note that values close to one correspond to periods of low inflation rates (regime 2) and values close zero correspond to periods of high inflation rates (regime 1).
Create and save all plots with function plot_nl
nl_plots <- plot_nl(results_nl)Show single impulse responses for each regime:
- The first plot shows the response of the first variable (GDP_gap) to a shock in the first variable (GDP_gap) in regime 1.
- The second plot shows the response of the first variable (GDP_gap) to a shock in the second variable (Inflation) in regime 2.
# Load packages
# lpirfs does not depend on those packages so they have to be additionally installed
library(ggpubr)
library(gridExtra)
# Save plots based on states
s1_plots <- sapply(nl_plots$gg_s1, ggplotGrob)
s2_plots <- sapply(nl_plots$gg_s2, ggplotGrob)
plot(s1_plots[[1]])
plot(s2_plots[[2]])
Show all impulse responses of regime 1 (high inflation rates):
marrangeGrob(s1_plots, nrow = ncol(endog_data), ncol = ncol(endog_data), top = NULL)
Show all impulse responses of regime 2 (low inflation rates):
marrangeGrob(s2_plots, nrow = ncol(endog_data), ncol = ncol(endog_data), top = NULL)
Example: Linear impulse responses with identified shock
This example replicates results from the Supplementary Appendix by Ramey and Zubairy (2018) (RZ-18). They use local projections to re-evaluate findings in Auerbach and Gorodnichenko (2012) (AG-12).
# Load data which is taken from https://www.journals.uchicago.edu/doi/10.1086/696277
ag_data <- ag_data
sample_start <- 8
sample_end <- dim(ag_data)[1]
# Endogenous variables
endog_data <- ag_data[sample_start:sample_end,3:5]
# Shock variable
shock <- endog_data[, 1]
# Estimate model with constant and 4 lags
results_lin_iv <- lp_lin_iv(endog_data,
shock = shock,
lags_endog_lin = 4,
trend = 0,
confint = 1.96,
hor = 20)
# Make and save plots
iv_lin_plots <- plot_lin(results_lin_iv)- The first element of ‘iv_lin_plots’ shows the response of the first variable (Gov) to the identified shock (Gov).
- The second element of ‘iv_lin_plots’ shows the response of the second variable (Tax) to the the identified shock (Gov).
- …
This plot equals the left plot in the mid-panel of Figure 12 in the Supplementary Appendix by RZ-18.
iv_lin_plots[[1]]
Multiply the log output response by a conversion factor of ~5.6 (AG-12; RZ-18).
GYsc <- mean(exp(endog_data$GDP)/exp(endog_data$Gov))
multiplier_mean <- results_lin_iv$irf_lin_mean*GYsc
multiplier_up <- results_lin_iv$irf_lin_up*GYsc
multiplier_low <- results_lin_iv$irf_lin_low*GYsc
results_lin_iv <- list(irf_lin_mean = multiplier_mean,
irf_lin_up = multiplier_up,
irf_lin_low = multiplier_low,
specs = results_lin_iv$specs)
# Make new plots
iv_lin_plots <- plot_lin(results_lin_iv)Compare with the right plot in the mid-panel of Figure 12 in the Supplementary Appendix by RZ-18.
iv_lin_plots[[3]]
Example: Nonlinear impulse responses with identified shock
# Load data
ag_data <- ag_data
sample_start <- 7
sample_end <- dim(ag_data)[1]
endog_data <- ag_data[sample_start:sample_end, 3:5]The shock variable is created by RZ-18 and available as supplementary data to their paper. The government spending shock is included in the package.
shock <- as.data.frame(ag_data$Gov_shock_mean[sample_start:sample_end])AG-12 include four lags of the 7-quarter moving average growth rate as exogenous regressors in their model (see RZ-18).
exog_data <- as.data.frame(ag_data$GDP_MA[sample_start:sample_end])AG-12 use the 7-quarter moving average growth rate as switching variable. They adjust it to have suffiently long recession periods.
switching_variable <- ag_data[sample_start:sample_end, 6] - 0.8# Estimate nonlinear model
results_nl_iv <- lp_nl_iv(endog_data,
lags_endog_nl = 3,
shock = shock,
exog_data = exog_data,
lags_exog = 4,
trend = 0,
confint = 1.96,
hor = 20,
switching = switching_variable,
use_logistic = TRUE,
lag_switching = TRUE,
use_hp = FALSE,
lambda = NaN,
gamma = 3)Multiply the log output responses by a conversion factor of ~5.6 (AG-12; RZ-18).
GYsc <- mean(exp(endog_data$GDP)/exp(endog_data$Gov))
irf_s1_mean <- results_nl_iv$irf_s1_mean*GYsc
irf_s1_up <- results_nl_iv$irf_s1_up*GYsc
irf_s1_low <- results_nl_iv$irf_s1_low*GYsc
irf_s2_mean <- results_nl_iv$irf_s2_mean*GYsc
irf_s2_up <- results_nl_iv$irf_s2_up*GYsc
irf_s2_low <- results_nl_iv$irf_s2_low*GYsc
results_nl_iv <- list(irf_s1_mean = irf_s1_mean,
irf_s1_up = irf_s1_up,
irf_s1_low = irf_s1_low,
irf_s2_mean = irf_s2_mean,
irf_s2_up = irf_s2_up,
irf_s2_low = irf_s2_low,
specs = results_nl_iv$specs)Make and save plots
plots_nl_iv <- plot_nl(results_nl_iv)
s1_plots <- sapply(plots_nl_iv$gg_s1, ggplotGrob)
s2_plots <- sapply(plots_nl_iv$gg_s2, ggplotGrob)# Show plot
plot(plots_nl_iv$gg_s1[[3]])
This plot corresponds to the red dotted line in the lower panel (right plot) of Figure 12 in the Supplementary Appendix of RZ-18. It shows the response of GDP to a shock in government spending equal to 1% of GDP during periods of economic expansions. It is the same data, identification scheme and threshold definition as in AG-12 (see RZ-18). ```
# Show plot
plot(plots_nl_iv$gg_s2[[3]])
This plot corresponds to the blue dotted line in the lower panel (right plot) of Figure 12 in the Supplementary Appendix of RZ-18. It shows the response of GDP to a shock in government spending equal to 1% of GDP during periods of economic slack. It is the same data, identification scheme and threshold definition as in AG-12 (see RZ-18).
Example: Linear impulse responses via 2SLS
# Set seed
set.seed(007)
# Load data
ag_data <- ag_data
sample_start <- 7
sample_end <- dim(ag_data)[1]
# Endogenous data
endog_data <- ag_data[sample_start:sample_end,3:5]
# Variable to shock with (government spending)
shock <- ag_data[sample_start:sample_end, 3]
# Generate instrument variable that is correlated with government spending
instrum <- 0.9*shock$Gov + rnorm(length(shock$Gov), 0, 0.02) %>%
as.data.frame()
# Estimate linear model via SLS
results_lin_iv <- lp_lin_iv(endog_data,
lags_endog_lin = 4,
shock = shock,
instrum = instrum,
use_twosls = TRUE,
trend = 0,
confint = 1.96,
hor = 20)
# Create all plots
iv_lin_plots <- plot_lin(results_lin_iv)# Show response of GDP
iv_lin_plots[[3]]
Example: Linear impulse responses for panel-data
The general formula for panel irfs is:
\[ y_{i,t+h} = \alpha_{i,h} + shock_{i,t}\beta_h + \boldsymbol{x}_{i,t}\boldsymbol{\gamma}_h + \varepsilon_{i, t + h}, \] where \(\alpha_{i, h}\) is a fixed effect, and \(\boldsymbol{x}_{i,t}\) is a vector of control variables that have to be specified. The variable \(shock\) can also be estimated by an instrument variable approach.
# This example is based on the STATA code 'LPs_basic_doall.do', provided on
# Òscar Jordà's website (https://sites.google.com/site/oscarjorda/home/local-projections)
# It estimates impulse reponses of the ratio of (mortgage lending/GDP) to a
# +1% change in the short term interest rate
# Load libraries to download and read excel file from the website
library(lpirfs)
library(httr)
library(readxl)
library(dplyr)
# Retrieve the JST Macrohistory Database
url_jst <-"http://www.macrohistory.net/JST/JSTdatasetR3.xlsx"
GET(url_jst, write_disk(jst_link <- tempfile(fileext = ".xlsx")))
#> Response [http://www.macrohistory.net/JST/JSTdatasetR3.xlsx]
#> Date: 2019-04-11 11:38
#> Status: 200
#> Content-Type: application/vnd.openxmlformats-officedocument.spreadsheetml.sheet
#> Size: 635 kB
#> <ON DISK> /tmp/RtmpLVnl0G/file6c846068b06f.xlsx
jst_data <- read_excel(jst_link, 2L)
# Swap the first two columns so that 'country' is the first (cross section) and 'year' the
# second (time section) column
jst_data <- jst_data %>%
dplyr::filter(year <= 2013) %>%
dplyr::select(country, year, everything())
# Prepare variables. This is based on the 'data.do' file
data_set <- jst_data %>%
mutate(stir = stir) %>%
mutate(mortgdp = 100*(tmort/gdp)) %>%
mutate(hpreal = hpnom/cpi) %>%
group_by(country) %>%
mutate(hpreal = hpreal/hpreal[year==1990][1]) %>%
mutate(lhpreal = log(hpreal)) %>%
mutate(lhpy = lhpreal - log(rgdppc)) %>%
mutate(lhpy = lhpy - lhpy[year == 1990][1]) %>%
mutate(lhpreal = 100*lhpreal) %>%
mutate(lhpy = 100*lhpy) %>%
ungroup() %>%
mutate(lrgdp = 100*log(rgdppc)) %>%
mutate(lcpi = 100*log(cpi)) %>%
mutate(lriy = 100*log(iy*rgdppc)) %>%
mutate(cay = 100*(ca/gdp)) %>%
mutate(tnmort = tloans - tmort) %>%
mutate(nmortgdp = 100*(tnmort/gdp)) %>%
dplyr::select(country, year, mortgdp, stir, ltrate, lhpy,
lrgdp, lcpi, lriy, cay, nmortgdp)
# Use data_sample from 1870 to 2013 and exclude WWI and WWII
data_sample <- seq(1870, 2016)[which(!(seq(1870, 2016) %in%
c(seq(1914, 1918),
seq(1939, 1947))))]
# Estimate panel model
results_panel <- lp_lin_panel(data_set = data_set,
data_sample = data_sample,
endog_data = "mortgdp",
cumul_mult = TRUE,
shock = "stir",
diff_shock = TRUE,
panel_model = "within",
panel_effect = "individual",
robust_cov = "vcovSCC",
c_exog_data = "cay",
c_fd_exog_data = colnames(data_set)[c(seq(4,9),11)],
l_fd_exog_data = colnames(data_set)[c(seq(3,9),11)],
lags_fd_exog_data = 2,
confint = 1.67,
hor = 5)# Create and plot irfs
plot_lin_panel <- plot_lin(results_panel)
plot(plot_lin_panel[[1]])
# Create and add instrument to data_set
set.seed(123)
data_set <- data_set %>%
group_by(country) %>%
mutate(instrument = 0.8*stir + rnorm(length(stir), 0, sd(na.omit(stir))/10)) %>%
ungroup()
# Estimate panel model with iv
results_panel <- lp_lin_panel(data_set = data_set,
data_sample = data_sample,
endog_data = "mortgdp",
cumul_mult = TRUE,
shock = "stir",
diff_shock = TRUE,
iv_reg = TRUE,
instrum = "instrument",
panel_model = "within",
panel_effect = "individual",
robust_cov = "vcovSCC",
c_exog_data = "cay",
c_fd_exog_data = colnames(data_set)[c(seq(4,9),11)],
l_fd_exog_data = colnames(data_set)[c(seq(3,9),11)],
lags_fd_exog_data = 2,
confint = 1.67,
hor = 5)# Create and plot irfs
plot_lin_panel <- plot_lin(results_panel)
plot(plot_lin_panel[[1]])
Use GMM to estimate panel model.
# Use a much smaller sample to have fewer T than N
data_sample <- seq(2000, 2012)
# Running this example gives a warning at each iteration as the data set is not well suited for
# GMM. GMM is based on N-asymptotics and the data set only contains 27 countries
results_panel <- lp_lin_panel(data_set = data_set,
data_sample = data_sample,
endog_data = "mortgdp",
cumul_mult = TRUE,
shock = "stir",
diff_shock = TRUE,
use_gmm = TRUE,
gmm_model = "onestep",
gmm_effect = "twoways",
gmm_transformation = "ld",
l_exog_data = "mortgdp",
lags_exog_data = 2,
l_fd_exog_data = colnames(data_set)[c(4, 6)],
lags_fd_exog_data = 1,
confint = 1.67,
hor = 5)
# Create and plot irfs
plot_lin_panel <- plot_lin(results_panel)
plot(plot_lin_panel[[1]])
Example: Nonlinear impulse responses for panel-data
# Use data_sample from 1870 to 2013 and exclude WWI and WWII
data_sample <- seq(1870, 2016)[which(!(seq(1870, 2016) %in%
c(seq(1914, 1918),
seq(1939, 1947))))]
# Estimate panel model
results_panel <- lp_nl_panel(data_set = data_set,
data_sample = data_sample,
endog_data = "mortgdp",
cumul_mult = TRUE,
shock = "stir",
diff_shock = TRUE,
panel_model = "within",
panel_effect = "individual",
robust_cov = "vcovSCC",
switching = "lrgdp",
lag_switching = TRUE,
use_hp = TRUE,
lambda = 6.25,
gamma = 10,
c_exog_data = "cay",
c_fd_exog_data = colnames(data_set)[c(seq(4,9),11)],
l_fd_exog_data = colnames(data_set)[c(seq(3,9),11)],
lags_fd_exog_data = 2,
confint = 1.67,
hor = 5)# Create irf plots
nl_plots <- plot_nl(results_panel)
# Compare states
library(ggpubr)
library(gridExtra)
combine_plots <- list(nl_plots$gg_s1[[1]], nl_plots$gg_s2[[1]])
marrangeGrob(combine_plots, nrow = 1, ncol = 2, top = NULL)
References
Ahmed, M.I. , Cassou, S.C. (2016) “Does Consumer Confidence Affect Durable Goods Spending During Bad and Good Economic Times Equally?” Journal of Macroeconomics, 50(1): 86-97. doi:10.1016/j.jmacro.2016.08.008
Auerbach, A. J., and Gorodnichenko Y. (2012). “Measuring the Output Responses to Fiscal Policy.” American Economic Journal: Economic Policy, 4 (2): 1-27. doi:10.1257/pol.4.2.1
Auerbach, A. J., and Gorodnichenko Y. (2013). “Fiscal Multipliers in Recession and Expansion” NBER Working Paper Series. Nr 17447. National Bureau of Economic Research
Croissant, Y., Millo, G. (2008). “Panel Data Econometrics in R: The plm Package.” , 27(2), 1-43. doi: 10.18637/jss.v027.i02.
Jordà, Ò. (2005) “Estimation and Inference of Impulse Responses by Local Projections.” American Economic Review, 95 (1): 161-182. doi:10.1257/0002828053828518
Jordà, Ò, Schularick, M., Taylor, A.M. (2015), “Betting the house”, Journal of International Economics, 96, S2-S18. doi:10.1016/j.jinteco.2014.12.011
Jordà, Ò., Schualrick, M., Taylor, A.M. (2018). “Large and State-Dependent Effects of Quasi-Random Monetary Experiments”, FRBSF working paper. Working Paper
Newey W.K., and West K.D. (1987). “A Simple, Positive-Definite, Heteroskedasticity and Autocorrelation Consistent Covariance Matrix.” Econometrica, 55: 703–708. doi:10.2307/1913610
Ramey, V.A., Zubairy, S. (2018). “Government Spending Multipliers in Good Times and in Bad: Evidence from US Historical Data.” Journal of Political Economy, 126(2): 850 - 901. doi:10.1086/696277
Ravn, M.O., Uhlig, H. (2002). “On Adjusting the Hodrick-Prescott Filter for the Frequency of Observations.” Review of Economics and Statistics, 84(2), 371-376. doi:10.1162/003465302317411604
License
GPL (>= 2)