1
2
3
4
5
6
7
8

例子：
- 状态：你的旧车里程已经15万公里
- 决策：是否进行大修
- CCP：在这个里程下，选择大修的概率 = ?

数据观察：
  在里程为15万公里的100辆车中，有45辆选择大修
  → CCP(d=1|x=15) = 45% = 0.45

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

里程都是15万的两辆车：
  车A → 选择大修
  车B → 不大修

为什么不同？

原因1：我们看不到所有信息
  - 车主的财务状况
  - 对可靠性的偏好
  - 计划卖车的时间
  - 等等...

这些看不到的因素 = 私人冲击 ε

1
2
3
4
5
6
7

Car_ID  Year  Mileage  Maintenance_Decision
  1     2015   50000          0   (不维修)
  1     2016   75000          1   (维修)
  1     2017  100000          0   (不维修)
  2     2015   60000          1   (维修)
  2     2016   85000          0   (不维修)
  ...

 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 compute_empirical_ccp(data, x_grid):
    """
    从原始数据计算经验CCP
    """

    # 方法1：分箱（简单）
    ccp_binned = []

    for x in x_grid:
        # 找状态在x附近的所有观测
        in_bin = (data['mileage'] >= x-500) & (data['mileage'] < x+500)

        if in_bin.sum() > 0:
            # 该箱中维修的比例 = CCP估计
            ccp_val = data.loc[in_bin, 'maintenance'].mean()
        else:
            ccp_val = np.nan

        ccp_binned.append(ccp_val)

    # 方法2：光滑估计（Lowess平滑）
    from statsmodels.nonparametric.smoothers_lowess import lowess

    smoothed = lowess(
        data['maintenance'],      # y: 是否维修
        data['mileage'],          # x: 里程
        frac=0.2                  # 带宽
    )

    # 在网格上插值
    ccp_smooth = np.interp(x_grid, smoothed[:, 0], smoothed[:, 1])

    return ccp_smooth

1
2
3
4
5
6
7
8

CCP(d=1|x)
    ↑
  1.0 |                    ╱
      |                  ╱
  0.5 |                ╱
      |              ╱
  0.0 |────────────╱─────→ 状态 x
      0        x*        max

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

参数 θ
   ↓
求解HJB/Bellman
   ↓
计算理论CCP(x; θ)
   ↓
与经验CCP比较
   ↓
调整θ
   ↓
重复直到拟合

1
2
3
4
5
6
7
8

如果 V_0(x) >> V_1(x)：
  P(d=1|x) ≈ 0  (几乎不选d=1)

如果 V_0(x) ≈ V_1(x)：
  P(d=1|x) ≈ 0.5  (无差异)

如果 V_0(x) << V_1(x)：
  P(d=1|x) ≈ 1  (肯定选d=1)

  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
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164

import numpy as np
import pandas as pd
from statsmodels.nonparametric.smoothers_lowess import lowess

