Keywords: Local projections, significance bands, monetary policy shocks, R.

Introduction

Inoue, Jordà, and Kuersteiner (2023) propose employing significance bands with confidence bands to enhance impulse response function visualization. This approach allows for a more nuanced interpretation of statistical estimates. Confidence bands provide information on the uncertainty of each impulse response coefficient. In contrast, significance bands are better suited for visually assessing the null hypothesis, which posits that the entire impulse response is zero.

Standard statistical software and the Lagrange multiplier principle make constructing significance bands straightforward under the null hypothesis. This ease of construction is a noteworthy advantage. The bands’ calculations do not depend on the impulse response horizon, and the paper provides approximate formulas. The authors also discuss bootstrap methods as a less assumption-reliant approach.

Monte Carlo simulations show that IJK’s proposed significance bands outperform traditional confidence bands. The simulations reveal that these bands have favorable size and power properties compared to solely considering whether confidence bands include zero. Thus, the authors argue that significance bands should be a standard graphical representation for impulse response functions.

In this paper, I replicate key findings from IJK’s study. Specifically, I focus on estimating local projections (LPs) for the Romer and Romer monetary policy shock (Romer and Romer, 2023.) I then test these estimates using the significance bands method.

Load packages

Start by loading the necessary packages

library(tidyverse) # for data processing and plotting
library(broom)     # for tidying regression output
library(readxl)    # for reading .xls files
library(dynlm)     # for running dynamic regressions with leads and lags
library(sandwich)  # for Newey-West standard errors
library(lmtest)    # for coefficient tests
library(car)       # for linear hypothesis tests
library(patchwork) # for combing figures

Basic setup

Set the basic parameters

H     <- 20    # LP horizon (in quarters)
pval  <- 0.05  # p-value
lags  <- 4     # lags to be used in LP regressions
nwlag <- H     # lags to be used in Newey-West
block <- H     # bootstrap block size

Set the sample

smpl_start <- ymd("1947-01-01")
smpl_end   <- ymd("2016-12-31")

Set seed for replication

set.seed(12345)

Data

Load the Romer and Romer (2023) dataset (link to source)

df_raw <- read_excel(
  "data/Romer&Romer_Data.xlsx",
  sheet = "Quarterly Data",
  skip = 13
)

Some preprocessing steps including computing the logarithm of the GDP price index (GDPPI) and renaming it as y, and renaming the monetary policy shock variable (RRDUMY) to s:

df <- df_raw %>% 
  rename(
    date = 1,
    s    = RRDUMMY
  ) %>% 
  mutate(
    year    = as.integer(substr(date, 1, 4)),
    quarter = as.integer(substr(date, 6, 6)),
    month   = (quarter - 1) * 3 + 1,
    date    = make_date(year, month, 1),
    y       = 100 * log(GDPPI),
    dy      = c(NA, 400 * diff(log(GDPPI))),
  ) %>% 
  filter(between(date, smpl_start, smpl_end)) %>% 
  select(date, s, y, dy) %>% 
  drop_na()

Compute the long differences of \(y_t\), i.e., \(\Delta^h y_{t+h} \equiv y_{t+h} - y_{t-1}\) for \(h=1,…,H\):

df_fwd <- df %>% 
  bind_cols(
    map_dfc(
      0:H, ~ df %>%
        transmute(!!paste0("fy", .x) := lead(y, .x) - lag(y, 1)))
  )

Orthogonalization

n this section, we orthogonalize the response and treatment (shock) variables with respect to the control variables, specifically the lags of both and a constant. This step offers a clear advantage by simplifying the regression process, making it easier to understand and run. The primary motivation behind this approach is the Frisch-Waugh-Lowell theorem, which suggests that the coefficient of interest remains numerically identical, whether we use the full model or the orthogonalized model for estimation. However, it’s essential to correct for the degrees of freedom when calculating the standard errors to ensure our results are accurate and robust.

Define the portion of the regression formula that includes both the constant and the control variables (4 lags of s and y):

formula_x = "L(y, 1:lags) + L(s, 1:lags)"

Orthogonalize \(z_t\) with respect to the controls and store the residual in a variable named e_s:

formula_s <- paste0("s ~ ", formula_x) %>% as.formula()

model_s <- dynlm(formula_s, data = df_fwd %>% ts())

