課題13：動学的離散選択 — Rust の NFXP と Hotz-Miller の CCP

計量経済学II

重要この課題のねらい

第13回でやった 動学的離散選択（Rust のバスエンジン交換モデル）を、自分の手で最初から最後まで通す。価値関数反復で DP を解き、パネルを生成し、NFXP（Rust）と Hotz-Miller の CCP（renewal action の finite dependence）の両方でパラメータを推定し、そして交換費補助の反実仮想を、正しい動学モデルと myopic（\(\beta=0\)）モデルとで比べる。今日の技術を全部使う、重めの課題だ。焦らず、まず DP を解いて S 字が出ることを確かめ、それから推定・反実仮想へ進もう。

• Part 1：\(K=15\)・\(\beta=0.95\) 固定の縮小版 Rust モデルを自作し、DP を解き、200台×100ヶ月のパネルを生成する（40点）。
• Part 2：NFXP と CCP で真値を回収し、交換費補助の反実仮想を動学と myopic で比べる（60点）。

配点は合計100点。手を動かせば2〜4時間で終わる分量。穴埋めコードの骨格を惜しみなく与えるので、まず埋めることから始めよ。

この課題は講義ノート（lecture13.qmd）と同じ構造だが、設定（状態数・パラメータ値・状態遷移）は変えてある。ノートのコードをコピペするだけでは通らないので、自分で理解して埋めること。とくに \(K=15\)、\(\theta_1=0.25\)、\(RC=5.0\)、状態遷移の確率が講義と違う点に注意。

ノート記法（第13回の再掲）

\(\beta\) は割引因子（\(=0.95\) に固定）。効用のパラメータは \(\theta\)。状態 \(x\)（走行距離、\(0 \ldots K-1\)）、行動 \(a \in \{0=\text{維持}, 1=\text{交換}\}\)。flow utility は \(u(x,0) = -\theta_1 x\)、\(u(x,1) = -RC\)。T1EV ショックのもとで \(V(x) = \log(e^{v(x,0)} + e^{v(x,1)}) + \gamma\)、選択確率は logit。

---

Part 1：縮小版 Rust モデルを自作する（40点）

舞台設定

社用車フリートの整備担当を考える。各車両のエンジンについて、毎月「維持（そのまま使う）」か「交換（載せ替える）」かを決める。走行距離を \(K=15\) 状態のグリッド \(x \in \{0, \ldots, 14\}\) で表す。

真のパラメータ（あなたはこれを後で推定で回収する）：

記号	意味	値
\(\beta\)	割引因子（固定、推定しない）	0.95
\(\theta_1\)	維持費の傾き（1状態あたり）	0.25
\(RC\)	エンジン交換費	5.0

状態遷移（講義とは確率が違う）：維持すると距離は \(+0/+1/+2\) 状態進み、その確率は \((0.25, 0.50, 0.25)\)。交換すると \(x=0\) にリセットしてから同じ分布で進む。

問1（15点）：状態遷移行列と価値関数反復

状態遷移行列 \(P_0\)（維持）・\(P_1\)（交換）を作り、価値関数反復で DP を解く関数 solve_vfi を完成させよ。\(V(x)\) の更新は Emax = logsum：\(V^{\text{new}}(x) = \log(e^{v(x,0)} + e^{v(x,1)}) + \gamma\)。数値安定化のため「最大値を引いてから exp」すること。

# ---- 基本設定 ----
K     <- 15
beta  <- 0.95
theta1_true <- 0.25
RC_true     <- 5.0
p_trans <- c(0.25, 0.50, 0.25)          # 維持で +0/+1/+2 進む確率（講義と違う！）
EULER   <- 0.5772156649015329
xs <- 0:(K - 1)

# 状態遷移行列（1-index に注意。状態値は x-1）
build_transition <- function(K, p) {
  nj <- length(p)
  P0 <- matrix(0, K, K); P1 <- matrix(0, K, K)
  for (x in 1:K) {
    xv <- x - 1
    for (j in 1:nj) {
      jj <- j - 1
      P0[x, min(xv + jj, K - 1) + 1] <- P0[x, min(xv + jj, K - 1) + 1] + p[j]
      P1[x, min(0  + jj, K - 1) + 1] <- P1[x, min(0  + jj, K - 1) + 1] + p[j]
    }
  }
  list(P0 = P0, P1 = P1)
}
tr <- build_transition(K, p_trans)
P0 <- tr$P0; P1 <- tr$P1

