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

## ----setup, message = FALSE---------------------------------------------------
library(ivreg2r)
library(dplyr)
data(mroz)
data(klein)
data(griliches)
data(nlswork)

## ----mroz-subset--------------------------------------------------------------
mroz_work <- mroz |> filter(inlf == 1)
nrow(mroz_work)

## ----mroz-baseline------------------------------------------------------------
fit_2sls <- ivreg2(lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6,
                   data = mroz_work)
summary(fit_2sls)

## ----liml---------------------------------------------------------------------
fit_liml <- ivreg2(lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6,
                   data = mroz_work, method = "liml")
summary(fit_liml)

## ----fuller-------------------------------------------------------------------
fit_fuller1 <- ivreg2(lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6,
                      data = mroz_work, fuller = 1)
fit_fuller1$kclass_value  # the k-class value, Stata's e(kclass)
tidy(fit_fuller1) |> filter(term == "educ")

## ----kclass-------------------------------------------------------------------
N <- nobs(fit_2sls)
L <- fit_2sls$rankzz
K <- fit_2sls$rank
k_nagar <- 1 + (L - K) / N
fit_nagar <- ivreg2(lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6,
                    data = mroz_work, kclass = k_nagar)
c(L = L, K = K, k = k_nagar)
tidy(fit_nagar) |> filter(term == "educ")

## ----coviv--------------------------------------------------------------------
fit_coviv <- ivreg2(lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6,
                    data = mroz_work, method = "liml", vcov = "robust",
                    small = TRUE, coviv = TRUE)
all.equal(coef(fit_coviv), coef(fit_liml))
tidy(fit_coviv) |> filter(term == "educ") |> select(term, estimate, std.error)

## ----klein-formula------------------------------------------------------------
klein_form <- consump ~ l(profits, 1) | profits + wagetot |
  govt + taxnetx + year + wagegovt + capital1 + l(totinc, 1)
fit_k_liml <- ivreg2(klein_form, data = klein, tvar = "yr", method = "liml")
summary(fit_k_liml)

## ----klein-suite--------------------------------------------------------------
fit_k_2sls   <- ivreg2(klein_form, data = klein, tvar = "yr")
fit_k_fuller <- ivreg2(klein_form, data = klein, tvar = "yr",
                       method = "liml", fuller = 1)
# Nagar's k = 1 + (L - K)/N = 1 + 4/21, which the help file rounds to 1.19
fit_k_nagar  <- ivreg2(klein_form, data = klein, tvar = "yr",
                       method = "kclass", kclass = 1.19)
klein_models <- list("2SLS" = fit_k_2sls, LIML = fit_k_liml,
                     "Fuller(1)" = fit_k_fuller, Nagar = fit_k_nagar)
bind_rows("2SLS" = tidy(fit_k_2sls), LIML = tidy(fit_k_liml),
          "Fuller(1)" = tidy(fit_k_fuller), Nagar = tidy(fit_k_nagar),
          .id = "estimator") |>
  filter(term %in% c("profits", "wagetot")) |>
  select(estimator, term, estimate)

## ----klein-coviv-cue----------------------------------------------------------
fit_k_coviv <- ivreg2(klein_form, data = klein, tvar = "yr",
                      method = "liml", coviv = TRUE)
fit_k_cue   <- ivreg2(klein_form, data = klein, tvar = "yr", method = "cue")
max(abs(coef(fit_k_coviv) - coef(fit_k_cue)))

## ----sargan-------------------------------------------------------------------
glance(fit_2sls)[, c("overid_stat", "overid_p")]

## ----orthog-------------------------------------------------------------------
fit_orthog <- ivreg2(
  lw ~ s + expr + tenure + rns + smsa + factor(year) | iq |
    med + kww + age + mrt,
  data = griliches, method = "gmm2s", orthog = c("age", "mrt")
)
diagnostics(fit_orthog) |> filter(test == "orthog")

## ----endog--------------------------------------------------------------------
fit_endog <- ivreg2(lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6,
                    data = mroz_work, endog = "educ")
diagnostics(fit_endog) |> filter(test == "endogeneity")

