住房与区位选择结构模型：Epple-Sieg 排序与 Bayer 离散选择

"""
代码块1：Epple-Sieg (1999) 排序模型 —— 连续社区排序与均衡房价求解
  - J 个社区按设施 a_j 升序；家庭收入对数正态；间接效用 V=(y-p)^alpha * a
  - 边界家庭 y_j* 在社区 j 与 j+1 间无差异，解析递推均衡房价
"""
import numpy as np
from scipy.stats import lognorm

np.random.seed(20240912)

# ---------- 1. 环境设定 ----------
J = 12                          # 社区数
alpha = 0.5                     # 偏好参数
a = 1.0 + 0.15 * np.arange(1, J + 1)       # 社区设施单调递增

# 收入分布：log y ~ N(mu, sigma^2)，中位收入约 50
mu, sigma = np.log(50.0), 0.40
Y_dist = lognorm(s=sigma, scale=np.exp(mu))
N = 40000
y_sim = Y_dist.rvs(size=N, random_state=7)

# ---------- 2. 由供给份额确定均衡边界收入 ----------
quantiles = np.arange(1, J) / J            # 1/J, ..., (J-1)/J
y_star = Y_dist.ppf(quantiles)             # J-1 个边界收入

# ---------- 3. 解析递推均衡房价 ----------
r = (a[:-1] / a[1:]) ** (1.0 / alpha)      # r_j = (a_j/a_{j+1})^(1/alpha)
p0 = 5.0                                   # 锚定最低社区房价
p_eq = np.zeros(J)
p_eq[0] = p0
for j in range(J - 1):
    # p_{j+1} = y_j*(1 - r_j) + r_j p_j
    p_eq[j + 1] = y_star[j] * (1.0 - r[j]) + r[j] * p_eq[j]

# ---------- 4. 模拟家庭选择，验证均衡 ----------
def choose_communities(p, a):
    resid = y_sim[:, None] - p[None, :]
    V = np.full((N, J), -np.inf)
    m = resid > 0
    V[m] = (resid[m] ** alpha) * np.tile(a, (N, 1))[m]
    return np.argmax(V, axis=1)
choice = choose_communities(p_eq, a)

# ---------- 5. 输出 ----------
print(f"{'社区':>4} {'设施a_j':>8} {'房价p_j':>9} {'平均收入':>9} {'份额':>7}")
for j in range(J):
    sel = choice == j
    inc = y_sim[sel].mean() if sel.any() else np.nan
    print(f"{j:>4} {a[j]:>8.2f} {p_eq[j]:>9.2f} {inc:>9.1f} {sel.mean():>7.3f}")

mean_inc = np.array([y_sim[choice==j].mean() if (choice==j).any() else np.nan
                     for j in range(J)])
print("平均收入随设施单调递增?", bool(np.all(np.diff(mean_inc[~np.isnan(mean_inc)])>=-1e-9)))

"""
代码块2：Bayer-McMillan-Rueben 离散选择 + IV 估计
  U_ij = X_j'beta - alpha p_j + gamma Xbar_j + xi_j + eps_ij
  p_j 与 xi_j 内生相关 -> OLS 偏误；用成本冲击 cost_shock 作 IV
"""
import numpy as np
from scipy.optimize import minimize

np.random.seed(20240912)
J, K = 100, 2
X = np.random.randn(J, K)
X[:,0] = (X[:,0]-X[:,0].mean())/X[:,0].std()   # 住房面积
X[:,1] = (X[:,1]-X[:,1].mean())/X[:,1].std()   # 到地铁距离

xi = 0.35 * np.random.randn(J)                  # 不可观测质量（与p相关）
cost_shock = np.random.randn(J)                 # 外生成本冲击（IV）
p = 0.5 + 0.4*X[:,0] - 0.3*X[:,1] + 0.6*xi + 1.5*cost_shock
p = np.abs(p) + 0.1
Xbar = np.random.randn(J)                       # 邻里构成
beta_true = np.array([1.0, -0.5]); alpha_true = -1.2; gamma_true = 0.3