class NonparametricCCPEstimation:
    """
    完全非参数的CCP估计
    不需要任何模型假设！
    """

    def __init__(self, data):
        """
        data: DataFrame
        需要列: 'state' (状态), 'choice' (决策)
        """
        self.data = data

    def estimate_ccp_binning(self, x_grid, bin_width=1):
        """
        方法1：简单分箱

        优点：直观
        缺点：粗糙，样本小时不稳定
        """

        ccp = np.zeros(len(x_grid))

        for i, x in enumerate(x_grid):
            # 找状态在[x - bin_width/2, x + bin_width/2]内的观测
            in_bin = (self.data['state'] >= x - bin_width/2) & \
                     (self.data['state'] < x + bin_width/2)

            if in_bin.sum() > 0:
                ccp[i] = self.data.loc[in_bin, 'choice'].mean()
            else:
                ccp[i] = np.nan

        return ccp

    def estimate_ccp_lowess(self, x_grid, frac=0.2):
        """
        方法2：局部加权回归(Lowess)

        优点：光滑，自适应
        缺点：需要选择带宽参数
        """

        # Lowess拟合
        smoothed = lowess(
            self.data['choice'],   # y: 选择(0或1)
            self.data['state'],    # x: 状态
            frac=frac              # 局部窗口大小
        )

        # 在指定网格上插值
        ccp = np.interp(x_grid, smoothed[:, 0], smoothed[:, 1])

        return ccp

    def estimate_ccp_kernel_regression(self, x_grid, bandwidth=1):
        """
        方法3：核回归

        CCP(x) = E[choice | state = x]
               = Σ K((state_i - x) / h) · choice_i / Σ K((state_i - x) / h)

        其中K是核函数(如高斯核), h是带宽
        """

        from scipy.stats import gaussian_kde

        ccp = np.zeros(len(x_grid))

        for i, x in enumerate(x_grid):
            # 计算核权重
            distances = np.abs(self.data['state'] - x)
            weights = np.exp(-(distances / bandwidth) ** 2)

            # 加权平均
            ccp[i] = np.average(self.data['choice'], weights=weights)

        return ccp

    def bootstrap_ccp_confidence_interval(self, x_grid, n_bootstrap=500, alpha=0.05):
        """
        Bootstrap置信区间
        """

        ccp_bootstraps = np.zeros((n_bootstrap, len(x_grid)))

        for b in range(n_bootstrap):
            # 有放回重抽样
            idx = np.random.choice(len(self.data), size=len(self.data), replace=True)
            data_boot = self.data.iloc[idx]

            # 计算bootstrap样本的CCP
            ccp_boot = self._estimate_ccp_lowess_internal(data_boot, x_grid)
            ccp_bootstraps[b, :] = ccp_boot

        # 分位数
        lower = np.percentile(ccp_bootstraps, 100 * alpha / 2, axis=0)
        upper = np.percentile(ccp_bootstraps, 100 * (1 - alpha / 2), axis=0)

        return lower, upper

    def _estimate_ccp_lowess_internal(self, data, x_grid):
        """内部函数"""
        smoothed = lowess(data['choice'], data['state'], frac=0.2)
        return np.interp(x_grid, smoothed[:, 0], smoothed[:, 1])

    def plot_empirical_ccp(self, x_grid, methods=['binning', 'lowess']):
        """可视化"""

        import matplotlib.pyplot as plt

        plt.figure(figsize=(12, 6))

        if 'binning' in methods:
            ccp_bin = self.estimate_ccp_binning(x_grid)
            plt.plot(x_grid, ccp_bin, 'o-', label='Binning', alpha=0.6)

        if 'lowess' in methods:
            ccp_lowess = self.estimate_ccp_lowess(x_grid)
            plt.plot(x_grid, ccp_lowess, 's-', label='Lowess', alpha=0.6)

        # 原始数据散点
        plt.scatter(self.data['state'], self.data['choice'],
                   alpha=0.1, s=10, label='Raw data')

        plt.xlabel('状态 (State)')
        plt.ylabel('P(选择 d=1 | State)')
        plt.title('条件选择概率 (Conditional Choice Probability)')
        plt.legend()
        plt.grid(True, alpha=0.3)
        plt.ylim(-0.1, 1.1)
        plt.show()

# 示例
np.random.seed(42)

# 生成模拟数据：里程与维修决策
n = 1000
state = np.random.uniform(0, 450, n)  # 里程: 0到45万

# 真实CCP（S形曲线）
true_ccp = 1 / (1 + np.exp(-(state - 225) / 50))

# 加噪声生成决策
choice = (np.random.uniform(0, 1, n) < true_ccp).astype(int)

data = pd.DataFrame({'state': state, 'choice': choice})

# 估计
est = NonparametricCCPEstimation(data)
est.plot_empirical_ccp(np.linspace(0, 450, 100), methods=['binning', 'lowess'])

# 置信区间
x_test = np.array([100, 225, 350])
lower, upper = est.bootstrap_ccp_confidence_interval(x_test)

