library(tidyverse)
library(tidyMacro)
library(zoo)
library(patchwork)
set_theme(fThemeTidyMacro())Local Projections with Shocks and IV
A guided tour via Jordà and Taylor (2025)
This chapter is a replication of Jorda and Taylor (2025) with fLP() and fLPIV(). Syntax conventions common to all four local-projection engines — formulas, ..macros, l()/f(), multi-band conf, thread rules, output shapes — live in the Local Projections syntax primer; this chapter focuses on the empirical exercise.
The empirical exercise projects log CPI (×100) on a unit Romer-Romer monetary policy shock over 1985Q1-2007Q4, using long differences CPI_{t+h} - CPI_{t-1} and four lags of real GDP growth, CPI inflation, and the short-term rate change as controls. Two features specific to reproducing Jorda and Taylor (2025) in MATLAB (LPmodel.m) matter here:
- The Newey-West bandwidth rule is Jordà’s classic
= h, i.e.nw_offset = 0(the syntax primer explains the default). Under this rule the standard errors matchOLSmodel.mbit-for-bit. - The LP-IV specification uses a fixed HAC bandwidth of 6 (
nw_lags_iv = 6), matching Stata’svce(hac nw 6).
1 Setup
2 Data
data("JordaTaylor2025")
jt_lp <- JordaTaylor2025 |>
mutate(Date = as.Date(as.yearqtr(Date, format = "%YQ%q"))) |>
dplyr::filter(Date >= "1985-01-01", Date <= "2007-10-01") |>
transmute(
Date,
CPI = 100 * logCPI,
RRShock,
dLogGDP,
dLogCPI,
dSTIR
)
jt_lp |> head()
#> # A tibble: 6 × 6
#> Date CPI RRShock dLogGDP dLogCPI dSTIR
#> <date> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1985-01-01 -56.7 0.0864 4.16 4.57 -0.790
#> 2 1985-04-01 -56.1 -0.146 3.48 2.47 -0.560
#> 3 1985-07-01 -55.5 0.251 6.14 2.31 -0.0200
#> 4 1985-10-01 -54.9 0.0732 2.48 2.31 0.200
#> 5 1986-01-01 -54.4 -0.0949 4.09 1.95 -0.270
#> 6 1986-04-01 -54.0 0.277 1.09 1.60 -0.910jt_lp |>
select(Date, CPI, RRShock, dLogGDP, dSTIR) |>
pivot_longer(-Date, names_to = "variable") |>
ggplot(aes(x = Date, y = value)) +
geom_line(linewidth = 0.7) +
facet_wrap(~ variable, scales = "free", ncol = 2) +
scale_x_date(date_labels = "%Y.%m") +
labs(x = NULL, y = NULL)
3 1 · Classic LP-OLS replication — one fit, two bands
The Jordà-Taylor (2025) result. A single call to fLP() delivers both the 68% and 95% bands: pass conf = c(68, 95) and the bands are rebuilt from the raw irfs_se in the fit object — no need to estimate twice.
controls <- c("dLogGDP", "dLogCPI", "dSTIR")
lp_jt <- fLP(
CPI ~ RRShock + l(..controls, 1:4),
data = jt_lp,
horizons = 17, # MATLAB nsteps = 18 means h = 0, ..., 17
shock = "RRShock",
conf = c(68, 95), # multi-band from one fit
nw_lags = 0,
nw_offset = 0, # classic Jordà rule; matches OLSmodel.m
cumulative = TRUE,
n_threads = 1
)3.1 Direct plot
fPlotLP() draws one ribbon per band automatically (widest first, so the tighter band sits visually on top).
p1 <- fPlotLP(lp_jt) +
labs(title = "CPI response to a unit Romer-Romer shock (Auto)",
y = "Percent")3.2 Customized plot from tidy data
Every plot can also be assembled by hand from the tidy plot data.
p2 <- fPlotLP(lp_jt, return_data = TRUE) |>
ggplot(aes(x = horizon, y = point)) +
geom_ribbon(aes(ymin = lower, ymax = upper, group = conf, alpha = conf, fill = NULL),
fill = "#407EC9") +
geom_line(linewidth = 0.9, colour = "#910048") +
geom_hline(yintercept = 0, linewidth = 0.4, linetype = "dashed") +
scale_alpha_manual(values = c(`0.95` = 0.18, `0.68` = 0.34)) +
scale_x_continuous(breaks = seq(0, 15, 5)) +
labs(title = "CPI response to Romer-Romer shock (Manual)", y = "Percent")
# Compare charts
p1 + p2
fPlotLP(..., return_data = TRUE)
4 2 · Many responses to one shock: Multi-equation LHS
fLP() accepts a vector-valued LHS: c(y1, y2, y3) ~ shock + .... Each LHS variable is estimated as its own local projection against the shared RHS — in parallel across horizons (OpenMP) — and returned as a single fLP object. Combined with conf = c(68, 95), one call yields a full multi-response × multi-band figure.
lp_multi <- fLP(
c(CPI, dSTIR) ~ RRShock + l(..controls, 1:4),
data = jt_lp,
horizons = 17,
shock = "RRShock",
conf = c(68, 95),
nw_lags = 0,
nw_offset = 0,
cumulative = TRUE,
n_threads = 1
)fPlotLP(lp_multi, facet_ncol = 2) +
labs(title = "Responses to a Romer-Romer shock (1 fit, 2 LHS × 2 bands)")
Behind the scenes each response has its own regression; only the RHS is shared. This is materially faster than looping in R because the horizon loop is parallelised in C++ and only the OLS cross-products are formed once per horizon. Note also that fLP() does not auto-augment the RHS — the RHS you write is exactly what enters each regression. For instance, above the dSTIR equation happens to have dSTIR_l1..l4 on the RHS only because dSTIR is a member of the ..controls macro; if you drop it from controls, its own lags disappear from the RHS.
5 3 · Level vs cumulative responses
cumulative = TRUE regresses y_{t+h} - y_{t-1} on the RHS; FALSE regresses the level y_{t+h}. Comparing them makes the “level” vs “long-difference” distinction explicit.
lp_cum <- fLP(
CPI ~ RRShock + l(..controls, 1:4),
data = jt_lp,
horizons = 17,
shock = "RRShock",
conf = c(68, 90),
nw_lags = 0,
nw_offset = 0,
cumulative = TRUE,
n_threads = 1
)
lp_lev <- fLP(
CPI ~ RRShock + l(..controls, 1:4),
data = jt_lp,
horizons = 17,
shock = "RRShock",
conf = c(68, 90),
nw_lags = 0,
nw_offset = 0,
cumulative = FALSE,
n_threads = 1
)
fPlotLP(lp_cum) + ggtitle('Cumulative Response') + fPlotLP(lp_lev) + ggtitle('Level Response')
6 4 · Newey-West bandwidth: matching MATLAB vs the MA&R default
Two bandwidth rules give identical point estimates but different SE bands: the classic Jordà = h rule (nw_offset = 0, matches LPmodel.m bit-for-bit) and the Miranda-Agrippino & Ricco = h + 1 rule (nw_offset = 1, the fLP() default). The syntax primer covers the general convention; here we visualise the resulting band gap on this replication.
lp_jorda <- fLP(CPI ~ RRShock + l(..controls, 1:4),
data = jt_lp, horizons = 17, shock = "RRShock",
conf = 90, nw_lags = 0, nw_offset = 0,
cumulative = TRUE, n_threads = 1)
lp_mar <- fLP(CPI ~ RRShock + l(..controls, 1:4),
data = jt_lp, horizons = 17, shock = "RRShock",
conf = 90, nw_lags = 0, nw_offset = 1,
cumulative = TRUE, n_threads = 1)
# IRFs are identical; only the NW bandwidth (and hence SEs) differ.
all.equal(as.numeric(lp_jorda$irfs), as.numeric(lp_mar$irfs))
#> [1] TRUEbind_rows(
fPlotLP(lp_jorda, return_data = TRUE) |> mutate(rule = "Jordà (h)"),
fPlotLP(lp_mar, return_data = TRUE) |> mutate(rule = "Miranda-Agrippino (h+1)")
) |>
ggplot(aes(x = horizon, y = point, colour = rule, fill = rule)) +
geom_ribbon(aes(ymin = lower, ymax = upper, color = rule), alpha = 0.5) +
geom_line(linewidth = 0.9) +
geom_hline(yintercept = 0, linewidth = 0.4, linetype = "dashed") +
labs(title = "CPI response — NW bandwidth comparison at 90%",
x = "Horizon", y = "Percent", colour = NULL, fill = NULL)
7 5 · LP-IV — Unemployment response to an FFR shock
fLPIV() estimates local projections with an external instrument, mirroring fLP()’s formula-based syntax. Replicates Jordà & Taylor (2025) Figure 6a: the unemployment response to a unit federal-funds-rate shock, instrumented by the Romer-Romer (RR) narrative shock.
The specification (Ex6 tab of the JT2025 workbook, monthly 1985M1-2000M1):
- LHS:
Unemployment(long difference,y_{t+h} - y_{t-1}) - Endogenous treatment:
FFRates(the fed funds rate) - Exogenous controls: six lags of
Unemployment,Inflation,FFRates - Instruments: contemporaneous RR shock plus six lags (over-identified 2SLS)
- HAC: Newey-West with a fixed bandwidth of 6 (Stata’s
vce(hac nw 6))
The result is delivered horizon-by-horizon via Frisch-Waugh-Lovell 2SLS with a delta-method HAC standard error — the same recipe as Cesa-Bianchi’s LPmodel.m in LP-IV mode.
data("JordaTaylor2025IV")
jt_iv <- JordaTaylor2025IV |>
filter(Date >= as.POSIXct("1985-01-01", tz = "UTC"),
Date <= as.POSIXct("2000-01-01", tz = "UTC")) |>
transmute(Date,
urate = Unemployment,
infl = Inflation,
ffr = FFRates,
rz = RRShock)jt_iv |>
pivot_longer(-Date, names_to = "variable") |>
ggplot(aes(x = Date, y = value)) +
geom_line(linewidth = 0.7) +
facet_wrap(~ variable, scales = "free", ncol = 2) +
labs(x = NULL, y = NULL)
iv_ctrl <- c("urate", "infl", "ffr")
lp_iv <- fLPIV(
urate ~ ffr + l(..iv_ctrl, 1:6),
instruments = ~ l(rz, 0:6),
data = jt_iv,
endog = "ffr",
horizons = 48,
conf = c(68, 95),
nw_lags_iv = 6, # fixed HAC bandwidth = Stata's vce(hac nw 6)
cumulative = TRUE,
n_threads = 1
)
lp_iv
#>
#> Local Projections — IV (fLPIV)
#> ---------------------------------------------
#> Original formula : urate ~ ffr + l(..iv_ctrl, 1:6)
#> Expanded formula : urate ~ ffr + urate_l1 + urate_l2 + urate_l3 + urate_l4 + urate_l5 + urate_l6 + infl_l1 + infl_l2 + infl_l3 + infl_l4 + infl_l5 + infl_l6 + ffr_l1 + ffr_l2 + ffr_l3 + ffr_l4 + ffr_l5 + ffr_l6
#> Instruments : ~rz_l0 + rz_l1 + rz_l2 + rz_l3 + rz_l4 + rz_l5 + rz_l6
#> Endogenous : ffr
#> Controls : urate_l1, urate_l2, urate_l3, urate_l4, urate_l5, urate_l6, infl_l1, infl_l2, infl_l3, infl_l4, infl_l5, infl_l6, ffr_l1, ffr_l2, ffr_l3, ffr_l4, ffr_l5, ffr_l6
#> LHS variables : urate
#> Horizons : 0 to 48
#> Cumulative : TRUE
#> Confidence : 95%, 68%
#> HAC bandwidth : fixed = 6 (vce(hac nw 6))
#> Observations : 175
#>
#> IRF (endog = 'ffr'):
#> urate
#> 0 -0.1006
#> 1 -0.2264
#> 2 -0.2281
#> 3 -0.3843
#> 4 -0.3026
#> 5 -0.0642
#> 6 -0.0791
#> 7 0.0232
#> 8 0.0965
#> 9 0.1766
#> 10 0.2850
#> 11 0.4309
#> 12 0.5384
#> 13 0.5911
#> 14 0.6026
#> 15 0.5316
#> 16 0.5267
#> 17 0.5582
#> 18 0.5159
#> 19 0.4835
#> 20 0.8204
#> 21 0.6397
#> 22 0.7474
#> 23 0.8150
#> 24 0.7392
#> 25 0.9504
#> 26 0.9749
#> 27 0.9934
#> 28 1.0497
#> 29 0.9379
#> 30 0.9229
#> 31 1.0563
#> 32 1.0085
#> 33 0.8723
#> 34 0.8391
#> 35 0.6559
#> 36 0.4724
#> 37 0.3058
#> 38 0.1818
#> 39 0.1458
#> 40 0.3319
#> 41 0.2168
#> 42 0.2054
#> 43 0.2100
#> 44 0.0409
#> 45 0.1004
#> 46 0.1305
#> 47 0.1604
#> 48 0.1674
#>
#> First-stage F (h = 0..H): 7.17, 7.3, 7.03, 6.95, 6.92, 6.81, 6.8, 6.59, 6.68, 6.66, 6.59, 6.68, 6.73, 6.68, 6.6, 6.5, 6.46, 6.41, 6.36, 6.33, 6.29, 6.25, 6.22, 6.17, 6.1, 6.18, 6.11, 6.08, 6.05, 5.95, 5.9, 5.8, 5.75, 5.71, 5.72, 5.67, 5.64, 5.58, 5.54, 5.49, 5.47, 5.43, 5.26, 5.24, 5.21, 5.18, 5.09, 5.08, 4.97The print method summarises the endogenous treatment, control set, instrument formula, and the Newey-West bandwidth rule (fixed vs horizon-varying).
7.1 Impulse response and first-stage strength
fPlotLP(lp_iv) +
labs(title = "Unemployment to a unit FFR shock (LP-IV)",
subtitle = "Instrument: Romer-Romer shock (0 + 6 lags), NW bandwidth = 6",
y = "Percentage points")
tibble::tibble(horizon = 0:48, F = as.numeric(lp_iv$Fstat_fs)) |>
ggplot(aes(x = horizon, y = F)) +
geom_col(fill = "#407EC9", alpha = 0.75) +
geom_hline(yintercept = 10, colour = "#910048", linetype = "dashed") +
labs(x = "Horizon", y = "First-stage F",
title = "First-stage F across horizons")
The Fstat_fs and rsqr_fs fields expose the per-horizon first-stage diagnostics for the (residualised) IV projection — worth checking whenever the point estimates look surprising. The df convention (Th - kc - nz) is documented in the syntax primer.
8 6 · Tidy output for downstream work
tidy.fLP() returns a long data frame. For single-band fits, the shape matches the historical broom-style output (estimate, se, lower, upper). For multi-band fits, one lower_X / upper_X column pair is produced per confidence level. It works on both fLP() and fLPIV() results.
tidy.fLP(lp_jt) |> head(6)
#> horizon lhs shock estimate se lower_95 upper_95 lower_68
#> 1 0 CPI RRShock 0.03926033 0.06578167 -0.08966937 0.1681900 -0.02615677
#> 2 1 CPI RRShock 0.08884440 0.11743912 -0.14133205 0.3190209 -0.02794386
#> 3 2 CPI RRShock -0.02728809 0.16816012 -0.35687588 0.3022997 -0.19451625
#> 4 3 CPI RRShock -0.06389951 0.20501135 -0.46571437 0.3379153 -0.26777466
#> 5 4 CPI RRShock -0.07363836 0.24941001 -0.56247300 0.4151963 -0.32166611
#> 6 5 CPI RRShock -0.22802839 0.31100676 -0.83759044 0.3815337 -0.53731152
#> upper_68
#> 1 0.10467743
#> 2 0.20563266
#> 3 0.13994007
#> 4 0.13997564
#> 5 0.17438939
#> 6 0.08125473tidy.fLP(lp_jorda) |> head()
#> horizon lhs shock estimate se lower upper
#> 1 0 CPI RRShock 0.03926033 0.06578167 -0.06894088 0.1474616
#> 2 1 CPI RRShock 0.08884440 0.11743912 -0.10432577 0.2820146
#> 3 2 CPI RRShock -0.02728809 0.16816012 -0.30388688 0.2493107
#> 4 3 CPI RRShock -0.06389951 0.20501135 -0.40111317 0.2733141
#> 5 4 CPI RRShock -0.07363836 0.24941001 -0.48388132 0.3366046
#> 6 5 CPI RRShock -0.22802839 0.31100676 -0.73958899 0.2835322tidy.fLP(lp_iv) |> head(6)
#> horizon lhs shock estimate se lower_95 upper_95 lower_68
#> 1 0 urate ffr -0.10059663 0.09992538 -0.2964468 0.095253515 -0.1999682
#> 2 1 urate ffr -0.22636130 0.12998155 -0.4811205 0.028397862 -0.3556225
#> 3 2 urate ffr -0.22805006 0.20592301 -0.6316518 0.175551633 -0.4328318
#> 4 3 urate ffr -0.38434557 0.19670393 -0.7698782 0.001187053 -0.5799593
#> 5 4 urate ffr -0.30262493 0.22765158 -0.7488138 0.143563973 -0.5290148
#> 6 5 urate ffr -0.06418874 0.26109180 -0.5759193 0.447541780 -0.3238335
#> upper_68
#> 1 -0.001225047
#> 2 -0.097100122
#> 3 -0.023268294
#> 4 -0.188731792
#> 5 -0.076235021
#> 6 0.1954560549 Summary
nw_offset = 0(classic Jordà= hrule) andnw_lags_iv = 6reproduce the Jorda and Taylor (2025) MATLAB / Stata standard errors.- Both
fLP()andfLPIV()handle multi-bandconfand, forfLP(), multi-equation LHS in a single fit — see the syntax primer for the shared conventions.
References
Jorda, Oscar, and Alan M. Taylor. 2025. “Local Projections.” Journal of Economic Literature 63 (1): 59–110. https://doi.org/10.1257/jel.20241521.