# 平均效用 + 份额反演所需的观测份额（简化，模拟市场）
delta = X @ beta_true + alpha_true*p + gamma_true*Xbar + xi
R = 300
v = np.random.randn(R)
mu = 0.1 * v[:,None] * X[None,:,0]
u = delta[None,:] + mu
P = np.exp(u - u.max(1, keepdims=True)); P /= P.sum(1, keepdims=True)
S = P.mean(0) * np.exp(0.02*np.random.randn(J)); S /= S.sum()
S0 = max(1.0 - S.sum(), 1e-3)
y_dep = np.log(S) - np.log(S0)                  # Berry 反演被解释变量

# (a) Naive OLS
Z_ols = np.column_stack([X, p, Xbar])
b_ols = np.linalg.lstsq(Z_ols, y_dep, rcond=None)[0]
print("OLS  alpha =", round(b_ols[2],3), "(真实 -1.2)")

# (b) IV/2SLS：cost_shock 作 p 的工具
W_iv = np.column_stack([X, Xbar, cost_shock])
p_hat = W_iv @ np.linalg.lstsq(W_iv, p, rcond=None)[0]
Z_iv = np.column_stack([X, p_hat, Xbar])
b_iv = np.linalg.lstsq(Z_iv, y_dep, rcond=None)[0]
print("IV   alpha =", round(b_iv[2],3), "(真实 -1.2)")

"""
代码块3：反事实 —— 设施改善 -> 房价资本化 -> 福利(CV)
  承接代码1：a1=优质学区+20%；a2=郊区通勤改善
  CV：在新(p,a)下维持旧效用 V_base 所需收入 y'，CV = y - y'
"""
import numpy as np
from scipy.stats import lognorm
np.random.seed(20240912)
J, alpha = 12, 0.5
a0 = 1.0 + 0.15*np.arange(1, J+1)
Y_dist = lognorm(s=0.40, scale=np.exp(np.log(50.0)))
N = 40000; y_sim = Y_dist.rvs(size=N, random_state=7)

def solve_prices(a, p0=5.0):
    q = np.arange(1, J)/J; ys = Y_dist.ppf(q)
    r = (a[:-1]/a[1:])**(1/alpha)
    p = np.zeros(J); p[0] = p0
    for j in range(J-1):
        p[j+1] = ys[j]*(1-r[j]) + r[j]*p[j]
    return p

def choose(p, a):
    res = y_sim[:,None] - p[None,:]
    V = np.full((N,J), -np.inf); m = res>0
    V[m] = (res[m]**alpha)*np.tile(a,(N,1))[m]
    return np.argmax(V, axis=1)

# 基线
p_base = solve_prices(a0); c0 = choose(p_base, a0)
V_base = (y_sim - p_base[c0])**alpha * a0[c0]

# 反事实A：优质学区(末4社区)设施+20%
a1 = a0.copy(); a1[-4:] *= 1.20
pA = solve_prices(a1); cA = choose(pA, a1)
y_eq = (V_base / a1[cA])**(1/alpha) + pA[cA]
CV_A = y_sim - y_eq
print("学区房价上涨%:", np.round(np.mean((pA[-4:]-p_base[-4:])/p_base[-4:])*100,1))
print("平均 CV:", round(CV_A.mean(),2))

# 反事实B：地铁（郊区通勤改善）
a2 = a0.copy(); a2[:6] *= 1.15; a2[6:] *= 1.03
pB = solve_prices(a2); cB = choose(pB, a2)
y_eqB = (V_base/a2[cB])**(1/alpha) + pB[cB]
CV_B = y_sim - y_eqB
print("市区房价变化%:", np.round(np.mean((pB[6:]-p_base[6:])/p_base[6:])*100,1))
print("平均 CV:", round(CV_B.mean(),2))