## ----endog-identity-----------------------------------------------------------
c(automatic = diagnostics(fit_2sls) |> filter(test == "endogeneity") |> pull(statistic),
  endog_arg = diagnostics(fit_endog) |> filter(test == "endogeneity") |> pull(statistic))

## ----hols-orthog--------------------------------------------------------------
fit_hols <- ivreg2(
  lwage ~ exper + expersq + educ | 0 | age + kidslt6 + kidsge6,
  data = mroz_work, orthog = "educ"
)
summary(fit_hols)

## ----ar-overid----------------------------------------------------------------
diagnostics(fit_liml) |>
  filter(test %in% c("anderson_rubin_overid_lr", "anderson_rubin_overid_lin"))

## ----stock-wright-------------------------------------------------------------
diagnostics(fit_2sls) |> filter(test == "stock_wright")

## ----weak-instruments---------------------------------------------------------
weak_form <- lw ~ s + expr + tenure + rns + smsa + factor(year) | iq | age + mrt
fit_weak <- ivreg2(weak_form, data = griliches, vcov = "robust")
summary(fit_weak)

sy10 <- diagnostics(fit_weak) |> filter(test == "sy_iv_size_10") |> pull(statistic)
sy10

## ----weak-values, include = FALSE---------------------------------------------
kp_f_weak <- diagnostics(fit_weak) |>
  filter(test == "weak_id_robust") |>
  pull(statistic)

## ----redundant----------------------------------------------------------------
fit_redund <- ivreg2(weak_form, data = griliches, vcov = "robust",
                     redundant = "mrt")
diagnostics(fit_redund) |> filter(test == "redundancy")

## ----redund-values, include = FALSE-------------------------------------------
redund_p <- diagnostics(fit_redund) |>
  filter(test == "redundancy") |>
  pull(p_value)

## ----b0-----------------------------------------------------------------------
ols_fit <- ivreg2(lwage ~ exper + expersq, data = mroz_work)
b0 <- c(educ = 0, coef(ols_fit))
fit_b0 <- ivreg2(
  lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6,
  data = mroz_work, vcov = "robust", b0 = b0
)
diagnostics(fit_b0) |> filter(test == "overid")  # the objective value, Stata's e(j)

## ----b0-identity--------------------------------------------------------------
fit_reg <- ivreg2(
  lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6,
  data = mroz_work, vcov = "robust"
)
b0_objective   <- diagnostics(fit_b0)  |> filter(test == "overid")       |> pull(statistic)
stock_wright_S <- diagnostics(fit_reg) |> filter(test == "stock_wright") |> pull(statistic)
all.equal(b0_objective, stock_wright_S)

diagnostics(fit_reg) |> filter(test == "stock_wright")

## ----b0-values, include = FALSE-----------------------------------------------
sw_p  <- diagnostics(fit_reg) |> filter(test == "stock_wright") |> pull(p_value)
sw_df <- diagnostics(fit_reg) |> filter(test == "stock_wright") |> pull(df)

## ----nlswork-overview---------------------------------------------------------
nrow(nlswork)
length(unique(nlswork$idcode))  # persons
length(unique(nlswork$year))    # interview years

## ----twoway-------------------------------------------------------------------
fit_1way <- ivreg2(
  ln_wage ~ grade + age + ttl_exp + tenure,
  data = nlswork, clusters = ~ idcode
)

fit_2way <- ivreg2(
  ln_wage ~ grade + age + ttl_exp + tenure,
  data = nlswork, clusters = ~ idcode + year
)

## ----twoway-compare-----------------------------------------------------------
tibble(
  term    = tidy(fit_1way)$term,
  se_1way = tidy(fit_1way)$std.error,
  se_2way = tidy(fit_2way)$std.error
) |>
  mutate(ratio = se_2way / se_1way)

## ----two-way-iv---------------------------------------------------------------
data(cigar)
cigar <- cigar |>
  mutate(
    lsales  = log(sales),
    lrprice = log(price / cpi),
    lrndi   = log(ndi / cpi),
    lrpimin = log(pimin / cpi)
  )

fit_iv_2way <- ivreg2(
  lsales ~ lrndi | lrprice | lrpimin + l(lrprice, 1),
  data = cigar, tvar = "year", ivar = "state",
  clusters = ~ state + year
)
summary(fit_iv_2way)

