动态优化结构估计大全

1
2

回归：工作小时数 = α + β × 当日价格补贴 + ε
发现：β = 0.5（补贴每增加10%，工作时间增加5%）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

┌─────────────────────────────────────────────────────────────────────────────┐
│                        结构估计工作流                                        │
├─────────────────────────────────────────────────────────────────────────────┤
│                                                                              │
│  ① 经济问题 → ② 模型设定 → ③ 数据对接 → ④ 识别分析                        │
│       │            │           │            │                               │
│       ▼            ▼           ▼            ▼                               │
│   "要回答      "状态/行动    "观测什么？  "哪些variation               │
│    什么？"     /收益/转移"   缺什么？"    识别什么？"                      │
│                                                                              │
│       ⑤ 估计方法选择 → ⑥ 数值实现 → ⑦ 推断与检验 → ⑧ 反事实              │
│            │              │             │              │                    │
│            ▼              ▼             ▼              ▼                    │
│       "MLE/GMM/       "优化器/       "标准误/       "政策模拟/           │
│        CCP/SMM？"      收敛？"        拟合？"        福利？"               │
│                                                                              │
└─────────────────────────────────────────────────────────────────────────────┘

1
2
3
4
5
6
7

                识别 = 找到从θ到数据统计量的可逆映射

        θ = (RC, θ_c)  ──映射M──>  可观测的行为变化 P(d|x)
              ↑                           ↑
          结构参数                    实际能估计的

        识别成功 ⟺ M 是单射（不同θ导致不同P）

1
2
3
4
5
6
7

外层优化（找θ使似然最大）
    │
    └──> 内层不动点（给定θ，解EV(x; θ)）
              │
              └──> 计算P(d|x; θ)
                       │
                       └──> 计算对数似然

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

NFXP算法：
─────────────────────────────────────────────────────
输入：数据{(x_it, d_it)}，初始参数θ^(0)，容差ε_in, ε_out

OUTER LOOP（优化θ）:
  while ||∇L(θ)|| > ε_out:

    INNER LOOP（解EV）:
      初始化 EV^(0)(x) = 0 for all x
      while ||EV^(k+1) - EV^(k)|| > ε_in:
        EV^(k+1)(x) = T(EV^(k); θ)  # 贝尔曼迭代
      end

    计算选择概率 P(d|x; θ, EV)
    计算似然 L(θ) 和梯度 ∇L(θ)
    更新 θ ← θ + step × direction

  end

输出：θ̂_MLE
─────────────────────────────────────────────────────

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

def solve_EV(theta, beta, G, tol=1e-10, max_iter=1000):
    """
    求解EV的不动点
    theta: (RC, theta_c1, theta_c2)
    G: 转移概率矩阵 [n_x, n_x, 2] (最后维度是action)
    """
    RC, tc1, tc2 = theta
    n_x = G.shape[0]
    x_grid = np.arange(n_x)

    # 当期效用（不含epsilon）
    u0 = -tc1 * x_grid - tc2 * x_grid**2  # 继续使用
    u1 = -RC - tc1 * 0 - tc2 * 0**2        # 更换（里程归零）
    u1 = np.full(n_x, u1)

    EV = np.zeros(n_x)

    for _ in range(max_iter):
        # 条件价值函数
        v0 = u0 + beta * G[:,:,0] @ EV
        v1 = u1 + beta * G[:,:,1] @ EV

        # Log-sum-exp（数值稳定版）
        max_v = np.maximum(v0, v1)
        EV_new = max_v + np.log(np.exp(v0 - max_v) + np.exp(v1 - max_v))

        if np.max(np.abs(EV_new - EV)) < tol:
            break
        EV = EV_new

    return EV

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

from scipy.optimize import minimize