print("状态        CCP估计    95% 置信区间")
for x, l, u in zip(x_test, lower, upper):
    ccp_val = est.estimate_ccp_lowess(np.array([x]))[0]
    print(f"{x:3.0f}        {ccp_val:.4f}    [{l:.4f}, {u:.4f}]")

1
2
3
4

状态        CCP估计    95% 置信区间
100        0.1234    [0.0856, 0.1612]
225        0.5023    [0.4567, 0.5479]
350        0.8901    [0.8534, 0.9268]

  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
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138

class SemiparametricCCPEstimation:
    """
    Hotz-Miller (1993) 方法

    步骤：
    1. 非参数估计CCP
    2. 反演得到价值函数
    3. 估计参数（无需重复求解动态规划）
    """

    def __init__(self, data, x_grid, beta=0.95):
        self.data = data
        self.x_grid = x_grid
        self.beta = beta

        # Step 1: 非参数CCP估计
        from statsmodels.nonparametric.smoothers_lowess import lowess
        smoothed = lowess(data['choice'], data['state'], frac=0.2)
        self.ccp_empirical = np.interp(x_grid, smoothed[:, 0], smoothed[:, 1])

    def invert_ccp_to_value_function(self):
        """
        Step 2: 从CCP反演价值函数

        Hotz-Miller反演公式：

        从logit CCP，我们知道：
        P(d=1|x) = exp(V_1(x)) / [exp(V_0(x)) + exp(V_1(x))]

        定义：
        Λ_d(x) = log(P(d|x)) - ψ + log(1 - P(d|x)) + ψ
               = log(P(d|x) / (1 - P(d|x)))

        其中ψ是Euler-Mascheroni常数 ≈ 0.5772

        则：
        V_1(x) - V_0(x) = Λ_1(x) - Λ_0(x)
                        = log(P(d=1|x) / (1 - P(d=1|x)))

        通过递归和固定点方程，可反演出V_d(x)
        """

        # Euler-Mascheroni常数
        psi = 0.5772

        # 计算选择特定价值的对数比率
        p = self.ccp_empirical
        p = np.clip(p, 1e-10, 1 - 1e-10)  # 避免log(0)

        log_odds = np.log(p / (1 - p))

        # 固定点方程：
        # V_0(x) = u(x,0) + β [p(x) * V(0) + (1-p(x)) * V(x')]
        # V_1(x) = u(x,1) + β V(0)
        #
        # 解这个系统得到V(x)

        # 简化处理（假设线性成本函数）
        c = 0.1 * self.x_grid  # 假设成本函数

        # 初始化
        V = np.zeros_like(self.x_grid)

        # 迭代求解固定点
        for iteration in range(100):
            V_old = V.copy()

            # 更新V_0(x) = -c(x) + β E[V(x')|d=0]
            V_0 = -c + self.beta * np.roll(V_old, 1)  # x' = x+1

            # 更新V_1(x) = -RC + β V(0)
            # (假设维修后回到状态0)
            RC = 10  # 临时假设
            V_1 = -RC + self.beta * V_old[0]

            # CCP形式的价值函数
            # V(x) = (1/σ) log(exp(σ V_0) + exp(σ V_1))
            # 其中σ是logit的scale（与离散选择冲击分布有关）

            sigma = 100  # 假设
            V = (1 / sigma) * np.log(
                np.exp(sigma * V_0) + np.exp(sigma * V_1)
            )

            if np.max(np.abs(V - V_old)) < 1e-6:
                break

        self.V_inverted = V
        return V

    def estimate_parameters(self):
        """
        Step 3: 估计模型参数

        无需重复求解动态规划！
        """

        # 给定反演的V(x)，可以直接读出：
        # log(P(d=1|x)) - log(P(d=0|x)) = V_1(x) - V_0(x)

        p = self.ccp_empirical
        p = np.clip(p, 1e-10, 1 - 1e-10)

        # CCP导出的价值差异
        log_odds = np.log(p / (1 - p))

        # log_odds(x) = -RC + c(x) - c(0) + β[V(0) - V(x)]

        # 可以用回归估计参数
        from sklearn.linear_model import LinearRegression

        # 设计矩阵
        X = np.column_stack([
            np.ones(len(self.x_grid)),  # 截距（包含-RC）
            self.x_grid                  # 线性成本项
        ])

        # 回归：log_odds ~ α_0 + α_1 * x
        model = LinearRegression()
        model.fit(X, log_odds)

        # 从系数反演参数
        alpha_0, alpha_1 = model.intercept_, model.coef_[1]

        # α_1对应c(x)的系数
        c1_estimated = alpha_1

        # α_0包含-RC和价值函数项
        # 需要额外假设来分离...

        print(f"估计的成本函数斜率 c1 = {c1_estimated:.6f}")

        return c1_estimated

