---
title: "Estimating Dynamic Causal Effects with hetiv"
author: "Daniel Kaufmann, Marc Burri, Valentin Grob"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Estimating Dynamic Causal Effects with hetiv}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 9,
  fig.height = 5,
  out.width = "100%"
)
```

## Introduction

**hetiv** provides tools for measuring and identifying multi-dimensional structural shocks in dynamic models using two complementary IV approaches:

- **Heteroskedasticity-IV** (Rigobon 2003, Rigobon and Sack 2004, Lewis 2022, Burri and Kaufmann, 2026b): exploits the higher
  variance of outcome variables on policy event days relative to control days to
  identify structural shocks without requiring external instruments.
- **Proxy-IV** (Mertens and Ravn 2013, Stock and Watson 2018): uses an external
  instrument (proxy) that is correlated with the shock of interest but
  uncorrelated with other shocks.

Both approaches are implemented as local projection IV estimators following
Jordà (2005), which directly produce impulse response functions (IRFs) across
multiple horizons. Inference is based on heteroskedasticity-robust standard
errors by default. Montiel Olea et al. (2025) show that these standard errors
suffice for local-projection impulse responses under weak conditions, even
though multi-step forecast errors are typically serially correlated. Newey-West
HAC standard errors are also available via `cov_type = "NW"` as an optional
robustness check.

The package allows for multiple shocks and endogenous variables. Weak-instrument HAR inference is implemented via the generalised minimum eigenvalue test of Lewis and Mertens (2025), which nests the classical Stock-Yogo (2005) test for the univariate homoskedastic case. `gweakivtest()` is a direct port of the Matlab files by Lewis and Mertens (2025) available on **https://karelmertens.com/research/**.

In addition, the package provides `kfpredict()` for predicting the underlying unobserved shocks based on the Kalman filter, as suggested by Burri and Kaufmann (2026a).

This vignette demonstrates the functionality using a simulated four-variable VAR
with two structural shocks.

## Simulated data

We simulate data from a VAR(2) with $N = 4$ variables, $E = 2$ event (policy)
shocks, $R = 2$ regular shocks, and $P = 2$ lags over $T = 500$ observations. The model reads:

$$y_t = \Psi \varepsilon_t + \Gamma v_t + \Phi(L) y_{t-1} + \beta X_t \ \ \text{for } t\in \mathcal{P}$$

$$y_t =  \Gamma v_t + \Phi(L) y_{t-1} + \beta X_t \ \ \text{for } t\in \mathcal{C}$$

where $\mathcal{P}, \mathcal{C}$ denote policy event and other days, respectively, $\Psi$ is the impact matrix of $E$ policy event shocks and $\Gamma$ the impact matrix of $R$ other shocks. $\Phi(L)$ is a conformable lag polynomial, and $X_t$ is a matrix of deterministic terms.

Regular shocks occur on all days. A policy event occurs every 10th period (approximately 10% of observations). The latter introduces heteroskedasticity in the variance-covariance matrix of the reduced-form residuals `hetiv()` exploits to identify the impulse response functions. Shocks are drawn from a standard normal distribution. The impact matrix is lower-triangular, corresponding to the identifying assumption by Burri and Kaufmann (2026b) that the first shock has a contemporaneous effect on all variables, while the second shock has no contemporaneous effect on the first variable. Deterministic weekday patterns are added to variables 3 and 4 to illustrate the role of controls. The event indicator `Ind` equals 1 on policy event days and 0 on control days.


```{r dgp}
library(hetiv)

# Dimensions
N <- 4 # variables
E <- 2 # event shocks
R <- 2 # regular shocks
P <- 2 # VAR lag order
H <- 20 # IRF horizons
Nevn <- 10 # Frequency of event shocks

# Impact matrix for event shocks (N x E); lower triangular for recursive ID
PsiE <- matrix(0, N, E)
PsiE[, 1] <- c(1.0, 0.5, 0.3, 0.2)
PsiE[, 2] <- c(0.0, 1.0, -0.4, 0.3)

SigE <- 4 # event shock variance (scalar, applies to all E shocks)

