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

## ----setup--------------------------------------------------------------------
library(formulaiv)

## -----------------------------------------------------------------------------
# Set up BH market access data
y <- ma$emp_growth # N x 1 vector
x <- ma$dma0 # N x 1 vector
z <- x # N x 1 vector
controls <- ma[, c("distance_B", "scaled_lat", "scaled_lon")] # N x J dataframe
f <- ma[, paste0("ma_nlink", 1:1999)] - ma$ma2007 # N x S dataframe
pbar <- rep(1 / 1999, 1999) # S x 1 vector
# In the BH market access example, N = 275, S = 1999, J = 3

# Set up BH high speed railway data
# Generate dummies of line opening status in 2016
for (i in 1:1999) {
  line[[paste0("open2016_sim", i)]] <- as.integer(
    line[[paste0("year_operate", i)]] <= 2016
  )
}

# Generate probability of line opening across S simulations
g <- line[, paste0("open2016_sim", 1:1999)] # L x S dataframe
qbar <- as.numeric(as.matrix(g) %*% pbar) # L x 1 vector

## ----eval=FALSE---------------------------------------------------------------
# BH_sens_joint_cons_no_controls <- formulaiv(
#   y = y,
#   x = x,
#   z = z,
#   f = f,
#   eps = seq(1, 20, 0.2),
#   cons = list(name = "joint", pbar = pbar)
# )$beta

## ----eval=FALSE---------------------------------------------------------------
# BH_sens_joint_cons_with_controls <- formulaiv(
#   y = y,
#   x = x,
#   z = z,
#   f = f,
#   eps = seq(1, 20, 0.2),
#   cons = list(name = "joint", pbar = pbar),
#   controls = controls
# )$beta

## ----eval=FALSE---------------------------------------------------------------
# BH_sens_marginal_cons_no_controls <- formulaiv(
#   y = y,
#   x = x,
#   z = z,
#   f = f,
#   eps = seq(1, 2.5, 0.25),
#   cons = list(name = "marginal", g = g, qbar = qbar)
# )$beta

## ----eval=FALSE---------------------------------------------------------------
# BH_sens_marginal_cons_with_controls <- formulaiv(
#   y = y,
#   x = x,
#   z = z,
#   f = f,
#   eps = seq(1, 2.5, 0.25),
#   cons = list(name = "marginal", g = g, qbar = qbar),
#   controls = controls
# )$beta