市场结构与进入退出：Bresnahan-Reiss 进入门槛模型

*==============================================================*
* Bresnahan-Reiss (1991) 市场进入模型 — Stata 实现
* 模拟数据: 200 个市场, n_firms=进入企业数(0..5),
*           pop=市场规模(对数人口), income=市场收入
*==============================================================*
clear all
set more off
set seed 20240912                    // 固定随机种子, 结果可复现
set obs 200

* --- 1. 模拟市场数据 ---
gen pop      = rnormal()              // 市场规模(标准化)
gen income   = rnormal()              // 另一市场特征
gen eps      = rnormal()              // 不可观测利润冲击
* 真实结构利润潜变量: pop 正向吸引进入, income 负向
gen latent   = 0.9*pop - 0.3*income + eps

* 按门槛把潜利润切成有序企业数 (0..5)
gen n_firms = .
replace n_firms = 0 if latent < -1.2
replace n_firms = 1 if latent >= -1.2 & latent < 0.0
replace n_firms = 2 if latent >=  0.0 & latent < 1.2
replace n_firms = 3 if latent >=  1.2 & latent < 2.4
replace n_firms = 4 if latent >=  2.4 & latent < 3.6
replace n_firms = 5 if latent >=  3.6

* --- 2. 有序 probit 估计 ---
oprobit n_firms pop income
estimates store br_probit

* 门槛/cutoff 解释: 第 n 个企业进入的利润门槛
estimates table br_probit, keep(/cut)

* --- 3. 有序 logit 稳健性 ---
ologit n_firms pop income
estimates store br_logit

* --- 4. 预测: 给定市场特征预测进入企业数 ---
predict p1 p2 p3 p4 p5, pr          // 各类别概率
* 期望企业数 = sum_n n * P(N=n)
gen E_n = 0*p0 + 1*p1 + 2*p2 + 3*p3 + 4*p4 + 5*p5
sum E_n

* --- 5. 反事实: 市场规模 pop 增加 1 个标准差 ---
replace pop = pop + 1.0
predict E_n_cf, pr
gen E_n_cf2 = 0*p0 + 1*p1 + 2*p2 + 3*p3 + 4*p4 + 5*p5
sum E_n_cf2

# -*- coding: utf-8 -*-
# 代码块2: Bresnahan-Reiss 有序 probit MLE (手写, 可运行)
import numpy as np
from scipy.stats import norm
from scipy.optimize import minimize

# ---------- 1. 模拟 M=200 个市场 ----------
rng = np.random.default_rng(7)
M = 200
S = rng.normal(0, 1, M)                 # 市场规模
x2 = rng.normal(0, 1, M)                # 另一市场特征
X = np.column_stack([np.ones(M), S, x2])
K = X.shape[1]
beta_true = np.array([0.0, 0.9, -0.3])  # 市场规模正向、收入负向
eps = rng.normal(0, 1, M)
latent = X @ beta_true + eps
tau_true = np.array([-1.2, 0.0, 1.2, 2.4, 3.6])
N_obs = np.digitize(latent, tau_true)   # 进入企业数 0..5
print("进入企业数分布:", np.bincount(N_obs, minlength=6))

# ---------- 2. 有序 probit 负对数似然 ----------
n_cat = 6
def neg_ll(params, X, y):
    b = params[:K]
    tau = np.concatenate([[-np.inf], params[K:K+n_cat-1], [np.inf]])
    xb = X @ b
    P = np.zeros(len(y))
    for n in range(n_cat):
        Pn = norm.cdf(tau[n+1]-xb) - norm.cdf(tau[n]-xb)
        P += (y == n)*Pn
    return -np.sum(np.log(np.clip(P, 1e-12, 1)))

tau0 = np.quantile(latent, np.linspace(0.1, 0.9, n_cat-1))
res = minimize(neg_ll, np.concatenate([[0,0.5,0.0], tau0]),
               args=(X, N_obs), method='L-BFGS-B')
b_hat, tau_hat = res.x[:K], res.x[K:]
print("beta 估计:", np.round(b_hat, 3), " (真实:", beta_true, ")")
print("门槛 tau 估计:", np.round(tau_hat, 3))

# ---------- 3. 门槛解释 ----------
for k in range(n_cat-1):
    print(f"  tau_{k+1}={tau_hat[k]:.3f}: 利润超此门槛才有第{k+2}家企业")

# -*- coding: utf-8 -*-
# 代码块3: 反事实预测 (接续代码块2 的 b_hat, tau_hat)
def predict_N(Xnew, b, tau):
    """给定市场特征 X, 返回期望进入企业数 E[N]"""
    xb = Xnew @ b
    edges = np.concatenate([[-np.inf], tau, [np.inf]])
    cum = norm.cdf(edges[:, None] - xb[None, :])  # (n_cat+1, M)
    Pn = np.diff(cum, axis=0)                       # (n_cat, M)
    return (Pn * np.arange(n_cat)[:, None]).sum(axis=0)

# (a) 基准: 用原始 X
N_base = predict_N(X, b_hat, tau_hat)

# (b) 反事实1: 市场规模 S 增加 1 个标准差
X_cf = X.copy(); X_cf[:, 1] = X_cf[:, 1] + 1.0
N_cf = predict_N(X_cf, b_hat, tau_hat)

# (c) 反事实2: 固定成本门槛整体上移 0.8 (进入更难)
N_cf2 = predict_N(X, b_hat, tau_hat + 0.8)

print(f"基准平均进入企业数        = {N_base.mean():.3f}")
print(f"市场规模+1σ 后平均企业数   = {N_cf.mean():.3f}  (+{N_cf.mean()-N_base.mean():.3f})")
print(f"固定成本门槛+0.8 后企业数  = {N_cf2.mean():.3f}  ({N_base.mean()-N_cf2.mean():.3f})")