flow_utils <- function(theta1, RC) list(u0 = -theta1 * xs, u1 = rep(-RC, K))

# ---- この関数を完成させる ----
solve_vfi <- function(theta1, RC, tol = 1e-10, maxit = 2000, bta = beta) {
  uu <- flow_utils(theta1, RC)
  V  <- rep(0, K)
  errs <- numeric(0)
  for (it in 1:maxit) {
    v0 <- uu$u0 + bta * as.vector(P0 %*% V)     # choice-specific value（維持）
    v1 <- ______                                # 【穴埋め】choice-specific value（交換）
    m  <- pmax(v0, v1)
    Vnew <- ______                              # 【穴埋め】logsum + Euler 定数（数値安定化込み）
    err  <- max(abs(Vnew - V))
    errs <- c(errs, err)
    V <- Vnew
    if (err < tol) break
  }
  v0 <- uu$u0 + bta * as.vector(P0 %*% V)
  v1 <- uu$u1 + bta * as.vector(P1 %*% V)
  list(V = V, v0 = v0, v1 = v1, errs = errs, niter = length(errs))
}

ccp_replace <- function(v0, v1) {
  m <- pmax(v0, v1)
  ______                                        # 【穴埋め】P(交換|x) の logit
}

# ---- 動作確認（穴埋め後に eval: true にして実行）----
sol <- solve_vfi(theta1_true, RC_true)
Pr  <- ccp_replace(sol$v0, sol$v1)
cat("収束反復数:", sol$niter, "\n")
cat("P(交換|x) が単調増加:", all(diff(Pr) >= -1e-12), "\n")
cat("P(交換|x):", round(Pr, 3), "\n")

# 収束（対数目盛）と S 字を1枚ずつ
library(patchwork)
p_conv <- ggplot(data.frame(iter = seq_along(sol$errs), err = sol$errs), aes(iter, err)) +
  geom_line(color = "#1f77b4") + scale_y_log10() +
  labs(title = "価値関数反復の収束", x = "反復", y = "更新幅（対数）")
p_scurve <- ggplot(data.frame(x = xs, P = Pr), aes(x, P)) +
  geom_line(color = "#d62728", linewidth = 1) + geom_point() +
  labs(title = "交換確率 P(交換|x)", x = "走行距離 x", y = "P(交換|x)")
p_conv + p_scurve

ヒントヒント：状態遷移行列の作り方と logsum

• 遷移行列：各行 \(x\) が「次期状態の確率分布」になる（行和=1）。維持は今の \(x\) から距離が進む、交換は \(x=0\) から進む。min(..., K-1) で上限に張り付かせるのを忘れずに。\(P_1\) は全行が同じになるはず（交換=リセットなので現在状態に依存しない）。これが renewal action の性質で、Part 2 の CCP で効いてくる。
• logsum の数値安定化：log(exp(v0) + exp(v1)) は、\(v\) が大きいと exp がオーバーフローする。各状態で m <- pmax(v0, v1) を引いてから m + log(exp(v0 - m) + exp(v1 - m)) とする（第2回の overflow 対策と同じ）。最後に + EULER（Euler 定数 \(\gamma\)）を足す。

問2（15点）：パネルを生成し、replacement cycle を可視化する

解いた交換方針に従って動く 200台×100ヶ月のパネルを生成せよ。各期、各車両は自分の状態 \(x\) の交換確率 \(P(\text{交換}\mid x)\) に従って交換を決め、状態が遷移する。数台ぶんの走行距離の軌跡（のこぎり歯 = replacement cycle）を描くこと。

