实验 72 寡头迭代：Cournot 与 Bertrand

import numpy as np

# 景区门口的烤肠摊（和第 16 讲一样）：游客的需求 P = 11 - 0.01 Q，每根成本 2 元，没有固定成本
A, B, MC = 11.0, 0.01, 2.0

def price(total):
    return max(A - B * total, 0.0)

def profit(q, others, c=MC):
    """我卖 q 根、别的摊一共卖 others 根时，我的利润"""
    return (price(q + others) - c) * q

def best_reply(others, c=MC):
    """令边际收益 = 边际成本，得到最佳卖量的公式；卖量不能是负的"""
    return max((A - c - B * others) / (2 * B), 0.0)

# 核对 1：公式，和"把 0 到 900 的卖量全试一遍"，答案一样吗？
grid = np.arange(0, 900.01, 0.5)
for others in (0, 150, 300, 450):
    brute = grid[np.argmax([profit(q, others) for q in grid])]
    print(f"别的摊卖 {others:>3} 根：公式 {best_reply(others):6.1f}，全部试一遍 {brute:6.1f}")

# 核对 2：N 家一样的摊，均衡里每家卖 (A - c) / (B (N + 1))
print("\n  N   每家卖    价格    每家利润")
for n in (1, 2, 3, 5, 10):
    q = (A - MC) / (B * (n + 1))
    print(f"{n:>3} {q:>8.1f} {price(n * q):>7.2f} {profit(q, (n - 1) * q):>10.2f}")

# 核对 3：两家轮流调整（第 16 讲那一行输出）。每一轮：甲按乙调整，然后乙按甲调整
q1, q2 = 100.0, 500.0
path = [(q1, q2)]
for _ in range(5):
    q1 = best_reply(q2); q2 = best_reply(q1); path.append((q1, q2))
print("\n轮流调整：" + " -> ".join(f"({a:.0f}, {b:.0f})" for a, b in path))

别的摊卖   0 根：公式  450.0，全部试一遍  450.0
别的摊卖 150 根：公式  375.0，全部试一遍  375.0
别的摊卖 300 根：公式  300.0，全部试一遍  300.0
别的摊卖 450 根：公式  225.0，全部试一遍  225.0

  N   每家卖    价格    每家利润
  1    450.0    6.50    2025.00
  2    300.0    5.00     900.00
  3    225.0    4.25     506.25
  5    150.0    3.50     225.00
 10     81.8    2.82      66.94

轮流调整：(100, 500) -> (200, 350) -> (275, 312) -> (294, 303) -> (298, 301) -> (300, 300)

def iterate(q0, rounds, mode="simul", lam=1.0, costs=None):
    """N 家摊反复调整卖量。
    mode="simul"：大家同时看上一轮对方的卖量，一起调整（lam < 1 表示只调整一部分，叫"惯性"）
    mode="alt"  ：轮流来，每家调整时，看到的是别人最新的卖量
    返回一个数组：第 r 行是第 r 轮之后每家的卖量"""
    q = np.array(q0, dtype=float); n = len(q)
    costs = np.full(n, MC) if costs is None else np.asarray(costs, dtype=float)
    hist = [q.copy()]
    for _ in range(rounds):
        if mode == "simul":
            new = np.array([best_reply(q.sum() - q[i], costs[i]) for i in range(n)])
            q = (1 - lam) * q + lam * new
        else:
            for i in range(n):
                q[i] = (1 - lam) * q[i] + lam * best_reply(q.sum() - q[i], costs[i])
        hist.append(q.copy())
    return np.array(hist)

# 先看两家：同时调整，和轮流调整，都能走到 (300, 300) 吗？
for mode in ("simul", "alt"):
    h = iterate([100, 500], 12, mode)
    print(f"{mode:>5}：" + " -> ".join(f"({a:.0f},{b:.0f})" for a, b in h[[0, 1, 2, 3, 6, 12]]))

# 再看 N 家，从均衡附近出发（每家比均衡多卖 20 根），同时调整 8 轮
print("\n同时调整，每家的卖量（均衡值 = 450 / (N + 1)）：")
for n in (2, 3, 4, 5):
    star = (A - MC) / (B * (n + 1))
    start = np.full(n, star); start[0] += 20
    h = iterate(start, 8, "simul")
    print(f"N = {n}（均衡 {star:6.1f}）：甲的卖量 " + " ".join(f"{x:6.1f}" for x in h[:9, 0]))

