动态离散选择模型估计方法的详细分析与对比

1
2
3
4

目标函数：max_θ Σᵢ Σₜ [dᵢₜ log P_θ(dᵢₜ|xᵢₜ) + (1-dᵢₜ)log(1-P_θ(dᵢₜ|xᵢₜ))]

其中P_θ(dᵢₜ|xᵢₜ)来自于求解贝尔曼方程：
V(x) = max_d [π(d,x) + β∫V(x')dF(x'|x,d)]

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

优化外层：
  初始化 θ
  for iter = 1, 2, ...:

    固定点内层：
      初始化V(x)
      repeat:
        V_new(x) = max_d [π(d,x) + β∫V_old(x')dF(x'|x,d)]
      until ||V_new - V_old|| < ε

    计算似然 L(θ)
    计算梯度 ∇L(θ)
    更新 θ = θ + step × ∇L(θ)

  return θ̂

1
2
3
4
5
6

第一阶段：非参数估计CCP
  P̂(x) = (非参数方法：核方法、局部多项式等)
  无需任何参数假设

第二阶段：参数反演
  使用固定的P̂(x)，反演结构参数θ

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

反演关系（Hotz-Miller映射）：
P̂(x) = ∫ exp(u(d,x) - v(x)) / Σ_d' exp(u(d',x) - v(x')) dε

其中：
  u(d,x) = 当期利润
  v(x) = 可以从CCP反演出的量
  ε = i.i.d极值噪声

反演过程：
  从P̂(x)和v(x)，反演θ的关键参数

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

目标函数：min_θ m̂ₙ(θ)' W m̂ₙ(θ)

其中：
  m̂ₙ(θ) = (1/n)Σᵢ g(Zᵢ;θ) （样本矩）
  W = 权重矩阵
  g(·;θ) = 矩条件（通常来自一阶条件或模型蕴含）

标准误：
  Var(θ̂) ≈ (Ĝ'WĜ)⁻¹ Ĝ'W Ŝ W Ĝ (Ĝ'WĜ)⁻¹
  其中Ĝ = 梯度，Ŝ = 矩的协方差

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

目标函数：min_θ [m̂_real - m_sim(θ)]' W [m̂_real - m_sim(θ)]

其中：
  m̂_real = 真实数据的矩
  m_sim(θ) = (1/N_sim T) Σₛ Σₜ m(Y_s,t;θ) （虚拟数据的矩）

虚拟数据生成：
  初始化x₀ ~ F₀
  for t = 1, ..., T:
    P_θ(x_t) = CCP（来自某处，如logit、原模型等）
    d_t ~ Bernoulli(P_θ(x_t))
    x_{t+1} ~ F_θ(x_{t+1}|x_t, d_t)

1
2
3
4
5
6
7
8

映射函数：β(θ) = arg min_β Q(β; Z_sim(θ))

反演目标：min_θ [β(θ) - β̂]' W [β(θ) - β̂]

其中：
  Z_sim(θ) = 用参数θ模拟的虚拟数据
  β(θ) = 虚拟数据上拟合辅助模型的参数
  β̂ = 真实数据上拟合辅助模型的参数

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

目标：p(θ|Y) ∝ p(Y|θ) × p(θ)

计算：MCMC（通常用Metropolis-Hastings）
  for iter = 1, ..., n:

    提议新参数：θ_new ~ q(θ_new|θ_old)

    接受概率：α = min(1, [p(Y|θ_new)p(θ_new)] / [p(Y|θ_old)p(θ_old)])

    计算p(Y|θ_new)需要：
      求解贝尔曼方程 → 得到CCP → 计算似然

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

第一步：非参数估计经验CCP
  P̂(x) = 从数据直接估计

第二步：贝叶斯反演
  p(θ|Y) ∝ L(Y|θ,P̂) × p(θ)

  其中似然使用固定的P̂（不随θ变化）

  MCMC：
    似然计算：L ∝ ∏ᵢₜ P̂(dᵢₜ|xᵢₜ)^{dᵢₜ} [1-P̂(dᵢₜ|xᵢₜ)]^{1-dᵢₜ}
    无需求解任何DP

1
2
3
4
5
6
7

外层优化：O(I_opt)              # I_opt ≈ 50-500次迭代
  每次迭代：
    内层固定点：O(I_bell × S²)  # I_bell ≈ 100-1000次迭代
                                # S = 离散状态数 ≈ 100-10000
    梯度计算：O(k × S²)          # k = 参数数

总复杂度：O(I_opt × I_bell × S²) ≈ O(100 × 500 × 10000²) = O(5×10¹²)

1
2
3
4
5
6
7
8
9

第一阶段：非参数CCP估计
  核密度估计：O(n × S)  # n = 观测数，S = 状态点数
  ≈ O(10000 × 100) = O(10⁶)

第二阶段：参数反演
  优化：O(I_opt)次迭代
  每次迭代：O(n) （计算似然）

总复杂度：O(n × S + I_opt × n) ≈ O(10⁶)

1
2
3
4
5
6
7
8
9

Method              Time per evaluation    Total optimization time
───────────────────────────────────────────────────────────
NFXP                2-5秒                  5-30分钟（50-300次评估）
CCP                 <0.1秒                <1分钟（100+次迭代）
GMM                 <0.1秒                <1分钟
SMM（N_sim=10000） 1-2秒                  10-50分钟
间接推断            0.5-1秒               5-20分钟
贝叶斯NFXP         2-5秒/迭代× 5000      10000-25000分钟（166-416小时！）
贝叶斯CCP          <0.1秒/迭代 × 5000    <50分钟

1
2

GMM ≈ CCP < 间接推断 ≈ SMM < NFXP << 贝叶斯NFXP
                                     但 贝叶斯NFXP >> 贝叶斯CCP

1
2
3
4
5
6
7
8

Var(θ̂_NFXP) = [I(θ₀)]^{-1}  （Fisher信息矩阵的逆）

其中I(θ) = -E[∂²logL/∂θ∂θ']

特性：
- 最小的渐近方差（Cramér-Rao下界）
- 当模型正确时，达到最有效
- 但只在大样本时有效

1
2
3
4
5
6
7
8
9

Var(θ̂_CCP) ≈ (G'G)^{-1} G' (Var(P̂) ⊗ I) G (G'G)^{-1}

其中：
- G = ∂θ/∂P̂ （参数对CCP的敏感性）
- Var(P̂) = 非参数CCP估计的方差

效率损失来自：
1. 第一步CCP估计的误差
2. 样本分割（某种意义上）

1
2
3
4
5
6
7
8
9

Var(θ̂_GMM) = (Ĝ'WĜ)^{-1} Ĝ'W Ŝ W Ĝ (Ĝ'WĜ)^{-1}

其中：
- Ĝ = 梯度
- Ŝ = 矩条件的协方差
- W = 权重矩阵

最优权重：W_opt = Ŝ^{-1}
此时与NFXP接近（如果矩条件来自一阶条件）

1
2
3
4
5

Var(θ̂_SMM) ≈ (Ĝ'WĜ)^{-1} Ĝ'W [Ŝ + var(虚拟数据)] W Ĝ (Ĝ'WĜ)^{-1}

额外方差来自：虚拟数据的随机性

解决：增加N_sim，或固定随机种子

1
2
3
4
5

Var(θ̂_II) ≈ [∂β/∂θ]^{-1} Var(β̂) [∂β/∂θ]^{-T}

两个来源：
1. 映射的陡峭度（∂β/∂θ小 → 方差大）
2. 辅助模型参数的精度（Var(β̂)）

1
2
3
4
5
6
7
8

Var(θ̂_Bayes) ≈ Var(θ|Y)  （后验方差）

取决于：
1. 似然的信息量
2. 先验的强弱
3. 样本量

通常 > MLE的方差（除非有非常强的好先验）

1
2
3
4
5
6
7
8
9

参数估计和标准误：

           NFXP        CCP        GMM        SMM      贝叶斯

θ_RC：    7.89(1.48)  7.92(2.15) 7.85(1.89) 7.88(1.95) 7.86(1.72)

θ_cost:   0.0001(.00002) .00009(.000035) .00008(.000028) .00007(.000031) .00008(.000029)

β:        0.9950(.0018) 0.9948(.0025) 0.9951(.0022) 0.9949(.0023) 0.9950(.0020)

1
2
3
4
5

J = n × m̂'(θ̂) W m̂(θ̂)

在零假设下（模型正确），J ~ χ²(# 超额条件)

p值 < 0.05：拒绝模型

1
2
3

对每个MCMC样本，模拟虚拟数据
计算虚拟数据和真实数据的某个统计量
看是否显著不同

  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

START: 有动态离散选择模型要估计

┌─────────────────────────────────────────────────────────┐
│ Step 1: 数据与模型的基本特征评估                         │
└─────────────────────────────────────────────────────────┘

问题A1: 样本量n有多大?
├─ n < 500（小样本）
│  └─> 路线A（倾向贝叶斯或CCP）
│
├─ 500 ≤ n < 5000（中等样本）
│  └─> 路线B（大多数方法可行）
│
└─ n ≥ 5000（大样本）
   └─> 路线C（所有方法都可考虑）

问题A2: 参数维度k有多大?
├─ k ≤ 5
│  └─> NFXP可考虑
│
├─ 5 < k ≤ 10
│  └─> CCP、GMM推荐
│
└─ k > 10
   └─> SMM、间接推断、GMM

问题A3: 状态维度d有多大?
├─ d = 1
│  └─> 任何方法都可以
│
├─ d = 2
│  └─> NFXP/CCP/GMM可行，SMM快
│
├─ d = 3
│  └─> CCP/GMM/SMM推荐
│
└─ d ≥ 4
   └─> 只有SMM、GMM、间接推断可行

问题A4: 有先验信息吗?
├─ 有强的先验（如从文献确定的折现因子）
│  └─> 贝叶斯方法有优势
│
├─ 有弱的先验（一般范围）
│  └─> 任何方法配合弱信息先验
│
└─ 无先验
   └─> 频率派方法更直接

┌─────────────────────────────────────────────────────────┐
│ Step 2: 计算资源与时间约束评估                           │
└─────────────────────────────────────────────────────────┘

问题B1: 有多少计算资源?
├─ 限制严格（笔记本电脑）
│  └─> GMM, CCP, 间接推断
│
├─ 中等（工作站）
│  └─> 任何方法都可以
│
└─ 充足（高性能集群）
   └─> 包括贝叶斯NFXP

问题B2: 对模型诊断/检验的需求?
├─ 需要自动的诊断检验
│  └─> GMM（J检验）
│
├─ 需要灵活的后验分析
│  └─> 贝叶斯方法
│
└─ 只要点估计和标准误
   └─> 任何方法都可以

┌─────────────────────────────────────────────────────────┐
│ Step 3: 模型specification的确信度                        │
└─────────────────────────────────────────────────────────┘

问题C1: 对模型specification有多确信?
├─ 非常确信（已发表的标准模型）
│  └─> NFXP（最有效）
│
├─ 相当确信（但有小调整）
│  └─> CCP或频率派方法
│
└─ 不太确信（探索性研究）
   └─> GMM（有诊断）或CCP（稳健性好）

┌─────────────────────────────────────────────────────────┐
│ Step 4: 实施经验与学习成本                               │
└─────────────────────────────────────────────────────────┘

问题D1: 有实施这些方法的经验吗?
├─ 丰富经验
│  └─> 可选任何方法
│
├─ 有基本经验
│  └─> CCP、GMM、SMM更安全
│
└─ 没有经验
   └─> GMM最好入门（最多教程和软件包）

END

1
2
3
4
5
6

if 有好的矩条件可设计:
    使用 GMM
else if 存在容易拟合的辅助模型:
    使用 间接推断
else:
    使用 SMM（虚拟数据的矩）

1
2
3
4
5
6
7

Step 1: 用GMM估计所有参数
Step 2: 检查J统计量，p值
        if p < 0.05:
            模型可能误设，需修改
        else:
            模型与数据一致
Step 3: 如果用多个矩条件集合，比较结果稳健性

1
2
3
4

Step 1: 获取后验样本 {θ^(s)}
Step 2: 对每个θ^(s)，求解修改后的DP，得到新的最优政策
Step 3: 根据新政策模拟虚拟数据
Step 4: 计算政策改变的影响（后验均值±标准差）

1
2
3
4
5
6
7
8
9

Step 1: 学习GMM基础
        推荐书：Cameron & Trivedi (2005) Ch.6

Step 2: 实施简单的GMM估计
        开源例子很多

Step 3: 学习动态模型中的GMM（如Arellano-Bond）

Step 4: 逐渐升级到CCP或SMM

1
2
3
4
5
6

w_速度 = 0.2   （计算时间重要程度）
w_精度 = 0.25  （渐近效率）
w_稳健 = 0.15  （对模型误设的稳健性）
w_诊断 = 0.15  （模型检验能力）
w_易用 = 0.15  （实施和学习的难度）
w_维数 = 0.1   （处理高维的能力）

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

# 步骤1：数据准备
# 输入：观测数据 {(x_it, d_it), i=1,...,n, t=1,...,T}
# 输出：P̂(x)

# 步骤2：离散化状态空间（如果连续）
# 将状态x分组到K个bin

# 步骤3：计算经验CCP
for each state bin k:
    n_k = 总观测数（state in bin k）
    d_k = 该bin中d=1的观测数
    P̂(x_k) = d_k / n_k

# 步骤4：光滑（可选）
# 用核方法或局部多项式光滑P̂

return P̂(·)

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

# 步骤1：设定目标函数
def objective(theta, P_hat):
    # 从固定的P̂和参数θ，反演结构参数
    # 使用Hotz-Miller映射：
    # u(d,x) - v(x) = log(P̂/(1-P̂)) + conditional expectation term

    # 最小化：Σ_it [u(d_it,x_it) - u(1-d_it,x_it) - predicted log odds]

    return loss

# 步骤2：优化
theta_hat = minimize(objective, theta_init, args=(P_hat,))

return theta_hat

 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

class CCPEstimator:

    def __init__(self, data, state_bins=20):
        self.data = data
        self.state_bins = state_bins

    def estimate_ccp_nonparametric(self):
        """第一阶段：非参数CCP估计"""

        # 离散化状态
        state_grid = np.linspace(self.data['x'].min(),
                                  self.data['x'].max(),
                                  self.state_bins)

        # 分配每个观测到最近的状态bin
        state_indices = np.digitize(self.data['x'], state_grid)

        # 计算经验CCP
        P_hat = {}
        for k in range(self.state_bins):
            mask = (state_indices == k)
            if mask.sum() > 0:
                P_hat[k] = self.data['d'][mask].mean()
            else:
                P_hat[k] = np.nan

        # 光滑（处理NaN）
        P_hat_smooth = self.smooth_ccp(P_hat)

        self.P_hat = P_hat_smooth
        self.state_grid = state_grid

        return P_hat_smooth

    def smooth_ccp(self, P_hat):
        """用局部多项式光滑"""
        from scipy.interpolate import UnivariateSpline

        indices = [k for k, v in P_hat.items() if not np.isnan(v)]
        x = np.array([self.state_grid[k] for k in indices])
        y = np.array([P_hat[k] for k in indices])

        spline = UnivariateSpline(x, y, s=0.01)

        return spline

    def likelihood(self, theta):
        """使用固定的P̂计算对数似然"""

        ll = 0
        for i in range(len(self.data)):
            x = self.data['x'][i]
            d = self.data['d'][i]

            # P̂(x)的值（用插值）
            p = self.P_hat(x)

            # 对数似然
            ll += d * np.log(p + 1e-10) + (1-d) * np.log(1-p + 1e-10)

        return ll

    def estimate_parameters(self, theta_init):
        """第二阶段：参数反演"""

        def objective(theta):
            return -self.likelihood(theta)

        result = minimize(objective, theta_init, method='BFGS')

        self.theta_hat = result.x

        return result.x

# 使用示例
estimator = CCPEstimator(data, state_bins=25)
P_hat = estimator.estimate_ccp_nonparametric()
theta_hat = estimator.estimate_parameters(theta_init=[15, 0.99])

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

矩条件的来源：

1. 一阶条件（Euler equation）：
   E[∂log L/∂θ] = 0

2. 模型预测的统计量：
   E[m(Z,θ)] = 0

3. 过度识别条件（如果矩数 > 参数数）

例子（Rust汽车维修）：
   m₁(θ) = E[d_t - P_θ(x_t)]  （CCP匹配）
   m₂(θ) = E[x_{t+1} - E_θ[x_{t+1}|x_t,d_t]]  （状态转移）
   m₃(θ) = E[d_t × x_t] - E_θ[d_t × x_t]  （交互项）

  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

class GMMEstimator:

    def __init__(self, data, moment_functions):
        self.data = data
        self.moment_functions = moment_functions  # list of functions

    def compute_moments(self, theta):
        """计算样本矩"""

        m = []
        for moment_fn in self.moment_functions:
            # 对每个观测计算矩条件
            moment_values = [moment_fn(self.data[i], theta)
                           for i in range(len(self.data))]

            # 取平均
            m.append(np.mean(moment_values))

        return np.array(m)

    def compute_moment_variance(self, theta):
        """计算矩条件的方差"""

        m_mat = []
        for i in range(len(self.data)):
            row = [moment_fn(self.data[i], theta)
                  for moment_fn in self.moment_functions]
            m_mat.append(row)

        m_mat = np.array(m_mat)

        # 样本协方差
        S = np.cov(m_mat.T)

        return S

    def gmm_objective(self, theta, W):
        """GMM目标函数"""

        m = self.compute_moments(theta)

        return m.T @ W @ m

    def estimate(self, theta_init, two_step=True):
        """两步GMM估计"""

        # 第一步：用单位矩阵作权重
        W1 = np.eye(len(self.moment_functions))

        result1 = minimize(lambda t: self.gmm_objective(t, W1),
                          theta_init, method='BFGS')
        theta1 = result1.x

        if two_step:
            # 计算最优权重（协方差矩阵的逆）
            S = self.compute_moment_variance(theta1)
            W2 = np.linalg.inv(S)

            # 第二步：用最优权重重新估计
            result2 = minimize(lambda t: self.gmm_objective(t, W2),
                             theta1, method='BFGS')

            theta_hat = result2.x
            J_stat = len(self.data) * self.gmm_objective(theta_hat, W2)

            return theta_hat, J_stat

        else:
            theta_hat = theta1
            J_stat = len(self.data) * result1.fun

            return theta_hat, J_stat

    def J_test(self, theta_hat):
        """过度识别检验"""

        from scipy.stats import chi2

        # ... (计算J统计量)

        p_value = 1 - chi2.cdf(J_stat, df=len(self.moment_functions) - len(theta_hat))

        return J_stat, p_value

# 使用示例
def moment_1(obs, theta):
    """CCP匹配矩"""
    x, d = obs
    P_theta_x = logistic(theta[0] + theta[1] * x)
    return d - P_theta_x

def moment_2(obs, theta):
    """状态转移矩"""
    x, d, x_next = obs
    expected_x = 10 + 0.5 * x + 5 * d  # 某个转移模式
    return x_next - expected_x

estimator = GMMEstimator(data, moment_functions=[moment_1, moment_2])
theta_hat, J = estimator.estimate(theta_init=[0, 0.1])
J_stat, p_value = estimator.J_test(theta_hat)

print(f"参数估计：{theta_hat}")
print(f"J统计量：{J_stat:.4f}")
print(f"p值：{p_value:.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
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

import pymc as pm
import arviz as az

class BayesianCCPEstimator:

    def __init__(self, data, priors):
        self.data = data
        self.priors = priors  # dict of prior specifications

        # 第一步：估计经验CCP（与CCP方法相同）
        self.P_hat = self.estimate_ccp_nonparametric()

    def estimate_ccp_nonparametric(self):
        """非参数CCP估计"""
        # （与CCP方法相同的实现）
        ...
        return P_hat_func

    def build_model(self):
        """构建PyMC模型"""

        with pm.Model() as model:

            # 先验定义
            theta_rc = pm.Gamma('theta_rc',
                               alpha=self.priors['rc_alpha'],
                               beta=self.priors['rc_beta'])

            theta_cost = pm.Normal('theta_cost',
                                  mu=self.priors['cost_mu'],
                                  sigma=self.priors['cost_sigma'])

            beta = pm.Beta('beta',
                          alpha=self.priors['beta_alpha'],
                          beta=self.priors['beta_beta'])

            # 似然（使用固定的P̂）
            d_pred = self.compute_likelihood(theta_rc, theta_cost, beta)

            d_obs = pm.Bernoulli('d_obs', p=d_pred,
                                observed=self.data['d'])

        return model

    def compute_likelihood(self, theta_rc, theta_cost, beta):
        """
        计算预测的选择概率
        使用固定的P̂
        """
        # 这里连接P̂和参数
        # （取决于具体的反演关系）
        ...
        return p_pred

    def sample(self, n_samples=2000, burn_in=1000):
        """MCMC抽样"""

        model = self.build_model()

        with model:
            trace = pm.sample(n_samples, tune=burn_in,
                            return_inferencedata=True,
                            cores=4)

        self.trace = trace

        return trace

    def posterior_summary(self):
        """后验统计"""

        return az.summary(self.trace)

    def plot_diagnostics(self):
        """诊断图"""

        az.plot_trace(self.trace)
        az.plot_forest(self.trace)

# 使用示例
priors = {
    'rc_alpha': 4,      # Gamma(4, 0.2) for RC
    'rc_beta': 0.2,
    'cost_mu': 0,
    'cost_sigma': 0.5,
    'beta_alpha': 10,   # Beta(10, 2) for discount factor
    'beta_beta': 2
}

estimator = BayesianCCPEstimator(data, priors)
trace = estimator.sample(n_samples=4000, burn_in=1000)
summary = estimator.posterior_summary()
estimator.plot_diagnostics()

1
2
3
4
5
6

用GMM：
  矩条件1：E[d_t - logistic(θ_rc + θ_β × x_t)] = 0
  矩条件2：E[Δx_{t+1} - E[Δx|x_t, d_t]] = 0

  结果：θ̂_rc ≈ 7.9, θ̂_β ≈ 0.995
  计算时间：< 5分钟

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

Step 1: 用CCP快速估计（初值）
        结果：θ̂ ≈ [7.9, 0.9945]
        时间：2分钟

Step 2: 用NFXP从CCP初值精化
        结果：θ̂ ≈ [7.88, 0.9950]
        时间：20分钟（从好初值收敛快）

Step 3: 用bootstrap计算标准误
        结果：SE ≈ [1.5, 0.002]
        时间：30分钟

总计：~50分钟，获得最优估计和信度

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

用贝叶斯CCP：
  先验：
    RC ~ Gamma(4, 0.2)     （均值20）
    β ~ Beta(10, 2)        （均值0.833）

  MCMC：4000迭代 + 1000 burn-in

  结果：
    RC: 7.85 ± 1.7  [95% CI: 4.8, 11.2]
    β:  0.9948 ± 0.002  [95% CI: 0.991, 0.999]

  计算时间：15分钟

  优势：
    - 完整的后验分布
    - 95%可信区间直接给出
    - 可用于后验预测和政策模拟

 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

# 步骤1：设计矩条件（5个条件匹配5个参数）

def moment_ccp_1(obs, theta):
    """CCP条件1：关于主设备的替换决策"""
    age1, age2, d_replace1 = obs
    P = logistic(theta[0] + theta[1] * age1 + theta[2] * age2)
    return d_replace1 - P

def moment_ccp_2(obs, theta):
    """CCP条件2：关于辅助设备的替换决策"""
    age1, age2, d_replace2 = obs
    P = logistic(theta[3] + theta[2] * age1 + theta[4] * age2)
    return d_replace2 - P

def moment_state_transition(obs, theta):
    """状态转移条件：年龄增长的预期"""
    age1, age2, d1, d2, age1_next, age2_next = obs
    expected_age1_next = age1 + 1 + (1 - d1) * depreciation_param
    return age1_next - expected_age1_next

# 步骤2：GMM估计
estimator = GMMEstimator(data, moment_functions=[
    moment_ccp_1, moment_ccp_2, moment_state_transition, ...
])

theta_hat, J = estimator.estimate(theta_init=[1, 0.05, -0.02, 0.5, 0.03])

# 步骤3：诊断
print(f"J统计量: {J:.4f}")
print(f"p值: {p_value:.4f}")

if p_value > 0.05:
    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

# 当矩条件难以设计时

from scipy.optimize import minimize

def simulate_data(theta, n_sim=1000):
    """模拟虚拟数据"""
    # ...
    return simulated_y

def compute_moments_real(data):
    """真实数据的矩"""
    return {
        'mean_d': data['d'].mean(),
        'var_d': data['d'].var(),
        'corr_age_decision': correlation(data['age'], data['d']),
        # ... 更多矩
    }

def smm_objective(theta, data):
    """SMM目标函数"""

    # 虚拟数据的矩
    sim_y = simulate_data(theta)
    m_sim = compute_moments(sim_y)

    # 真实数据的矩
    m_real = compute_moments_real(data)

    # 距离
    diff = np.array(list(m_sim.values())) - np.array(list(m_real.values()))

    return diff @ diff

# 估计
theta_hat = minimize(smm_objective, theta_init, args=(data,)).x

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

问1：样本量 > 5000？
   是 → Q2
   否 → Q4

问2：参数数 ≤ 5？
   是 → CCP 或 NFXP（有充足计算资源时）
   否 → Q3

问3：能设计好的矩条件？
   是 → GMM
   否 → SMM 或 间接推断

问4（小样本）：有先验吗？
   是 → 贝叶斯CCP
   否 → CCP + 稳健标准误

问5（所有情况）：需要模型诊断？
   是 → GMM（J检验）
   否 → 任何方法都可以

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

# ❌ 错误做法
while iter < max_iter:
    V_new = max_d [π(d) + β*E_V]
    if iter % 10 == 0 and iter > 10:
        break  # 不足以收敛！

# ✓ 正确做法
while True:
    V_new = max_d [π(d) + β*E_V]
    if ||V_new - V_old|| < 1e-8:  # 严格的收敛标准
        break
    V_old = V_new

1
2
3
4
5
6

# ❌ 错误：直接求导
dL/dθ = numerical_gradient(likelihood, theta, eps=1e-6)

# ✓ 正确：自适应eps
eps = 1e-6 * (1 + |theta|)
dL/dθ = (likelihood(theta+eps) - likelihood(theta-eps)) / (2*eps)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

# ❌ 错误做法
state_bins = 50  # 太多bins，每个bin样本少
P_hat = observed_frequency_in_bins

# 结果：很多bin的P̂ = 0 或 1，导致反演问题

# ✓ 正确做法
# 选择bin数使得每个bin平均有10-20个观测
optimal_bins = int(np.sqrt(n))  # 或用交叉验证
P_hat = smooth_kernal_or_localpolynomial(state, d)

1
2
3
4
5
6
7
8
9

# ❌ 过度光滑：P̂接近常数，无信息
P_hat = heavy_smoothing_kernel(state, d)

# ❌ 过度拟合：P̂有噪声，反演不稳定
P_hat = nearest_neighbor(state, d)

# ✓ 正确：交叉验证选择带宽
bandwidth = cross_validation_bandwidth(state, d)
P_hat = kernel_density(state, d, bandwidth)

1
2
3
4
5
6
7

# ❌ 选择的矩对参数不敏感
moment_1 = E[d]  # 这是标量，无法识别所有参数

# ✓ 正确：选择有区分性的矩
moment_1 = E[d]  # 识别平均替换率
moment_2 = E[d × x]  # 识别成本对年龄的敏感性
moment_3 = E[d × x²]  # 加强识别

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

# ❌ 错误：总是用单位矩阵
W = I

# ✓ 正确：
# 第一步：W = I，粗估θ
theta_1 = minimize(m'*I*m, ...)

# 第二步：计算最优权重
S = cov(moments)
W_opt = inv(S)

# 第三步：W = W_opt，精化估计
theta_2 = minimize(m'*W_opt*m, ...)

1
2
3
4
5
6
7

# ❌ 虚拟样本太小，矩估计有噪声
N_sim = 500, T = 10
# 导致m_sim(θ)非常不平滑，优化困难

# ✓ 虚拟样本足够大
N_sim = 10000, T = 50
# 或者固定随机种子使m_sim(θ)确定

1
2
3
4
5
6
7
8
9

# ❌ 错误的决策过程
d_t ~ Bernoulli(P_true)  # 用真实分布
x_{t+1} ~ F_true | d_t   # 用真实转移

# 这不对！应该用参数化的

# ✓ 正确的虚拟数据生成
d_t ~ Bernoulli(P_θ(x_t))      # 用模型的CCP
x_{t+1} ~ F_θ | x_t, d_t       # 用参数化的转移

1
2
3
4
5

# ❌ 辅助模型太简单，无法区分参数
aux_model: d ~ 常数  # 无法识别任何参数

# ✓ 辅助模型有充足复杂度
aux_model: logit(d) ~ α + β*x + γ*x²

1
2
3
4
5
6
7
8
9

# ❌ 直接反演，假设映射可逆
beta_hat_grid = compute_beta_on_grid(theta_grid)
theta_hat = invert(beta_hat)  # 可能无唯一解！

# ✓ 先检验可逆性
# 绘制β(θ)的曲线，看是否单调或可逆
plot_mapping(theta_grid, beta_hat_grid)
# 检查Jacobian的秩
rank_J = rank(jacobian_beta_wrt_theta)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

# ❌ 运行MCMC后直接用样本
for i in range(n_samples):
    theta_i = theta_samples[i]
    # ...

# ✓ 正确的诊断步骤
az.plot_trace(trace)  # 检查trace图
az.summary(trace)     # 检查Rhat和ESS
# 丢弃burn-in样本
post_samples = trace.posterior[burn_in:]

1
2
3
4
5

# ❌ 无信息先验导致问题
beta ~ Uniform(0, 1)  # 在高维中这是很弱的信息！

# ✓ 弱信息先验
beta ~ Beta(10, 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
27
28

□ 数据质量检查
  ☐ 是否有异常值或数据输入错误？
  ☐ 样本量是否足够？
  ☐ 数据是否包含足够的变异？

□ 模型specification检查
  ☐ 模型假设是否清晰？
  ☐ 参数数与样本量的比例合理？
  ☐ 是否有识别问题？

□ 初值检查
  ☐ 初值是否在合理范围？
  ☐ 从多个初值开始是否得到相同结果？

□ 收敛检查
  ☐ 优化是否收敛（梯度接近零）？
  ☐ Hessian是否为负定（最大化）？
  ☐ 目标函数值是否合理？

□ 结果合理性检查
  ☐ 估计的参数在理论范围内？
  ☐ 符号和大小符合经济直觉？
  ☐ 与文献中的结果可比较？

□ 稳健性检查
  ☐ 结果对样本选择敏感吗？
  ☐ 删除某些观测后结果是否稳定？
  ☐ 对模型specification的微小改变敏感吗？

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

START

问题1：你对方法学有深入理解吗？
├─ 否 → 用GMM（最容易学习）
└─ 是 → 问题2

问题2：模型多复杂？
├─ 简单（k≤5, d≤1） → 问题3
├─ 中等（k≤10, d≤2） → 问题3
└─ 复杂（k>10 或 d>2） → 问题4

问题3：你有充足的计算资源吗？
├─ 是 → NFXP或CCP都可以
└─ 否 → CCP（最快）

问题4（复杂模型）：能设计好矩条件吗？
├─ 能 → GMM
├─ 不能，但能想到辅助模型 → 间接推断
└─ 都不能 → SMM

END

1
2
3
4
5
6

NFXP：  最优但最慢 (Optimal but slowest)
CCP：   平衡的首选 (Balanced choice)
GMM：   诊断最好 (Best diagnostics)
SMM：   万能但贵 (Universal but expensive)
间接：  灵活但需好模型 (Flexible if good auxiliary model)
贝叶斯：最完整但最慢 (Most complete but slowest)

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

第一步：选择主要方法
  推荐：CCP + 健壮性检查
  或：  NFXP（如果计算可行）
  或：  贝叶斯CCP（如果要完整分布）

第二步：选择替选/对比方法
  推荐：至少再用一个方法验证
  例如：CCP主估计 + GMM对比 + 贝叶斯稳健性

第三步：诊断与报告
  NFXP/CCP：报告标准误、初值敏感性
  GMM：报告J检验、过度识别检验
  贝叶斯：报告先验、MCMC诊断、敏感性分析