def neg_log_likelihood(theta, data, beta, G):
    """
    负对数似然（用于最小化）
    data: DataFrame with columns ['x', 'd']
    """
    EV = solve_EV(theta, beta, G)

    RC, tc1, tc2 = theta
    x = data['x'].values
    d = data['d'].values

    # 计算选择概率
    u0 = -tc1 * x - tc2 * x**2 + beta * EV[x]  # 注意：需要插值或索引
    u1 = -RC + beta * EV[0]  # 更换后里程归零

    # 数值稳定的softmax
    max_u = np.maximum(u0, u1)
    prob1 = np.exp(u1 - max_u) / (np.exp(u0 - max_u) + np.exp(u1 - max_u))
    prob1 = np.clip(prob1, 1e-10, 1-1e-10)  # 避免log(0)

    # 对数似然
    ll = d * np.log(prob1) + (1-d) * np.log(1-prob1)
    return -np.sum(ll)

# 多初值优化
best_result = None
for _ in range(20):
    theta0 = np.random.uniform([5, 0.001, 0], [15, 0.01, 0.001])
    result = minimize(neg_log_likelihood, theta0, args=(data, beta, G),
                      method='L-BFGS-B',
                      bounds=[(0.1, 50), (1e-6, 0.1), (0, 0.01)])
    if best_result is None or result.fun < best_result.fun:
        best_result = result

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