# N = 5 时，市场上整体发生了什么？看总卖量和价格
n = 5; star = (A - MC) / (B * (n + 1))
start = np.full(n, star); start[0] += 20
h = iterate(start, 12, "simul")
print("\nN = 5：总卖量和价格（均衡：总卖量 750，价格 3.50）")
for r in range(0, 13, 2):
    print(f"  第 {r:>2} 轮：总卖量 {h[r].sum():7.1f}，价格 {price(h[r].sum()):.2f}，各家卖量 {np.round(h[r], 0)}")

simul：(100,500) -> (200,400) -> (250,350) -> (275,325) -> (297,303) -> (300,300)
  alt：(100,500) -> (200,350) -> (275,312) -> (294,303) -> (300,300) -> (300,300)

同时调整，每家的卖量（均衡值 = 450 / (N + 1)）：
N = 2（均衡  300.0）：甲的卖量  320.0  300.0  305.0  300.0  301.2  300.0  300.3  300.0  300.1
N = 3（均衡  225.0）：甲的卖量  245.0  225.0  235.0  220.0  232.5  218.7  231.9  218.4  231.7
N = 4（均衡  180.0）：甲的卖量  200.0  180.0  195.0  165.0  206.3  142.5  237.2   94.7  308.2
N = 5（均衡  150.0）：甲的卖量  170.0  150.0  170.0  120.0  215.0   22.5  406.3    0.0  450.0

N = 5：总卖量和价格（均衡：总卖量 750，价格 3.50）
  第  0 轮：总卖量   770.0，价格 3.30，各家卖量 [170. 150. 150. 150. 150.]
  第  2 轮：总卖量   830.0，价格 2.70，各家卖量 [170. 165. 165. 165. 165.]
  第  4 轮：总卖量  1070.0，价格 0.30，各家卖量 [215. 214. 214. 214. 214.]
  第  6 轮：总卖量  2030.0，价格 0.00，各家卖量 [406. 406. 406. 406. 406.]
  第  8 轮：总卖量  2250.0，价格 0.00，各家卖量 [450. 450. 450. 450. 450.]
  第 10 轮：总卖量  2250.0，价格 0.00，各家卖量 [450. 450. 450. 450. 450.]
  第 12 轮：总卖量  2250.0，价格 0.00，各家卖量 [450. 450. 450. 450. 450.]

# 为什么：均衡附近，"同时调整"每一轮把偏差乘以一个数。这个数，就是调整的"雅可比矩阵"的特征值
def one_round_matrix(n, mode):
    """在均衡附近，一轮调整把「偏差」变成多少：一个 n × n 的矩阵（对均衡处的微小偏差做一轮，逐列读出来）"""
    star = np.full(n, (A - MC) / (B * (n + 1))); G = np.zeros((n, n))
    for j in range(n):
        start = star.copy(); start[j] += 1e-3
        G[:, j] = (iterate(start, 1, mode)[1] - star) / 1e-3
    return G

print("一轮调整之后，偏差被放大（>1）还是缩小（<1）？最大特征值的绝对值：")
print("  N     同时调整   轮流调整")
for n in (2, 3, 4, 5, 10, 20):
    rho = lambda mode: np.abs(np.linalg.eigvals(one_round_matrix(n, mode))).max()
    print(f"{n:>3} {rho('simul'):>11.3f} {rho('alt'):>10.3f}   （同时调整：(N - 1) / 2 = {(n - 1) / 2:.1f}）")

# 药方一：加"惯性"，每一轮只调整一部分 lam。用二分法找：最大的、还能收敛的 lam
def converges(n, lam, rounds=4000):
    star = (A - MC) / (B * (n + 1))
    start = np.full(n, star); start[0] += 20
    h = iterate(start, rounds, "simul", lam)
    return abs(h[-1] - star).max() < 1e-6

print("\n药方一：惯性。最大的可收敛调整比例 lam（对照公式 4 / (N + 1)）：")
for n in (4, 5, 10, 20):
    lo, hi = 0.0, 1.0
    for _ in range(25):
        mid = (lo + hi) / 2
        if converges(n, mid): lo = mid
        else: hi = mid
    print(f"N = {n:>2}：{lo:.4f}（公式 {4 / (n + 1):.4f}）")