# Impact matrix for regular shocks (N x R)
PsiR <- matrix(0, N, R)
PsiR[, 1] <- c(1.0, -0.3, 0.2, 0.1)
PsiR[, 2] <- c(0.2, 1.0, -0.1, 0.4)

# VAR coefficient matrices at lags 1 and 2
Phi <- array(0, dim = c(N, N, P))
Phi[, , 1] <- matrix(c(
  0.60, 0.05, -0.04, 0.03,
  0.06, 0.40, 0.05, -0.03,
  -0.05, 0.04, 0.50, 0.06,
  0.03, -0.03, 0.05, 0.70
), N, N, byrow = TRUE)
Phi[, , 2] <- matrix(c(
  0.10, 0.03, -0.02, 0.02,
  0.04, 0.10, 0.03, -0.02,
  -0.03, 0.02, 0.10, 0.03,
  0.02, -0.02, 0.03, 0.10
), N, N, byrow = TRUE)

# Simulate — seed is set internally by simulatedata()
Nobs <- 500
Nbin <- 100
sim <- simulatedata(
  Phi = Phi, SigE = SigE, PsiE = PsiE, PsiR = PsiR,
  Nobs = Nobs, Nbin = Nbin, N = N, R = R, E = E,
  Nevn = Nevn, P = P, eDist = 0, seed = 42
)

# Extract simulated data and event indicator
y_data <- sim$y
Ind <- as.integer(sim$IndE[, 1])

# Add deterministic weekday variation into variables 3 and 4
y_data[, 3] <- y_data[, 3] + 0.5 * seq_len(Nobs) %% 5
y_data[, 4] <- y_data[, 4] + 0.3 * seq_len(Nobs) %% 5
```

To illustrate the proxy-IV approach, we construct a noisy external instrument for the two event shocks by adding Gaussian noise to the true shocks. The proxy is set to zero on control days. This mirrors the typical situation in high-frequency identification of monetary policy shocks, where high-frequency surprises are recorded around policy announcement windows and are zero otherwise in the event-study sample.

```{r proxy}
# RNG state after simulatedata(seed = 42) is deterministic, so no extra seed needed
e_proxy <- sim$eE
e_proxy[, 1] <- e_proxy[, 1] + rnorm(Nobs, sd = 1)
e_proxy[, 2] <- e_proxy[, 2] + rnorm(Nobs, sd = 1)
e_proxy[Ind == 0, ] <- 0
```

## Estimation

### Heteroskedasticity-IV without controls

First, we estimate a misspecified model, failing to control for lagged dependent variables and weekday-variation. Note that we impose a normalization on the impact on the endogenous variable for every shock. The value of the normalization can be chosen by `norm`. Here, we choose the normalization that corresponds to the true impact matrix.

```{r het-no-controls}
res_het <- hetiv(
  y       = y_data,
  O       = y_data,
  Ind     = Ind,
  P       = 0,
  H       = H,
  E       = E,
  norm    = 1,
  details = TRUE
)
```

### Heteroskedasticity-IV with controls

Second, we estimate the correctly specified model, with $P = 2$ lags and the weekday dummies. Note that the weekday dummies absorb the deterministic pattern that we added to the data. The regressions generally include a constant term. Therefore, we only add four weekday dummies. 

```{r het-with-controls}
# Weekday dummies: four indicator variables (one left out as reference)
X_data <- matrix(0, nrow = Nobs, ncol = 4)
for (i in 1:4) X_data[, i] <- as.integer((seq_len(Nobs)) %% 5 == (i - 1))