simulate_panel <- function(n_bus, n_month, theta1, RC, burn = 50) {
  sol <- solve_vfi(theta1, RC)
  Pr  <- ccp_replace(sol$v0, sol$v1)
  total_T <- n_month + burn
  X <- matrix(0L, n_bus, total_T); A <- matrix(0L, n_bus, total_T)
  x <- sample(0:(K - 1), n_bus, replace = TRUE)
  for (t in 1:total_T) {
    u <- runif(n_bus)
    a <- ______                                 # 【穴埋め】状態 x の交換確率で交換判定（0/1）
    X[, t] <- x
    A[, t] <- a
    jump <- sample(0:2, n_bus, replace = TRUE, prob = p_trans)
    x <- ______                                 # 【穴埋め】交換なら 0+jump にリセット、維持なら x+jump（上限 K-1）
  }
  list(X = X[, (burn + 1):total_T], A = A[, (burn + 1):total_T])
}

n_bus <- 200; n_month <- 100
panel <- simulate_panel(n_bus, n_month, theta1_true, RC_true)
X <- panel$X; A <- panel$A
cat("全体の交換率:", round(mean(A), 4), "  平均走行距離:", round(mean(X), 3), "\n")

# 数台の軌跡（赤点 = 交換した月）
show_bus <- c(1, 2, 3, 4)
df_traj <- do.call(rbind, lapply(show_bus, function(b)
  data.frame(month = 1:n_month, mileage = X[b, ], bus = paste0("車", b), replace = A[b, ])))
ggplot(df_traj, aes(month, mileage)) +
  geom_line(color = "grey50") +
  geom_point(data = subset(df_traj, replace == 1), color = "#d62728", size = 1.6) +
  facet_wrap(~ bus, ncol = 2) +
  labs(title = "走行距離の軌跡と交換タイミング", x = "月", y = "走行距離 x")

ヒントヒント：交換判定と状態リセット

• 交換判定：u < Pr[x + 1] で TRUE/FALSE が出るので as.integer(...) で 0/1 に。Pr は 1-index なので状態 \(x\) の確率は Pr[x + 1]。
• 状態更新：ifelse(a == 1L, pmin(0 + jump, K - 1), pmin(x + jump, K - 1))。交換したら 0 から、維持なら今の \(x\) から距離が進む。
• burn-in：初期状態をランダムに置いたので、最初の数十期は定常分布に落ち着くまでの過渡期。burn = 50 期を捨てて、定常に近い部分だけ残す。

問3（10点）：状態ごとの交換頻度を確認する

状態ごとの経験交換頻度 \(\hat P(\text{交換}\mid x)\) を計算し、モデルの交換確率 \(P(\text{交換}\mid x)\) と重ねて描け。観測数の少ない（薄い）状態がどこか、\(x=0\) の交換頻度がどうなっているかを確認し、2〜3文でコメントせよ（この観察が Part 2 の CCP 推定で効いてくる）。

freq_by_state  <- sapply(0:(K - 1), function(xv) { m <- (X == xv); if (sum(m) > 0) mean(A[m]) else NA })
count_by_state <- sapply(0:(K - 1), function(xv) sum(X == xv))

df_freq <- data.frame(x = xs, emp = freq_by_state, model = Pr)
ggplot(df_freq, aes(x)) +
  geom_line(aes(y = model), color = "#1f77b4", linewidth = 1) +
  geom_point(aes(y = emp), color = "#d62728", size = 1.8) +
  labs(title = "経験頻度（赤）vs モデル（青）", x = "走行距離 x", y = "P(交換|x)")
cat("状態ごとの観測数:", count_by_state, "\n")

ノートこの課題での記述統計の意味

高走行距離の状態は観測が薄く、経験頻度がガタつく（交換でリセットされるので、そもそも高い \(x\) に長居しない）。逆に \(x=0\)（交換直後）は観測は多いが交換はまず起きず、経験頻度がほぼ0になる。この「\(\hat P = 0\)（\(\log\) が \(-\infty\)）」と「薄い状態のノイズ」が、CCP 推定で対処すべき2つの問題だ。

---

Part 2：推定と反実仮想（60点）

問4（15点）：NFXP（Rust の方法）で推定する

NFXP の負の対数尤度を完成させ、optim で \((\theta_1, RC)\) を推定せよ。内側で価値関数反復を解き切るのがキモ。\(\beta = 0.95\) に固定すること（推定しない）。真値回収表と Hessian からの標準誤差を出し、実行時間を測れ。