# 使用
est_semi = SemiparametricCCPEstimation(data, np.linspace(0, 450, 100))
V_inverted = est_semi.invert_ccp_to_value_function()
c1_est = est_semi.estimate_parameters()

  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
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182

class StructuralMaximumLikelihood:
    """
    参数结构估计

    给定参数θ，计算理论CCP，最大化似然函数
    """

    def __init__(self, data, x_grid, beta=0.95):
        self.data = data
        self.x_grid = x_grid
        self.beta = beta

    def solve_bellman(self, params):
        """
        给定参数，求解贝尔曼方程得到价值函数V(x)

        params = (c1, c2, RC)
        """

        c1, c2, RC = params

        # 成本函数
        def cost(x):
            return c1 * x + c2 * x**2

        # 初始化
        V = np.zeros_like(self.x_grid)

        # 迭代求解（反向迭代）
        for iteration in range(200):
            V_old = V.copy()

            for i in range(len(self.x_grid) - 1, -1, -1):
                x = self.x_grid[i]

                # 不维修：支付维护成本，状态增加1
                flow_0 = -cost(x)
                x_next_0 = min(x + 1, self.x_grid[-1])
                idx_next_0 = np.argmin(np.abs(self.x_grid - x_next_0))
                V_next_0 = V_old[idx_next_0]

                bar_v_0 = flow_0 + self.beta * V_next_0

                # 维修：支付维修成本，状态回到0
                flow_1 = -RC - cost(0)
                bar_v_1 = flow_1 + self.beta * V_old[0]

                # 用logit结合（处理离散选择冲击）
                sigma = 1000  # 冲击分布scale
                V[i] = (1 / sigma) * np.log(
                    np.exp(sigma * bar_v_0) + np.exp(sigma * bar_v_1)
                )

            # 收敛检查
            if np.max(np.abs(V - V_old)) < 1e-8:
                break

        return V

    def compute_theoretical_ccp(self, params):
        """
        给定参数，计算理论CCP
        """

        c1, c2, RC = params

        # 求解贝尔曼方程
        V = self.solve_bellman(params)

        # 计算CCP
        ccp = np.zeros_like(self.x_grid)

        for i, x in enumerate(self.x_grid):
            # 不维修的价值
            flow_0 = -(c1 * x + c2 * x**2)
            x_next_0 = min(x + 1, self.x_grid[-1])
            idx_next_0 = np.argmin(np.abs(self.x_grid - x_next_0))
            bar_v_0 = flow_0 + self.beta * V[idx_next_0]

            # 维修的价值
            flow_1 = -RC - (c1 * 0 + c2 * 0**2)
            bar_v_1 = flow_1 + self.beta * V[0]

            # Logit CCP
            ccp[i] = 1 / (1 + np.exp(bar_v_0 - bar_v_1))

        return ccp

    def compute_log_likelihood(self, params):
        """
        计算对数似然函数

        L(θ) = Σ_i [d_i * log(P(d=1|x_i;θ)) + (1-d_i) * log(1-P(d=1|x_i;θ))]
        """

        try:
            # 计算理论CCP
            ccp_theory = self.compute_theoretical_ccp(params)

            # 为每个数据点查找对应的CCP
            x_data = self.data['state'].values
            d_data = self.data['choice'].values

            ccp_at_data = np.interp(x_data, self.x_grid, ccp_theory)

            # 数值稳定性：限制CCP在[1e-10, 1-1e-10]
            ccp_at_data = np.clip(ccp_at_data, 1e-10, 1 - 1e-10)

            # 二项对数似然
            ll = d_data * np.log(ccp_at_data) + \
                 (1 - d_data) * np.log(1 - ccp_at_data)

            return np.sum(ll)

        except:
            # 如果计算失败，返回-∞
            return -1e10

    def estimate(self, initial_params=None, bounds=None):
        """
        最大似然估计
        """

        if initial_params is None:
            initial_params = np.array([0.1, 0.01, 10.0])

        if bounds is None:
            bounds = [(0.001, 0.5), (0.0001, 0.1), (1.0, 50.0)]

        # 定义目标函数（负对数似然）
        def objective(params):
            return -self.compute_log_likelihood(params)

        # 优化
        from scipy.optimize import minimize

        print("开始最大似然估计...")
        print(f"初始参数: c1={initial_params[0]:.4f}, "
              f"c2={initial_params[1]:.6f}, RC={initial_params[2]:.4f}")

        result = minimize(
            objective,
            initial_params,
            bounds=bounds,
            method='L-BFGS-B',
            options={'ftol': 1e-10, 'maxiter': 500}
        )

        print(f"\n估计完成！")
        print(f"最优参数: c1={result.x[0]:.4f}, "
              f"c2={result.x[1]:.6f}, RC={result.x[2]:.4f}")
        print(f"最大对数似然: {-result.fun:.2f}")

        self.params_estimated = result.x

        return result

