# 個人の価格係数は正。nu はどの計算でも同じ積分点を使う。
rc_probs <- function(p, quality, a, sigma, nu) {
alpha <- exp(a + sigma * nu)
v <- sweep(-outer(alpha, p), 2, quality, "+")
vmax <- pmax(0, apply(v, 1, max))
ev <- exp(v - vmax)
ev / (exp(-vmax) + rowSums(ev))
}
rc_market <- function(p, quality, a, sigma, nu) {
prob <- rc_probs(p, quality, a, sigma, nu)
alpha <- exp(a + sigma * nu)
jac <- crossprod(prob, prob * alpha) / length(nu)
diag(jac) <- diag(jac) - colMeans(prob * alpha)
list(s = colMeans(prob), jac = jac)
}
owner_matrix <- function(groups) outer(groups, groups, "==") * 1
markup_at <- function(p, groups, quality, a, sigma, nu) {
demand <- rc_market(p, quality, a, sigma, nu)
Delta <- -owner_matrix(groups) * t(demand$jac)
as.vector(solve(Delta, demand$s))
}
solve_prices <- function(mc, groups, quality, a, sigma, nu,
p_init = mc + 1, tol = 1e-10, maxit = 2000) {
p <- p_init
for (it in seq_len(maxit)) {
target <- mc + markup_at(p, groups, quality, a, sigma, nu)
err <- max(abs(target - p))
if (err < tol) return(list(p = p, iterations = it, error = err))
p <- 0.5 * p + 0.5 * target
}
stop("価格の反復が収束しなかった")
}
# 複数市場の確率をまとめて計算。各積分点について市場を並べ、列は商品。
rc_all <- function(delta, price, a, sigma, nu, slopes = TRUE) {
nt <- nrow(price); nj <- ncol(price); nr <- length(nu)
mean_alpha <- exp(a + sigma^2 / 2)
deviation <- exp(a + sigma * nu) - mean_alpha
v <- matrix(0, nt * nr, nj)
for (j in seq_len(nj)) {
v[, j] <- rep(delta[, j], nr) -
rep(price[, j], nr) * rep(deviation, each = nt)
}
vmax <- pmax(0, v[, 1])
for (j in 2:nj) vmax <- pmax(vmax, v[, j])
ev <- exp(v - vmax)
prob <- ev / (exp(-vmax) + rowSums(ev))
share <- matrix(0, nt, nj)
own_slope <- matrix(0, nt, nj)
alpha <- rep(exp(a + sigma * nu), each = nt)
for (j in seq_len(nj)) {
share[, j] <- rowMeans(matrix(prob[, j], nt, nr))
if (slopes) own_slope[, j] <- -rowMeans(matrix(
alpha * prob[, j] * (1 - prob[, j]), nt, nr))
}
list(share = share, own_slope = own_slope)
}
invert_shares <- function(S_obs, price, a, sigma, nu,
tol = 1e-9, maxit = 2000) {
delta <- log(S_obs) - log(1 - rowSums(S_obs))
for (it in seq_len(maxit)) {
predicted <- rc_all(delta, price, a, sigma, nu, slopes = FALSE)$share
correction <- log(S_obs) - log(predicted)
delta <- delta + correction
if (max(abs(correction)) < tol) {
return(list(delta = delta, iterations = it,
error = max(abs(correction))))
}
}
stop("シェアの反転が収束しなかった")
}
vec <- function(m) as.vector(t(m))
# その市場の需要を固定したときの1人あたり消費者余剰。
cs_value <- function(p, quality, a, sigma, nu) {
alpha <- exp(a + sigma * nu)
v <- sweep(-outer(alpha, p), 2, quality, "+")
vmax <- pmax(0, apply(v, 1, max))
inclusive <- vmax + log(exp(-vmax) + rowSums(exp(v - vmax)))
mean(inclusive / alpha)
}
market_outcomes <- function(p, mc, quality, a, sigma, nu, market_size = 10000) {
s <- rc_market(p, quality, a, sigma, nu)$s
c(quantity = market_size * sum(s),
profit = market_size * sum((p - mc) * s),
consumer_surplus = market_size * cs_value(p, quality, a, sigma, nu))
}