## ----rf-----------------------------------------------------------------------
# Same robust fit as fit_reg above, re-estimated with the reduced form stored
fit_rf <- ivreg2(
  lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6,
  data = mroz_work, vcov = "robust", reduced_form = "rf"
)
fit_rf$reduced_form$coefficients

## ----rf-ar-identity-----------------------------------------------------------
rf_F <- fit_rf$reduced_form$f_stat
ar_F <- diagnostics(fit_rf) |> filter(test == "anderson_rubin_f") |> pull(statistic)
all.equal(rf_F, ar_F)
c(rf_F = rf_F, ar_F = ar_F)

## ----rf-system----------------------------------------------------------------
fit_system <- ivreg2(
  lwage ~ exper + expersq | educ | age + kidslt6 + kidsge6,
  data = mroz_work, reduced_form = "system"
)

# Equation names
fit_system$reduced_form$depvar

# Cross-equation VCV dimensions
dim(fit_system$reduced_form$vcov)

## ----partial------------------------------------------------------------------
fit_full <- ivreg2(
  lw ~ s + expr + tenure + rns + smsa + factor(year) | iq | med + kww + age,
  data = griliches, small = TRUE
)

fit_partial <- ivreg2(
  lw ~ s + expr + tenure + rns + smsa + factor(year) | iq | med + kww + age,
  data = griliches, small = TRUE, partial = "factor(year)"
)

shared <- c("s", "expr", "tenure", "rns", "smsa", "iq")
all.equal(coef(fit_full)[shared], coef(fit_partial)[shared])

## ----partial-se---------------------------------------------------------------
all.equal(
  tidy(fit_full)    |> filter(term %in% shared) |> pull(std.error),
  tidy(fit_partial) |> filter(term %in% shared) |> pull(std.error)
)

## ----dofminus-fe--------------------------------------------------------------
gril_within <- griliches |>
  group_by(year) |>
  mutate(
    lw     = lw     - mean(lw)     + mean(griliches$lw),
    s      = s      - mean(s)      + mean(griliches$s),
    expr   = expr   - mean(expr)   + mean(griliches$expr),
    tenure = tenure - mean(tenure) + mean(griliches$tenure)
  ) |>
  ungroup()
G <- length(unique(griliches$year))

fit_lsdv <- ivreg2(lw ~ s + expr + tenure + factor(year),
                   data = griliches, small = TRUE)
fit_fe <- ivreg2(lw ~ s + expr + tenure, data = gril_within,
                 small = TRUE, dofminus = G - 1)

slopes <- c("s", "expr", "tenure")
all.equal(coef(fit_lsdv)[slopes], coef(fit_fe)[slopes])
all.equal(
  tidy(fit_lsdv) |> filter(term %in% slopes) |> pull(std.error),
  tidy(fit_fe)   |> filter(term %in% slopes) |> pull(std.error)
)
c(sigma_lsdv = fit_lsdv$sigma, sigma_fe = fit_fe$sigma,
  df_lsdv = fit_lsdv$df.residual, df_fe = fit_fe$df.residual)

## ----modelsummary-setup, eval = requireNamespace("modelsummary", quietly = TRUE)----
library(modelsummary)

## ----ms-basic, eval = requireNamespace("modelsummary", quietly = TRUE)--------
modelsummary(fit_2sls, output = "markdown")

## ----ms-compare, eval = requireNamespace("modelsummary", quietly = TRUE)------
modelsummary(klein_models, output = "markdown",
             gof_map = c("nobs", "r.squared", "adj.r.squared", "sigma"))

## ----ms-exponentiate, eval = requireNamespace("modelsummary", quietly = TRUE)----
modelsummary(fit_2sls, output = "markdown", exponentiate = TRUE)

## ----ms-custom, eval = requireNamespace("modelsummary", quietly = TRUE)-------
modelsummary(
  list("2SLS" = fit_2sls, "LIML" = fit_liml),
  output = "markdown",
  coef_rename = c(educ = "Education", exper = "Experience",
                  expersq = "Experience²"),
  gof_map = c("nobs", "r.squared", "sigma", "weak_id_stat")
)