res_het_X <- hetiv(
  y       = y_data,
  O       = y_data,
  X       = X_data,
  Ind     = Ind,
  P       = P,
  H       = H,
  E       = E,
  norm    = 1,
  details = TRUE
)
```

### Proxy-IV

Third, we use the external instrument to estimate the model via proxy-IV. The same lags are included as in the second model. However, the weekday dummies are dropped due to collinearity, as only the event days are used to identify the responses. Note that the proxy-IV function still requires `Ind` as an input, because, in case the instrument is missing, but `Ind = 1`, we can later on predict the unobserved shock using the Kalman filter with `kfpredict()` (see Burri and Kaufmann, 2026a).

```{r proxy-iv}
res_proxy <- proxyiv(
  y         = y_data,
  O         = y_data,
  Z         = e_proxy,
  Ind       = Ind,
  P         = P,
  H         = H,
  E         = E,
  norm      = 1,
  recursive = FALSE,
  details   = TRUE
)
```

### Proxy-IV with recursive zero restriction

If the proxies are valid, we do not need additional restrictions to identify the impulse responses. However, the function allows us to additionally impose zero restrictions mirroring the assumptions used by `hetiv()`. This is useful in a situation where there is some doubt about the validity of the external instruments, and the researcher is convinced that a zero restriction is valid.

```{r proxy-iv-recursive}
res_proxy_rec <- proxyiv(
  y         = y_data,
  O         = y_data,
  Z         = e_proxy,
  Ind       = Ind,
  P         = P,
  H         = H,
  E         = E,
  norm      = 1,
  recursive = TRUE,
  details   = TRUE
)
```

## Impact matrices

The estimated impact matrix $\Psi$ gives the contemporaneous responses of all $N$
variables to each structural shock. We can compare the estimates across the four
specifications against the true $\Psi$. At first sight, the differences are small. As we will see below, there are still some differences in terms of the accuracy of the estimates, the predicted shocks, and the strength of the instruments.

```{r impact-matrices}
tab_1 <- round(data.frame(
  True = PsiE[, 1],
  HET_IV = res_het$Psi[, 1],
  HET_IV_X = res_het_X$Psi[, 1],
  Proxy_IV = res_proxy$Psi[, 1],
  Proxy_IV_rec = res_proxy_rec$Psi[, 1]
), 2)

knitr::kable(
  tab_1,
  col.names = c(
    "True", "HET-IV without controls", "HET-IV with controls",
    "Proxy-IV with controls", "Proxy-IV with recursive restriction"
  ),
  caption = paste(
    "Impact matrix estimates ($\\Psi$) for shock 1 across different",
    "specifications"
  )
)

tab_2 <- round(data.frame(
  True = PsiE[, 2],
  HET_IV = res_het$Psi[, 2],
  HET_IV_X = res_het_X$Psi[, 2],
  Proxy_IV = res_proxy$Psi[, 2],
  Proxy_IV_rec = res_proxy_rec$Psi[, 2]
), 2)

knitr::kable(
  tab_2,
  col.names = c(
    "True", "HET-IV without controls", "HET-IV with controls",
    "Proxy-IV with controls", "Proxy-IV with recursive restriction"
  ),
  caption = paste(
    "Impact matrix estimates ($\\Psi$) for shock 2 across different",
    "specifications"
  )
)
```

## Impulse response functions

### Impulse responses with confidence intervals

We can assess the accuracy of the estimates, and the dynamic causal effects, using `plotirf()`. The function plots IRFs with shaded confidence bands for one estimation approach. The confidence levels can be chosen freely. Here, we show the 90%, 95%, and 99% confidence intervals for the HET-IV estimates with controls.

```{r plotirf, fig.height = 7}
var_labels <- paste0("Variable ", 1:N)

plots_het_X <- plotirf(
  IRFest = res_het_X$irf,
  IRFse  = res_het_X$se,
  HTick  = 5,
  Labels = var_labels,
  ci     = c(0.90, 0.95, 0.99)
)

