Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Wilke Humphreys Bargaining Model
# Replication code for Wilke, Humphreys 2020, Field Experiments, Theory, and External Validity
# The analysis was run using
## R version 3.5.3 (2019-03-11)
## DeclareDesign version 0.18.0
## tidyverse version 1.2.1
## bbmle version 1.0.20
## knitr version 1.22
# Install and load packages -----------------------------------------------
packages <- c(
"DeclareDesign",
"tidyverse",
"bbmle",
"knitr")
invisible(lapply(packages, function(x) if (!require(x, character.only=T)){install.packages(x);library(x, character.only = T)}))
rm(packages)
# Set number of simulations -----------------------------------------------
sims = 1000
# Define structural estimator ---------------------------------------------
# Maximum Likelihood Estimator to estimate kappa (k), delta (d) and q:
structural_estimator <- function(data,
Y = "Y_2_obs",
# function to calculate
# predicted price for rational
# customers given z_i and delta,
# defaults to n=2:
y = function(Z, d)
Z * d + (1 - Z) * (1 - d))
{
# Define negative log likelihood as a function of kappa, delta and q
LL <- function(k, d, q) {
m <- with(data, y(Z, d))
R <- q * dbeta(data[Y][[1]], k * 3 / 4, k * 1 / 4) +
(1 - q) * dbeta(data[Y][[1]], k * m, k * (1 - m))
- sum(log(R))
}
# Estimation
M <- mle2(
LL,
method = "L-BFGS-B",
start = list(k = 2, d = 0.50, q = 0.50),
lower = c(k = 1, d = 0.01, q = 0.01),
upper = c(k = 1000, d = 0.99, q = 0.99)
)
# Format output from estimation
out <- data.frame(coef(summary(M)), outcome = Y)
names(out) <- c("estimate", "std.error", "statistic",
"p.value", "outcome")
# Use estimates of q and delta to predict average treatment effects (ATEs)
# Predicted ATE for n=2
out[4, 1] <- (1 - out["q", "estimate"]) * (2 * out["d", "estimate"] - 1)
# Predicted ATE for n=infinity
out[5, 1] <- (1 - out["q", "estimate"]) * (2 * out["d", "estimate"] /
(1 + out["d", "estimate"]) - 1)
out
}
# Declare the design ----------------------------------------------------------
# Define parameter values:
N = 500 # Sample size
d = 0.8 # True delta (unknown)
k = 6 # Nuisance parameter that affects variance (unknown)
q = 0.5 # Share of behavioral types in the population (unknown)
e = 3/4 # Price paid by behavioral customers (known)
random = TRUE # Switch to control whether
# first movers are randomly assigned
# Design declaration:
design <-
# Define the population
declare_population(N = N,
# indicator for behavioral type (norm = 1)
norm = rbinom(N, 1, q),
# probablity of going first without random assignment
p = ifelse(norm == 1, .8, .5)
) +
# Define mean potential outcomes for n = 2
declare_potential_outcomes(Y_2 ~ norm*e + (1-norm)*(Z*d + (1-Z)*(1-d))) +
# Define mean potential outcomes for n = infinity
declare_potential_outcomes(Y_inf ~ norm*e +
(1-norm)*(Z*d/(1+d) +
(1-Z)*(1-d/(1+d)))) +
# Define estimands (quantities we want to estimate)
declare_estimand(ATE_2 = mean(Y_2_Z_1 - Y_2_Z_0), # ATE n = 2
ATE_inf = mean(Y_inf_Z_1 - Y_inf_Z_0), # ATE n = infinity
k = k, # kappa
d = d, # delta
q = q) + # q
# Declare assignment process (random assignment if random = TRUE)
declare_assignment(prob_unit = if(random) rep(0.5, N) else p, simple = TRUE) +
# Declare revealed potential outcomes
declare_reveal(Y_2, Z) +
declare_reveal(Y_inf, Z) +
# Get draws from beta distribution given means for n = 2 and n = infinity
declare_step(fabricate, Y_2_obs = rbeta(N, Y_2*k, (1-Y_2)*k),
Y_inf_obs = rbeta(N, Y_inf*k, (1-Y_inf)*k)
) +
# Declare estimators
# Difference-in-means for n = 2
declare_estimator(Y_2_obs ~ Z,
estimand = "ATE_2",
label = "DIM_2") +
# Difference-in-means for n = infinity
declare_estimator(Y_inf_obs ~ Z,
estimand = "ATE_inf",
label = "DIM_inf") +
# MLE for n = 2
declare_estimator(handler = tidy_estimator(structural_estimator),
estimand = c("k","d", "q", "ATE_2", "ATE_inf"),
label = "Struc_2") +
# MLE for n = infinity
declare_estimator(handler = tidy_estimator(structural_estimator),
Y = "Y_inf_obs",
y = function(Z, d) Z*d/(1+d) + (1-Z)*(1-d/(1+d)),
estimand = c("k","d","q","ATE_2", "ATE_inf"),
label = "Struc_inf")
# Diagnose the design -----------------------------------------------------
diagnosis_1 <- diagnose_design(design, sims = sims)
# Change design to remove random assignment -------------------------------
design_2 <- redesign(design, random = FALSE)
# Diagnose design without random assignment -------------------------------
diagnosis_2 <- diagnose_design(design_2, sims = sims)
# Change design to assume q = 0 -------------------------------------------
# Turn random assignment on again
design_3 <- redesign(design_2, random = TRUE)
# New estimator that assumes q = 0
structural_estimator_2 <- function(data,
Y = "Y_2_obs",
# function to calculate
# predicted price for rational
# customers given z_i and delta,
# defaults to n=2:
y = function(Z, d)
Z * d + (1 - Z) * (1 - d)) {
# Define negative log likelihood as a function of kappa and delta
LL <- function(k, d) {
m <- with(data, y(Z, d))
R <- dbeta(data[Y][[1]], k * m, k * (1 - m))
- sum(log(R))
}
# Estimation
M <- mle2(LL, start = list(k = 5, d = 0.5))
# Format output from estimation
out <- data.frame(coef(summary(M)), outcome = Y)
names(out) <- c("estimate", "std.error", "statistic",
"p.value", "outcome")
# Use estimates of delta to predict average treatment effects (ATEs)
# Predicted ATE for n=2
out[3, "estimate"] <- 2 * out["d", "estimate"] - 1
# Predicted ATE for n=infinity
out[4, "estimate"] <- 2 * (out["d", "estimate"] / (1 + out["d", "estimate"])) - 1
out
}
# Update MLE for n = 2
design_3 <- replace_step(
design_3,
11,
declare_estimator(
handler =
tidy_estimator(structural_estimator_2),
estimand =
c("k", "d", "ATE_2", "ATE_inf"),
label = "Struc_2_norm"
)
)
# Update MLE for n = infinity
design_3 <- replace_step(
design_3,
12,
declare_estimator(
handler =
tidy_estimator(structural_estimator_2),
Y = "Y_inf_obs",
y = function(Z, d)
Z * d / (1 + d) + (1 - Z) *
(1 - d / (1 + d)),
estimand = c("k", "d", "ATE_2", "ATE_inf"),
label = "Struc_inf_norm"
)
)
# Diagnose design with wrong model ----------------------------------------
diagnosis_3 <- diagnose_design(design_3 , sims = sims)
# Change data generating process so that behavioral customers pay 1/2 -----
design_4 <- redesign(design_3, e = 1/2)
# Diagnose design with wrong but observationally equivalent model ---------
diagnosis_4 <- diagnose_design(design_4 , sims = sims)
# Helper function to clean diagnoses --------------------------------------
cleanDiagn <- function(diagnosis) {
out <- reshape_diagnosis(diagnosis)
out[out=="NA"] <- ""
out <- out[!out$`Design Label` == "",]
names(out)[which(names(out) == "Estimator Label")] <- "Estimator"
names(out)[which(names(out) == "Mean Estimand")] <- "Estimand"
return(out)
}
# Result Tables -----------------------------------------------------------
# Tables in the main text refer to maximum likelihood estimators as "MLE_n"
# The code here reproduces the tables in the appendix which refer to the same estimators as "Struc_n"
# Specifically:
# MLE_2 = Struc_2
# MLE_{\infty} = Struc_inf
# MLE'_2 = Struc_2_norm
# MLE'_{\infty} = Struc_inf_norm
# Diagnosis of design with correct model (tables 2, 3 and 4 in text and table 8 in appendix)
result1 <- cleanDiagn(diagnosis_1)
kable(result1[, - c(1, 4, 5, 8, 9, 11:13)], row.names = F, caption = "Detailed diagnosis of design with correct model")
# Diagnosis of design with correct model but no randomization (table 5 in text and table 9 in appendix)
result2 <- cleanDiagn(diagnosis_2)
kable(result2[, - c(1, 2, 5, 6, 9, 10, 12:14)], row.names = F,
caption = "Detailed diagnosis of design with correct model but no randomization")
# Diagnosis of design with incorrect model (6 in text and table 10 in appendix)
diagnosis_3$diagnosands_df
result3 <- cleanDiagn(diagnosis_3)
kable(result3[, - c(1:2, 5, 6, 9, 10, 12:14)], row.names = F,
caption = "Detailed diagnosis of design with incorrect model (assume $q=0$)")
# Diagnosis of design with incorrect yet observationally equivalent model (table 7 in text and table 11 in appendix)
result4 <- cleanDiagn(diagnosis_4)
kable(result4[, - c(1:2, 5, 6, 9, 10, 12:14)], row.names = F,
caption = "Detailed Diagnosis of design with incorrect yet observationally equivalent model")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment