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

## ----pkg, message = FALSE-----------------------------------------------------
library(infometrics)

## ----dgp----------------------------------------------------------------------
set.seed(10)
n  <- 200
x1 <- runif(n, 0, 20)                 # clean regressor
x2 <- runif(n, 0, 20)
x3 <- x2 + rnorm(n, 0, 0.05)          # near-collinear with x2 (extreme)

b  <- c(1, -2, 3, 0)                  # (Intercept), x1, x2, x3
y  <- b[1] + b[2] * x1 + b[3] * x2 + b[4] * x3 + rnorm(n, 0, sqrt(2))
dat <- data.frame(y, x1, x2, x3)

truth <- c(`(Intercept)` = 1, x1 = -2, x2 = 3, x3 = 0)
cor(x2, x3)                           # near-perfect collinearity

## ----fit1---------------------------------------------------------------------
ols <- lm(y ~ x1 + x2 + x3, data = dat)

Zvec <- seq(-6, 6, length.out = 5)   # coefficient support (shared)
sig  <- sd(dat$y)                    # standard deviation of y
k    <- max(3, sqrt(2 * log(n)))     # sample-size-aware multiple
vvec <- seq(-k * sig, k * sig, length.out = 5)

gme <- linreg(y ~ x1 + x2 + x3, data = dat,
              Z = Zvec, v = vvec)             # nu = 0.5, uniform priors => GME

## ----table1-------------------------------------------------------------------
se <- function(fit) sqrt(diag(vcov(fit)))
tab1 <- data.frame(
  Truth    = truth,
  OLS      = coef(ols),
  `OLS SE` = se(ols),
  GME      = coef(gme),
  `GME SE` = se(gme),
  check.names = FALSE
)
knitr::kable(tab1, digits = 3,
             caption = "OLS versus GME under near-perfect collinearity of x2 and x3.")

## ----prior--------------------------------------------------------------------
# we can recover probabilities associated with the  
# prior information, using the inverse_ce() function  
temp <- data.frame(y=2.8, s1 = Zvec[1], s2 = Zvec[2], 
                   s3= Zvec[3], s4= Zvec[4], s5= Zvec[5])
p0_x2 <- inverse_ce(y ~ s1 + s2 + s3 + s4 + s5 - 1, data=temp)$p_hat
p0_other_x <- rep(1/5, 5)
p0 <- rbind(p0_other_x, p0_other_x, p0_x2, p0_other_x)

gce <- linreg(y ~ x1 + x2 + x3, data = dat,
              Z = Zvec, v = vvec, p0 = p0)    # GCE: informative signal prior

## ----table2-------------------------------------------------------------------
tab2 <- data.frame(
  Truth          = truth,
  `OLS`      = coef(ols),
  `GME` = coef(gme),
  `GCE`   = coef(gce),
  check.names     = FALSE
)
knitr::kable(tab2, digits = 3,
             caption = "OLS, GME (uniform prior), and GCE (informative prior).")

## ----nu-----------------------------------------------------------------------
nus <- c(0.1, 0.5, 0.9)
fit_nu <- lapply(nus, function(nu)
  linreg(y ~ x1 + x2 + x3, data = dat, Z = Zvec, v = vvec, p0=p0, nu = nu))

tab3 <- rbind(sapply(fit_nu, coef),
              `R-squared` = sapply(fit_nu, `[[`, "r.squared"),
              `Signal S`  = sapply(fit_nu, `[[`, "S"))
tab3 <- cbind(Truth = c(truth, `R-squared` = NA, `Signal S` = NA), round(tab3, 4))
colnames(tab3) <- c("Truth", paste0("nu = ", nus))
knitr::kable(tab3, digits = 3,
             caption = "Coefficients, ordinary R-squared, and signal normalized entropy across nu.")

## ----foc-ok-------------------------------------------------------------------
c(foc_residual = gme$foc_residual, converged = gme$converged)

## ----foc-bad------------------------------------------------------------------
bad <- linreg(y ~ x1 + x2 + x3, data = dat,
              Z = Zvec, v = c(-0.01, 0, 0.01))   # far too narrow
round(coef(bad), 3)
c(foc_residual = bad$foc_residual, converged = bad$converged)

## ----ent----------------------------------------------------------------------
c(signal_S    = gme$S,
  pseudo_R2   = 1 - gme$S,
  ordinary_R2 = gme$r.squared)

## ----er-----------------------------------------------------------------------
summary(gme)

## ----wald---------------------------------------------------------------------
wald_test <- function(fit, R, r = NULL) {
  R <- rbind(R)                                 # allow a plain vector
  if (is.null(r)) r <- rep(0, nrow(R))
  d <- as.vector(R %*% coef(fit)) - r
  W <- as.numeric(t(d) %*% solve(R %*% vcov(fit) %*% t(R), d))
  c(W = W, df = nrow(R), p.value = pchisq(W, df = nrow(R), lower.tail = FALSE))
}

# H0: beta_x3 = 0  (single restriction; equals the squared t-ratio)
wald_test(gme, c(0, 0, 0, 1))

# H0: beta_x2 = beta_x3 = 0  (joint test of the collinear pair)
wald_test(gme, rbind(c(0, 0, 1, 0),
                     c(0, 0, 0, 1)))

## ----iv-dgp-------------------------------------------------------------------
set.seed(125)
n  <- 200
z  <- rnorm(n)                       # instrument
u  <- rnorm(n)                       # shared shock => endogeneity
xe <- 0.7 * z + u + rnorm(n)         # endogenous regressor
y  <- 1 + 1.5 * xe + u               # true slope = 1.5

X   <- cbind(`(Intercept)` = 1, xe = xe)   # design: intercept + endogenous xe
IV  <- cbind(`(Intercept)` = 1, z  = z)    # instruments: intercept (self) + z
Ziv <- matrix(c(-10, 0, 10), nrow = 2, ncol = 3, byrow = TRUE)  # signal support

## ----iv-fit-------------------------------------------------------------------
ols_iv <- lm(y ~ xe)
b_2sls <- solve(crossprod(IV, X), crossprod(IV, y))   # just-identified 2SLS
gme_iv <- linreg_iv(y, X, IV, Ziv)

tab_iv <- data.frame(
  Truth    = c(`(Intercept)` = 1, xe = 1.5),
  OLS      = coef(ols_iv),
  `2SLS`   = as.vector(b_2sls),
  `GME-IV` = coef(gme_iv),
  check.names = FALSE
)
knitr::kable(tab_iv, digits = 3,
             caption = "OLS is biased by endogeneity; 2SLS and GME-IV recover the true slope.")

## ----iv-summary---------------------------------------------------------------
summary(gme_iv)