# ---- NFXP の負の対数尤度（内側で DP を解き切る）----
neg_loglik_nfxp <- function(theta, X, A) {
  theta1 <- theta[1]; RC <- theta[2]
  if (theta1 <= 0 || RC <= 0) return(1e10)
  sol <- solve_vfi(theta1, RC, tol = 1e-8, maxit = 1000)   # ← 内側の DP
  m   <- pmax(sol$v0, sol$v1)
  lae <- m + log(exp(sol$v0 - m) + exp(sol$v1 - m))        # log(exp v0 + exp v1)
  logp1 <- ______                                          # 【穴埋め】log P(交換|x)
  logp0 <- ______                                          # 【穴埋め】log P(維持|x)
  xr <- as.vector(X) + 1; ar <- as.vector(A)
  ll <- sum(ifelse(ar == 1L, logp1[xr], logp0[xr]))
  -ll
}

t0 <- Sys.time()
fit_nfxp <- optim(par = c(0.15, 3.0), fn = neg_loglik_nfxp, X = X, A = A,
                  method = "Nelder-Mead", control = list(reltol = 1e-9, maxit = 2000),
                  hessian = TRUE)
t1 <- Sys.time()
nfxp_time <- as.numeric(difftime(t1, t0, units = "secs"))

se_nfxp <- sqrt(diag(solve(fit_nfxp$hessian)))
recovery <- data.frame(
  parameter = c("theta1", "RC"),
  truth     = c(theta1_true, RC_true),
  estimate  = round(fit_nfxp$par, 4),
  se        = round(se_nfxp, 4),
  ci_lower  = round(fit_nfxp$par - 1.96 * se_nfxp, 4),
  ci_upper  = round(fit_nfxp$par + 1.96 * se_nfxp, 4)
)
kable(recovery, caption = "NFXP：真値回収と標準誤差")
cat("NFXP 実行時間:", round(nfxp_time, 2), "秒\n")

警告識別上の注意：\(\beta\) を選択データだけで推定する難しさ

この課題では \(\beta\) を推定対象に加えない。選択データだけでは割引因子とflow utilityを一般に分離できず、追加の除外制約・外部情報・正規化が要る。ここの尤度へ無制約に\(\beta\)を足すと目的関数が平坦・不安定になり、探索が\(\beta>1\)や負の値へ飛ぶこともある。\(\beta=0.95\)を固定して\((\theta_1,RC)\)だけを推定し、必要なら別の固定値で感度分析する。

問5（15点）：Hotz-Miller の CCP で推定する

CCP アプローチを実装せよ。renewal action（交換=リセット）の finite dependence で、価値関数を解かずに線形回帰で \(\theta\) を得る。恒等式は