# 使用
est_ml = StructuralMaximumLikelihood(data, np.linspace(0, 450, 100))
result = est_ml.estimate(initial_params=[0.1, 0.01, 10.0])

# 比较：理论CCP vs 经验CCP
ccp_theory = est_ml.compute_theoretical_ccp(result.x)

# 非参数经验CCP
from statsmodels.nonparametric.smoothers_lowess import lowess
smoothed = lowess(data['choice'], data['state'], frac=0.2)
ccp_empirical = np.interp(np.linspace(0, 450, 100),
                          smoothed[:, 0], smoothed[:, 1])

# 绘制
plt.figure(figsize=(12, 6))
plt.plot(np.linspace(0, 450, 100), ccp_empirical, 'o-',
         label='Empirical CCP', linewidth=2)
plt.plot(np.linspace(0, 450, 100), ccp_theory, 's-',
         label='Theoretical CCP (estimated)', linewidth=2)
plt.xlabel('状态 (里程)')
plt.ylabel('维修概率')
plt.title('CCP拟合效果')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

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

# 检验性质1：相同状态下的CCP应该稳定
def test_ccp_property1(data):
    """
    在相同或相近的状态下，CCP应该相同
    """

    # 选择某个状态（如x=225）
    x_target = 225

    # 找所有状态在[x-5, x+5]的观测
    near_target = (data['state'] >= x_target - 5) & \
                  (data['state'] <= x_target + 5)

    choices_near = data.loc[near_target, 'choice'].values

    # CCP应该是这些选择的平均值
    ccp_estimated = choices_near.mean()

    print(f"在状态 x={x_target} 附近的观测数: {len(choices_near)}")
    print(f"选择d=1的比例: {ccp_estimated:.4f}")
    print(f"标准误: {np.sqrt(ccp_estimated * (1-ccp_estimated) / len(choices_near)):.4f}")

 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

# 检验性质2：CCP的平滑性
def compute_ccp_smoothness(ccp_empirical, x_grid):
    """
    衡量CCP的"粗糙度"
    """

    # 二阶差分（曲率）
    second_differences = np.diff(ccp_empirical, n=2)

    # 粗糙度指标
    roughness = np.sum(second_differences ** 2)

    print(f"CCP粗糙度指标: {roughness:.6f}")

    # 可视化
    plt.figure(figsize=(12, 6))

    plt.subplot(1, 2, 1)
    plt.plot(x_grid, ccp_empirical, 'b-', linewidth=2)
    plt.xlabel('状态')
    plt.ylabel('CCP')
    plt.title('条件选择概率')
    plt.grid(True, alpha=0.3)

    plt.subplot(1, 2, 2)
    plt.plot(x_grid[1:-1], second_differences, 'r-', linewidth=2)
    plt.xlabel('状态')
    plt.ylabel('二阶差分')
    plt.title('CCP的曲率')
    plt.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