CCP两步法：
─────────────────────────────────────────────────────
第一步：估计"可约化对象"（不涉及结构参数）
  (a) 估计CCP：p̂₁(x) = #{d=1 且 x=x} / #{x=x}  # 或用Logit平滑
  (b) 估计转移概率：Ĝ(x'|x,d) = 频率计数 / 非参数估计

第二步：估计结构参数（单次优化）
  构造伪似然：
    L̃(θ) = Σ [d·log p̃₁(x;θ,p̂) + (1-d)·log(1-p̃₁(x;θ,p̂))]

  其中 p̃₁ 是用θ和p̂构造的"伪CCP"

  最大化 L̃(θ) 得到 θ̂

─────────────────────────────────────────────────────

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

def ccp_estimator(data, beta, n_forward=100):
    """
    CCP两步估计器（前向模拟版本）
    """
    # 第一步：估计CCP和转移概率
    ccp = data.groupby('x')['d'].mean().values  # 简单频率估计
    ccp = np.clip(ccp, 0.001, 0.999)  # 避免边界问题

    # 估计里程增量分布
    delta_x_dist = estimate_transition(data)  # 假设已实现

    # 第二步：构造前向价值并估计
    def pseudo_likelihood(theta):
        RC, tc1, tc2 = theta
        n_x = len(ccp)

        # 模拟前向价值
        V_forward = np.zeros(n_x)
        for x in range(n_x):
            # 模拟从x出发、使用当前CCP的期望贴现效用
            V_forward[x] = simulate_forward_value(x, ccp, theta, beta,
                                                   delta_x_dist, n_forward)

        # 构造伪CCP
        u0 = -tc1 * np.arange(n_x) - tc2 * np.arange(n_x)**2
        u1 = -RC + np.zeros(n_x)

        v0 = u0 + beta * expected_continuation(V_forward, 0, delta_x_dist)
        v1 = u1 + beta * V_forward[0]  # 更换后从0开始

        pseudo_ccp = 1 / (1 + np.exp(v0 - v1))

        # 伪似然
        x = data['x'].values
        d = data['d'].values
        ll = d * np.log(pseudo_ccp[x]) + (1-d) * np.log(1-pseudo_ccp[x])
        return -np.sum(ll)

    result = minimize(pseudo_likelihood, [10, 0.005, 0.0001], method='Nelder-Mead')
    return result.x

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

SMM算法：
─────────────────────────────────────────────────────
for each θ in optimization:
    1. 解模型（可选：解不动点或用CCP）
    2. 模拟S条路径，每条长度T
    3. 计算模拟矩：m̂_sim(θ) = (1/S)Σ m(x^s, d^s)
    4. 计算矩差：g(θ) = m̂_data - m̂_sim(θ)
    5. 计算目标：Q(θ) = g(θ)' W g(θ)

最小化 Q(θ) 得到 θ̂_SMM
─────────────────────────────────────────────────────

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

def smm_objective(theta, data_moments, beta, G, n_sim=100, T_sim=200, seed=42):
    """
    SMM目标函数
    data_moments: 从真实数据计算的矩
    """
    np.random.seed(seed)  # 固定模拟随机性（CUE估计时重要）

    # 解模型
    EV = solve_EV(theta, beta, G)
    ccp = compute_ccp(EV, theta, beta)

    # 模拟数据
    sim_d = []
    sim_x = []
    for s in range(n_sim):
        x = 0  # 初始里程
        for t in range(T_sim):
            p1 = ccp[min(x, len(ccp)-1)]
            d = np.random.binomial(1, p1)
            sim_d.append(d)
            sim_x.append(x)
            # 转移
            if d == 1:
                x = sample_delta_x(G)  # 更换后从增量开始
            else:
                x = min(x + sample_delta_x(G), len(ccp)-1)

    sim_x = np.array(sim_x)
    sim_d = np.array(sim_d)

    # 计算模拟矩
    sim_moments = np.array([
        sim_d.mean(),                    # 更换率
        (sim_d * sim_x).mean(),          # 更换×里程
        sim_x.mean(),                    # 平均里程
        ((sim_x > 50) * sim_d).mean(),   # 高里程时的更换率
    ])

    # 矩差（需要匹配维度和类型）
    g = data_moments - sim_moments

    # 目标函数（这里用单位矩阵作为权重）
    return g @ g

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

from sklearn.linear_model import LogisticRegression

def indirect_inference(theta, real_data, beta, G, n_sim=50):
    """
    间接推断目标函数
    """
    # 在真实数据上估计辅助模型
    aux_model_real = LogisticRegression()
    X_real = real_data[['x', 'x_sq']].values  # x 和 x^2
    y_real = real_data['d'].values
    aux_model_real.fit(X_real, y_real)
    psi_real = np.concatenate([aux_model_real.intercept_, aux_model_real.coef_.flatten()])

    # 模拟数据
    sim_data = simulate_data(theta, beta, G, n_obs=len(real_data)*n_sim)

    # 在模拟数据上估计辅助模型
    aux_model_sim = LogisticRegression()
    X_sim = sim_data[['x', 'x_sq']].values
    y_sim = sim_data['d'].values
    aux_model_sim.fit(X_sim, y_sim)
    psi_sim = np.concatenate([aux_model_sim.intercept_, aux_model_sim.coef_.flatten()])

    # 目标函数
    return np.sum((psi_real - psi_sim)**2)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

import pymc as pm
import pytensor.tensor as pt

def bayesian_rust_estimation(data, beta, G, n_draws=5000, n_tune=1000):
    """
    贝叶斯NFXP（使用PyMC）
    注意：这是简化版本，真实应用需要更细致的处理
    """
    with pm.Model() as model:
        # 先验
        RC = pm.HalfNormal('RC', sigma=20)
        tc1 = pm.HalfNormal('tc1', sigma=0.01)
        tc2 = pm.HalfNormal('tc2', sigma=0.001)

        # 似然（需要用pytensor实现EV求解器——这里简化处理）
        # 实际中，可能需要用pm.Potential + 自定义似然
        theta = pt.stack([RC, tc1, tc2])

        # 假设有函数 compute_ll_pytensor 用pytensor实现
        ll = compute_ll_pytensor(theta, data, beta, G)
        pm.Potential('likelihood', ll)

        # 采样
        trace = pm.sample(n_draws, tune=n_tune, cores=4)

    return trace

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

def discretize_state(n_points=100, x_max=450, method='uniform'):
    """
    状态空间离散化（里程）
    """
    if method == 'uniform':
        return np.linspace(0, x_max, n_points)
    elif method == 'chebyshev':
        # 切比雪夫节点，适合插值
        k = np.arange(1, n_points+1)
        nodes = 0.5 * (x_max) * (1 - np.cos((2*k - 1) * np.pi / (2*n_points)))
        return nodes

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

def construct_transition_matrix(x_grid, delta_x_probs, x_max):
    """
    构造转移概率矩阵
    delta_x_probs: 里程增量的概率分布（如 [p_0, p_1, p_2]）
    """
    n_x = len(x_grid)
    n_delta = len(delta_x_probs)

    # 不更换时的转移矩阵
    G0 = np.zeros((n_x, n_x))
    for i, x in enumerate(x_grid):
        for delta, p in enumerate(delta_x_probs):
            new_x = x + delta
            j = np.searchsorted(x_grid, new_x)  # 找最近的网格点
            j = min(j, n_x - 1)  # 边界处理
            G0[i, j] += p

    # 更换时的转移矩阵（从0开始）
    G1 = np.zeros((n_x, n_x))
    for delta, p in enumerate(delta_x_probs):
        j = min(delta, n_x - 1)
        G1[:, j] = p  # 所有状态转移到相同的里程增量

    return G0, G1

1
2
3
4
5
6
7

def logsumexp_stable(a, b):
    """数值稳定的log-sum-exp"""
    max_val = np.maximum(a, b)
    return max_val + np.log1p(np.exp(-np.abs(a - b)))

# 或使用scipy
from scipy.special import logsumexp

1
2
3

def safe_prob(p, eps=1e-10):
    """避免log(0)和log(1)"""
    return np.clip(p, eps, 1 - eps)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

def multi_start_optimization(objective, bounds, n_starts=20, method='L-BFGS-B'):
    """
    多初值优化
    """
    results = []

    # 拉丁超立方采样初值
    from scipy.stats.qmc import LatinHypercube
    sampler = LatinHypercube(d=len(bounds))
    samples = sampler.random(n=n_starts)

    for i in range(n_starts):
        # 将[0,1]映射到bounds
        x0 = np.array([b[0] + samples[i,j]*(b[1]-b[0])
                       for j, b in enumerate(bounds)])

        try:
            res = minimize(objective, x0, method=method, bounds=bounds)
            results.append(res)
        except:
            continue

    # 选择最优
    best = min(results, key=lambda r: r.fun)

    # 检查收敛一致性
    top_3 = sorted(results, key=lambda r: r.fun)[:3]
    print("Top 3 function values:", [r.fun for r in top_3])
    print("Top 3 parameters:", [r.x for r in top_3])

    return best

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

def diagnose_convergence(result, objective, bounds, n_restart=5):
    """
    诊断优化收敛质量
    """
    diagnostics = {}

    # 1. 检查梯度是否接近零
    from scipy.optimize import approx_fprime
    grad = approx_fprime(result.x, objective, 1e-8)
    diagnostics['grad_norm'] = np.linalg.norm(grad)

    # 2. 检查Hessian正定性
    from numdifftools import Hessian
    hess_func = Hessian(objective)
    H = hess_func(result.x)
    eigenvalues = np.linalg.eigvalsh(H)
    diagnostics['hessian_positive_definite'] = np.all(eigenvalues > 0)
    diagnostics['min_eigenvalue'] = eigenvalues.min()

    # 3. 从最优解附近重新优化
    restart_results = []
    for _ in range(n_restart):
        x0 = result.x + np.random.normal(0, 0.01, size=len(result.x))
        x0 = np.clip(x0, [b[0] for b in bounds], [b[1] for b in bounds])
        res = minimize(objective, x0, method='L-BFGS-B', bounds=bounds)
        restart_results.append(res.fun)

    diagnostics['restart_variation'] = np.std(restart_results)

    return diagnostics

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

from joblib import Parallel, delayed

def parallel_smm(thetas, data_moments, beta, G, n_jobs=-1):
    """
    并行计算多个参数点的SMM目标值
    """
    def eval_one(theta):
        return smm_objective(theta, data_moments, beta, G)

    results = Parallel(n_jobs=n_jobs)(
        delayed(eval_one)(theta) for theta in thetas
    )

    return np.array(results)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22

def compute_mle_se(theta_hat, neg_loglik, data):
    """
    MLE标准误（基于Hessian）
    """
    from numdifftools import Hessian

    # 数值Hessian
    hess_func = Hessian(lambda t: neg_loglik(t, data))
    H = hess_func(theta_hat)

    # 信息矩阵 = -Hessian / N
    I_hat = H / len(data)

    # 方差 = 信息矩阵逆 / N
    try:
        var_theta = np.linalg.inv(I_hat) / len(data)
        se = np.sqrt(np.diag(var_theta))
    except np.linalg.LinAlgError:
        print("Warning: Hessian is singular!")
        se = np.full(len(theta_hat), np.nan)

    return se

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

def bootstrap_se(data, estimator_func, n_bootstrap=200, n_jobs=-1):
    """
    Bootstrap标准误
    """
    def one_bootstrap(seed):
        np.random.seed(seed)
        # 重抽样
        idx = np.random.choice(len(data), size=len(data), replace=True)
        boot_data = data.iloc[idx]
        # 估计
        try:
            return estimator_func(boot_data)
        except:
            return np.full(len(data.columns), np.nan)

    results = Parallel(n_jobs=n_jobs)(
        delayed(one_bootstrap)(s) for s in range(n_bootstrap)
    )

    results = np.array([r for r in results if not np.any(np.isnan(r))])

    se = np.std(results, axis=0)
    ci_low = np.percentile(results, 2.5, axis=0)
    ci_high = np.percentile(results, 97.5, axis=0)

    return se, ci_low, ci_high

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

def plot_model_fit(data, theta_hat, beta, G):
    """
    可视化模型拟合
    """
    import matplotlib.pyplot as plt

    # 计算模型预测的CCP
    EV = solve_EV(theta_hat, beta, G)
    ccp_model = compute_ccp(EV, theta_hat, beta)

    # 数据中的CCP
    ccp_data = data.groupby('x')['d'].mean()

    fig, ax = plt.subplots(figsize=(10, 6))

    ax.scatter(ccp_data.index, ccp_data.values, alpha=0.5, label='Data')
    ax.plot(np.arange(len(ccp_model)), ccp_model, 'r-', linewidth=2, label='Model')

    ax.set_xlabel('Mileage')
    ax.set_ylabel('Replacement Probability')
    ax.legend()
    ax.set_title('Model Fit: Observed vs. Predicted CCP')

    return fig

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

def goodness_of_fit(data, theta_hat, beta, G):
    """
    拟合优度检验
    """
    # 模型预测
    EV = solve_EV(theta_hat, beta, G)
    ccp_model = compute_ccp(EV, theta_hat, beta)

    # 计算各类统计量
    x = data['x'].values
    d = data['d'].values

    # 预测概率
    pred_prob = ccp_model[x]

    # Log-likelihood
    ll = np.sum(d * np.log(pred_prob + 1e-10) + (1-d) * np.log(1 - pred_prob + 1e-10))

    # Pseudo R-squared (McFadden)
    ll_null = len(data) * (d.mean() * np.log(d.mean()) + (1-d.mean()) * np.log(1-d.mean()))
    pseudo_r2 = 1 - ll / ll_null

    # 按里程分组的拟合误差
    data['pred'] = pred_prob
    data['x_bin'] = pd.cut(data['x'], bins=10)
    fit_by_bin = data.groupby('x_bin').agg({'d': 'mean', 'pred': 'mean'})
    rmse = np.sqrt(((fit_by_bin['d'] - fit_by_bin['pred'])**2).mean())

    return {
        'log_likelihood': ll,
        'pseudo_r2': pseudo_r2,
        'rmse_by_bin': rmse
    }

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

def j_test(data, theta_hat, moments_func, W):
    """
    J统计量过度识别检验
    """
    N = len(data)
    g = moments_func(theta_hat, data)  # 矩条件向量
    n_moments = len(g)
    n_params = len(theta_hat)

    # J统计量
    J = N * g @ W @ g

    # 自由度 = 过度识别的数量
    df = n_moments - n_params

    # p值
    from scipy.stats import chi2
    p_value = 1 - chi2.cdf(J, df)

    return {'J': J, 'df': df, 'p_value': p_value}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

def sensitivity_analysis(data, baseline_theta, variations):
    """
    对关键假设的敏感性分析
    """
    results = []

    # 基准结果
    baseline_result = estimate_model(data, beta=0.9999)
    results.append({'spec': 'baseline', **dict(zip(['RC', 'tc1', 'tc2'], baseline_result))})

    for var_name, var_value in variations.items():
        if var_name == 'beta':
            result = estimate_model(data, beta=var_value)
        elif var_name == 'cost_functional_form':
            result = estimate_model_alt_cost(data, form=var_value)
        # ... 其他变体

        results.append({'spec': f'{var_name}={var_value}',
                       **dict(zip(['RC', 'tc1', 'tc2'], result))})

    return pd.DataFrame(results)

# 使用示例
variations = {
    'beta': [0.99, 0.999, 0.9999],
    'cost_functional_form': ['linear', 'quadratic', 'cubic'],
    'error_distribution': ['logit', 'probit']
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

## 稳健性检查清单

### 模型规格

-   [ ] 成本函数形式：线性 vs 二次 vs 非参数
-   [ ] 折现因子：β ∈ {0.95, 0.99, 0.999, 0.9999}
-   [ ] 误差分布：Logit vs Probit

### 数据处理

-   [ ] 去除极端值（<1%, >99%分位）
-   [ ] 不同的状态离散化粒度
-   [ ] 不同的样本期间

### 估计方法

-   [ ] NFXP vs CCP vs GMM 结果比较
-   [ ] 不同初值的收敛性
-   [ ] Bootstrap vs Hessian 标准误

### 拟合质量

-   [ ] 整体更换率拟合
-   [ ] 分里程段的更换率拟合
-   [ ] 样本外预测（如适用）

1
2
3
4
5
6
7
8
9

反事实模拟流程：
─────────────────────────────────────────────────────
1. 用估计的 θ̂ 作为真实偏好参数
2. 修改政策环境（如成本函数、转移概率）
3. 在新环境下重新解模型（解新的贝尔曼方程）
4. 计算新均衡下的行为/福利
5. 比较 counterfactual - baseline
6. 报告不确定性（如用 bootstrap）
─────────────────────────────────────────────────────

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

def counterfactual_cost_subsidy(theta_hat, beta, G, subsidy_rate=0.5):
    """
    反事实：更换成本补贴
    """
    RC, tc1, tc2 = theta_hat

    # 基准情况
    EV_base = solve_EV(theta_hat, beta, G)
    ccp_base = compute_ccp(EV_base, theta_hat, beta)
    avg_replacement_rate_base = ccp_base.mean()  # 简化：均匀分布状态

    # 反事实：RC降低
    RC_cf = RC * (1 - subsidy_rate)
    theta_cf = (RC_cf, tc1, tc2)

    EV_cf = solve_EV(theta_cf, beta, G)
    ccp_cf = compute_ccp(EV_cf, theta_cf, beta)
    avg_replacement_rate_cf = ccp_cf.mean()

    # 结果
    return {
        'baseline_replacement_rate': avg_replacement_rate_base,
        'counterfactual_replacement_rate': avg_replacement_rate_cf,
        'change': avg_replacement_rate_cf - avg_replacement_rate_base,
        'change_percent': (avg_replacement_rate_cf - avg_replacement_rate_base) / avg_replacement_rate_base * 100
    }

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

def counterfactual_tech_improvement(theta_hat, beta, G, cost_reduction=0.3):
    """
    反事实：维护成本降低
    """
    RC, tc1, tc2 = theta_hat

    # 反事实：维护成本参数降低
    tc1_cf = tc1 * (1 - cost_reduction)
    tc2_cf = tc2 * (1 - cost_reduction)
    theta_cf = (RC, tc1_cf, tc2_cf)

    # 解反事实模型
    EV_cf = solve_EV(theta_cf, beta, G)
    ccp_cf = compute_ccp(EV_cf, theta_cf, beta)

    # 模拟长期平稳分布
    stationary_dist = simulate_stationary_distribution(ccp_cf, G, n_sim=100000)

    # 计算平均里程和更换频率
    avg_mileage_cf = np.dot(stationary_dist, np.arange(len(stationary_dist)))
    avg_replacement_rate_cf = np.dot(stationary_dist, ccp_cf)

    return {
        'avg_mileage_counterfactual': avg_mileage_cf,
        'avg_replacement_rate_counterfactual': avg_replacement_rate_cf,
        'stationary_distribution': stationary_dist
    }

def simulate_stationary_distribution(ccp, G, n_sim=100000, T=500):
    """
    模拟长期平稳分布
    """
    G0, G1 = G
    n_x = len(ccp)

    # 构造联合转移矩阵
    # P(x'|x) = p(d=0|x) * G0(x,x') + p(d=1|x) * G1(x,x')
    P = np.zeros((n_x, n_x))
    for x in range(n_x):
        P[x, :] = (1 - ccp[x]) * G0[x, :] + ccp[x] * G1[x, :]

    # 找平稳分布（特征向量法）
    eigenvalues, eigenvectors = np.linalg.eig(P.T)
    idx = np.argmin(np.abs(eigenvalues - 1))
    stationary = np.real(eigenvectors[:, idx])
    stationary = stationary / stationary.sum()

    return stationary

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

def counterfactual_with_uncertainty(data, estimator_func, cf_func, n_bootstrap=200):
    """
    带不确定性的反事实分析
    """
    # 点估计
    theta_hat = estimator_func(data)
    cf_point = cf_func(theta_hat)

    # Bootstrap
    cf_bootstrap = []
    for b in range(n_bootstrap):
        # 重抽样
        boot_idx = np.random.choice(len(data), size=len(data), replace=True)
        boot_data = data.iloc[boot_idx]

        try:
            theta_b = estimator_func(boot_data)
            cf_b = cf_func(theta_b)
            cf_bootstrap.append(cf_b)
        except:
            continue

    # 整理结果
    cf_array = np.array([list(cf.values()) for cf in cf_bootstrap])

    results = {}
    for i, key in enumerate(cf_point.keys()):
        results[key] = {
            'point': cf_point[key],
            'se': np.std(cf_array[:, i]),
            'ci_low': np.percentile(cf_array[:, i], 2.5),
            'ci_high': np.percentile(cf_array[:, i], 97.5)
        }

    return results

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

## 反事实分析结果

### 情景：更换成本补贴 50%

| 指标 | 基准 | 反事实 | 变化 | 95% CI |
| --- | --- | --- | --- | --- |
| 更换率 | 3.2% | 5.8% | +2.6pp | [2.1, 3.2] |
| 平均里程 | 28.5 | 21.3 | -7.2 | [-9.1, -5.4] |
| 年均更换次数 | 1.2 | 2.0 | +0.8 | [0.6, 1.0] |

### 主要发现

1. **更换率大幅上升**：50%的成本补贴使更换率几乎翻倍（从 3.2%到 5.8%）

2. **车队里程分布左移**：由于更频繁更换，平均里程下降 25%

3. **政策成本**：假设车队规模 100 辆，年补贴支出约为...

### 注意事项

-   反事实假设政策对维护成本无溢出效应
-   未考虑一般均衡效应（如发动机价格变化）
-   置信区间基于 200 次 bootstrap