# 药方二：轮流调整。N = 10，每家从 100 根出发，看离均衡（81.8）最远的那家差多少
n = 10; star = (A - MC) / (B * (n + 1))
h = iterate(np.full(n, 100.0), 60, "alt")
print("\n药方二：轮流调整，N = 10，离均衡最远的那家差多少：")
print("  第 0、1、2、5、10、20、40 轮：" + "、".join(f"{abs(h[r] - star).max():.2f}" for r in (0, 1, 2, 5, 10, 20, 40)))
print("轮流调整，要几轮才能让每家离均衡都不超过 0.5 根：")
for n in (2, 5, 10, 20, 40):
    star = (A - MC) / (B * (n + 1)); h = iterate(np.full(n, 100.0), 4000, "alt")
    r = int(np.argmax(abs(h - star).max(axis=1) < 0.5))
    print(f"  N = {n:>2}：{r} 轮")

一轮调整之后，偏差被放大（>1）还是缩小（<1）？最大特征值的绝对值：
  N     同时调整   轮流调整
  2       0.500      0.250   （同时调整：(N - 1) / 2 = 0.5）
  3       1.000      0.354   （同时调整：(N - 1) / 2 = 1.0）
  4       1.500      0.458   （同时调整：(N - 1) / 2 = 1.5）
  5       2.000      0.544   （同时调整：(N - 1) / 2 = 2.0）
 10       4.500      0.784   （同时调整：(N - 1) / 2 = 4.5）
 20       9.500      0.924   （同时调整：(N - 1) / 2 = 9.5）

药方一：惯性。最大的可收敛调整比例 lam（对照公式 4 / (N + 1)）：
N =  4：0.7985（公式 0.8000）
N =  5：0.6654（公式 0.6667）
N = 10：0.3630（公式 0.3636）
N = 20：0.1901（公式 0.1905）

药方二：轮流调整，N = 10，离均衡最远的那家差多少：
  第 0、1、2、5、10、20、40 轮：18.18、81.82、31.92、10.02、3.19、0.28、0.00
轮流调整，要几轮才能让每家离均衡都不超过 0.5 根：
  N =  2：5 轮
  N =  5：9 轮
  N = 10：18 轮
  N = 20：62 轮
  N = 40：198 轮

def cournot(costs):
    """成本各不相同的 Cournot 均衡。公式：q_i = (A - (n+1) c_i + 所有人的成本之和) / (B (n+1))；
    如果某家算出来是负的，它就不卖了，把它去掉再算一遍。"""
    costs = np.asarray(costs, dtype=float)
    active = np.ones(len(costs), dtype=bool)
    while True:
        n = active.sum()
        q = np.where(active, (A - (n + 1) * costs + costs[active].sum()) / (B * (n + 1)), 0.0)
        if (q >= 0).all(): break
        active[np.argmax(np.where(q < 0, costs, -1))] = False     # 去掉成本最高的、卖量为负的那家
    p = price(q.sum())
    return q, p, (p - costs) * q

# 三家摊：甲成本 1.5（自家有便宜的货源），乙 2.0，丙 2.5
costs = [1.5, 2.0, 2.5]
q, p, pi = cournot(costs)
h = iterate(np.full(3, 100.0), 300, "alt", costs=costs)
print("三家摊，成本 1.5 / 2.0 / 2.5")
print(f"  公式：卖量 {np.round(q, 1)}，价格 {p:.3f}，利润 {np.round(pi, 1)}")
print(f"  轮流调整 300 轮：卖量 {np.round(h[-1], 1)}（和公式一致）")
share = q / q.sum()
hhi = (share ** 2).sum()
elasticity = (1 / B) * p / q.sum()
lerner = (share * (p - np.array(costs)) / p).sum()
print(f"  份额 {np.round(share, 3)}，HHI = {hhi:.4f}，市场弹性 = {elasticity:.3f}")
print(f"  按份额加权的加价率 = {lerner:.4f}；HHI ÷ 弹性 = {hhi / elasticity:.4f}（相等）")

# 成本差距太大时，贵的那家就退出了
print("\n丙的成本一路涨：")
for c3 in (2.5, 4.0, 6.0, 8.0):
    q, p, pi = cournot([1.5, 2.0, c3])
    print(f"  丙成本 {c3:.1f}：卖量 {np.round(q, 1)}，价格 {p:.3f}" + ("（丙退出）" if q[2] == 0 else ""))