df_e_s <- model_s$residuals %>% 
  as_tibble() %>% 
  rename(e_s = 1)

Orthogonalize \(\Delta_h y_{t+h}\) for \(h=0,...,H\) with respect to the controls and store the residuals in e_fy0,…,e_fyH:

mat_e_y <- matrix(NA, length(df_e_s$e_s), H+1)

for (i in 0:H) {
  
  formula_y <- paste0("fy", i, "~ ", formula_x) %>%  as.formula()
  model_y   <- dynlm(formula_y, data = df_fwd %>% ts())
  
  if (i == 0) {
    mat_e_y[ , i+1] <- model_y$residuals %>% as.vector()
  } else{
    mat_e_y[ , i+1] <- c(model_y$residuals %>% as.vector(), rep(NA, i))
  }

}

df_e_y <- mat_e_y %>% 
  as_tibble()

colnames(df_e_y) <- paste0("e_fy", 0:H)

Combine the orthogonalized residuals into a single dataframe:

df_e_ys <- df_e_y %>% 
  bind_cols(df_e_s) %>% 
  mutate(date = 1:n())

Local projections

Run local projections using the orthogonalized series. Specifically, estimate the given regression model3:

\[ e^y_{t+h}=e^s_t \beta_h+u_{t+h} \quad \text { for } h=0,1, \ldots, H-1 ; \quad t=1, \ldots, T \]

and store the \(\beta_h\) coefficients and their corresponding standard errors for each horizon.

mat_lp <- matrix(data = NA, ncol = 7, nrow = H+1)  # matrix to store estimation results

for (i in 0:H) {
  
  formula_e_y <- paste0("e_fy", i, " ~ e_s -1") %>% as.formula()
  lp_fit      <- dynlm(formula_e_y, data = df_e_ys %>% ts())
  
  nw_se     <- sqrt(diag(NeweyWest(lp_fit, lag = i, prewhite = FALSE)))
  n         <- length(lp_fit$residuals)
  nw_se_adj <- nw_se * sqrt((n - 1) / (n - 1 - 8)) # adjust for degrees of freedom

  beta_i      <- lp_fit$coefficients["e_s"]
  se_beta_i   <-  nw_se_adj["e_s"]
  
  mat_lp[i+1, 1] <- i
  mat_lp[i+1, 2] <- beta_i
  mat_lp[i+1, 3] <- se_beta_i
  mat_lp[i+1, 4] <- beta_i + 1 * se_beta_i
  mat_lp[i+1, 5] <- beta_i - 1 * se_beta_i
  mat_lp[i+1, 6] <- beta_i + 2 * se_beta_i
  mat_lp[i+1, 7] <- beta_i - 2 * se_beta_i

}

df_lp <- mat_lp %>% 
  as_tibble()

colnames(df_lp) <- 
  c("horizon", "mid", "se", "up1sd", "down1sd", "up2sd", "down2sd")

Significance bands

Significance bands using asymptotic approximation:

  1. Calculate the sample average of the product (the orthogonalized) \(s_t z_t\). Call this \(\hat{\gamma}_{s z}\).
  2. Construct the auxiliary variable \(\eta_t=y_t z_t\) and regress \(\eta_t\) on a constant. The Newey-West estimate of the standard error of the intercept coefficient is an estimate of \(s_{\hat{\eta}}\).
  3. An estimate of \(\sigma / \sqrt{T-h}\), call it \(\hat{s}_{\beta_h}\), is therefore: \[ \hat{s}_{\beta_h}=\frac{\hat{s}_\eta}{\hat{\gamma}_{s z}} \]
  4. Construct the significance bands as: \[ \left[\zeta_{\alpha / 2 H} \hat{\mathrm{s}}_{\beta_h}, \zeta_{1-\alpha / 2 H} \hat{\mathrm{s}}_{\beta_h}\right] \]
gamma      <- df_e_ys$e_s^2
gamma_mean <- mean(gamma)

eta        <- df_e_ys$e_fy0 * df_e_ys$e_s^2
model_eta  <- lm(eta ~ 1)
nw_se_eta  <- sqrt(diag(NeweyWest(model_eta, lag = nwlag)))
seta       <- nw_se_eta["(Intercept)"]

sbeta      <- seta/gamma_mean

bju <- qnorm(1 - (pval / (2 * (H + 1)))) * sbeta
bjd <- -bju