for (j in seq_len(E)) {
  idx <- ((j - 1) * N + 1):(j * N)
  panel <- cowplot::plot_grid(plotlist = plots_het_X[idx], ncol = 2)
  title <- cowplot::ggdraw() +
    cowplot::draw_label(
      paste0("HET-IV with controls — Shock ", j, " (90%, 95%, 99% CI)"),
      size = 11
    )
  print(cowplot::plot_grid(title, panel, ncol = 1, rel_heights = c(0.05, 1)))
}
```

### HET-IV estimates versus true IRF

We can compare the estimates to the true impulse responses. `plot2irf()` overlays two sets of IRFs. We compare the HET-IV estimate (blue) against the theoretical IRF computed from the known VAR parameters (red). The function `computeirf()` computes the theoretical IRF. The standard errors for the true IRF are set to zero, so that no confidence bands are plotted for the true IRF.

```{r plot2irf-true, fig.height = 7}
irf_true <- computeirf(PsiE, Phi, H, cum = FALSE)
irf_true_se <- array(0, dim = dim(irf_true), dimnames = dimnames(irf_true))

plots_vs_true <- plot2irf(
  IRF1   = res_het_X$irf,
  IRF1se = res_het_X$se,
  IRF2   = irf_true,
  IRF2se = irf_true_se,
  HTick  = 5,
  Labels = var_labels,
  ci     = 0.90
)

for (j in seq_len(E)) {
  idx <- ((j - 1) * N + 1):(j * N)
  panel <- cowplot::plot_grid(plotlist = plots_vs_true[idx], ncol = 2)
  title <- cowplot::ggdraw() +
    cowplot::draw_label(
      paste0("HET-IV with controls (blue) vs True (red) — Shock ", j, " (90% CI)"),
      size = 11
    )
  print(cowplot::plot_grid(title, panel, ncol = 1, rel_heights = c(0.05, 1)))
}
```



### HET-IV estimates versus misspecified model

We can also compare the accuracy of the estimates using the correct and misspecified models.

```{r plot2irf-misspecified, fig.height = 7}
plots_vs_true <- plot2irf(
  IRF1   = res_het_X$irf,
  IRF1se = res_het_X$se,
  IRF2   = res_het$irf,
  IRF2se = res_het$se,
  HTick  = 5,
  Labels = var_labels,
  ci     = 0.90
)

for (j in seq_len(E)) {
  idx <- ((j - 1) * N + 1):(j * N)
  panel <- cowplot::plot_grid(plotlist = plots_vs_true[idx], ncol = 2)
  title <- cowplot::ggdraw() +
    cowplot::draw_label(
      paste0("HET-IV with controls (blue) vs without controls (red) — Shock ", j, " (90% CI)"),
      size = 11
    )
  print(cowplot::plot_grid(title, panel, ncol = 1, rel_heights = c(0.05, 1)))
}
```

We see that the confidence intervals are wider using the misspecified model (red). This is especially true for the third and fourth variables that are affected by the deterministic weekday pattern. The misspecified model does not account for this pattern, which leads to more uncertainty and biased estimates. The intervals are also somewhat wider for variables 1 and 2 because we fail to control for lagged dependent variables.


### HET-IV versus Proxy-IV
Next, we compare the HET-IV estimates to the Proxy-IV estimates. We obtain relatively similar point estimates across the two approaches.

```{r plot2irf-het-proxy, fig.height = 7}
plots_vs_proxy <- plot2irf(
  IRF1   = res_het_X$irf,
  IRF1se = res_het_X$se,
  IRF2   = res_proxy$irf,
  IRF2se = res_proxy$se,
  HTick  = 5,
  Labels = var_labels,
  ci     = 0.90
)

for (j in seq_len(E)) {
  idx <- ((j - 1) * N + 1):(j * N)
  panel <- cowplot::plot_grid(plotlist = plots_vs_proxy[idx], ncol = 2)
  title <- cowplot::ggdraw() +
    cowplot::draw_label(
      paste0("HET-IV (blue) vs Proxy-IV (red) — Shock ", j, " (90% CI)"),
      size = 11
    )
  print(cowplot::plot_grid(title, panel, ncol = 1, rel_heights = c(0.05, 1)))
}
```


## Shock extraction

`kfpredict()` recovers structural shocks from reduced-form residuals via a
Kalman filter prediction. It requires the estimated impact matrix `Psi`, the
residual covariance matrices on event and control days (`Sig`, `SigR`), and the
residuals `et`.

```{r kfpredict}
shocks_het <- kfpredict(
  Sig = res_het$Sig, SigR = res_het$SigR,
  Psi = res_het$Psi, et = res_het$et
)

