结构计量经济学：从直觉到代码的完整指南

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

def check_local_identification(moment_func, theta0, eps=1e-6):
    """检查矩条件的Jacobian是否满秩"""
    k = len(theta0)
    jacobian = np.zeros((len(moment_func(theta0)), k))

    for i in range(k):
        theta_plus = theta0.copy()
        theta_plus[i] += eps
        theta_minus = theta0.copy()
        theta_minus[i] -= eps
        jacobian[:, i] = (moment_func(theta_plus) - moment_func(theta_minus)) / (2 * eps)

    rank = np.linalg.matrix_rank(jacobian)
    print(f"Jacobian rank: {rank}, # of parameters: {k}")
    return rank >= k

 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

* 输入：效用向量 V（存储为向量）
* 计算Logit选择概率
* 输入：效用向量 V（存储为向量）
* 输出：概率向量

program define logit_choice_prob, rclass
    args V

    matrix V_matrix = `V'
    local J = rowsof(V_matrix)

    * 计算 exp(V) - max(V)，以数值稳定
    matrix V_centered = J(`J', 1, 0)
    local max_V = 0

    forval j = 1/`J' {
        local v_j = V_matrix[`j', 1]
        if `j' == 1 | `v_j' > `max_V' {
            local max_V = `v_j'
        }
    }

    * 计算指数
    matrix exp_V = J(`J', 1, 0)
    local sum_exp = 0

    forval j = 1/`J' {
        local v_j = V_matrix[`j', 1]
        local exp_v = exp(`v_j' - `max_V')
        matrix exp_V[`j', 1] = `exp_v'
        local sum_exp = `sum_exp' + `exp_v'
    }

    * 计算概率
    matrix prob = J(`J', 1, 0)
    forval j = 1/`J' {
        local exp_v_j = exp_V[`j', 1]
        local p_j = `exp_v_j' / (1 + `sum_exp')
        matrix prob[`j', 1] = `p_j'
    }

    return matrix prob = prob
end

 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

* 计算GMM渐近方差
* 输入：矩条件 g_matrix (N x L), 权重矩阵 W (L x L), N
* 输出：方差矩阵 var_matrix

program define gmm_variance, rclass
    args g_matrix W_matrix

    local N = rowsof(g_matrix)
    local L = colsof(g_matrix)

    * 计算均值梯度 G = mean(g)
    matrix G = J(`L', 1, 0)
    forval l = 1/`L' {
        local sum_g = 0
        forval i = 1/`N' {
            local sum_g = `sum_g' + `g_matrix'[`i', `l']
        }
        matrix G[`l', 1] = `sum_g' / `N'
    }

    * 计算 Sigma = (1/N) g'g
    matrix Sigma = `g_matrix'' * `g_matrix' / `N'

    * 计算方差：(G'WG)^{-1} G'W Σ WG (G'WG)^{-1} / N
    matrix GWG = G' * `W_matrix' * G
    matrix GWG_inv = inv(GWG)
    matrix temp1 = GWG_inv * G' * `W_matrix'
    matrix temp2 = temp1 * Sigma * `W_matrix' * G * GWG_inv
    matrix var_matrix = temp2 / `N'

    return matrix variance = var_matrix
end

 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 gmm_objective(theta, data, instruments, W, penalty=0.0):
    """
    GMM目标函数

    Args:
        theta: 参数向量
        data: 数据字典
        instruments: 工具变量矩阵 Z (N x L)
        W: 权重矩阵 (L x L)
        penalty: 正则化系数
    """
    # 计算结构残差
    xi = compute_structural_residual(theta, data)  # 你的模型特定函数

    # 样本矩
    g = (instruments.T @ xi) / len(xi)  # (L,)

    # GMM准则
    Q = g.T @ W @ g

    # 可选：加入正则化避免极端参数
    Q += penalty * np.sum(theta**2)

    return Q

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

def blp_contraction(delta_init, s_obs, theta, X, nu_draws, tol=1e-12, max_iter=1000):
    """
    BLP收缩映射：求解市场份额反演

    δ^{t+1} = δ^t + log(s_obs) - log(s_pred(δ^t))
    """
    delta = delta_init.copy()

    for iteration in range(max_iter):
        # 计算预测份额
        s_pred = predict_share(delta, theta, X, nu_draws)

        # 收缩映射更新
        delta_new = delta + np.log(s_obs) - np.log(s_pred + 1e-15)

        # 检查收敛
        if np.max(np.abs(delta_new - delta)) < tol:
            return delta_new, iteration

        delta = delta_new

    raise ValueError(f"Contraction did not converge after {max_iter} iterations")

 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

from scipy.optimize import minimize, differential_evolution

def robust_optimization(objective, theta_init, bounds):
    """
    多阶段优化策略
    """
    # Stage 1: 全局搜索找到好的起始点
    result_global = differential_evolution(
        objective,
        bounds,
        maxiter=100,
        polish=False  # 不在最后做局部优化
    )

    # Stage 2: 从全局最优附近做局部精细搜索
    result_local = minimize(
        objective,
        result_global.x,
        method='L-BFGS-B',
        bounds=bounds,
        options={'ftol': 1e-10}
    )

    # Stage 3: 多起始点验证
    best_result = result_local
    for _ in range(5):
        x0 = result_local.x + np.random.randn(len(theta_init)) * 0.1
        x0 = np.clip(x0, [b[0] for b in bounds], [b[1] for b in bounds])

        result = minimize(objective, x0, method='L-BFGS-B', bounds=bounds)
        if result.fun < best_result.fun:
            best_result = result

    return best_result

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

def get_smart_starting_values(data):
    """
    从简约式估计获取起始值
    """
    # 简单Logit作为起点
    from sklearn.linear_model import LogisticRegression

    lr = LogisticRegression()
    lr.fit(data['X'], data['y'])

    theta_init = {
        'beta': lr.coef_.flatten(),
        'alpha': -lr.coef_[0, -1],  # 假设最后一列是价格
        'sigma': np.ones(lr.coef_.shape[1]) * 0.5  # 随机系数初始化为小值
    }

    return theta_init

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

import jax
import jax.numpy as jnp

@jax.jit
def gmm_objective_jax(theta, data, instruments, W):
    """JAX版本，支持自动微分"""
    xi = compute_structural_residual_jax(theta, data)
    g = (instruments.T @ xi) / len(xi)
    return g.T @ W @ g

# 自动获取梯度函数
gmm_gradient = jax.grad(gmm_objective_jax)

# 同时计算值和梯度
gmm_value_and_grad = jax.value_and_grad(gmm_objective_jax)

  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
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276

import numpy as np
from scipy.optimize import minimize
from scipy.linalg import inv

class BLPEstimator:
    """
    Berry-Levinsohn-Pakes (1995) 需求估计器
    """

    def __init__(self, X, p, s_obs, Z, ns=500, seed=42):
        """
        Args:
            X: 产品特征 (J x K)
            p: 价格 (J,)
            s_obs: 观测市场份额 (J,)
            Z: 工具变量 (J x L)
            ns: 模拟抽样数
            seed: 随机种子
        """
        self.X = X
        self.p = p
        self.s_obs = s_obs
        self.Z = Z
        self.J, self.K = X.shape
        self.L = Z.shape[1]
        self.ns = ns

        # 固定随机抽样
        np.random.seed(seed)
        self.nu = np.random.randn(ns, self.K + 1)  # K个属性 + 1个价格

        # 预计算
        self.W = inv(Z.T @ Z)  # 初始权重矩阵
        self.PZ = Z @ inv(Z.T @ Z) @ Z.T  # 投影矩阵

    def predict_shares(self, delta, sigma):
        """
        预测市场份额：通过模拟积分

        s_j = (1/ns) Σ_r exp(δ_j + μ_jr) / (1 + Σ_k exp(δ_k + μ_kr))
        """
        # 构建个体偏离项 μ
        # mu[r, j] = X[j,:] @ (sigma_beta * nu[r, :K]) - sigma_alpha * nu[r, K] * p[j]
        sigma_beta = sigma[:-1]
        sigma_alpha = sigma[-1]

        mu = self.X @ (sigma_beta * self.nu[:, :-1].T).T  # (ns, J)
        mu -= sigma_alpha * np.outer(self.nu[:, -1], self.p)  # (ns, J)

        # 效用
        V = delta + mu  # (ns, J)

        # 选择概率（包含outside option）
        exp_V = np.exp(V - V.max(axis=1, keepdims=True))  # 数值稳定
        denom = 1 + exp_V.sum(axis=1, keepdims=True)
        prob = exp_V / denom  # (ns, J)

        # 平均份额
        shares = prob.mean(axis=0)

        return shares

    def contraction_mapping(self, sigma, tol=1e-12, max_iter=1000):
        """
        BLP收缩映射：从观测份额反推 δ
        """
        # 初始化
        delta = np.log(self.s_obs) - np.log(1 - self.s_obs.sum())  # Logit式初始值

        for it in range(max_iter):
            s_pred = self.predict_shares(delta, sigma)

            # 避免数值问题
            s_pred = np.maximum(s_pred, 1e-15)

            # 收缩更新
            delta_new = delta + np.log(self.s_obs) - np.log(s_pred)

            # 收敛检验
            if np.max(np.abs(delta_new - delta)) < tol:
                return delta_new, True

            delta = delta_new

        return delta, False

    def linear_parameters(self, delta):
        """
        给定 δ，用IV回归恢复线性参数 (β̄, ᾱ)

        δ = X β̄ - ᾱ p + ξ
        """
        # 将价格合并到设计矩阵
        Xp = np.column_stack([self.X, -self.p])  # (J, K+1)

        # 2SLS
        # θ = (X'PZ X)^{-1} X'PZ δ
        XpPZ = Xp.T @ self.PZ
        theta_linear = inv(XpPZ @ Xp) @ XpPZ @ delta

        beta_bar = theta_linear[:-1]
        alpha_bar = theta_linear[-1]

        return beta_bar, alpha_bar

    def gmm_objective(self, sigma):
        """
        GMM目标函数
        """
        # Step 1: 收缩映射
        delta, converged = self.contraction_mapping(sigma)

        if not converged:
            return 1e10  # 惩罚未收敛

        # Step 2: 恢复线性参数
        beta_bar, alpha_bar = self.linear_parameters(delta)

        # Step 3: 计算结构残差
        xi = delta - self.X @ beta_bar + alpha_bar * self.p

        # Step 4: 矩条件
        g = self.Z.T @ xi / self.J  # (L,)

        # Step 5: GMM准则
        Q = g.T @ self.W @ g

        return Q

    def estimate(self, sigma_init=None, method='L-BFGS-B'):
        """
        主估计函数
        """
        if sigma_init is None:
            sigma_init = np.ones(self.K + 1) * 0.5

        # 设置参数边界（标准差必须为正）
        bounds = [(0.01, 5.0)] * (self.K + 1)

        # 优化
        result = minimize(
            self.gmm_objective,
            sigma_init,
            method=method,
            bounds=bounds,
            options={'disp': True, 'maxiter': 1000}
        )

        # 提取最终估计
        sigma_hat = result.x
        delta_hat, _ = self.contraction_mapping(sigma_hat)
        beta_bar_hat, alpha_bar_hat = self.linear_parameters(delta_hat)

        # 存储结果
        self.results = {
            'sigma': sigma_hat,
            'beta_bar': beta_bar_hat,
            'alpha_bar': alpha_bar_hat,
            'delta': delta_hat,
            'gmm_obj': result.fun,
            'converged': result.success
        }

        return self.results

    def compute_elasticities(self):
        """
        计算自价格弹性和交叉价格弹性
        """
        if not hasattr(self, 'results'):
            raise ValueError("需要先运行 estimate()")

        sigma = self.results['sigma']
        delta = self.results['delta']
        alpha_bar = self.results['alpha_bar']
        sigma_alpha = sigma[-1]

        # 计算弹性矩阵
        elasticities = np.zeros((self.J, self.J))

        for r in range(self.ns):
            # 个体r的价格系数
            alpha_r = alpha_bar + sigma_alpha * self.nu[r, -1]

            # 个体选择概率
            mu_r = self.X @ (sigma[:-1] * self.nu[r, :-1])
            mu_r -= sigma_alpha * self.nu[r, -1] * self.p
            V_r = delta + mu_r
            exp_V = np.exp(V_r - V_r.max())
            prob_r = exp_V / (1 + exp_V.sum())

            # 弹性公式
            for j in range(self.J):
                for k in range(self.J):
                    if j == k:
                        elasticities[j, j] -= alpha_r * self.p[j] * (1 - prob_r[j]) * prob_r[j]
                    else:
                        elasticities[j, k] += alpha_r * self.p[k] * prob_r[k] * prob_r[j]

        elasticities /= self.ns

        # 归一化为弹性
        for j in range(self.J):
            elasticities[j, :] /= self.predict_shares(delta, sigma)[j]

        return elasticities


# ========== 使用示例 ==========
if __name__ == "__main__":
    # 生成模拟数据
    np.random.seed(123)
    J = 50  # 产品数
    K = 3   # 产品特征数

    # 真实参数
    TRUE_BETA_BAR = np.array([1.0, 0.5, 0.3])
    TRUE_ALPHA_BAR = 2.0
    TRUE_SIGMA = np.array([0.5, 0.3, 0.2, 1.0])  # 最后一个是价格系数的标准差

    # 产品特征
    X = np.random.randn(J, K)

    # 价格（包含内生性：与未观测质量相关）
    xi = np.random.randn(J) * 0.5  # 未观测质量
    cost = np.random.randn(J) * 0.3  # 成本冲击（工具变量来源）
    p = 2 + cost + 0.5 * xi + np.random.randn(J) * 0.1  # 价格内生

    # 工具变量：成本冲击 + 竞争者特征
    Z = np.column_stack([
        cost,
        X.sum(axis=0) - X,  # 其他产品的特征和
        np.random.randn(J)  # 额外的外生变量
    ])

    # 模拟真实市场份额
    def simulate_true_shares(X, p, xi, beta_bar, alpha_bar, sigma, ns=10000):
        np.random.seed(999)
        nu = np.random.randn(ns, K + 1)
        delta = X @ beta_bar - alpha_bar * p + xi

        shares = np.zeros(J)
        for r in range(ns):
            mu = X @ (sigma[:-1] * nu[r, :-1]) - sigma[-1] * nu[r, -1] * p
            V = delta + mu
            exp_V = np.exp(V - V.max())
            prob = exp_V / (1 + exp_V.sum())
            shares += prob
        return shares / ns

    s_obs = simulate_true_shares(X, p, xi, TRUE_BETA_BAR, TRUE_ALPHA_BAR, TRUE_SIGMA)

    # 估计
    estimator = BLPEstimator(X, p, s_obs, Z, ns=500)
    results = estimator.estimate()

    print("\n" + "="*50)
    print("估计结果:")
    print("="*50)
    print(f"σ (随机系数标准差): {results['sigma']}")
    print(f"真实值:              {TRUE_SIGMA}")
    print()
    print(f"β̄ (均值特征系数):   {results['beta_bar']}")
    print(f"真实值:              {TRUE_BETA_BAR}")
    print()
    print(f"ᾱ (均值价格系数):   {results['alpha_bar']:.4f}")
    print(f"真实值:              {TRUE_ALPHA_BAR}")
    print()
    print(f"GMM目标函数值:       {results['gmm_obj']:.6f}")
    print(f"收敛状态:            {results['converged']}")

    # 计算弹性
    elasticities = estimator.compute_elasticities()
    print()
    print("自价格弹性（前5个产品）:")
    print(np.diag(elasticities)[:5])

 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

def post_estimation_diagnostics(estimator):
    """
    估计后诊断
    """
    results = estimator.results

    # 1. 预测份额 vs 实际份额
    s_pred = estimator.predict_shares(results['delta'], results['sigma'])
    corr = np.corrcoef(estimator.s_obs, s_pred)[0, 1]
    print(f"份额相关性: {corr:.4f} (应该 > 0.95)")

    # 2. 残差分析
    xi = results['delta'] - estimator.X @ results['beta_bar'] + results['alpha_bar'] * estimator.p
    print(f"残差均值: {xi.mean():.6f} (应该接近0)")
    print(f"残差标准差: {xi.std():.4f}")

    # 3. 弹性合理性
    elasticities = estimator.compute_elasticities()
    own_elasticities = np.diag(elasticities)
    print(f"自价格弹性范围: [{own_elasticities.min():.2f}, {own_elasticities.max():.2f}]")
    print(f"弹性为负的比例: {(own_elasticities < 0).mean():.2%} (应该是100%)")

    # 4. 过度识别检验 (J-test)
    g = estimator.Z.T @ xi / len(xi)
    J_stat = len(xi) * g.T @ estimator.W @ g
    df = estimator.L - (estimator.K + 1 + len(results['sigma']))
    if df > 0:
        from scipy.stats import chi2
        p_value = 1 - chi2.cdf(J_stat, df)
        print(f"过度识别J统计量: {J_stat:.2f}, df={df}, p-value={p_value:.4f}")

    return {
        'share_correlation': corr,
        'xi_mean': xi.mean(),
        'xi_std': xi.std(),
        'elasticity_range': (own_elasticities.min(), own_elasticities.max())
    }

 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 robustness_checks(X, p, s_obs, Z):
    """
    稳健性检验套件
    """
    results_list = []

    # 1. 多起始点
    print("检验1: 起始值敏感性")
    for seed in range(10):
        np.random.seed(seed)
        sigma_init = np.abs(np.random.randn(X.shape[1] + 1)) * 0.5 + 0.1

        estimator = BLPEstimator(X, p, s_obs, Z, ns=500, seed=42)
        results = estimator.estimate(sigma_init=sigma_init)

        results_list.append({
            'seed': seed,
            'sigma': results['sigma'].copy(),
            'alpha': results['alpha_bar'],
            'obj': results['gmm_obj']
        })

    # 检查收敛一致性
    alphas = [r['alpha'] for r in results_list]
    print(f"  α估计值范围: [{min(alphas):.3f}, {max(alphas):.3f}]")
    print(f"  标准差: {np.std(alphas):.4f}")

    # 2. 模拟数敏感性
    print("\n检验2: 模拟数敏感性")
    for ns in [200, 500, 1000]:
        estimator = BLPEstimator(X, p, s_obs, Z, ns=ns, seed=42)
        results = estimator.estimate()
        print(f"  ns={ns}: α = {results['alpha_bar']:.4f}")

    return results_list