Significance bands using the Wold-Block Bootstrap:

  1. Calculate the sample average of \(s_t z_t\). Call this \(\hat{\gamma}_{s z}\).
  2. Construct the auxiliary variable \(\eta_t=y_t z_t\) and regress \(\eta_t\) on a constant. The Wild Block bootstrap estimate of the standard error of the intercept coefficient is an estimate of \(s_{\hat{\eta}}\).
  3. An estimate of \(\sigma / \sqrt{T-h}\), call it \(\hat{s}_{\beta_h}^b\), is therefore: \[ \hat{s}_{\beta_h}^b=\frac{\hat{s}_{\eta \eta}^b}{\hat{\gamma}_{s z}} \]
  4. Construct the significance bands as: \[ \left[\zeta_{\alpha / 2 H} \hat{s}_{\beta_{h^{\prime}}}^b \zeta_{1-\alpha / 2 H} \hat{s}_{\beta_h}^b\right] \]
set.seed(12345)   # for replication

n         <- length(df_e_ys$e_fy0)
block_id  <- floor((1:n + block - 1) / block)

model_eta <- lm(eta ~ 1)
nw_se_eta <- sqrt(diag(vcovBS(model_eta, cluster = block_id, R = 1000)))
seta      <- nw_se_eta["(Intercept)"]

sbeta     <- seta/gamma_mean

bootju <- qnorm(1 - (pval / (2 * (H + 1)))) * sbeta
bootjd <- -bootju

Combine the estimation results of significance bands into a single dataframe:

df_sb <- tibble(
  bju    = rep(bju, H+1),
  bjd    = rep(bjd, H+1),
  bootju = rep(bootju, H+1),
  bootjd = rep(bootjd, H+1)
)

Next, plot the LP with the confidence bands (one and two standard deviations) and significance bands (Newey-West and bootstrap):

p1 <- df_lp %>% 
  bind_cols(df_sb) %>% 
  ggplot(aes(x = horizon)) +
  geom_line(aes(y = mid), color = 'darkgreen', linewidth = 1.5) +
  geom_line(aes(y = bju), color = "red", linetype = 2, size = 1) +
  geom_line(aes(y = bjd), color = "red", linetype = 2, size = 1) +
  geom_line(aes(y = bootju), color = "blue", linetype = 2, size = 1) +
  geom_line(aes(y = bootjd), color = "blue", linetype = 2, size = 1) +
  geom_ribbon(aes(ymin = down1sd, ymax = up1sd),
              fill = "darkgreen", linetype = 0, alpha = 0.2) +
  geom_ribbon(aes(ymin = down2sd, ymax = up2sd),
              fill = "darkgreen", linetype = 0, alpha = 0.1) +
  scale_x_continuous(expand = c(0,0), breaks = seq(0, H, 1)) +
  # scale_y_continuous(expand = c(0,0), limits = c(-4, 4)) +
  labs(
    x = "Quarters after the shock",
    y = "Percentage points",
    title = "Cumulative response of GDP Price Index to a monetary policy shock",
  ) +
  geom_hline(aes(yintercept = 0)) +
  theme_light(12) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )

p1

The figure illustrates the impulse response of the long difference, calculated from 100 times the logarithm of the GDP Price Index, after a monetary policy shock as described by Romer and Romer (2023). The model incorporates four lagged variables both for the GDP Price Index and the monetary policy shock. Dark and light green shaded areas represent confidence bands corresponding to one and two standard errors, respectively. Additionally, dashed blue lines indicate classical Newey-West significance bands, while dashed red lines represent bootstrapped significance bands.

In examining the graphical representation, it is evident that the pointwise null hypothesis positing no effect remains unchallenged within the bounds of 2 standard deviations, persisting even at an extended horizon of 20 quarters. Conversely, when assessing a broader perspective, the overarching hypothesis of an absolute absence of effect finds little support in the observed data. The impulse response distinctly falls beneath the significance bands, a conclusion consistent with both the Newey-West and the bootstrapped versions. This delineation underscores the nuanced distinctions between point-specific and overarching effects in the context of economic implications.

Stacked LPs

