Welfare Calculation
Using our demand estimates, we evaluate the impact that the monitoring program has on consumer welfare and on firm profits in the status quo. To do this, we simulate a counterfactual scenario in which monitoring is not available, and compare it to the status quo baseline with monitoring.
- For each driver \(i\), simulate random coefficients (private risk) \(L \in \mathbb{N}^{+}\) times
- For each draw \(l \in \{ 1, ..., L \}\), calculate ex-ante utility directly and the corresponding certainty equivalent. First-period choice probabilities are also calculated, which gives us the monitoring share. Expected cost of the first semester can be calculated directly. But we also need to form an expectation of the second-period cost (and prices) in order to calculate total surplus (and profit):
- Simulate \(K \in \mathbb{N}^{+}\) draws of first-period claim occurrence and monitoring score based on private risk . Each draw pins down the renewal price change that driver \(i\) would face in the second period. All other prices remain constant. For each first-period choice \(d\), we can then calculate the second period choice probability and the corresponding expected cost
- To do this, we need to generate choice sequence. First we calculate choice probability in \(t = 0\) for each draw of \(\lambda\)
| 1 |
1 |
0 |
3 |
100 |
0.04 |
… |
0.02 |
… |
0.25 |
… |
0.20 |
| 1 |
2 |
0 |
1 |
110 |
0.05 |
… |
0.03 |
… |
0.05 |
… |
0.10 |
| 1 |
3 |
0 |
2 |
90 |
0.05 |
… |
0.03 |
… |
0.35 |
… |
0.40 |
| 1 |
4 |
0 |
2 |
95 |
0.05 |
… |
0.03 |
… |
0.35 |
… |
0.30 |
| 2 |
1 |
0 |
1 |
130 |
0.02 |
… |
0.01 |
… |
0.15 |
… |
0.20 |
# Number of Draws
L <- 100
K <- 50
## df for t = 0, 1
df_t0 <- df %>%
dplyr::filter(t == 0)
df_t1 <- df %>%
dplyr::filter(t == 1)
## Simulate Private Risk Lambda
sim_lambda <- gen_lambda(df, cost_par, L) %>%
dplyr::select(starts_with("lambda"))
sim_lambda_expand <- df %>%
dplyr::bind_cols(sim_lambda)
for(k in 1:(K - 1)){
sim_lambda_expand <- dplyr::bind_cols(sim_lambda_expand, sim_lambda)
}
column_order_lambda <- sim_lambda_expand %>%
dplyr::select(starts_with("lambda")) %>%
names() %>%
sort()
sim_lambda_expand <- sim_lambda_expand %>%
dplyr::select("i", "d", "t", "m", "f",
"x_1", "x_2", "x_3", "x_4",
"price", "prior_firm", "choice",
"claims", "score", "R_s", "R_C",
all_of(column_order_lambda))
colnames(sim_lambda_expand) <- c("i", "d", "t", "m", "f",
"x_1", "x_2", "x_3", "x_4",
"price", "prior_firm", "choice",
"claims", "score", "R_s", "R_C",
paste(rep("lambda", L * K),
rep(1:L, K, each = K),
rep(1:K, L),
sep = "_"))
## Lambda with L columns
sim_lambda_t0 <- df %>%
dplyr::bind_cols(sim_lambda) %>%
dplyr::filter(t == 0)
sim_lambda_t1 <- df %>%
dplyr::bind_cols(sim_lambda) %>%
dplyr::filter(t == 1)
## Lambda with L * K columns
sim_lambda_expand_t0 <- sim_lambda_expand %>%
dplyr::filter(t == 0)
sim_lambda_expand_t1 <- sim_lambda_expand %>%
dplyr::filter(t == 1)
## Calculate Choice Probability in t = 0
sim_choice_probability_t0 <- function(df, theta, cost_par, lambda, Col){
# Parameters
eta_f <- theta[1:4]
xi <- theta[5:6]
alpha_m0 <- theta[7:9]
alpha_m1 <- theta[10:13]
beta <- theta[14]
risk_aversion <- theta[15]
# Cost Parameters
theta_score <- cost_par[8:11]
sigma_score <- cost_par[12]
l0 <- cost_par[13]
a_l <- cost_par[14]
# Price
price <- df$price * matrix(rep(1, dim(df)[1] * Col), nrow = dim(df)[1])
colnames(price) <- c(paste("price", 1:Col, sep = "_"))
price <- data.frame(price) %>%
tibble::as_tibble()
# Generate Demand Friction
demand_friction <- lambda %>%
dplyr::mutate_at(vars(starts_with("lambda")),
funs(ifelse(f == prior_firm,
0,
eta_f[1] + eta_f[2]*x_1 + eta_f[3]*x_2 + eta_f[4]*x_3)
+ ifelse(d == 1 & t == 0,
xi[1] + xi[2]*log(.),
0))) %>%
dplyr::rename_with(\(x) str_replace(x, "lambda", "demand_friction"),
starts_with("lambda")) %>%
dplyr::select(matches("demand_friction"))
# Generate E_score
E_score <- lambda %>%
dplyr::mutate_at(vars(starts_with("lambda")),
funs(ifelse(t == 0 & m == 1,
exp(theta_score[1] + theta_score[2]*log(.)
+ theta_score[3]*x_1 + theta_score[4]*x_2
+ sigma_score^2 / 2),
0))) %>%
dplyr::rename_with(\(x) str_replace(x, "lambda", "E_score"),
starts_with("lambda"))
# Generate E_oop
E_oop <- lambda %>%
dplyr::mutate_at(vars(starts_with("lambda")),
funs(.*exp(-.)*(l0^(a_l)/(a_l - 1))*policy_limit^(-a_l + 1))) %>%
dplyr::rename_with(\(x) str_replace(x, "lambda", "E_oop"),
starts_with("lambda")) %>%
dplyr::select(matches("E_oop"))
# Generate E_oop_sq
E_oop_sq <- lambda %>%
dplyr::mutate_at(vars(starts_with("lambda")),
funs(.*exp(-.)*((2 * l0^(a_l)) / ((a_l - 1) * (a_l - 1))) * policy_limit^(-a_l + 2))) %>%
dplyr::rename_with(\(x) str_replace(x, "lambda", "E_oop_sq"),
starts_with("lambda")) %>%
dplyr::select(matches("E_oop_sq"))
# Generate E_R_s
E_R_s <- E_score %>%
dplyr::mutate_at(vars(starts_with("E_score")),
funs(ifelse(t == 0 & m == 1,
(alpha_m1[1] + alpha_m1[2]*x_1 + alpha_m1[3]*x_2 + alpha_m1[4]*.) / beta,
(alpha_m0[1] + alpha_m0[2]*x_1 + alpha_m0[3]*x_2) / beta))) %>%
dplyr::rename_with(\(x) str_replace(x, "E_score", "E_R_s"),
starts_with("E_score")) %>%
dplyr::select(matches("E_R_s"))
# Generate E_R_s_sq
E_R_s_sq <- E_score %>%
dplyr::mutate_at(vars(starts_with("E_score")),
funs(ifelse(t == 0 & m == 1,
(alpha_m1[1] + alpha_m1[2]*x_1 + alpha_m1[3]*x_2 + alpha_m1[4]*.) / beta^2,
(alpha_m0[1] + alpha_m0[2]*x_1 + alpha_m0[3]*x_2) / beta^2))) %>%
dplyr::rename_with(\(x) str_replace(x, "E_score", "E_R_s_sq"),
starts_with("E_score")) %>%
dplyr::select(matches("E_R_s_sq"))
# Generate E_R_C
E_R_C <- lambda %>%
dplyr::mutate_at(vars(starts_with("lambda")),
funs(0.98 * exp(-.) + 1.2 * (1 - exp(-.)))) %>%
dplyr::rename_with(\(x) str_replace(x, "lambda", "E_R_C"),
starts_with("lambda")) %>%
dplyr::select(matches("E_R_C"))
# Generate E_R_C_sq
E_R_C_sq <- lambda %>%
dplyr::mutate_at(vars(starts_with("lambda")),
funs(0.98^2 * exp(-.) + 1.2^2 * (1 - exp(-.)))) %>%
dplyr::rename_with(\(x) str_replace(x, "lambda", "E_R_C_sq"),
starts_with("lambda")) %>%
dplyr::select(matches("E_R_C_sq"))
# Calculate E_h
E_h <- - price - demand_friction - E_oop - E_R_s * E_R_C * price
E_h <- E_h %>%
dplyr::rename_with(\(x) str_replace(x, "price", "E_h"),
starts_with("price"))
# Calculate E_h_sq
E_h_sq <- price^2 + demand_friction^2 + E_oop_sq + E_R_s_sq*E_R_C_sq*price^2
+ 2*price*demand_friction + 2*price*E_oop + 2*price*(E_R_s*E_R_C*price)
+ 2*demand_friction*E_oop + 2*demand_friction*(E_R_s*E_R_C*price)
+ 2*E_oop*(E_R_s*E_R_C*price)
E_h_sq <- E_h_sq %>%
dplyr::rename_with(\(x) str_replace(x, "price", "E_h_sq"),
starts_with("price"))
# Calculate v
v <- (E_h - (risk_aversion / 2) * E_h_sq) %>%
dplyr::rename_with(\(x) str_replace(x, "E_h", "v"),
starts_with("E_h"))
## Calculate Choice Probability
output <- df %>%
dplyr::bind_cols(v) %>%
dplyr::group_by(i, t) %>%
dplyr::mutate_at(vars(starts_with("v")),
funs(exp(.) / sum(exp(.)))) %>%
dplyr::ungroup() %>%
dplyr::rename_with(\(x) str_replace(x, "v", "choice_probability"),
starts_with("v"))
return(output)
}
sim_choice_probability_t0(df_t0, theta, cost_par, sim_lambda_t0, L)
- Next, we simulate \(K \in \mathbb{N}^{+}\) draws of first-period claim occurrence and monitoring score based on private risk for each \(\lambda\) (thus we simulate \(L \times K\) claims and monitoiring score in total)
| 1 |
1 |
0 |
0.05 |
… |
0.03 |
0 |
… |
1 |
-3.0 |
… |
-1.0 |
| 1 |
2 |
0 |
0.05 |
… |
0.03 |
0 |
… |
0 |
0.0 |
… |
0.0 |
| 1 |
3 |
0 |
0.05 |
… |
0.03 |
1 |
… |
0 |
0.0 |
… |
0.0 |
| 1 |
4 |
0 |
0.05 |
… |
0.03 |
0 |
… |
0 |
0.0 |
… |
0.0 |
| 2 |
1 |
0 |
0.02 |
… |
0.01 |
2 |
… |
0 |
-2.5 |
… |
-1.2 |
- Based on each claim and monitoring score, renewal price multiplier \(R_{0}\) is calculated
| 1 |
1 |
0 |
0.98 |
… |
1.20 |
0.99 |
… |
0.95 |
0.97 |
… |
1.15 |
| 1 |
2 |
0 |
0.98 |
… |
0.98 |
1.02 |
… |
0.98 |
1.00 |
… |
0.96 |
| 1 |
3 |
0 |
1.20 |
… |
0.98 |
1.02 |
… |
0.98 |
1.22 |
… |
0.96 |
| 1 |
4 |
0 |
0.98 |
… |
0.98 |
1.02 |
… |
0.98 |
1.00 |
… |
0.96 |
| 2 |
1 |
0 |
1.20 |
… |
0.98 |
0.87 |
… |
0.90 |
1.07 |
… |
0.88 |
sim_renewal_price_multiplier <- function(df, theta, cost_par, lambda, Col){
## Generate Claims
lambda_mat <- lambda %>%
dplyr::select(starts_with("lambda")) %>%
as.matrix()
claims <- matrix(rpois(dim(df)[1]*Col, lambda = lambda_mat), nrow = dim(df)[1])
colnames(claims) <- paste("claims", 1:Col, sep = "_")
claims <- data.frame(claims) %>%
tibble::as_tibble()
## Generate Monitoring Score
theta_score <- cost_par[8:11]
sigma_score <- cost_par[12]
mu_score <- lambda %>%
dplyr::mutate_at(vars(starts_with("lambda")),
funs(theta_score[1] + theta_score[2]*log(.)
+ theta_score[3]*x_1 + theta_score[4]*x_2)) %>%
dplyr::rename_with(\(x) str_replace(x, "lambda", "mu_score"),
starts_with("lambda")) %>%
dplyr::select(starts_with("mu_score")) %>%
as.matrix()
score <- matrix(rlnorm(dim(df)[1]*Col, mu_score, sigma_score), nrow = dim(df)[1])
colnames(score) <- paste("score", 1:Col, sep = "_")
score <- data.frame(score) %>%
tibble::as_tibble()
## Calculate Claim Surcharge R_C
R_C <- claims %>%
dplyr::mutate_at(vars(starts_with("claims")),
funs(ifelse(. == 0, 0.98, 1.2))) %>%
dplyr::rename_with(\(x) str_replace(x, "claims", "R_C"),
starts_with("claims"))
## Calculate Base R_s
### Note that, for simplicity, assume that R0 is deterministic conditional on s.
### In reality, the spread of baseline R0 without monitoring may have subtle
### nonlinear effects on consumer choice, which we assume away
alpha_m0 <- theta[7:9]
alpha_m1 <- theta[10:13]
beta <- theta[14]
R_s <- df %>%
dplyr::select(m, x_1, x_2) %>%
dplyr::bind_cols(score) %>%
dplyr::mutate_at(vars(starts_with("score")),
funs(ifelse(m == 0,
(alpha_m0[1] + alpha_m0[2]*x_1 + alpha_m0[3]*x_2) / beta,
(alpha_m1[1] + alpha_m1[2]*x_1 + alpha_m1[3]*x_2 + alpha_m1[4]*.) / beta))) %>%
dplyr::rename_with(\(x) str_replace(x, "score", "R_s"),
starts_with("score")) %>%
dplyr::select(-c(m, x_1, x_2))
# ## Calculate Base Renewal Price
# price <- df$price * matrix(rep(1, dim(df)[1] * L * K), nrow = dim(df)[1])
# colnames(price) <- c(paste("price", 1:(L * K), sep = "_"))
#
# price <- data.frame(price) %>%
# tibble::as_tibble()
R_0 <- (R_C * R_s) %>%
dplyr::rename_with(\(x) str_replace(x, "R_C", "R_0"),
starts_with("R_C"))
output <- df %>%
dplyr::bind_cols(R_0) %>%
dplyr::select(c("i", "d", matches("R_0")))
return(output)
}
sim_renewal_price_multiplier(df_t0, theta, cost_par, sim_lambda_expand_t0, L*K)
- Renewal price itself is expanded to the following table:
| 1 |
1 |
1 |
1 |
100 |
0.97 |
… |
1.15 |
97 |
… |
115 |
| 1 |
2 |
1 |
1 |
90 |
0.97 |
… |
1.15 |
87 |
… |
105 |
| 1 |
3 |
1 |
1 |
95 |
0.97 |
… |
1.15 |
92 |
… |
105 |
| 1 |
1 |
1 |
2 |
110 |
1.00 |
… |
0.96 |
110 |
… |
105 |
| 1 |
2 |
1 |
2 |
90 |
1.00 |
… |
0.96 |
90 |
… |
85 |
| 1 |
3 |
1 |
2 |
95 |
1.00 |
… |
0.96 |
95 |
… |
90 |
| 1 |
1 |
1 |
3 |
110 |
1.22 |
… |
0.96 |
130 |
… |
105 |
| 1 |
2 |
1 |
3 |
90 |
1.22 |
… |
0.96 |
110 |
… |
85 |
| 1 |
3 |
1 |
3 |
95 |
1.22 |
… |
0.96 |
115 |
… |
90 |
| 1 |
1 |
1 |
4 |
110 |
1.00 |
… |
0.96 |
110 |
… |
105 |
| 1 |
2 |
1 |
4 |
90 |
1.00 |
… |
0.96 |
90 |
… |
85 |
| 1 |
3 |
1 |
4 |
95 |
1.00 |
… |
0.96 |
95 |
… |
90 |
| 2 |
1 |
1 |
1 |
130 |
1.07 |
… |
0.88 |
137 |
… |
105 |
prior_price <- df_t0 %>%
dplyr::select(i, d, price) %>%
dplyr::rename(prior_price = price)
renewal_price_multiplier <- sim_renewal_price_multiplier(df_t0, theta, cost_par, sim_lambda_expand_t0, L*K)
renewal_price <- expand.grid(i = 1:N, d = 1:3, t = 1, prior_choice = 1:4) %>%
tibble::as_tibble() %>%
dplyr::arrange(i, prior_choice) %>%
dplyr::mutate(key = ifelse(d == 1 & prior_choice == 1, 1, d + 1)) %>%
dplyr::left_join(prior_price, by = c("i", "key" = "d")) %>%
dplyr::left_join(renewal_price_multiplier, by = c("i", "prior_choice" = "d"))
- Finally we calculate choice probability in \(t = 1\) based on each renewal price and \(\lambda\) in \(t = 1\). Note that \(\lambda\) in \(t = 1\) is different from that of \(t = 0\), and hence we need to draw again
| 1 |
1 |
1 |
1 |
97 |
… |
115 |
0.06 |
… |
0.04 |
… |
0.45 |
… |
0.40 |
| 1 |
2 |
1 |
1 |
87 |
… |
105 |
0.06 |
… |
0.04 |
… |
0.25 |
… |
0.30 |
| 1 |
3 |
1 |
1 |
92 |
… |
105 |
0.06 |
… |
0.04 |
… |
0.30 |
… |
0.30 |
| 1 |
1 |
1 |
2 |
110 |
… |
105 |
0.06 |
… |
0.04 |
… |
0.25 |
… |
0.20 |
| 1 |
2 |
1 |
2 |
90 |
… |
85 |
0.06 |
… |
0.04 |
… |
0.25 |
… |
0.40 |
| 1 |
3 |
1 |
2 |
95 |
… |
90 |
0.06 |
… |
0.04 |
… |
0.50 |
… |
0.40 |
| 1 |
1 |
1 |
3 |
130 |
… |
105 |
0.06 |
… |
0.04 |
… |
0.20 |
… |
0.10 |
| 1 |
2 |
1 |
3 |
110 |
… |
85 |
0.06 |
… |
0.04 |
… |
0.35 |
… |
0.50 |
| 1 |
3 |
1 |
3 |
115 |
… |
90 |
0.06 |
… |
0.04 |
… |
0.45 |
… |
0.40 |
| 1 |
1 |
1 |
4 |
110 |
… |
105 |
0.06 |
… |
0.04 |
… |
0.55 |
… |
0.30 |
| 1 |
2 |
1 |
4 |
90 |
… |
85 |
0.06 |
… |
0.04 |
… |
0.15 |
… |
0.30 |
| 1 |
3 |
1 |
4 |
95 |
… |
90 |
0.06 |
… |
0.04 |
… |
0.35 |
… |
0.40 |
| 2 |
1 |
1 |
1 |
137 |
… |
105 |
0.03 |
… |
0.07 |
… |
0.25 |
… |
0.20 |
- Thus, while choice probability in \(t = 0\) has \((4N \times L)\) matrix form, choice probability in \(t = 1\) is \((12N \times LK)\) matrix, due to \(K\) times draw of claims and monitoring score as well as path dependence of choice
## Calculate Certainty Equivalent for Each Draw
### theta_estを入れるのが正しいが挙動確認のためthetaを入れている
cal_certainty_equivalent <- function(df, theta, cost_par, R, lambda){
## Risk Aversion
risk_aversion <- theta[15]
## Calculate CE
### We apply `cal_v` to obtain v
output <- cal_v(df, theta, cost_par, R, lambda) %>%
dplyr::mutate_at(vars(starts_with("v")),
funs(1 / risk_aversion - sqrt((1 / risk_aversion^2) - 2*. / risk_aversion))) %>%
dplyr::rename_with(\(x) str_replace(x, "v", "certatinty_equivalent"),
starts_with("v"))
return(output)
}
#tmp_price <- sim_renewal_price(df_t0, theta, cost_par, R, sim_lambda_expand)
# Calculate Choice Probability in t = 1
cal_v_modified <- function(df, theta, cost_par, R, lambda_t0, lambda_t1){
# Parameters
eta_f <- theta[1:4]
xi <- theta[5:6]
alpha_m0 <- theta[7:9]
alpha_m1 <- theta[10:13]
beta <- theta[14]
risk_aversion <- theta[15]
# Cost Parameters
theta_score <- cost_par[8:11]
sigma_score <- cost_par[12]
l0 <- cost_par[13]
a_l <- cost_par[14]
# Price
price <- sim_renewal_price(df, theta, cost_par, R, lambda_t0)
# Generate Demand Friction
demand_friction <- df %>%
dplyr::left_join(lambda_t1, by = c("i", "f")) %>%
dplyr::mutate_at(vars(starts_with("lambda")),
funs(~ifelse(f == prior_firm,
0,
eta_f[1] + eta_f[2]*x_1 + eta_f[3]*x_2 + eta_f[4]*x_3)
+ ifelse(d == 1 & t == 0,
xi[1] + xi[2]*log(.),
0))) %>%
dplyr::rename_with(\(x) str_replace(x, "lambda", "demand_friction"),
starts_with("lambda")) %>%
dplyr::select(matches("demand_friction"))
# Generate E_score
E_score <- df %>%
dplyr::left_join(lambda_t1, by = c("i", "f")) %>%
dplyr::mutate_at(vars(starts_with("lambda")),
funs(ifelse(t == 0 & m == 1,
exp(theta_score[1] + theta_score[2]*log(.)
+ theta_score[3]*x_1 + theta_score[4]*x_2
+ sigma_score^2 / 2),
0))) %>%
dplyr::rename_with(\(x) str_replace(x, "lambda", "E_score"),
starts_with("lambda"))
# Generate E_oop
E_oop <- df %>%
dplyr::left_join(lambda_t1, by = c("i", "f")) %>%
dplyr::mutate_at(vars(starts_with("lambda")),
funs(.*exp(-.)*(l0^(a_l)/(a_l - 1))*policy_limit^(-a_l + 1))) %>%
dplyr::rename_with(\(x) str_replace(x, "lambda", "E_oop"),
starts_with("lambda")) %>%
dplyr::select(matches("E_oop"))
# Generate E_oop_sq
E_oop_sq <- df %>%
dplyr::left_join(lambda_t1, by = c("i", "f")) %>%
dplyr::mutate_at(vars(starts_with("lambda")),
funs(.*exp(-.)*((2 * l0^(a_l)) / ((a_l - 1) * (a_l - 1))) * policy_limit^(-a_l + 2))) %>%
dplyr::rename_with(\(x) str_replace(x, "lambda", "E_oop_sq"),
starts_with("lambda")) %>%
dplyr::select(matches("E_oop_sq"))
# Generate E_R_s
E_R_s <- E_score %>%
dplyr::mutate_at(vars(starts_with("E_score")),
funs(ifelse(t == 0 & m == 1,
(alpha_m1[1] + alpha_m1[2]*x_1 + alpha_m1[3]*x_2 + alpha_m1[4]*.) / beta,
(alpha_m0[1] + alpha_m0[2]*x_1 + alpha_m0[3]*x_2) / beta))) %>%
dplyr::rename_with(\(x) str_replace(x, "E_score", "E_R_s"),
starts_with("E_score")) %>%
dplyr::select(matches("E_R_s"))
# Generate E_R_s_sq
E_R_s_sq <- E_score %>%
dplyr::mutate_at(vars(starts_with("E_score")),
funs(ifelse(t == 0 & m == 1,
(alpha_m1[1] + alpha_m1[2]*x_1 + alpha_m1[3]*x_2 + alpha_m1[4]*.) / beta^2,
(alpha_m0[1] + alpha_m0[2]*x_1 + alpha_m0[3]*x_2) / beta^2))) %>%
dplyr::rename_with(\(x) str_replace(x, "E_score", "E_R_s_sq"),
starts_with("E_score")) %>%
dplyr::select(matches("E_R_s_sq"))
# Generate E_R_C
E_R_C <- df %>%
dplyr::left_join(lambda_t1, by = c("i", "f")) %>%
dplyr::mutate_at(vars(starts_with("lambda")),
funs(0.98 * exp(-.) + 1.2 * (1 - exp(-.)))) %>%
dplyr::rename_with(\(x) str_replace(x, "lambda", "E_R_C"),
starts_with("lambda")) %>%
dplyr::select(matches("E_R_C"))
# Generate E_R_C_sq
E_R_C_sq <- df %>%
dplyr::left_join(lambda_t1, by = c("i", "f")) %>%
dplyr::mutate_at(vars(starts_with("lambda")),
funs(0.98^2 * exp(-.) + 1.2^2 * (1 - exp(-.)))) %>%
dplyr::rename_with(\(x) str_replace(x, "lambda", "E_R_C_sq"),
starts_with("lambda")) %>%
dplyr::select(matches("E_R_C_sq"))
# Calculate E_h
E_h <- - price - demand_friction - E_oop - E_R_s * E_R_C * price
E_h <- E_h %>%
dplyr::rename_with(\(x) str_replace(x, "price", "E_h"),
starts_with("price"))
# Calculate E_h_sq
E_h_sq <- price^2 + demand_friction^2 + E_oop_sq + E_R_s_sq*E_R_C_sq*price^2
+ 2*price*demand_friction + 2*price*E_oop + 2*price*(E_R_s*E_R_C*price)
+ 2*demand_friction*E_oop + 2*demand_friction*(E_R_s*E_R_C*price)
+ 2*E_oop*(E_R_s*E_R_C*price)
E_h_sq <- E_h_sq %>%
dplyr::rename_with(\(x) str_replace(x, "price", "E_h_sq"),
starts_with("price"))
# Calculate v
v <- (E_h - (risk_aversion / 2) * E_h_sq) %>%
dplyr::rename_with(\(x) str_replace(x, "E_h", "v"),
starts_with("E_h"))
output <- df %>%
dplyr::bind_cols(v)
return(output)
}
cal_choice_probability_modified <- function(df, theta, cost_par, R, lambda_t0, lambda_t1){
## Calculate v
v <- cal_v_modified(df, theta, cost_par, R, lambda_t0, lambda_t1)
## Calculate Choice Probability
output <- v %>%
dplyr::group_by(i, t) %>%
dplyr::mutate_at(vars(starts_with("v")),
funs(exp(.) / sum(exp(.)))) %>%
dplyr::ungroup() %>%
dplyr::rename_with(\(x) str_replace(x, "v", "choice_probability"),
starts_with("v"))
}
tmp_choice_prob <- cal_choice_probability_modified(df_t0, theta, cost_par, R, sim_lambda_expand_t0, sim_lambda_expand_t1)