--- title: "Time Series IV, Panel VCE, and GMM" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Time Series IV, Panel VCE, and GMM} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ## Introduction Classical and heteroskedasticity-robust standard errors both assume that the errors are serially uncorrelated. In time-series and panel instrumental-variables models that assumption rarely holds: shocks persist across periods, and in panels they may also be correlated across units in the same period. When it fails, the usual standard errors are inconsistent, the reported weak-identification and overidentification statistics use the wrong variance, and two-stage least squares is no longer the efficient use of the moment conditions. This vignette covers the tools `ivreg2r` provides for such settings: heteroskedasticity- and autocorrelation-consistent (HAC) standard errors and their kernel machinery, the panel variance estimators (Kiefer, Driscoll-Kraay, Stock-Watson, and cluster-plus-kernel), and generalized method of moments (GMM) estimation, including the continuously-updated estimator. Each dataset is introduced where it is first used. ```{r setup, message = FALSE} library(ivreg2r) library(dplyr) data(phillips) data(stockwatson) data(griliches) data(abdata) data(grunfeld) data(nlswork) data(wagepan) ``` ## HAC and AC standard errors The `phillips` dataset contains `r nrow(phillips)` annual U.S. observations (`r min(phillips$year)`--`r max(phillips$year)`) on inflation and unemployment. A simple expectations-augmented Phillips curve regresses the change in inflation on the unemployment rate. We pass the data frame as-is: `ivreg2()` drops incomplete cases per specification, exactly as Stata does, so each fit uses the largest estimation sample its own variables allow. The dataset ships several pre-built lag columns that most fits below do not use, so running `na.omit()` on the whole data frame first would drop rows whose only missing values are in those unused columns, silently changing the estimation sample. Kernel-based inference is requested with `kernel`, `bw` (bandwidth), and `tvar` (the time variable). The type of variance estimator follows the `vcov` setting. With `vcov = "iid"` a kernel gives autocorrelation-consistent (AC) standard errors, which assume homoskedasticity but allow serial correlation. With `vcov = "robust"` the same kernel gives HAC standard errors, which allow both. This matches Stata's `ivreg2`, where `bw`/`kernel` alone produces AC and combining `robust` with a bandwidth produces HAC. The explicit values `vcov = "AC"` and `vcov = "HAC"` request the same two estimators by name, which is how the AC fit below is written. The following three fits use the same Phillips curve and differ only in the variance estimator: ```{r hac-compare} f_ph <- cinf ~ unem fit_iid <- ivreg2(f_ph, data = phillips, small = TRUE) fit_ac <- ivreg2(f_ph, data = phillips, vcov = "AC", kernel = "bartlett", bw = 3, tvar = "year", small = TRUE) fit_hac <- ivreg2(f_ph, data = phillips, vcov = "robust", kernel = "bartlett", bw = 3, tvar = "year", small = TRUE) bind_rows(iid = tidy(fit_iid), AC = tidy(fit_ac), HAC = tidy(fit_hac), .id = "vce") |> select(vce, term, estimate, std.error) ``` The estimate column repeats each coefficient across the three `vce` rows: the point estimates are identical, because the variance estimator does not affect them, while the standard error on `unem` shifts as the error assumptions change. The HAC fit is the specification from Stata's `ivreg2` help file, run on the full N = `r nobs(fit_hac)` sample. The Bartlett kernel applies linearly declining weights to autocovariances up to lag `bw - 1`; combined with `vcov = "robust"` this is the Newey-West (1987) HAC estimator, and the help file pairs it with the equivalent `newey cinf unem, lag(2)`. ### Choosing a kernel `ivreg2r` supports eight kernels, matching Stata's `ivreg2`: | Kernel | Weights | Guaranteed PSD? | Typical use | |--------|---------|-----------------|-------------| | Bartlett | Linear decline | Yes | Default choice; Newey-West HAC | | Parzen | Piecewise cubic | Yes | Smoother than Bartlett | | Quadratic Spectral | $\sin/\cos$ based | Yes | Andrews (1991) | | Truncated | Uniform up to lag | No | All lags equally weighted; Kiefer equivalence | | Tukey-Hanning | Cosine taper | No | Standard spectral window | | Tukey-Hamming | Cosine taper | No | Variant of Hanning | | Daniell | $\sin(x)/x$ | No | Spectral analysis | | Tent | Tent-shaped spectral window | No | Spectral analysis | Bartlett, Parzen, and Quadratic Spectral guarantee a positive semi-definite covariance estimate (Andrews 1991). For the other kernels, use the `psd` option if negative eigenvalues arise (see the PSD corrections section below). ### Automatic bandwidth Use `bw = "auto"` for the Newey-West (1994) plug-in bandwidth selector, available for the Bartlett, Parzen, and Quadratic Spectral kernels. Pairing it with the Quadratic Spectral kernel reproduces a help-file example: ```{r hac-auto} fit_auto <- ivreg2(cinf ~ unem, data = phillips, kernel = "qs", bw = "auto", tvar = "year") fit_auto$bw ``` ### Kernel-based VCE with IV Kernel-based standard errors work with every estimation method. Instrumenting unemployment with its own first three lags gives an overidentified model with two overidentifying restrictions. The formula's first part is `~ 1`, which means no exogenous regressors beyond the constant, so the whole specification is Stata's `ivreg2 cinf (unem = L.unem L2.unem L3.unem), ...`: ```{r hac-iv} f_ph_iv <- cinf ~ 1 | unem | l(unem, 1) + l(unem, 2) + l(unem, 3) fit_iv_hac <- ivreg2(f_ph_iv, data = phillips, tvar = "year", vcov = "robust", kernel = "bartlett", bw = 3) summary(fit_iv_hac) ``` ```{r hac-iv-values, include = FALSE} overid_iv_hac <- diagnostics(fit_iv_hac) |> filter(test == "overid") ``` The diagnostics adapt automatically to the kernel-based variance: the identification tests are Kleibergen-Paap rk statistics and the overidentification test is a kernel-based Hansen J, here `r round(overid_iv_hac$statistic, 3)` on `r overid_iv_hac$df` degrees of freedom (p = `r round(overid_iv_hac$p_value, 3)`) — a rejection at the 5% level, and a useful caution: a variable's own lags are not automatically valid instruments when the errors are serially correlated, since the same persistence that motivates the HAC variance can reach the moment conditions themselves. The Stock-Yogo critical values printed alongside the weak-identification statistic are tabulated for the Cragg-Donald statistic under i.i.d. errors; Baum, Schaffer & Stillman (2007, p. 490) advise applying them to the Kleibergen-Paap rk Wald F only with caution. The Kleibergen-Paap identification tests always use the Bartlett kernel internally regardless of the kernel chosen for the variance, matching Stata's `ivreg2`. ### Lags and differences in the formula The instruments above were written with `l()`, one of two formula operators that build lags and differences on the fly against the time variable, mirroring Stata's `tsset`-aware `L.` and `D.` operators. `l(x, k)` is the k-th lag (Stata's `Lk.x`), and `d(x, k)` is the k-th order iterated difference $(1-L)^k x$ (Stata's `Dk.x`, so `d(x, 2)` is $x - 2\,l(x,1) + l(x,2)$). A range such as `l(unem, 1:3)` expands to three separate terms. See `?"ts-operators"` for the full reference. Lags are computed by time *value*, not row position: `l(x, k)` at time `t` is the value of `x` at time `t - k` within the same panel unit, and if no observation exists at `t - k` — the start of a series, or a gap in the time variable — the lag is `NA` and the row leaves the estimation sample. This is why `ivreg2()` warns when the time variable has gaps: the warning is informational, not an error, and the lags remain correct at their true temporal distance. You can skip this next step: the `l()` operator above already handles gaps, and nothing later in the vignette depends on it. When it is more convenient to build lag columns before fitting, the gap-safe idiom is to fill the calendar with `tidyr::complete()` first, so that `dplyr::lag()` steps by one period rather than by one row: ```{r ts-complete, message = FALSE, eval = requireNamespace("tidyr", quietly = TRUE)} library(tidyr) phillips_lagged <- phillips |> complete(year = full_seq(year, 1)) |> arrange(year) |> mutate(unem_lag1 = lag(unem)) head(phillips_lagged[, c("year", "unem", "unem_lag1")], 3) ``` On these annual data with no gaps, `l(unem, 1)` and this pre-computed `unem_lag1` coincide; on panels with gaps, filling the calendar first is what keeps a positional `lag()` aligned to calendar time. (`tidyr` is an optional Suggests, needed only for the `complete()` call in this pre-computation idiom; the `l()` and `d()` formula operators require neither it nor `dplyr`.) ## Panel variance estimators Panel data introduce richer correlation structures than a single time series. Beyond one-way and two-way clustering (covered in `vignette("introduction")` and `vignette("advanced-iv")`), `ivreg2r` provides three panel variance estimators, organized here by the asymptotic regime each one suits. The help file groups its own panel examples the same way. ### Large N, small T: the Arellano-Bond firm panel The `abdata` panel follows `r length(unique(abdata$id))` UK firms over `r min(abdata$year)`--`r max(abdata$year)` — many units, few periods. The help file's illustration for this regime estimates an employment equation by two-step GMM with the regressors instrumented by their own first and second differences (built here with `d()`), clustering on the firm with `clusters = ~ id`, which is Stata's `cluster(id)`: ```{r abdata-gmm} fit_ab <- ivreg2( n ~ 1 | w + k + ys | d(w, 1) + d(k, 1) + d(ys, 1) + d(w, 2) + d(k, 2) + d(ys, 2), data = abdata, tvar = "year", ivar = "id", method = "gmm2s", clusters = ~ id ) summary(fit_ab) ``` ```{r abdata-gmm-values, include = FALSE} overid_ab <- diagnostics(fit_ab) |> filter(test == "overid") ``` Clustering on the firm makes the moment covariance robust to arbitrary serial correlation within a firm's history, which is the correlation of concern when the panel is short. The differenced instruments cost the first two periods of each firm, so the estimation sample is N = `r nobs(fit_ab)` across `r fit_ab$n_clusters` firms, and the Hansen J of `r round(overid_ab$statistic, 3)` (p = `r round(overid_ab$p_value, 3)`) does not reject the overidentifying restrictions. Two-step GMM is treated in full in the GMM section below. Kiefer (1980) offers a homoskedastic alternative for the same regime. It assumes conditional homoskedasticity across all observations but allows arbitrary within-unit serial correlation, and is equivalent to AC inference with a truncated kernel spanning the full panel length $T$, so that every available within-panel lag receives full weight. The specification here drops down to the simpler uninstrumented `n ~ w + k` to isolate the variance estimator on a plain regression; that simplification is for exposition, not a requirement of `kiefer`, which works with any specification: ```{r kiefer} fit_kiefer <- ivreg2(n ~ w + k, data = abdata, kiefer = TRUE, tvar = "year", ivar = "id") summary(fit_kiefer) ``` The equivalence is exact. `kiefer = TRUE` sets the bandwidth to the full span $T$ internally, just as Stata's `kiefer` option does, while an *explicitly* supplied bandwidth must be less than $T$ in both programs, making $T - 1 = `r length(unique(abdata$year)) - 1`$ the largest explicit choice here. The two nonetheless coincide, because the truncated kernel gives every lag up to the bandwidth full weight and the largest lag a panel of length $T$ contains is $T - 1$ — so bandwidths $T$ and $T - 1$ select exactly the same lags: ```{r kiefer-equiv} bw_max <- length(unique(abdata$year)) - 1 fit_tru8 <- ivreg2(n ~ w + k, data = abdata, kernel = "truncated", bw = bw_max, tvar = "year", ivar = "id") all.equal(vcov(fit_kiefer), vcov(fit_tru8)) ``` Use Kiefer when conditional homoskedasticity is plausible but within-unit serial correlation of arbitrary form cannot be ruled out; if heteroskedasticity is also a concern, use Driscoll-Kraay or cluster-plus-kernel instead. The heteroskedastic counterpart of that equivalence also holds: clustering on the unit with robust standard errors is the same as a robust truncated kernel at the maximum bandwidth, since both weight every within-unit lag equally: ```{r cluster-kernel-equiv} fit_cl_rob <- ivreg2(n ~ w + k, data = abdata, clusters = ~ id, vcov = "robust") fit_tru_rob <- ivreg2(n ~ w + k, data = abdata, kernel = "truncated", bw = bw_max, vcov = "robust", tvar = "year", ivar = "id") all.equal(vcov(fit_cl_rob), vcov(fit_tru_rob)) ``` ### Small N, large T: the Grunfeld investment panel Driscoll and Kraay (1998) handle cross-sectional dependence together with persistence. The estimator clusters on the time variable and applies a kernel to scores summed within each period, so it is robust to common shocks — business cycles, policy changes — that move all units at once. It is a large-$T$ estimator and needs a reasonably long time-series dimension, so we switch to the `grunfeld` investment panel (`r length(unique(grunfeld$company))` firms over `r min(grunfeld$year)`--`r max(grunfeld$year)`): ```{r dk} fit_dk <- ivreg2(invest ~ mvalue + kstock, data = grunfeld, dkraay = 2, small = TRUE, tvar = "year", ivar = "company") summary(fit_dk) ``` The `dkraay` argument sets the bandwidth: here Bartlett with bandwidth 2, i.e. autocovariances up to one lag. Driscoll-Kraay is equivalent to clustering on the time variable with a kernel (Bartlett by default, overridable through `kernel`), an equivalence the help file demonstrates for exactly this model: ```{r dk-equiv} fit_dk_ck <- ivreg2(invest ~ mvalue + kstock, data = grunfeld, clusters = ~ year, kernel = "bartlett", bw = 2, small = TRUE, tvar = "year", ivar = "company") all.equal(vcov(fit_dk), vcov(fit_dk_ck)) ``` Even $T = `r length(unique(grunfeld$year))`$ is not very large, so Driscoll-Kraay results on panels of this length should be read with caution. ### Large N, moderate T: the NLS young-women panel Two-way cluster-plus-kernel accounts for both within-panel correlation (clustering on the unit) and autocorrelated cross-panel shocks (a kernel on time-aggregated scores). Specify both `clusters` and `kernel`. We use the `nlswork` panel of `r length(unique(nlswork$idcode))` women over `r length(unique(nlswork$year))` survey years: ```{r ck} fit_ck <- ivreg2( ln_wage ~ grade + age + ttl_exp + tenure, data = nlswork, clusters = ~ idcode + year, kernel = "truncated", bw = 2, tvar = "year", ivar = "idcode" ) summary(fit_ck) ``` The gap warning is expected: the NLS interviews were not conducted every calendar year, so most women have gaps in their `year` sequence, and the kernel steps by calendar time rather than by row. The variance decomposition is $V_{CK} = V_{cluster,i} + V_{kernel,t} - V_{HAC}$, where the third term is the observation-level HAC variance at the individual-time intersection, following the two-way structure of Cameron, Gelbach & Miller (2011). With the truncated kernel shown here, this construction is exactly Thompson's (2011) estimator: two-way cluster-plus-kernel with a truncated kernel whose bandwidth is the number of periods after which the common cross-panel shocks can be ignored. With only `r length(unique(nlswork$year))` survey years the time dimension is modest for the large-$T$ side of this asymptotics, so the kernel term should be read with the same caution as Driscoll-Kraay above. ### Fixed-effects panels: the Stock-Watson correction For a within-transformed fixed-effects panel regression, Stock & Watson (2008) derive a bias-adjusted heteroskedasticity-robust variance that corrects the degrees-of-freedom bias of the usual cluster-robust estimator when $T$ is fixed and $N$ is large. Request it with `sw = TRUE` and an `ivar`. Because this regime is $N \to \infty$ with fixed $T$, a large-$N$, short-$T$ panel is the appropriate home; we use `wagepan` (`r length(unique(wagepan$nr))` men over `r length(unique(wagepan$year))` years, 1980--1987). The correction is defined for the within-transformed regression, so we first subtract each man's own mean from the dependent variable and from every regressor — the within transformation that Stata's `xtreg, fe` and `xtivreg2, fe` apply internally — and then regress the demeaned outcome on the demeaned regressors: ```{r sw} wagepan_fe <- wagepan |> group_by(nr) |> mutate(lwage_dm = lwage - mean(lwage), expersq_dm = expersq - mean(expersq), married_dm = married - mean(married), union_dm = union - mean(union)) |> ungroup() fit_sw <- ivreg2(lwage_dm ~ 0 + expersq_dm + married_dm + union_dm, data = wagepan_fe, sw = TRUE, ivar = "nr", tvar = "year", dofminus = 545, vcov = "robust") summary(fit_sw) ``` The `~ 0 +` in the formula is Stata's `noconstant`: the within transformation has already removed every man's mean, so no constant remains to estimate. The `dofminus = 545` argument is Stata's `dofminus()`, which charges the `r length(unique(wagepan$nr))` absorbed person means to the residual degrees of freedom; without it the variance would treat the demeaned data as though no means had been estimated and would understate the standard errors. Experience itself (`exper`) is deliberately absent. Within a man it advances one year per calendar year, so after demeaning its variation is a common time trend shared by everyone observed in the same calendar year. The classic `wagepan` equation (Wooldridge) includes year dummies alongside the person effects, and the person effects plus year dummies absorb `exper` exactly, so that specification must drop it; we keep its conventional regressor set — `expersq`, `married`, and `union` — here as well. The `sw` option forces robust standard errors and is incompatible with clustering, kernels, Kiefer, and Driscoll-Kraay; we pass `vcov = "robust"` explicitly above rather than relying on the automatic override (leaving `vcov` at its default `"iid"` would trigger the same forced override, with a warning noting it). It appears in neither the help file nor the Baum, Schaffer & Stillman papers, and Stata's `ivreg2` labels it a beta feature; the correction itself is standard (Stock & Watson 2008). It is unrelated to the Sanderson-Windmeijer "SW" first-stage F statistic. ### When to use which | VCE | Cross-sectional heteroskedasticity | Within-panel autocorrelation | Cross-panel dependence | |-----|-------------------------------------|-------------------------------|------------------------| | Kiefer | No* | Yes | No | | Stock-Watson | Yes | No | No | | Driscoll-Kraay | Yes | Yes (kernel) | Yes | | Cluster-plus-kernel | Yes | Yes (clustering) | Yes (kernel) | | Two-way cluster | Yes | Yes (within groups) | Yes (within groups) | *Kiefer assumes conditional homoskedasticity across all observations, both across individuals and across time. ## Generalized method of moments ### When GMM helps Two-step efficient GMM (`method = "gmm2s"`) uses an optimal weighting matrix to exploit overidentifying restrictions more efficiently than 2SLS under non-i.i.d. errors. When the model is exactly identified, GMM and 2SLS produce the same coefficients. Two finite-sample cautions apply before reaching for it by default. The optimal weighting matrix is a function of fourth moments, which are hard to estimate well, so efficient-GMM Wald tests tend to over-reject in finite samples (Baum, Schaffer & Stillman 2003, §2.7, citing Hayashi 2000, p. 215). Performance also degrades as the instrument count grows relative to the sample (Baum, Schaffer & Stillman 2007, p. 488). In small samples with modest overidentification, robust 2SLS may be preferable. ### A HAC-GMM arc on the Stock-Watson data We follow the worked example of Baum, Schaffer & Stillman (2007, p. 476): a quarterly U.S. Phillips-curve equation from Stock & Watson's textbook dataset (`r nrow(stockwatson)` observations, `r min(stockwatson$year)`--`r max(stockwatson$year)`), regressing the change in inflation on the unemployment rate, instrumented by lags of GDP growth, the T-bill rate, the exchange rate, and the T-bond rate. The errors are serially correlated, so the GMM weighting matrix is HAC (Bartlett, bandwidth 5). In Stata terms the specification is `ivreg2 dinf (UR = ggdp_2 TBILL_1 ER_1 TBON_1)`, with `~ 1` again meaning no exogenous regressors beyond the constant. First the 2SLS baseline, then two-step GMM: ```{r sw-2sls} sw_formula <- dinf ~ 1 | UR | ggdp_2 + TBILL_1 + ER_1 + TBON_1 fit_sw_2sls <- ivreg2(sw_formula, data = stockwatson) ur_2sls <- tidy(fit_sw_2sls) |> filter(term == "UR") ur_2sls ``` ```{r gmm2s} fit_gmm <- ivreg2(sw_formula, data = stockwatson, method = "gmm2s", vcov = "robust", kernel = "bartlett", bw = 5, tvar = "date") ur_gmm <- tidy(fit_gmm) |> filter(term == "UR") summary(fit_gmm) ``` The first step estimates the model by 2SLS and forms the optimal weighting matrix from the residuals; the second step re-estimates with that weighting, and the Hansen (1982) J statistic tests the overidentifying restrictions. Following the published narrative, the i.i.d. 2SLS estimate of the unemployment coefficient is `r sprintf("%.3f", ur_2sls$estimate)` and clearly significant, while the HAC-GMM estimate falls to `r sprintf("%.3f", ur_gmm$estimate)` (p = `r round(ur_gmm$p.value, 3)`): once the serial correlation is accounted for, the apparent significance disappears. The continuously-updated estimator (CUE) completes the arc. CUE jointly optimizes the GMM objective over the parameter vector and the weighting matrix, and is the GMM generalization of LIML to non-i.i.d. errors. Continuing the same equation: ```{r cue} fit_cue <- ivreg2(sw_formula, data = stockwatson, method = "cue", vcov = "robust", kernel = "bartlett", bw = 5, tvar = "date") ur_cue <- tidy(fit_cue) |> filter(term == "UR") summary(fit_cue) ``` ```{r cue-values, include = FALSE} overid_cue <- diagnostics(fit_cue) |> filter(test == "overid") ``` CUE moves the coefficient further in the same direction, to `r sprintf("%.3f", ur_cue$estimate)` (p = `r round(ur_cue$p.value, 3)`), with the Hansen J at `r round(overid_cue$statistic, 3)` (p = `r round(overid_cue$p_value, 3)`) leaving the overidentifying restrictions comfortably unrejected. Like LIML, CUE has no finite moments of any order regardless of the degree of overidentification, and its finite-sample distribution can have very heavy tails (Hansen, Heaton & Yaron 1996). Its objective can have multiple local optima, and Baum, Schaffer & Stillman (2007, p. 479) warn that CUE estimation can be slow and problematic when the number of parameters is substantial. `ivreg2r` uses BFGS optimization with Nelder-Mead refinement, starting from the 2SLS estimates under i.i.d. weighting and from the two-step GMM estimates under a robust, cluster, or kernel variance, matching Stata's documented behavior. On ill-conditioned models the objective can possess economically degenerate minima with negative R-squared, a documented property of CUE itself (Hausman, Lewis, Menzel & Newey 2011); when it does, both Stata and `ivreg2r` converge to the same degenerate solution, so inspect the J statistic and R-squared before trusting a CUE fit. ### Reusing the weighting matrix GMM works with every variance type. On the Griliches (1976) wage equation — `r nrow(griliches)` young men from the U.S. National Longitudinal Survey, the help file's own heteroskedastic-GMM illustration following Hayashi (2000, p. 255) — two-step GMM with heteroskedasticity-robust weighting estimates the return to measured ability (`iq`) instrumented by test scores and family background. Here `factor(year)` expands to a set of year dummies, R's spelling of Stata's `i.year`: ```{r gril-gmm} gril_formula <- lw ~ s + expr + tenure + rns + smsa + factor(year) | iq | med + kww + age + mrt fit_g_gmm <- ivreg2(gril_formula, data = griliches, method = "gmm2s", vcov = "robust") tidy(fit_g_gmm) |> filter(term == "iq") ``` Every fit stores the estimated moment-condition covariance as `fit$S`, and every fit with a defined weighting matrix — all but LIML and k-class, exactly as in Stata — stores it as `fit$W`, matching `e(S)` and `e(W)`. This makes the help file's weighting-matrix reuse workflow direct. A robust 2SLS fit of the same equation computes the very moment covariance the two-step estimator uses in its second step, so passing that `S` back through `smatrix` reproduces the GMM fit exactly: ```{r gril-reuse} fit_g_2sls <- ivreg2(gril_formula, data = griliches, vcov = "robust") refit <- ivreg2(gril_formula, data = griliches, method = "gmm2s", vcov = "robust", smatrix = fit_g_2sls$S) all.equal(fit_g_2sls$S, fit_g_gmm$S) all.equal(coef(refit), coef(fit_g_gmm)) ``` Supplying `wmatrix` instead sets the *first-step* weighting matrix rather than fixing the second-step `S`, so the two-step estimator still recomputes `S` from the first-step residuals and the result need not match. For the reuse workflow, `smatrix` is the argument to use. For two-step GMM and CUE fits, `fit$W` is simply the inverse of the defining `fit$S`. ### Centering the moment conditions The `center` option subtracts the mean of the moment conditions before forming the covariance matrix `S` that feeds the variance and the diagnostics, which can improve the finite-sample behavior of the standard errors and the J test. No canonical example exists — the option postdates the help-file examples and both Baum-Schaffer-Stillman papers — so we show its effect directly on the Griliches GMM fit, reading the Hansen J from each fit's `diagnostics()` table (the package's one-row-per-test accessor, the diagnostic counterpart to `tidy()` for coefficients): ```{r gril-center} fit_g_cen <- ivreg2(gril_formula, data = griliches, method = "gmm2s", vcov = "robust", center = TRUE) bind_rows(uncentered = tidy(fit_g_gmm), centered = tidy(fit_g_cen), .id = "fit") |> filter(term == "iq") |> select(fit, term, estimate, std.error) overid_g_gmm <- diagnostics(fit_g_gmm) |> filter(test == "overid") overid_g_cen <- diagnostics(fit_g_cen) |> filter(test == "overid") overid_g_gmm overid_g_cen coef_shift <- max(abs(coef(fit_g_cen) - coef(fit_g_gmm))) coef_shift ``` Because two-step GMM uses `S` as its second-step weighting matrix, centering changes the coefficients as well as the inference — here by at most `r signif(coef_shift, 2)` — and shifts the J statistic from `r round(overid_g_gmm$statistic, 1)` to `r round(overid_g_cen$statistic, 1)`, with a slight change in the standard errors on the `iq` row shown above. For estimators whose coefficients do not depend on `S`, such as 2SLS with a centered robust variance, only the standard errors and the J statistic move. ## Partialling out regressors The `partial` argument removes specified exogenous regressors from every variable — dependent, endogenous, and excluded instruments — by Frisch-Waugh-Lovell projection before estimation, which absorbs high-dimensional fixed effects without carrying them as reported coefficients. The coefficients on partialled-out variables are not reported, but the remaining coefficients are numerically identical to including those variables directly, for 2SLS, LIML, and two-step GMM (the identity does not extend to CUE or to iterated GMM beyond two steps). Partialling also has a substantive role that the Griliches wage equation illustrates. Clustering that equation on `year` gives far fewer clusters than moment conditions, so the cluster-robust moment covariance is singular and the overidentification statistics cannot be computed. The warnings are the point, not a malfunction: they report exactly the rank deficiency that Baum, Schaffer & Stillman (2007, pp. 484--485) describe. This equation instruments `iq` with `med`, `kww`, and `age`, without the `mrt` term used in the GMM example above. Partialling out the year dummies leaves the reported coefficients unchanged but, as the counts after the fit show, does not by itself repair the deficiency: ```{r partial-fwl, warning = TRUE} gril_cl <- lw ~ s + expr + tenure + rns + smsa + factor(year) | iq | med + kww + age fit_full <- ivreg2(gril_cl, data = griliches, clusters = ~ year) fit_partial <- ivreg2(gril_cl, data = griliches, clusters = ~ year, partial = "factor(year)") shared <- c("s", "expr", "tenure", "rns", "smsa", "iq") all.equal(coef(fit_full)[shared], coef(fit_partial)[shared]) ``` The two fits agree on the shared coefficients to floating-point precision, confirming the projection identity even though both leave the Hansen J uncomputed. The counts explain why partialling alone does not help: the full fit has `r ncol(fit_full$S)` moment conditions against `r fit_full$n_clusters` clusters, and the projection — which absorbs the year dummies and the constant from both the regressor and instrument sets — still leaves `r ncol(fit_partial$S)` moment conditions against the same `r fit_partial$n_clusters` clusters. Two-step GMM is stricter than 2SLS here: it must invert the moment covariance, so with more moment conditions than clusters it fails outright, reproducing Stata's `r(506)` error rather than proceeding with a singular weighting matrix: ```{r partial-gmm-fail, error = TRUE} ivreg2(gril_cl, data = griliches, clusters = ~ year, partial = "factor(year)", method = "gmm2s") ``` Partialling out one further regressor brings the instrument count down to the number of clusters, and two-step GMM becomes feasible: ```{r partial-gmm-ok} fit_g_feas <- ivreg2(gril_cl, data = griliches, clusters = ~ year, partial = c("factor(year)", "rns"), method = "gmm2s") tidy(fit_g_feas) |> filter(term == "iq") diagnostics(fit_g_feas) |> filter(test == "overid") ``` With the instrument count now matched to the `r fit_g_feas$n_clusters` year clusters, the Hansen J is computed and the `iq` coefficient is reported. The special tokens `"_cons"` (partial out the constant, equivalent to demeaning every variable) and `"_all"` (partial out all exogenous regressors) are also accepted. Partialling adjusts `sdofminus` for the degrees of freedom it consumes; `nopartialsmall = TRUE` suppresses that adjustment. When fixed effects have been partialled out *manually*, outside `ivreg2()`, the `dofminus` and `sdofminus` arguments make the same adjustment, with a self-verifying demonstration in the degrees-of-freedom section of `vignette("advanced-iv")`. ## PSD corrections Kernels outside the Bartlett/Parzen/Quadratic-Spectral family do not guarantee a positive semi-definite estimate of `S`, the covariance of the orthogonality conditions that feeds both the variance and the chi-squared diagnostics. When `S` has negative eigenvalues, `psd` applies an eigenvalue correction: `psd = "psd0"` sets negative eigenvalues to zero (Politis 2007), and `psd = "psda"` replaces them with their absolute values (Stock & Watson 2008, Remark 8). Driscoll-Kraay with its default Bartlett kernel is positive semi-definite by construction, but the truncated kernel carries no such guarantee. On `wagepan`, Driscoll-Kraay with a truncated kernel produces an indefinite `S` — a genuinely binding case, and here `S` is also singular, so the Hansen J and Stock-Wright (2000) S statistics cannot be computed and are reported as `NA`, exactly as Stata leaves them blank: ```{r psd-binding, warning = TRUE} psd_formula <- lwage ~ exper + expersq + married + union | hours | educ + black fit_nopsd <- ivreg2(psd_formula, data = wagepan, dkraay = 2, kernel = "truncated", tvar = "year", ivar = "nr") min(eigen(fit_nopsd$S, symmetric = TRUE)$values) ``` Refitting with `psd = "psda"` repairs `S`; a warning reports each corrected matrix, where Stata performs the same correction silently: ```{r psd-corrected, warning = TRUE} fit_psda <- ivreg2(psd_formula, data = wagepan, dkraay = 2, kernel = "truncated", tvar = "year", ivar = "nr", psd = "psda") min(eigen(fit_psda$S, symmetric = TRUE)$values) ``` The correction is applied to `S` before the coefficient covariance is assembled from it, exactly as in Stata, so the standard errors change while the coefficients do not: ```{r psd-compare} bind_rows(uncorrected = tidy(fit_nopsd), psda = tidy(fit_psda), .id = "fit") |> select(fit, term, estimate, std.error) ``` Two quantities are deliberately untouched, again matching Stata: the conventional `vcov = "iid"` variance, in which `S` plays no role, and the Kleibergen-Paap identification statistics, which never receive the `psd` option. ## Suppressing identification statistics For simulation loops where the identification diagnostics are not needed, `noid = TRUE` suppresses the underidentification and weak-identification statistics (and the redundancy test), saving the rank computations they require: ```{r noid} fit_noid <- ivreg2(f_ph_iv, data = phillips, tvar = "year", noid = TRUE) diagnostics(fit_noid) ``` The `underid`, `weak_id`, and Stock-Yogo rows are absent from the returned tibble — the whole identification block is dropped from `summary()` and `glance()` as well — while the `overid` test and the other non-identification diagnostics are still reported. ## Estimator equivalences Several estimators coincide under specific conditions, which is useful for verification. | Equivalence | Conditions | Where shown | |-------------|------------|-------------| | Driscoll-Kraay = cluster(time) + kernel | Same kernel and bandwidth | In code above (grunfeld) | | Kiefer = truncated kernel | Panel data, AC, `bw = T - 1` (all within-panel lags weighted one) | In code above (abdata) | | Cluster-robust = truncated kernel, robust | Panel data, `bw = T - 1`, `vcov = "robust"` | In code above (abdata) | | LIML + COVIV = CUE (COVIV: `coviv = TRUE`, the 2SLS-bread covariance for LIML) | Conditional-homoskedasticity (i.i.d.) weighting; agreement up to optimizer tolerance | Demonstrated in `vignette("advanced-iv")` | | LIML = 2SLS | Exactly identified | Stated in `vignette("advanced-iv")` | | GMM2S = 2SLS | Exactly identified | Stated in the GMM section above | ## References - Andrews, D.W.K. (1991). "Heteroskedasticity and Autocorrelation Consistent Covariance Matrix Estimation." *Econometrica*, 59(3), 817--858. - Baum, C.F., Schaffer, M.E., and Stillman, S. (2003). "Instrumental Variables and GMM: Estimation and Testing." *Stata Journal*, 3(1), 1--31. - Baum, C.F., Schaffer, M.E., and Stillman, S. (2007). "Enhanced Routines for Instrumental Variables/Generalized Method of Moments Estimation and Testing." *Stata Journal*, 7(4), 465--506. - Cameron, A.C., Gelbach, J.B., and Miller, D.L. (2011). "Robust Inference with Multiway Clustering." *Journal of Business & Economic Statistics*, 29(2), 238--249. - Driscoll, J.C. and Kraay, A.C. (1998). "Consistent Covariance Matrix Estimation with Spatially Dependent Panel Data." *Review of Economics and Statistics*, 80(4), 549--560. - Griliches, Z. (1976). "Wages of Very Young Men." *Journal of Political Economy*, 84(4, Part 2), S69--S85. - Hansen, L.P. (1982). "Large Sample Properties of Generalized Method of Moments Estimators." *Econometrica*, 50(4), 1029--1054. - Hansen, L.P., Heaton, J., and Yaron, A. (1996). "Finite-Sample Properties of Some Alternative GMM Estimators." *Journal of Business & Economic Statistics*, 14(3), 262--280. - Hausman, J., Lewis, R., Menzel, K., and Newey, W. (2011). "Properties of the CUE Estimator and a Modification with Moments." *Journal of Econometrics*, 165(1), 45--57. - Hayashi, F. (2000). *Econometrics*. Princeton University Press. - Kiefer, N.M. (1980). "Estimation of Fixed Effect Models for Time Series of Cross-Sections with Arbitrary Intertemporal Covariance." *Journal of Econometrics*, 14(2), 195--202. - Kleibergen, F. and Paap, R. (2006). "Generalized Reduced Rank Tests Using the Singular Value Decomposition." *Journal of Econometrics*, 133(1), 97--126. - Newey, W.K. and West, K.D. (1987). "A Simple, Positive Semi-Definite, Heteroskedasticity and Autocorrelation Consistent Covariance Matrix." *Econometrica*, 55(3), 703--708. - Newey, W.K. and West, K.D. (1994). "Automatic Lag Selection in Covariance Matrix Estimation." *Review of Economic Studies*, 61(4), 631--653. - Politis, D.N. (2007). "Higher-Order Accurate, Positive Semi-Definite Estimation of Large-Sample Covariance and Spectral Density Matrices." Working paper, UC San Diego; published (2011) in *Econometric Theory*, 27(4), 703--744. - Stock, J.H. and Watson, M.W. (2008). "Heteroskedasticity-Robust Standard Errors for Fixed Effects Panel Data Regression." *Econometrica*, 76(1), 155--174. - Stock, J.H. and Wright, J.H. (2000). "GMM with Weak Identification." *Econometrica*, 68(5), 1055--1096. - 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: Essays in Honor of Thomas Rothenberg*. Cambridge University Press. - Thompson, S.B. (2011). "Simple Formulas for Standard Errors that Cluster by Both Firm and Time." *Journal of Financial Economics*, 99(1), 1--10.