In this section, we employ a technique proposed by Jordà, where we estimate the \(h\) local projection regression simultaneously as a panel anduse Driscoll and Kraay (1998) HAC robust standard errors for inference. The advantage of this estimation method is that it provides us with the coefficients’ variance-covariance matrix, thereby enabling us to test the hypothesis that \(\beta_h=0\) for all \(h\).

We start by pivoting the data to long format

df_e_ys_long <- 
  df_e_ys %>%
  pivot_longer(-c(date, e_s), names_to = "horizon", values_to = "e_fy") %>%
  mutate(
    horizon = as.factor(gsub("^e_fy", "", horizon))
  ) %>% 
  select(date, horizon, e_fy, e_s)

Next, we estimate the data stacked as a panel using the plm package:

plm_model <- plm::plm(e_fy ~ factor(horizon):e_s, data = df_e_ys_long, index = c("horizon", "date"), model = "within")

Finally, we conduct the joint hypothesis test, namely: \[ H_0: \beta_h=0, \text{for } h=0,\ldots,H \]

hypothesis <- sapply(0:H, function(i) paste0("factor(horizon)", i, ":e_s = 0"))

linearHypothesis(plm_model, hypothesis, vcov. = function(x) plm::vcovSCC(x, maxlag=nwlag))
## Linear hypothesis test
## 
## Hypothesis:
## factor(horizon)0:e_s = 0
## factor(horizon)1:e_s = 0
## factor(horizon)2:e_s = 0
## factor(horizon)3:e_s = 0
## factor(horizon)4:e_s = 0
## factor(horizon)5:e_s = 0
## factor(horizon)6:e_s = 0
## factor(horizon)7:e_s = 0
## factor(horizon)8:e_s = 0
## factor(horizon)9:e_s = 0
## factor(horizon)10:e_s = 0
## factor(horizon)11:e_s = 0
## factor(horizon)12:e_s = 0
## factor(horizon)13:e_s = 0
## factor(horizon)14:e_s = 0
## factor(horizon)15:e_s = 0
## factor(horizon)16:e_s = 0
## factor(horizon)17:e_s = 0
## factor(horizon)18:e_s = 0
## factor(horizon)19:e_s = 0
## factor(horizon)20:e_s = 0
## 
## Model 1: restricted model
## Model 2: e_fy ~ factor(horizon):e_s
## 
## Note: Coefficient covariance matrix supplied.
## 
##   Res.Df Df  Chisq Pr(>Chisq)    
## 1   5544                         
## 2   5523 21 704.75  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comparison to Romer and Romer (2023)

In Romer and Romer (2023), the researchers did not explore the effects on the price level as we did. Instead, they focused on directly assessing the implications for inflation. Specifically, they quantified their response as 400 times the log difference in the quarterly GDP Price Index. For context, we replicated their results.4

The left-hand side of the figure presents the LP and bands related to the price level. In contrast, the right-hand side highlights the LP and bands relevant to inflation. Significantly, we cannot reject the null hypothesis that there is no effect, based on the significance bands approach.

This seeming “discrepancy” arises from the distinct nature of the measurements. Analyzing the cumulative price level version raises the issue of whether the price level, over time, deviates either above or below its initial treatment point. On the other hand, individual quarterly inflation rates might seem negligible, but their combined effect over time can lead to a notable shift.

References

Driscoll, J. C., & Kraay, A. C. (1998). Consistent covariance matrix estimation with spatially dependent panel data. Review of economics and statistics, 80(4), 549-560.

Inoue, A., Jordà, Ò., & Kuersteiner, G. M. (2023). Significance Bands for Local Projections. arXiv preprint arXiv:2306.03073.

Romer, C. D., & Romer, D. H. (2023). Presidential Address: Does Monetary Policy Matter? The Narrative Approach after 35 Years. American Economic Review, 113(6), 1395-1423.

  1. Here’s a link to the GitHub repository that contains both the code and dataset utilized in this notebook.↩︎

  2. The views expressed in this notebook are my own and do not represent the views of the Bank of Israel or any of its staff members. I would like to thank Oscar Jordà for his most helpful feedback and for sharing his code. Any remaining errors are mine.↩︎

  3. The lpirf R package is excellent for estimating and visualizing both linear and nonlinear impulse responses through local projections. While I haven’t utilized lpirf in this analysis due to the need for greater flexibility, I highly recommend it!↩︎

  4. For those curious about the details of our replication process, the code is accessible in the notebook’s GitHub repository.↩︎