\[ \underbrace{\log \hat P(1\mid x) - \log \hat P(0\mid x) + \beta \sum_{x'}\big[P_1 - P_0\big](x'\mid x)\,\log \hat P(1\mid x')}_{=:\ y(x)} \;=\; \theta_1\, x - RC. \]

\(y(x)\) を経験 CCP から計算し、\(x\) に回帰する（傾き=\(\theta_1\)、切片=\(-RC\)）。\(\hat P=0\)対策にはJeffreys/add-half平滑化 \((n_1+0.5)/(n+1)\) を使う（Laplace/add-oneは \((n_1+1)/(n+2)\)）。薄い状態対策として観測数 \(\geq30\) の状態をWLSに使う。

diffP <- P1 - P0

# 状態ごとの交換回数 n1 と観測回数 n
n_x  <- sapply(0:(K - 1), function(xv) sum(X == xv))
n1_x <- sapply(0:(K - 1), function(xv) sum(A[X == xv]))

# Jeffreys/add-half 平滑化した CCP
Psm    <- (n1_x + 0.5) / (n_x + 1.0)
logP1h <- log(Psm); logP0h <- log(1 - Psm)

# 将来項の logP1（未観測状態は最も近い観測状態で代用。境界も安全に処理）
logP1h_full <- logP1h
observed_idx <- which(n_x > 0)
if (length(observed_idx) == 0) stop("CCPを計算できる観測状態がない")
for (idx in which(n_x == 0)) {
  nearest <- observed_idx[which.min(abs(observed_idx - idx))]
  logP1h_full[idx] <- logP1h[nearest]
}

# y(x) を計算 → 観測数>=30 の状態で重み付き回帰
y_hat <- ______                                 # 【穴埋め】logP1h - logP0h + beta * (P1-P0) 行 x logP1h_full
good  <- n_x >= 30
fit_ccp <- lm(y_hat ~ xs, weights = n_x, subset = good)

theta1_ccp <- coef(fit_ccp)[2]
RC_ccp     <- -coef(fit_ccp)[1]                 # 切片が -RC
cat("CCP：theta1 =", round(theta1_ccp, 4), " RC =", round(RC_ccp, 4),
    "  （使った状態数:", sum(good), "）\n")

# NFXP と CCP を並べる
compare_tbl <- data.frame(
  method = c("NFXP", "CCP"),
  theta1 = round(c(fit_nfxp$par[1], theta1_ccp), 4),
  RC     = round(c(fit_nfxp$par[2], RC_ccp), 4)
)
kable(compare_tbl, caption = "NFXP vs CCP", row.names = FALSE)

ヒントヒント：finite dependence がなぜ効くか（講義の復習）

交換すると状態が \(x=0\) にリセットされる（\(P_1\) が全行同じ）。だから「交換後の将来価値」は現在状態によらない定数になり、\(v(x,1)-v(x,0)\) の将来項の差をとると打ち消える。残るのは遷移確率（既知）と CCP（データから推定可能）だけ。だから価値関数反復なしで \(\theta\) が線形回帰で出る。まず真の CCP（ccp_replace(sol$v0, sol$v1)）で \(y(x) = \theta_1 x - RC\) が厳密に成り立つことを確認すると、恒等式の実装が正しいか検算できる（講義ノート参照）。

警告よくあるバグ：\(\gamma\) 定数の扱いと将来項の符号

• \(\gamma\)（Euler 定数）：CCP の恒等式には \(\gamma\) が現れない（\(v(x,1)-v(x,0)\) の差で \(\gamma\) が消えるから）。y_hat に \(\gamma\) を足さないこと。\(\gamma\) が要るのは価値関数反復（logsum）の方だけ。
• 将来項の符号：\(y(x)\) の第3項は \(+\beta (P_1 - P_0) \log \hat P(1\mid x')\) で、プラス（導出で \(-\log P(1\mid x')\) が2回符号を変える）。ここを間違えると \(\theta_1\)・\(RC\) が真値から大きく外れる。真の CCP での検算（\(y = \theta_1 x - RC\)）が一致するかで検知できる。

問6（15点）：反実仮想を動学と myopic で比べる

交換費 \(RC\) を20%補助したときの効果を、(a) 正しい動学モデル（\(\beta=0.95\)）と (b) myopic モデル（\(\beta=0\) で推定したパラメータを \(\beta=0\) の DP で使う）とで比べよ。長期交換率・平均走行距離は、方針が誘導するマルコフ連鎖の定常分布から計算する。

# 方針が誘導する状態遷移 M と定常分布
policy_transition <- function(Pr) (1 - Pr) * P0 + Pr * P1
ergodic_dist <- function(M) {
  A <- rbind(t(M) - diag(K), rep(1, K)); b <- c(rep(0, K), 1)
  d <- qr.solve(A, b); d / sum(d)
}
summarize_policy <- function(theta1, RC, label, bta = beta) {
  s  <- solve_vfi(theta1, RC, bta = bta)
  Pr <- ccp_replace(s$v0, s$v1)
  d  <- ergodic_dist(policy_transition(Pr))
  rep_rate <- sum(d * Pr); mean_mi <- sum(d * xs)
  cat(sprintf("  %s：交換率=%.4f, 平均走行距離=%.3f\n", label, rep_rate, mean_mi))
  list(rep_rate = rep_rate, mean_mi = mean_mi, d = d)
}

# (a) 正しい動学の反実仮想（真値 or NFXP 推定値を使ってよい）
base_dyn <- summarize_policy(theta1_true, RC_true,       "動学 補助前")
sub_dyn  <- summarize_policy(theta1_true, RC_true * 0.8, "動学 補助後")

# (b) myopic 推定：大きめの標本で beta=0 の logit MLE
panel_big <- simulate_panel(2000, 300, theta1_true, RC_true)
Xb <- panel_big$X; Ab <- panel_big$A
neg_loglik_myopic <- function(theta, X, A) {
  theta1 <- theta[1]; RC <- theta[2]
  if (theta1 <= 0 || RC <= 0) return(1e10)
  v0 <- -theta1 * xs; v1 <- rep(-RC, K)
  m  <- pmax(v0, v1); lae <- m + log(exp(v0 - m) + exp(v1 - m))
  xr <- as.vector(X) + 1; ar <- as.vector(A)
  -sum(ifelse(ar == 1L, (v1 - lae)[xr], (v0 - lae)[xr]))
}
fit_myopic <- optim(c(0.15, 3.0), neg_loglik_myopic, X = Xb, A = Ab, method = "Nelder-Mead")
cat("myopic 推定 theta1 =", round(fit_myopic$par[1], 4),
    " RC =", round(fit_myopic$par[2], 4), "（真値 0.25, 5.0）\n")

# myopic 分析者は beta=0 のモデルで反実仮想を予測する
base_myo <- summarize_policy(fit_myopic$par[1], fit_myopic$par[2],       "myopic 補助前", bta = 0)
sub_myo  <- summarize_policy(fit_myopic$par[1], fit_myopic$par[2] * 0.8, "myopic 補助後", bta = 0)

cf_tbl <- data.frame(
  analysis   = c("正しい動学", "myopic (誤り)"),
  rep_before = round(c(base_dyn$rep_rate, base_myo$rep_rate), 4),
  rep_after  = round(c(sub_dyn$rep_rate,  sub_myo$rep_rate), 4),
  rep_change_pct = round(c((sub_dyn$rep_rate / base_dyn$rep_rate - 1) * 100,
                           (sub_myo$rep_rate / base_myo$rep_rate - 1) * 100), 1)
)
kable(cf_tbl, caption = "RC 20%補助の効果予測：動学 vs myopic", row.names = FALSE)

警告よくあるバグ：交換後の状態リセットのバグ

反実仮想でも DP を解き直す（solve_vfi を新しい \(RC\) で呼ぶ）ことを忘れずに。\(RC\) を下げても価値関数を再計算しなければ、交換方針が変わらず、反実仮想が意味をなさない。また policy_transition の (1 - Pr) * P0 + Pr * P1 は、R の列方向リサイクルで「行 \(x\) を \(Pr[x]\) 倍」になる（行和=1）。ここは R の癖に助けられている箇所なので、rowSums(policy_transition(Pr)) が全部1になるか確認するとよい。

問7（10点）：設備保全担当への提言を3行で

あなたは社用車フリートの設備保全を統括する立場だとする。問6の結果（動学 vs myopic の反実仮想の食い違い、NFXP・CCP の推定値）に基づいて、経営・政策への含意を3行で書け。とくに「交換費補助という施策の効果を、なぜ動学モデルで評価すべきか（myopic だと何を誤るか）」に触れること。数値を1つは根拠に引くこと。

# ここに、判断材料を1つの表にまとめておくと書きやすい
# （交換率の変化・平均走行距離の変化・myopic とのズレ）

---

提出方法

• この .qmd の穴埋め（______ と #| eval: false）をすべて埋め、#| eval: false を外して最後まで render できる状態にする。
• 各問の出力（表・図）とコメントを含めること。
• render した HTML を提出。

重要配点（合計100点）

• 問1（15）：状態遷移行列と価値関数反復（S 字・幾何収束）
• 問2（15）：パネル生成と replacement cycle の可視化
• 問3（10）：状態ごとの交換頻度の確認とコメント
• 問4（15）：NFXP の実装・真値回収・SE（\(\beta\) 固定）
• 問5（15）：CCP（finite dependence）の実装・NFXP との比較
• 問6（15）：反実仮想を動学 vs myopic で比較
• 問7（10）：設備保全担当への提言（数値に基づく3行）