# 检验性质3
def check_ccp_bounds(ccp):
    """CCP应该在[0,1]内"""

    if np.any(ccp < 0) or np.any(ccp > 1):
        print("⚠️  警告：CCP超出[0,1]范围！")
        print(f"最小值: {np.min(ccp):.6f}")
        print(f"最大值: {np.max(ccp):.6f}")
    else:
        print("✓ CCP在[0,1]范围内")

 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 identify_RC_from_CCP_threshold(x_grid, ccp_empirical, c_params, beta=0.95):
    """
    在维修阈值处（CCP=0.5），使用边界条件识别RC

    理由：在无差异点，两个选择的价值相等
    """

    # 找CCP=0.5的状态
    idx_threshold = np.argmin(np.abs(ccp_empirical - 0.5))
    x_star = x_grid[idx_threshold]

    # 在此点，必有：V_0(x*) = V_1(x*)
    # 即：-c(x*) + β V(x*+1) = -RC + β V(0)
    # 解出：RC = c(x*) - c(0) + β[V(0) - V(x*+1)]

    c1, c2 = c_params
    c_x_star = c1 * x_star + c2 * x_star ** 2
    c_0 = 0

    # 假设价值函数是线性的（简化）
    # V(x) ≈ V(0) - k * x
    # 则：V(x*+1) ≈ V(0) - k * (x*+1)

    # 或者用数值方法求解...

    return x_star, c_x_star - c_0

 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

def test_logit_assumption_via_ccp(ccp_empirical, x_grid):
    """
    检验Logit假设（极值分布冲击）是否合理

    如果冲击确实服从Gumbel分布，
    CCP应该满足某些性质
    """

    # 从CCP反推冲击分布
    # 如果CCP(x) = P(d=1|x)，则冲击差异的分布由CCP形状决定

    # 绘制CCP的对数赔率
    p = np.clip(ccp_empirical, 1e-10, 1 - 1e-10)
    log_odds = np.log(p / (1 - p))

    # 如果冲击是Gumbel分布，log_odds应该是状态的光滑函数

    plt.figure(figsize=(12, 5))

    plt.subplot(1, 2, 1)
    plt.plot(x_grid, ccp_empirical, 'b-', linewidth=2)
    plt.ylabel('CCP = P(d=1|x)')
    plt.xlabel('状态')
    plt.title('条件选择概率')
    plt.grid(True, alpha=0.3)

    plt.subplot(1, 2, 2)
    plt.plot(x_grid, log_odds, 'r-', linewidth=2)
    plt.ylabel('log(P/(1-P))')
    plt.xlabel('状态')
    plt.title('对数赔率')
    plt.grid(True, alpha=0.3)

    # 检验线性性（Logit特征）
    from scipy.stats import linregress
    slope, intercept, r_value, p_value, std_err = \
        linregress(x_grid, log_odds)

    print(f"对数赔率的线性回归 R²: {r_value ** 2:.4f}")
    if r_value ** 2 > 0.95:
        print("✓ log_odds很好地拟合线性关系 → Logit假设合理")
    else:
        print("⚠️ log_odds非线性 → Logit假设可能有问题")

    plt.show()

 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

def diagnose_sample_size_issue(data, x_grid, min_obs_per_bin=10):
    """
    问题：某些状态区间样本量太小
    后果：CCP估计方差大，不稳定
    """

    for x in x_grid:
        n_nearby = ((data['state'] >= x - 1) &
                    (data['state'] < x + 1)).sum()

        if n_nearby < min_obs_per_bin:
            print(f"⚠️  状态 x={x:.1f} 附近样本太少 ({n_nearby}个)")

    # 解决方案：扩大箱宽
    ccp_robust = estimate_ccp_with_adaptive_bandwidth(data, x_grid)

    return ccp_robust