# 合并：十家一样的摊，其中 k 家合并成一家。合并后只是"少了 k - 1 个对手"，成本不变
print("\n十家一样的摊（每家成本 2）。k 家合并后，合并的这家摊，利润够不够补回合并前各自的利润之和？")
n = 10
pi_each = cournot([MC] * n)[2][0]
print(f"  合并前每家利润 {pi_each:.2f}，价格 {cournot([MC] * n)[1]:.3f}")
for k in range(2, 11):
    q, p, pi = cournot([MC] * (n - k + 1))
    gain = pi[0] - k * pi_each
    print(f"  k = {k:>2} 家合并：合并后价格 {p:.3f}，合并的这家利润 {pi[0]:7.2f}，合并前这 k 家的利润之和 {k * pi_each:7.2f}，差 {gain:+7.2f}")

# 如果合并能让成本下降（效率），降多少才够？二分法找：k = 2 的合并，合并后成本 c，使得合并这家的利润 = 2 × 合并前
def merged_profit(k, c_new):
    q, p, pi = cournot([c_new] + [MC] * (n - k)); return pi[0], p
print("\n两家合并，合并后成本降到多少，这次合并才有赚头？")
lo, hi = 0.0, MC
for _ in range(40):
    mid = (lo + hi) / 2
    if merged_profit(2, mid)[0] >= 2 * pi_each: lo = mid
    else: hi = mid
c_need = hi; pi_m, p_m = merged_profit(2, c_need)
print(f"  成本要从 2.00 降到 {c_need:.3f}（降 {100 * (MC - c_need) / MC:.1f}%）；这时价格 {p_m:.3f}，合并前是 {cournot([MC] * n)[1]:.3f}")

# 再问一个更严的问题：成本要降多少，才能让这次合并之后，价格不比合并前高？
lo, hi = 0.0, MC
for _ in range(40):
    mid = (lo + hi) / 2
    if merged_profit(2, mid)[1] <= cournot([MC] * n)[1]: lo = mid
    else: hi = mid
print(f"想让价格不涨，成本要降到 {lo:.3f}（降 {100 * (MC - lo) / MC:.1f}%），比「合并有赚头」要求的 {100 * (MC - c_need) / MC:.1f}% 多得多")

三家摊，成本 1.5 / 2.0 / 2.5
  公式：卖量 [275. 225. 175.]，价格 4.250，利润 [756.2 506.2 306.2]
  轮流调整 300 轮：卖量 [275. 225. 175.]（和公式一致）
  份额 [0.407 0.333 0.259]，HHI = 0.3443，市场弹性 = 0.630
  按份额加权的加价率 = 0.5468；HHI ÷ 弹性 = 0.5468（相等）

丙的成本一路涨：
  丙成本 2.5：卖量 [275. 225. 175.]，价格 4.250
  丙成本 4.0：卖量 [312.5 262.5  62.5]，价格 4.625
  丙成本 6.0：卖量 [333.3 283.3   0. ]，价格 4.833（丙退出）
  丙成本 8.0：卖量 [333.3 283.3   0. ]，价格 4.833（丙退出）

十家一样的摊（每家成本 2）。k 家合并后，合并的这家摊，利润够不够补回合并前各自的利润之和？
  合并前每家利润 66.94，价格 2.818
  k =  2 家合并：合并后价格 2.900，合并的这家利润   81.00，合并前这 k 家的利润之和  133.88，差  -52.88
  k =  3 家合并：合并后价格 3.000，合并的这家利润  100.00，合并前这 k 家的利润之和  200.83，差 -100.83
  k =  4 家合并：合并后价格 3.125，合并的这家利润  126.56，合并前这 k 家的利润之和  267.77，差 -141.21
  k =  5 家合并：合并后价格 3.286，合并的这家利润  165.31，合并前这 k 家的利润之和  334.71，差 -169.40
  k =  6 家合并：合并后价格 3.500，合并的这家利润  225.00，合并前这 k 家的利润之和  401.65，差 -176.65
  k =  7 家合并：合并后价格 3.800，合并的这家利润  324.00，合并前这 k 家的利润之和  468.60，差 -144.60
  k =  8 家合并：合并后价格 4.250，合并的这家利润  506.25，合并前这 k 家的利润之和  535.54，差  -29.29
  k =  9 家合并：合并后价格 5.000，合并的这家利润  900.00，合并前这 k 家的利润之和  602.48，差 +297.52
  k = 10 家合并：合并后价格 6.500，合并的这家利润 2025.00，合并前这 k 家的利润之和  669.42，差 +1355.58

