## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup, message = FALSE---------------------------------------------------
library(ivreg2r)
library(dplyr)
data(phillips)
data(stockwatson)
data(griliches)
data(abdata)
data(grunfeld)
data(nlswork)
data(wagepan)

## ----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)

## ----hac-auto-----------------------------------------------------------------
fit_auto <- ivreg2(cinf ~ unem, data = phillips,
                   kernel = "qs", bw = "auto", tvar = "year")
fit_auto$bw

## ----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)

## ----hac-iv-values, include = FALSE-------------------------------------------
overid_iv_hac <- diagnostics(fit_iv_hac) |> filter(test == "overid")

## ----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)

## ----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)

## ----abdata-gmm-values, include = FALSE---------------------------------------
overid_ab <- diagnostics(fit_ab) |> filter(test == "overid")

## ----kiefer-------------------------------------------------------------------
fit_kiefer <- ivreg2(n ~ w + k, data = abdata, kiefer = TRUE,
                     tvar = "year", ivar = "id")
summary(fit_kiefer)

## ----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))

## ----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))

## ----dk-----------------------------------------------------------------------
fit_dk <- ivreg2(invest ~ mvalue + kstock, data = grunfeld,
                 dkraay = 2, small = TRUE, tvar = "year", ivar = "company")
summary(fit_dk)

## ----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))

## ----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)

## ----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)

## ----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

## ----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)

## ----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)

## ----cue-values, include = FALSE----------------------------------------------
overid_cue <- diagnostics(fit_cue) |> filter(test == "overid")

## ----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")

## ----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))

## ----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

## ----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])

## ----partial-gmm-fail, error = TRUE-------------------------------------------
try({
ivreg2(gril_cl, data = griliches, clusters = ~ year,
       partial = "factor(year)", method = "gmm2s")
})

## ----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")

## ----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)

## ----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)

## ----psd-compare--------------------------------------------------------------
bind_rows(uncorrected = tidy(fit_nopsd), psda = tidy(fit_psda),
          .id = "fit") |>
  select(fit, term, estimate, std.error)

## ----noid---------------------------------------------------------------------
fit_noid <- ivreg2(f_ph_iv, data = phillips, tvar = "year", noid = TRUE)
diagnostics(fit_noid)