################################################################### # # Script for Replication of: Are Consumers Willing to Pay to Let # Cars Drive for Them? Analyzing Response to Autonomous Vehicles # Authors: Ricardo Daziano, Mauricio Sarrias, Benjamin Leard ################################################################### # # ======================================= # # I. Load Data and Packages---- # ======================================= # library("mlogit") library("devtools") #install_bitbucket("mauricio1986/gmnl") # Install the lastest version of gmnl library("gmnl") setwd("~/Documents/Work/Work with Daziano/Automation/Rfiles") # Set directory data.autom <- read.table("longinc.csv", sep = ",", header = T) # Load data autom.mlogit <- mlogit.data(data.autom, choice = "chosen", shape = "long", alt.var = "alternative", id = "id") # mlogit format # ======================================== # # II. Creating Dataset and Variables ---- # ======================================== # # Create ASCs autom.mlogit$asc2 <- as.numeric(autom.mlogit$alternative == "e") autom.mlogit$asc3 <- as.numeric(autom.mlogit$alternative == "h") autom.mlogit$asc4 <- as.numeric(autom.mlogit$alternative == "p") # Scaling variables: This helps for the convergence of complex models autom.mlogit$pvfc1r <- autom.mlogit$pvfc1r / 10 autom.mlogit$nage <- autom.mlogit$age / 10 autom.mlogit$nyrsdrivng <- autom.mlogit$yrsdrivng / 10 # Create new database with -1 * price autom.mlogitO <- mlogit.data(autom.mlogit, choice = "chosen", shape = "long", alt.var = "alternative", id = "id", opposite = c("pricen")) autom.mlogitO$ocharget <- -1 * autom.mlogitO$charget autom.mlogitO$opvfc1r <- -1 * autom.mlogitO$pvfc1r # Create interactions autom.mlogitO$autom1_gc <- autom.mlogitO$autom1 * autom.mlogitO$googlecar autom.mlogitO$autom1_ma <- autom.mlogitO$autom1 * autom.mlogitO$male autom.mlogitO$autom1_in <- autom.mlogitO$autom1 * autom.mlogitO$lninc autom.mlogitO$autom1_yd <- autom.mlogitO$autom1 * autom.mlogitO$nyrsdrivng autom.mlogitO$autom1_ws <- autom.mlogitO$autom1 * autom.mlogitO$west autom.mlogitO$autom1_md <- autom.mlogitO$autom1 * autom.mlogitO$midwest autom.mlogitO$autom1_ne <- autom.mlogitO$autom1 * autom.mlogitO$northeast autom.mlogitO$autom2_gc <- autom.mlogitO$autom2 * autom.mlogitO$googlecar autom.mlogitO$autom2_ma <- autom.mlogitO$autom2 * autom.mlogitO$male autom.mlogitO$autom2_in <- autom.mlogitO$autom2 * autom.mlogitO$lninc autom.mlogitO$autom2_yd <- autom.mlogitO$autom2 * autom.mlogitO$nyrsdrivng autom.mlogitO$autom2_ws <- autom.mlogitO$autom2 * autom.mlogitO$west autom.mlogitO$autom2_md <- autom.mlogitO$autom2 * autom.mlogitO$midwest autom.mlogitO$autom2_ne <- autom.mlogitO$autom2 * autom.mlogitO$northeast # ======================================== # # III. Table 4: Conditional Logits ---- # ======================================== # # Column 1: Montly cost mnl_1 <- gmnl(chosen ~ price + asc2 + asc3 + asc4 + ocost + lnrange + charget + autom1 + autom2 | 0, data = autom.mlogitO) summary(mnl_1) # Column 2: Cost per mile mnl_2 <- gmnl(chosen ~ price + asc2 + asc3 + asc4 + cost + lnrange + charget + autom1 + autom2 | 0, data = autom.mlogitO) summary(mnl_2) # Column 3: PVFC mnl_ps <- gmnl(chosen ~ price + asc2 + asc3 + asc4 + pvfc1r + lnrange + charget + autom1 + autom2 | 0 , data = autom.mlogitO) summary(mnl_ps) AIC(mnl_1) AIC(mnl_2) AIC(mnl_ps) BIC(mnl_1) BIC(mnl_2) BIC(mnl_ps) # Implied WTP measures wtp.gmnl(mnl_1, wrt = "price") wtp.gmnl(mnl_2, wrt = "price") wtp.gmnl(mnl_ps, wrt = "price") # ======================================== # # IV. Table 5: Models with parametric heterogeneity ---- # ======================================== # # MIXL-N mixl_1 <- gmnl(chosen ~ price + asc2 + asc3 + asc4 + pvfc1r + lnrange + charget + autom1 + autom2 | 0 , data = autom.mlogitO, model = "mixl", R = 500, ranp = c(lnrange = "n", charget = "n", autom1 = "n", autom2 = "n"), panel = TRUE, iterlim = 500, print.level = 2, method = "bhhh", print.init = TRUE) summary(mixl_1) AIC(mixl_1) BIC(mixl_1) # MIXL-N-OH mixl_hier <- gmnl(chosen ~ price + asc2 + asc3 + asc4 + pvfc1r + lnrange + charget + autom1 + autom2 | 0 | 0 | googlecar + male + lninc + nyrsdrivng + west + midwest + northeast, data = autom.mlogitO, model = "mixl", R = 500, ranp = c(lnrange = "n", charget = "n", autom1 = "n", autom2 = "n"), mvar = list(autom1 = c("googlecar", "male", "lninc", "nyrsdrivng", "west", "midwest", "northeast"), autom2 = c("googlecar", "male", "lninc", "nyrsdrivng", "west", "midwest", "northeast")), panel = TRUE, iterlim = 500, print.level = 3, method = "bhhh", print.init = TRUE) summary(mixl_hier) AIC(mixl_hier) BIC(mixl_hier) # MIXL-LN-I start_ln2 <- c(coef(mixl_1)[c("price", "asc2", "asc3", "asc4", "pvfc1r")], coef(mixl_1)["lnrange"], coef(mixl_1)["charget"], log(coef(mixl_1)[c("autom1", "autom1")]), rep(0.1, 4)) mixlLn_2 <- gmnl(chosen ~ price + asc2 + asc3 + asc4 + pvfc1r + lnrange + charget + autom1 + autom2 | 0 , data = autom.mlogitO, model = "mixl", R = 100, ranp = c(lnrange = "n", charget = "n", autom1 = "ln", autom2 = "ln"), panel = TRUE, iterlim = 500, print.level = 3, method = "bhhh", start = start_ln2, print.init = TRUE) summary(mixlLn_2) AIC(mixlLn_2) BIC(mixlLn_2) # Compute mean and standard deviation of log-normal distribution library("msm") mean_autom1 <- exp(coef(mixlLn_2)["autom1"] + (coef(mixlLn_2)["sd.autom1"]^2)/2) mean_autom1 # Automation 1's mean se.mean_autom1 <- deltamethod( ~exp(x8 + (x12 ^ 2) / 2), mean = coef(mixlLn_2), cov = vcov(mixlLn_2), ses = TRUE) 2 * (1 - pnorm(abs(mean_autom1 / se.mean_autom1))) # Automation 2's mean mean_autom2 <- exp(coef(mixlLn_2)["autom2"] + (coef(mixlLn_2)["sd.autom2"]^2)/2) mean_autom2 se.mean_autom2 <- deltamethod( ~exp(x9 + (x13 ^ 2)/2), mean = coef(mixlLn_2), cov = vcov(mixlLn_2), ses = TRUE) se.mean_autom2 2 * (1 - pnorm(abs(mean_autom2 / se.mean_autom2))) # Automation 1's standard deviation sd_autom1 <- exp(coef(mixlLn_2)["autom1"] + (coef(mixlLn_2)["sd.autom1"]^2)/2) * sqrt(exp(coef(mixlLn_2)["sd.autom1"] ^ 2 ) - 1) sd_autom1 se.sd_autom <- deltamethod( ~exp(x8 + (x12 ^ 2)/2) * sqrt(exp(x12 ^2) - 1), mean = coef(mixlLn_2), cov = vcov(mixlLn_2), ses = TRUE) se.sd_autom 2 * (1 - pnorm(abs(sd_autom1 / se.sd_autom))) # Automation 2's standard deviation sd_autom2 <- exp(coef(mixlLn_2)["autom2"] + (coef(mixlLn_2)["sd.autom2"]^2)/2) * sqrt(exp(coef(mixlLn_2)["sd.autom2"] ^ 2 ) - 1) sd_autom2 se.sd_autom2 <- deltamethod( ~exp(x9 + (x13 ^ 2)/2) * sqrt(exp(x13 ^2) - 1), mean = coef(mixlLn_2), cov = vcov(mixlLn_2), ses = TRUE) se.sd_autom2 2 * (1 - pnorm(abs(sd_autom2 / se.sd_autom2))) # MIXL-LN-II start_ln3 <- c(coef(mixl_1)[c("price", "asc2", "asc3", "asc4", "pvfc1r")], log(coef(mixl_1)["lnrange"]), coef(mixl_1)["charget"], log(coef(mixl_1)[c("autom1", "autom1")]), rep(0.1, 4)) mixlLn_3 <- gmnl(chosen ~ price + asc2 + asc3 + asc4 + pvfc1r + lnrange + charget + autom1 + autom2 | 0 , data = autom.mlogitO, model = "mixl", R = 100, ranp = c(lnrange = "ln", charget = "n", autom1 = "ln", autom2 = "ln"), panel = TRUE, iterlim = 500, print.level = 3, method = "bhhh", start = start_ln3, print.init = TRUE) summary(mixlLn_3) AIC(mixlLn_3) BIC(mixlLn_3) # Range's mean mean_range <- exp(coef(mixlLn_3)["lnrange"] + (coef(mixlLn_3)["sd.lnrange"] ^ 2) / 2) mean_range se.mean_range <- deltamethod( ~exp(x6 + (x10 ^ 2) / 2), mean = coef(mixlLn_3), cov = vcov(mixlLn_3), ses = TRUE) se.mean_range 2 * (1 - pnorm(abs(mean_range / se.mean_range))) # Range's standard deviation sd_range <- exp(coef(mixlLn_3)["lnrange"] + (coef(mixlLn_3)["sd.lnrange"] ^ 2) / 2) * sqrt(exp(coef(mixlLn_3)["sd.lnrange"] ^ 2 ) - 1) sd_range se.sd_range <- deltamethod( ~exp(x6 + (x10 ^ 2)/2) * sqrt(exp(x10 ^2) - 1), mean = coef(mixlLn_3), cov = vcov(mixlLn_3), ses = TRUE) se.sd_range 2 * (1 - pnorm(abs(sd_range / se.sd_range))) # Automation 1's mean mean_autom1 <- exp(coef(mixlLn_3)["autom1"] + (coef(mixlLn_3)["sd.autom1"] ^ 2) / 2) mean_autom1 se.mean_autom1 <- deltamethod( ~exp(x8 + (x12 ^ 2) / 2), mean = coef(mixlLn_3), cov = vcov(mixlLn_3), ses = TRUE) se.mean_autom1 2 * (1 - pnorm(abs(mean_autom1 / se.mean_autom1))) # Automation 2's mean mean_autom2 <- exp(coef(mixlLn_3)["autom2"] + (coef(mixlLn_3)["sd.autom2"] ^ 2) / 2) mean_autom2 se.mean_autom2 <- deltamethod( ~exp(x9 + (x13 ^ 2) / 2), mean = coef(mixlLn_3), cov = vcov(mixlLn_3), ses = TRUE) se.mean_autom2 2 * (1 - pnorm(abs(mean_autom2 / se.mean_autom2))) # Automation 1's standard deviation sd_autom1 <- exp(coef(mixlLn_3)["autom1"] + (coef(mixlLn_3)["sd.autom1"] ^ 2) / 2) * sqrt(exp(coef(mixlLn_3)["sd.autom1"] ^ 2 ) - 1) sd_autom1 se.sd_autom <- deltamethod( ~exp(x8 + (x12 ^ 2)/2) * sqrt(exp(x12 ^ 2) - 1), mean = coef(mixlLn_3), cov = vcov(mixlLn_3), ses = TRUE) se.sd_autom 2 * (1 - pnorm(abs(sd_autom1 / se.sd_autom))) # Automation 2's standard deviation sd_autom2 <- exp(coef(mixlLn_3)["autom2"] + (coef(mixlLn_3)["sd.autom2"] ^ 2) / 2) * sqrt(exp(coef(mixlLn_3)["sd.autom2"] ^ 2 ) - 1) sd_autom2 se.sd_autom2 <- deltamethod( ~exp(x9 + (x13 ^ 2) / 2) * sqrt(exp(x13 ^ 2) - 1), mean = coef(mixlLn_3), cov = vcov(mixlLn_3), ses = TRUE) se.sd_autom2 2 * (1 - pnorm(abs(sd_autom2 / se.sd_autom2))) # ======================================== # # V. Table 6: Models with no parametric heterogeneity---- # ======================================== # # Estimate LC model for initial values lc3_C <- gmnl(chosen ~ price + asc2 + asc3 + asc4 + pvfc1r + lnrange + autom1 + autom2 | 0 | 0 | 0 | 1 + days80miles + ownercar + nvehadd + googlecar + nage + male + married + childcnt + compcollege + hschorless + singfamh + apartment + ownhouse + nyrsdrivng + accident + preferdriving + fulltime + parttime + hmaker + white + conserv + lib + west + midwest + northeast + urban + lninc, data = autom.mlogitO, model = "lc", panel = TRUE, iterlim = 500, method = "bhhh", Q = 3) summary(lc3_C) # Starting values start_mm_3C <- c(coef(lc3_C)[1:(8 * 3)], rep(0, 6), rep(0.1, 2), rep(0, 6), rep(0.1, 2), rep(0, 6), rep(0.1, 2), coef(lc3_C)[((8 * 3) + 1):length(coef(lc3_C))]) # MM-MNL Model: Warning, it takes around 5 hours! mm_3_C <- gmnl(chosen ~ price + asc2 + asc3 + asc4 + pvfc1r + lnrange + autom1 + autom2 | 0 | 0 | 0 | 1 + days80miles + ownercar + nvehadd + googlecar + nage + male + married + childcnt + compcollege + hschorless + singfamh + apartment + ownhouse + nyrsdrivng + accident + preferdriving + fulltime + parttime + hmaker + white + conserv + lib + west + midwest + northeast + urban + lninc, data = autom.mlogitO, model = "mm", R = 500, ranp = c(price = "n", asc2 = "n", asc3 = "n", asc4 = "n", pvfc1r = "n", lnrange = "n", autom1 = "n", autom2 = "n"), panel = TRUE, iterlim = 500, print.level = 3, method = "bhhh", print.init = TRUE, fixed = c(rep(FALSE, 8 * 3), rep(TRUE, 6), rep(FALSE, 2), rep(TRUE, 6), rep(FALSE, 2), rep(TRUE, 6), rep(FALSE, 2), rep(FALSE, length(coef(lc3_C)[((8 * 3) + 1):length(coef(lc3_C))]))), start = start_mm_3C, Q = 3) summary(mm_3_C) AIC(mm_3_C) BIC(mm_3_C) # Probabilities for each class colMeans(mm_3_C$Qir$wnq) # ======================================== # # VI. Computing WTP ---- # ======================================== # # Plot Unconditional Distribution (WTP) set.seed(123) S <- 10000 # Automation 1 autom1_mixl <- -(coef(mixl_1)["autom1"] + abs(coef(mixl_1)["sd.autom1"]) * rnorm(S)) / coef(mixl_1)["price"] autom1_oh <- -(coef(mixl_hier)["autom1"] + mean(autom.mlogitO$autom1_gc) * coef(mixl_hier)["autom1.googlecar"] + mean(autom.mlogitO$autom1_ma) * coef(mixl_hier)["autom1.male"] + mean(autom.mlogitO$autom1_in) * coef(mixl_hier)["autom1.lninc"] + mean(autom.mlogitO$autom1_yd) * coef(mixl_hier)["autom1.nyrsdrivng"] + mean(autom.mlogitO$autom1_ws) * coef(mixl_hier)["autom1.west"] + mean(autom.mlogitO$autom1_md) * coef(mixl_hier)["autom1.midwest"] + mean(autom.mlogitO$autom1_ne) * coef(mixl_hier)["autom1.northeast"] + abs(coef(mixl_hier)["sd.autom1"]) * rnorm(S)) / coef(mixl_hier)["price"] autom1_log <- -(exp(coef(mixlLn_2)["autom1"] + abs(coef(mixlLn_2)["sd.autom1"]) * rnorm(S))) / coef(mixlLn_2)["price"] autom1_log2 <- -(exp(coef(mixlLn_3)["autom1"] + abs(coef(mixlLn_3)["sd.autom1"]) * rnorm(S))) / coef(mixlLn_3)["price"] autom1_mm_class_1 <- (-(coef(mm_3_C)["class.1.autom1"] + abs(coef(mm_3_C)["class.1.sd.autom1"]) * rnorm(S)) / coef(mm_3_C)["class.1.price"]) autom1_mm_class_2 <- (-(coef(mm_3_C)["class.2.autom1"] + abs(coef(mm_3_C)["class.2.sd.autom1"]) * rnorm(S)) / coef(mm_3_C)["class.2.price"]) autom1_mm_class_3 <- (-(coef(mm_3_C)["class.3.autom1"] + abs(coef(mm_3_C)["class.3.sd.autom1"]) * rnorm(S)) / coef(mm_3_C)["class.3.price"]) # Automation 2 autom2_mixl <- -(coef(mixl_1)["autom2"] + abs(coef(mixl_1)["sd.autom2"]) * rnorm(S)) / coef(mixl_1)["price"] autom2_oh <- -(coef(mixl_hier)["autom2"] + mean(autom.mlogitO$autom2_gc) * coef(mixl_hier)["autom2.googlecar"] + mean(autom.mlogitO$autom2_ma) * coef(mixl_hier)["autom2.male"] + mean(autom.mlogitO$autom2_in) * coef(mixl_hier)["autom2.lninc"] + mean(autom.mlogitO$autom2_yd) * coef(mixl_hier)["autom2.nyrsdrivng"] + mean(autom.mlogitO$autom2_ws) * coef(mixl_hier)["autom2.west"] + mean(autom.mlogitO$autom2_md) * coef(mixl_hier)["autom2.midwest"] + mean(autom.mlogitO$autom2_ne) * coef(mixl_hier)["autom2.northeast"] + abs(coef(mixl_hier)["sd.autom2"]) * rnorm(S)) / coef(mixl_hier)["price"] autom2_log <- -(exp(coef(mixlLn_2)["autom2"] + abs(coef(mixlLn_2)["sd.autom2"]) * rnorm(S))) / coef(mixlLn_2)["price"] autom2_log2 <- -(exp(coef(mixlLn_3)["autom2"] + abs(coef(mixlLn_3)["sd.autom2"]) * rnorm(S))) / coef(mixlLn_3)["price"] autom2_mm_class_1 <- (-(coef(mm_3_C)["class.1.autom2"] + abs(coef(mm_3_C)["class.1.sd.autom2"]) * rnorm(S)) / coef(mm_3_C)["class.1.price"]) autom2_mm_class_2 <- (-(coef(mm_3_C)["class.2.autom2"] + abs(coef(mm_3_C)["class.2.sd.autom2"]) * rnorm(S)) / coef(mm_3_C)["class.2.price"]) autom2_mm_class_3 <- (-(coef(mm_3_C)["class.3.autom2"] + abs(coef(mm_3_C)["class.3.sd.autom2"]) * rnorm(S)) / coef(mm_3_C)["class.3.price"]) # Figure 2 par(mfrow = c(1, 2), mar = c(1, 1, 3, 1), oma = c(4, 4, 0, 0)) par(cex = 0.7) plot(density(autom1_mixl), xlab = "WTP / 1000", lty = 1, lwd = 1, ylim = c(0, 0.17), main = "Automation 1") lines(density(autom1_oh), lty = 2, lwd = 1) lines(density(autom1_log), lty = 3, lwd = 1) lines(density(autom1_log2), lty = 4, lwd = 1) legend("topleft", c("MIXL-N", "MXIL-N-OH", "MIXL-LN-I", "MIXL-LN-II"), lty = c(1, 2, 3, 4, 5) ) plot(density(autom2_mixl), xlab = "WTP / 1000", lty = 1, lwd = 1, ylim = c(0, 0.17), main = "Automation 2") lines(density(autom2_oh), lty = 2, lwd = 1) lines(density(autom2_log), lty = 3, lwd = 1) lines(density(autom2_log2), lty = 4, lwd = 1) legend("topleft", c("MIXL-N", "MIXL-N-OH", "MIXL-LN-I", "MIXL-LN-II"), lty = c(1, 2, 3, 4, 5) ) # Figure 3 par(mfrow = c(1, 2), mar = c(1, 1, 3, 1), oma = c(4, 4, 0, 0)) par(cex = 0.7) plot(density(autom1_mm_class_1), xlab = "WTP / 1000", lty = 1, lwd = 1, ylim = c(0, 0.8), xlim = c(-20, 20), main = "Automation 1") lines(density(autom1_mm_class_2), lty = 2, lwd = 1) lines(density(autom1_mm_class_3), lty = 3, lwd = 1) legend("topleft", c("Class 1", "Class 2", "Class 3"), lty = c(1, 2, 3) ) plot(density(autom2_mm_class_1), xlab = "WTP / 1000", lty = 1, lwd = 1, ylim = c(0, 0.07), xlim = c(-40, 40), main = "Automation 2") lines(density(autom2_mm_class_2), lty = 2, lwd = 1) lines(density(autom2_mm_class_3), lty = 3, lwd = 1) legend("topleft", c("Class 1", "Class 2", "Class 3"), lty = c(1, 2, 3) )