两家合并，合并后成本降到多少，这次合并才有赚头？
  成本要从 2.00 降到 1.714（降 14.3%）；这时价格 2.871，合并前是 2.818
想让价格不涨，成本要降到 1.182（降 40.9%），比「合并有赚头」要求的 14.3% 多得多

demand = lambda p: np.maximum(1100.0 - 100.0 * p, 0.0)          # 把 P = 11 - 0.01 Q 反过来写：Q = 1100 - 100 P

def bertrand_br(p_rival, c, step):
    """产品一模一样、产能不限。我的价格只能是 step 的整数倍。在这些价格里，选利润最大的；
    如果有好几个一样好，选最接近我自己成本的那个（没有赚头，就别乱降价）。"""
    prices = np.arange(0.0, 11.0 + step / 2, step)
    d = demand(prices)
    sold = np.where(prices < p_rival - 1e-9, d, np.where(prices > p_rival + 1e-9, 0.0, d / 2))   # 便宜的全拿，一样平分
    pi = (prices - c) * sold
    best = np.flatnonzero(pi >= pi.max() - 1e-9)
    return prices[best[np.argmin(abs(prices[best] - c))]]

def price_war(p1, p2, c1=2.0, c2=2.0, step=0.01, max_moves=5000):
    """两家轮流定价（甲先），每次都是对对方当前价格的最佳反应，直到谁也不想动"""
    path = [(p1, p2)]; still = 0
    for t in range(max_moves):
        if t % 2 == 0: new = bertrand_br(p2, c1, step); changed = abs(new - p1) > 1e-9; p1 = new
        else:          new = bertrand_br(p1, c2, step); changed = abs(new - p2) > 1e-9; p2 = new
        path.append((p1, p2)); still = 0 if changed else still + 1
        if still == 2: break
    return np.array(path)

# 从垄断价 6.5 出发，成本一样
w = price_war(6.5, 6.5)
print(f"一分钱一档：从 (6.50, 6.50) 出发，共调整 {int((abs(np.diff(w, axis=0)).sum(axis=1) > 1e-9).sum())} 次，停在 ({w[-1][0]:.2f}, {w[-1][1]:.2f})；每家利润 {(w[-1][0] - MC) * demand(w[-1][0]) / 2:.2f}")
print("前几步：" + " -> ".join(f"({a:.2f},{b:.2f})" for a, b in w[:6]))

# 价格能报到多细，会不会改变结局？
print("\n价格的最小单位 step 不同：")
for step in (0.01, 0.05, 0.25, 1.00):
    w = price_war(6.5, 6.5, step=step)
    print(f"  step = {step:>4.2f}：停在 ({w[-1][0]:.2f}, {w[-1][1]:.2f})，每家利润 {(w[-1][0] - MC) * demand(w[-1][0]) / 2:7.2f}")

# 成本不一样：甲 2 元，乙 3 元（第 16 讲思考题 3）
w = price_war(6.5, 6.5, c1=2.0, c2=3.0)
d = demand(w[-1][0])
print(f"\n成本 2 vs 3：停在 ({w[-1][0]:.2f}, {w[-1][1]:.2f})；甲独占，卖 {d:.0f} 根，利润 {(w[-1][0] - 2.0) * d:.2f}；乙利润 0")

一分钱一档：从 (6.50, 6.50) 出发，共调整 450 次，停在 (2.01, 2.01)；每家利润 4.50
前几步：(6.50,6.50) -> (6.49,6.50) -> (6.49,6.48) -> (6.47,6.48) -> (6.47,6.46) -> (6.45,6.46)

价格的最小单位 step 不同：
  step = 0.01：停在 (2.01, 2.01)，每家利润    4.50
  step = 0.05：停在 (2.05, 2.05)，每家利润   22.38
  step = 0.25：停在 (2.25, 2.25)，每家利润  109.38
  step = 1.00：停在 (3.00, 3.00)，每家利润  400.00

成本 2 vs 3：停在 (2.99, 3.00)；甲独占，卖 801 根，利润 792.99；乙利润 0