shocks_het_X <- kfpredict(
  Sig = res_het_X$Sig, SigR = res_het_X$SigR,
  Psi = res_het_X$Psi, et = res_het_X$et
)

shocks_proxy <- kfpredict(
  Sig = res_proxy$Sig, SigR = res_proxy$SigR,
  Psi = res_proxy$Psi, et = res_proxy$et
)

shocks_proxy_rec <- kfpredict(
  Sig = res_proxy_rec$Sig, SigR = res_proxy_rec$SigR,
  Psi = res_proxy_rec$Psi, et = res_proxy_rec$et
)
```

We assess the accuracy of the predictions by comparing them to the true event shocks from the underlying model.

```{r shock-correlation}
true_shocks <- sim$eE

cor_df1 <- data.frame(
  True = true_shocks[, 1],
  Proxy = e_proxy[, 1],
  HET_IV = shocks_het[, 1],
  HET_IV_X = shocks_het_X[, 1],
  Proxy_IV_X = shocks_proxy[, 1],
  Proxy_IV_Rec = shocks_proxy_rec[, 1]
)

cor_df2 <- data.frame(
  True = true_shocks[, 2],
  Proxy = e_proxy[, 2],
  HET_IV = shocks_het[, 2],
  HET_IV_X = shocks_het_X[, 2],
  Proxy_IV_X = shocks_proxy[, 2],
  Proxy_IV_Rec = shocks_proxy_rec[, 2]
)

tab_all <- data.frame(
  rbind(
    round(cor(cor_df1, use = "complete.obs"), 2)[1, ],
    round(cor(cor_df2, use = "complete.obs"), 2)[1, ]
  )
)
rownames(tab_all) <- c("Shock 1", "Shock 2")

knitr::kable(
  tab_all,
  col.names = c(
    "True", "Proxy", "HET-IV without controls", "HET-IV with controls",
    "Proxy-IV with controls", "Proxy-IV with recursive restriction"
  ),
  row.names = TRUE,
  caption = "Correlation of predicted shocks with true shocks"
)
```

The correlation between the true shocks and the proxy is slightly lower than one due to an attenuation bias (see Burri and Kaufmann, 2026a), which depends on the variance of the noise term affecting the proxy. The correlation of the prediction is also relatively low if we use a misspecified model. However, if we use the correct model, the correlation is close to unity for Shock 1 and 0.94 for Shock 2. Proxy-IV with controls yields a high correlation as well, at least for Shock 1. Note that the specific values depend on the exact nature of the simulated data.

## Weak instrument test

`gweakivtest()` implements the generalised minimum eigenvalue test of Lewis and
Mertens (2025), which is robust to heteroskedasticity and autocorrelation and applicable to multiple endogenous regressors and multiple instruments.
Note that classical Stock-Yogo (2005) test is applicable only for the homoskedastic case.

Both `hetiv()` and `proxyiv()` return a `WeakData` object when `details =
TRUE`. This data frame contains the data expected by `gweakivtest()`. The code below extracts the relevant columns from the `WeakData` data frame and runs the weak instrument test for each of the four specifications. For illustration, we use the heteroskedasticity-robust version (`EHW`). For the HAR (Newey West) version, use the option `NW`. By default, the function uses a bias tolerance of 10% at a significance level of 5%.

```{r gweakivtest}
# Helper: extract y, Y, X, Z from WeakData and run gweakivtest
run_weaktest <- function(weakdata, E) {
  # y: outcome variable E+1 (not used as endogenous regressor)
  y <- weakdata[, paste0("y", E + 1)]
  # Y: first E outcome variables (endogenous regressors)
  Y <- weakdata[, paste0("y", 1:E)]
  # Z: the E instruments
  Z <- weakdata[, paste0("Z", 1:E), ]
  # X: lagged ("o*"), deterministic ("x*"), and indicator ("i*") controls.
  # gweakivtest() adds a constant term if one is missing.
  ctrl <- startsWith(colnames(weakdata), "o") |
    startsWith(colnames(weakdata), "x") |
    startsWith(colnames(weakdata), "i")
  X <- if (any(ctrl)) {
    weakdata[, ctrl, drop = FALSE]
  } else {
    matrix(numeric(0), nrow(weakdata), 0)
  }

  gweakivtest(y, Y, X, Z, cov_type = "EHW")
}

