3 Estimation
- Now we estimate the choice model, using data we generate
- Again, the choice model in the paper is as follows:
\[ \begin{align} u_{idt}(C, s) &= u_{\gamma}(w_{it} + h_{idt}(C, s)) \\ h_{idt}(C, s) &= -p_{idt} - \mathbf{1}_{d, t-1} \cdot \psi_{idt} - e(C, y_{d}) - p_{idt} \cdot R_{idt}(C, s) \\ \text{where } \psi_{idt} &= \mathbf{1}_{d, t-1} \cdot \eta_{0} + \mathbf{1}_{f_{d}, t-1} \cdot \eta_{it} + \mathbf{1}_{m_{d}} \cdot \mathbf{1}_{t=0} \cdot \xi_{it} \end{align} \]
- At each period \(t\), consumer \(i\) chooses \(d\) so as to maximize her expected utility:
\[ \begin{align} d_{it} &= \mathop{\rm argmax}\limits_{d \in D_{it}} \left\{ v_{idt} + \varepsilon_{idt} \right\} \\ \text{where } v_{idt} &= \mathbb{E}_{C, s}[u_{idt}(C, s)] = \mathbb{E}[h_{idt}] -\frac{\gamma}{2}\mathbb{E}[h_{idt}^{2}] \end{align} \]
We estimate the model by simulated maximum likelihood method. In doing so, we adopt a two-step estimation procedure since the cost model is easier to estimate but requires a large amount of data, while the demand model is making use of rich variations in choice but computationally demanding.
Therefore, we estimate risk and monitoring score parameters \((\theta_{\lambda}, \sigma_{\lambda}, \theta_{s}, \sigma_{s})\) first.
Then, we feed the estimates into the demand models as truth.
The overall log likelihood function is expressed as follows:
\[ \begin{align} \mathscr{L}_{i} = \sum_{t \leq T_{i}} \int_{\lambda} \underbrace{\mathscr{L}(R_{it}, s_{i}, C_{it}, d_{it} \mid \lambda, \psi, x_{it}, p_{it}, D_{it}, d_{i, t-1} ; \Theta)}_{(A): \text{obs. stoc outcome}} \cdot \underbrace{g_{\lambda}(\lambda \mid x_{it} ; \theta_{\lambda}, \sigma_{\lambda})}_{(B): \text{latent var.}} d\lambda \end{align} \]
- We can decompose the likelihood component A as follows:
\[ \begin{align} (A) = & \ln \Pr(d_{it} \mid \lambda, x_{it}, p_{it}, D_{it}, d_{i, t-1}; a, \psi_{0}, \psi_{1}, \theta_{\eta}, \theta_{\xi}, \alpha, \theta_{\beta}) \notag \\ & + \ln \Pr(C_{it} \mid \lambda, x_{it}) + \ln g(l_{it} \mid d_{it}, x_{it}; \alpha, \theta_{\beta}) \notag \\ & + \ln g_{s}(s_{i} \mid \lambda, x_{it}; \theta_{s}, \sigma_{s}) + \ln g_{R}(R_{idt} \mid C_{it}, s_{i}, \lambda, x_{it}; \theta_{R}, \theta_{R, m}, \sigma_{R}) \end{align} \]
3.1 Fuctions
- To estimate the choice model, we write the following functions:
gen_lambdais a function that draws unobserved consumer risk \(\lambda\). In this function, we draw \(\lambda\) for \(R\) timesLL_choiceis a function that calculates choice probability that a consumer chooses each plan for each draw of \(\lambda\)LL_renewal_priceis a function that calculates choice probability that a consumer chooses each plan for each draw of \(\lambda\)SLLis a function that calculates the simulated log-likelihood of consumer’s choice. We maximize this fucntion to get the estimates of choice model parameters
- Note that when optimizing an objective function, it is important to keep the realizations of the shocks the same across the evaluations of the objective function. If the realization of the shocks differ across the objective function evaluations, the optimization algorithm will not converge because it cannot distinguish the change in the value of the objective function due to the difference in the parameters and the difference in the realized shocks
- To avoid this problem, we generate \(\lambda\) outside the optimization algorithm
## Generate Lambda
gen_lambda <- function(df, cost_par, R){
# Parameters
theta_lambda <- cost_par[1:6]
sigma_lambda <- cost_par[7]
# Generate mu_lambda
output <- df %>%
dplyr::mutate(mu_lambda = theta_lambda[1]
+ theta_lambda[2]*x_1 + theta_lambda[3]*x_2
+ theta_lambda[4]*x_3 + theta_lambda[5]*x_4
+ theta_lambda[6]*ifelse(m == 1 & t == 0, 1, 0))
# Generate epsilon_lambda (draw R times)
epsilon_lambda <- matrix(rnorm(N * R, 0, sigma_lambda), nrow = N)
colnames(epsilon_lambda) <- c(paste("lambda", 1:R, sep = "_"))
epsilon_lambda <- data.frame(expand_grid(i = 1:N), epsilon_lambda) %>%
tibble::as_tibble()
## lambda
output <- output %>%
dplyr::left_join(epsilon_lambda, by = c("i")) %>%
dplyr::mutate_at(.vars = vars(starts_with("lambda")),
.funs = ~ exp(mu_lambda + .)) %>%
dplyr::select(-mu_lambda)
return(output)
}
# Calculate Choice Probability
cal_v <- function(df, theta, cost_par, R, lambda){
# 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] * R), nrow = dim(df)[1])
colnames(price) <- c(paste("price", 1:R, 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"))
output <- df %>%
dplyr::bind_cols(v)
return(output)
}
cal_choice_probability <- function(df, theta, cost_par, R, lambda){
## Calculate v
v <- cal_v(df, theta, cost_par, R, lambda)
## 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"))
}
# Log-likelihood for Choice Probability
LL_choice <- function(df, theta, cost_par, R, lambda){
## Calculate Choice Probability
choice_probability <- cal_choice_probability(df, theta, cost_par, R, lambda)
## Calculate Log-likelihood
LL_choice <- choice_probability %>%
dplyr::mutate_at(vars(starts_with("choice_probability")),
funs(choice * log(.))) %>%
dplyr::summarise_at(vars(starts_with("choice_probability")),
funs(sum(.)))
LL_choice <- as.matrix(LL_choice) %>%
apply(1, mean)
# Output
output <- LL_choice
return(output)
}
# Log-likelihood for Renewal Price
LL_renewal_price <- function(df, theta, cost_par, lambda){
# Parameters
alpha_m0 <- theta[7:9]
alpha_m1 <- theta[10:13]
beta <- theta[14]
# Generate log-likelihood of R_s
## In the paper, we assume away uncertainty in s
R_s <- df %>%
dplyr::filter(t == 1) %>%
dplyr::mutate(gamma_alpha = ifelse(score == 1,
alpha_m0[1] + alpha_m0[2]*x_1 + alpha_m0[3]*x_2,
alpha_m1[1] + alpha_m1[2]*x_1 + alpha_m1[3]*x_2 + alpha_m1[4]*score),
LL_gamma = log(((beta^gamma_alpha)/gamma(gamma_alpha)) * R_s^(gamma_alpha - 1) * exp(-beta * R_s)))
LL_R_s <- sum(R_s$LL_gamma)
## 作ったが、Choice Modelの推定においては不要
# Generate log-likelihood of R_C
## We first match lambda in t = 0, since claim surcharge is resulted from lambda in t = 0
# lambda_t0 <- lambda %>%
# dplyr::filter(t == 0, f == 1) %>%
# dplyr::select(i, m, matches("lambda"))
#
# R_C <- df %>%
# dplyr::filter(t == 1) %>%
# dplyr::mutate(m = ifelse(score == 1, 0, 1)) %>%
# dplyr::left_join(lambda_t0, by = c("i", "m")) %>%
# dplyr::mutate_at(vars(starts_with("lambda")),
# funs((exp(-.)^(R_C == 0.98)) * (1 - exp(-.)^(R_C == 1.2)))) %>%
# dplyr::rename_with(\(x) str_replace(x, "lambda", "LL_surcharge"),
# starts_with("lambda")) %>%
# dplyr::summarise_at(vars(starts_with("LL_surcharge")),
# funs(sum(.)))
#
# LL_R_C <- as.matrix(R_C) %>%
# apply(1, mean)
output <- LL_R_s
return(output)
}
## Simulated Log-likelihood
SLL <- function(df, theta, cost_par, R, lambda){
## Choice Probability
LL_choice <- LL_choice(df, theta, cost_par, R, lambda)
## Renewal Price
LL_renewal_price <- LL_renewal_price(df, theta, cost_par, lambda)
## Output
output <- (-1) * (LL_choice + LL_renewal_price)
print(output)
return(output)
}3.2 Optimization
# Cost Parameters
cost_par <- 0
cost_par[1:6] <- theta_lambda
cost_par[7] <- sigma_lambda
cost_par[8:11] <- theta_score
cost_par[12] <- sigma_score
cost_par[13] <- l0
cost_par[14] <- a_l
# Choice Parameters
theta <- 0
theta[1:4] <- eta_f
theta[5:6] <- xi
theta[7:9] <- alpha_m0
theta[10:13] <- alpha_m1
theta[14] <- beta
theta[15] <- risk_aversion
# Generate lambda
R <- 50
lambda <- gen_lambda(df, cost_par, R)
## Initial Value
# x0 <- rep(1, length(theta))
# x0[7] <- 20
# x0[10] <- 20
# x0[14] <- 20
# x0[15] <- 0
x0 <- numeric(length(theta))
x0[7] <- 1
x0[10] <- 1
x0[14] <- 1
x0 <- theta
## Optimization
res <- optim(par = x0,
fn = SLL,
control = list(maxit = 10^6),
method = "Nelder-Mead",
df = df,
cost_par = cost_par,
R = R,
lambda = lambda)
theta_est <- res$par