# 产能有限（每家最多 k 根）。便宜的先卖，卖完了顾客才去另一家；愿付价最高的顾客，先被便宜的那家拿走
def sales(p_me, p_rival, k_me, k_rival):
    d_me, d_rival = demand(p_me), demand(p_rival)
    cheaper = np.minimum(k_me, d_me)                                                           # 我更便宜：顾客都来，但我最多卖 k_me
    dearer = np.minimum(k_me, np.maximum(d_me - np.minimum(k_rival, d_rival), 0.0))            # 我更贵：对方先卖，剩下的才是我的
    same = np.minimum(k_me, np.maximum(d_me / 2, d_me - k_rival))                              # 一样价：平分，对方卖不完的归我
    return np.where(p_me < p_rival - 1e-9, cheaper, np.where(p_me > p_rival + 1e-9, dearer, same))

def capacity_br(p_rival, k_me, k_rival, step):
    prices = np.arange(0.0, 11.0 + step / 2, step)
    pi = (prices - MC) * sales(prices, p_rival, k_me, k_rival)
    best = np.flatnonzero(pi >= pi.max() - 1e-9)
    return prices[best[-1]]                                      # 一样好的，选最高的价格

def pure_equilibria(k, step=0.05):
    """网格上所有「谁也不想单独改价」的价格组合（两家产能一样，所以只要查一家的利润）"""
    prices = np.arange(MC, 11.0 + step / 2, step)
    best_vs = np.array([((prices - MC) * sales(prices, r, k, k)).max() for r in prices])       # 对方定 r，我最多能赚多少
    out = []
    for i, p1 in enumerate(prices):
        for j, p2 in enumerate(prices):
            u1 = (p1 - MC) * sales(p1, p2, k, k); u2 = (p2 - MC) * sales(p2, p1, k, k)
            if u1 >= best_vs[j] - 1e-9 and u2 >= best_vs[i] - 1e-9: out.append((round(float(p1), 2), round(float(p2), 2)))
    return out

print("产能各 k 根，有没有「谁也不想单独改价」的价格组合？（价格取到 5 分钱，只列出最低的一个）")
for k in (100, 200, 250, 300, 350, 500, 700, 800, 900, 1000):
    eq = pure_equilibria(k)
    print(f"  k = {k:>4}：" + (f"有，价格 {eq[0][0]:.2f}（网格上共 {len(eq)} 个组合）" if eq else "没有"))

# 没有稳定的价格，会怎样？让两家轮流、各自对对方当前的价格做最佳反应，看价格怎么走
def edgeworth(k, start=(6.0, 6.0), moves=3000, step=0.01):
    p = list(start); path = [tuple(p)]
    for t in range(moves):
        i = t % 2
        p[i] = capacity_br(p[1 - i], k, k, step); path.append(tuple(p))
    return np.array(path)

print("\n轮流最佳反应（价格一分钱一档，从 (6.00, 6.00) 出发，丢掉前 1000 次，看后面的）：")
for k in (300, 500, 700):
    w = edgeworth(k)[1000:]
    avg_profit = np.mean([((a - MC) * sales(a, b, k, k) + (b - MC) * sales(b, a, k, k)) / 2 for a, b in w])
    jumps = [t for t in range(1, len(w)) if w[t, 0] - w[t - 1, 0] > 0.5]
    period = f"，每 {jumps[1] - jumps[0]} 次调整循环一圈" if len(jumps) > 1 else "，不循环"
    print(f"  k = {k}：价格在 {w.min():.2f} 到 {w.max():.2f} 之间，平均价格 {w.mean():.2f}，每家平均利润 {avg_profit:.1f}{period}")

产能各 k 根，有没有「谁也不想单独改价」的价格组合？（价格取到 5 分钱，只列出最低的一个）
  k =  100：有，价格 9.00（网格上共 1 个组合）
  k =  200：有，价格 7.00（网格上共 1 个组合）
  k =  250：有，价格 6.00（网格上共 2 个组合）
  k =  300：有，价格 5.00（网格上共 2 个组合）
  k =  350：没有
  k =  500：没有
  k =  700：没有
  k =  800：有，价格 2.10（网格上共 1 个组合）
  k =  900：有，价格 2.00（网格上共 2 个组合）
  k = 1000：有，价格 2.00（网格上共 2 个组合）

轮流最佳反应（价格一分钱一档，从 (6.00, 6.00) 出发，丢掉前 1000 次，看后面的）：
  k = 300：价格在 5.01 到 5.01 之间，平均价格 5.01，每家平均利润 901.5，不循环
  k = 500：价格在 2.81 到 4.00 之间，平均价格 3.40，每家平均利润 525.5，每 120 次调整循环一圈
  k = 700：价格在 2.15 到 3.00 之间，平均价格 2.58，每家平均利润 237.0，每 86 次调整循环一圈