specs <- list(
  "HET-IV, no controls" = res_het$WeakData,
  "HET-IV, with controls" = res_het_X$WeakData,
  "Proxy-IV, with controls" = res_proxy$WeakData,
  "Proxy-IV, with controls + recursive restriction" = res_proxy_rec$WeakData
)

wt_results <- lapply(specs, run_weaktest, E = E)
```

```{r gweakivtest-table}
tab <- do.call(rbind, lapply(names(wt_results), function(nm) {
  r <- wt_results[[nm]]
  data.frame(
    Specification = nm,
    Statistic = round(r$gmin_generalized, 2),
    LM_CV = round(r$gmin_generalized_critical_value, 2),
    Strong = ifelse(r$gmin_generalized > r$gmin_generalized_critical_value,
      "Yes", "No"
    ),
    stringsAsFactors = FALSE
  )
}))

knitr::kable(
  tab,
  col.names = c(
    "Specification", "Statistic", "LM critical value",
    "Strong instruments?"
  ),
  caption = "Weak instrument test results (Lewis-Mertens generalised minimum eigenvalue test)"
)
```

The specification without controls is affected by a weak instrument problem. The correctly specified model passes the weak instrument tests. Proxy-IV fails the test. However, note that this depends on the degree of noise used to construct the proxy and does not suggest that `hetiv()` is generally superior to `proxyiv()`. Interestingly, at least for HET-IV, the critical values are much higher than the common rule of thumb for the Stock and Yogo (2005) test and higher than in the HAR test for the univariate case by Montiel Olea and Pflueger (2013) and Lewis (2022), which is around 23.


## References

Burri, M. and Kaufmann, D. (2026a). Measuring monetary policy shocks.
*IRENE Working Papers* 24-03, IRENE Institute of Economic Research,
University of Neuchâtel.

Burri, M. and Kaufmann, D. (2026b). Multiple monetary policy shocks from daily
data: A heteroskedasticity IV approach. *Economics Letters*, 268, 113091.
<https://doi.org/10.1016/j.econlet.2026.113091>

Jordà, Ò. (2005). Estimation and inference of impulse responses by local
projections. *American Economic Review*, 95(1), 161–182.

Lewis, D. J. (2022). Robust inference in models identified via
heteroskedasticity. *Review of Economics and Statistics*, 104(3), 510–524.

Lewis, D. J. and Mertens, K. (2025). A robust test for weak instruments for
2SLS with multiple endogenous regressors. *Review of Economic Studies*,
DOI: 10.1093/restud/rdaf103.

Mertens, K. and Ravn, M. O. (2013). The dynamic effects of personal and
corporate income tax changes in the United States. *American Economic Review*,
103(4), 1212–1247.

Montiel Olea, J. L. and Pflueger, C. E. (2013). A robust test for weak instruments. 
*Journal of Business & Economic Statistics*, 31(3), 358-369.

Montiel Olea, J. L., M. Plagborg-Møller, E. Qian, C. K. Wolf (2025). *Local projections or VARs? A primer for macroeconomists*. NBER Macroeconomics Annual 2025, vol. 40, pp. 1-64, National Bureau of Economic Research.

Rigobon, R. (2003). Identification through heteroskedasticity. *Review of
Economics and Statistics*, 85(4), 777–792.

Stock, J. H. and Watson, M. W. (2018). Identification and estimation of dynamic
causal effects in macroeconomics using external instruments. *Economic Journal*,
128(610), 917–948.

Stock, J. H. and Yogo, M. (2005). Testing for weak instruments in linear IV
regression. In D. W. K. Andrews and J. H. Stock (Eds.), *Identification and
inference for econometric models*, pp. 80–108. Cambridge University Press.