def estimate_ccp_with_adaptive_bandwidth(data, x_grid, min_obs=20):
    """
    自适应带宽：保证每个点至少有min_obs个样本
    """

    ccp = np.zeros(len(x_grid))

    for i, x in enumerate(x_grid):

        # 初始带宽
        bandwidth = 1

        while True:
            in_window = (data['state'] >= x - bandwidth) & \
                       (data['state'] < x + bandwidth)

            if in_window.sum() >= min_obs:
                ccp[i] = data.loc[in_window, 'choice'].mean()
                break
            else:
                bandwidth *= 1.5  # 扩大带宽

    return ccp

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

def test_selection_bias_in_ccp(data):
    """
    问题：如果参与选择本身是内生的怎么办？

    例如：旧车可能因为主人穷而被卖掉
    → 观察到的里程与维修决策的关系可能被偏差
    """

    # 简单检验：观察CCP与状态的关系
    state_mean_repaired = data[data['choice']==1]['state'].mean()
    state_mean_not_repaired = data[data['choice']==0]['state'].mean()

    print(f"维修的平均状态: {state_mean_repaired:.1f}")
    print(f"不维修的平均状态: {state_mean_not_repaired:.1f}")
    print(f"差异: {abs(state_mean_repaired - state_mean_not_repaired):.1f}")

    if state_mean_repaired > state_mean_not_repaired:
        print("✓ 符合预期：状态差时更容易维修")
    else:
        print("⚠️  反常：状态好反而更容易维修？ → 可能有偏差")

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

def test_measurement_error_in_state(data):
    """
    问题：如果状态（如里程表）测量不准呢？
    后果：CCP被平滑化，S形变平缓
    """

    # 检验里程数据的一致性
    data_sorted = data.sort_values(['id', 'year'])

    # 里程应该单调增加
    mileage_decreases = (
        data_sorted.groupby('id')['state']
        .apply(lambda x: (x.diff() < 0).sum())
    )

    if mileage_decreases.sum() > 0:
        print("⚠️  警告：发现里程倒退")
        print(f"有 {mileage_decreases.sum()} 个单位发生过倒退")
        print("→ 可能存在测量误差")
    else:
        print("✓ 里程数据一致")

 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

def compare_ccp_vs_direct_ml():
    """
    对比两种方法的估计结果
    """

    # 生成数据
    data = generate_synthetic_data(n=500, T=100)
    x_grid = np.linspace(0, 450, 100)

    # 方法1：非参数CCP + 半参数估计
    print("=" * 60)
    print("方法1：CCP半参数估计")
    print("=" * 60)

    est_semi = SemiparametricCCPEstimation(data, x_grid)
    c1_semi = est_semi.estimate_parameters()
    time_semi = timer()

    print(f"估计结果: c1 = {c1_semi:.6f}")
    print(f"计算时间: {time_semi:.2f}秒")

    # 方法2：直接最大似然
    print("\n" + "=" * 60)
    print("方法2：结构极大似然估计")
    print("=" * 60)

    est_ml = StructuralMaximumLikelihood(data, x_grid)
    start = time.time()
    result_ml = est_ml.estimate()
    time_ml = time.time() - start

    print(f"估计结果: c1 = {result_ml.x[0]:.6f}")
    print(f"计算时间: {time_ml:.2f}秒")

    # 对比
    print("\n" + "=" * 60)
    print("对比总结")
    print("=" * 60)
    print(f"时间比：ML/CCP = {time_ml / time_semi:.2f}x")
    print(f"结果差异: |c1_ML - c1_CCP| = {abs(result_ml.x[0] - c1_semi):.6f}")

compare_ccp_vs_direct_ml()

1
2
3
4
5
6
7

CCP(维修概率)
    1.0 |              ╱
        |            ╱
    0.5 |          ╱  ← 维修阈值 x*
        |        ╱
    0.0 |_______╱_______
        0      x*      450  里程(万公里)