def collusion(n):
    """N 家摊约好：每家卖垄断量的 1/N。返回：守约的利润、偷偷多卖一次的利润、永远回到 Cournot 的利润"""
    q_m = (A - MC) / (2 * B) / n
    q_dev = best_reply((n - 1) * q_m)
    q_c = (A - MC) / (B * (n + 1))
    return profit(q_m, (n - 1) * q_m), profit(q_dev, (n - 1) * q_m), profit(q_c, (n - 1) * q_c)

def delta_star(stay, dev, punish):
    """偷偷多卖一次，今天多得 dev - stay；之后每期少得 stay - punish。两边打平的耐心程度"""
    return (dev - stay) / (dev - punish)

# 在"无限期"里（这里取 4000 期，之后的折现小到可以忽略）直接比一比两条路
def value(first, later, delta, T=4000):
    return first + sum(later * delta ** t for t in range(1, T))

print("N 家摊合谋，偷偷多卖一次，被发现后永远回到 Cournot 的卖量：")
print("  N   守约    偷偷多卖   惩罚期的   门槛 δ*   δ*±0.01 时，守约的总收益 − 偷偷多卖的")
for n in (2, 3, 5, 10):
    stay, dev, pun = collusion(n); d = delta_star(stay, dev, pun)
    gaps = [value(stay, stay, x) - value(dev, pun, x) for x in (d - 0.01, d + 0.01)]
    print(f"{n:>3} {stay:>7.1f} {dev:>9.1f} {pun:>9.1f} {d:>9.3f}   {gaps[0]:>+8.1f}, {gaps[1]:>+8.1f}")

print("\n换成更狠的惩罚，或者换成比价格的市场，门槛 δ* 怎么变？")
print("  N   数量竞争      数量竞争        价格竞争（偷偷的人在垄断价上降一分钱，")
print("      Cournot 惩罚  价格战惩罚      吃下全部游客）+ 价格战惩罚")
p_m = (A + MC) / 2
for n in (2, 3, 5, 10):
    stay, dev, pun = collusion(n)
    steal = (p_m - 0.01 - MC) * demand(p_m - 0.01)
    print(f"{n:>3} {delta_star(stay, dev, pun):>10.3f} {delta_star(stay, dev, 0.0):>14.3f} {delta_star((p_m - MC) * demand(p_m) / n, steal, 0.0):>14.3f}（1 - 1/N = {1 - 1 / n:.3f}）")

delta = 0.9
print(f"\nδ = {delta} 时，最多几家还能维持合谋：")
for name, f in (("数量竞争 + Cournot 惩罚", lambda n: delta_star(*collusion(n))),
                ("数量竞争 + 价格战惩罚", lambda n: delta_star(*collusion(n)[:2], 0.0)),
                ("价格竞争 + 价格战惩罚", lambda n: 1 - 1 / n)):
    print(f"  {name}：{max(n for n in range(2, 400) if f(n) <= delta)} 家")

N 家摊合谋，偷偷多卖一次，被发现后永远回到 Cournot 的卖量：
  N   守约    偷偷多卖   惩罚期的   门槛 δ*   δ*±0.01 时，守约的总收益 − 偷偷多卖的
  2  1012.5    1139.1     900.0     0.529       -5.0,     +5.2
  3   675.0     900.0     506.2     0.571       -9.0,     +9.4
  5   405.0     729.0     225.0     0.643      -13.7,    +14.5
 10   202.5     612.6      66.9     0.752      -21.1,    +22.9

换成更狠的惩罚，或者换成比价格的市场，门槛 δ* 怎么变？
  N   数量竞争      数量竞争        价格竞争（偷偷的人在垄断价上降一分钱，
      Cournot 惩罚  价格战惩罚      吃下全部游客）+ 价格战惩罚
  2      0.529          0.111          0.500（1 - 1/N = 0.500）
  3      0.571          0.250          0.667（1 - 1/N = 0.667）
  5      0.643          0.444          0.800（1 - 1/N = 0.800）
 10      0.752          0.669          0.900（1 - 1/N = 0.900）

δ = 0.9 时，最多几家还能维持合谋：
  数量竞争 + Cournot 惩罚：33 家
  数量竞争 + 价格战惩罚：37 家
  价格竞争 + 价格战惩罚：10 家