5D Library › 信贷周期、自组织临界性与市场相变的统计物理 · English
第二章 沙堆与大缓和：自组织临界性
一粒沙引起的雪崩
1987 年，巴克、汤超与维森菲尔德提出了一个极其简单的模型 (Bak et al. 1987)。在一个方格网上，每个格子里放着一些沙粒；每次在随机一格加一粒沙；任何格子的沙粒达到 4 粒，就「倒塌」，把 4 粒沙分给上下左右四个邻居，落到边缘外的沙粒就丢失了。邻居得到沙粒后可能也达到 4 粒，于是继续倒塌，形成一次「雪崩」。
这个模型的惊人之处在于：没有人调节任何参数，沙堆会自己演化到一种特殊的状态，在这种状态下，加一粒沙可能什么都不引起，也可能引起波及整个网格的雪崩，雪崩的规模分布近似一个幂律：
巴克把这种现象称为自组织临界性，并在后来的著作中设想它可以解释地震、物种灭绝乃至经济波动中「大事件」的频繁出现 (Bak 1996)。
推导
48 × 48 的网格，先加 10,000 粒沙让沙堆进入稳定状态，然后再加 40,000 粒，记录每一粒引起的倒塌次数。最后用最大似然法估计幂律指数 (Clauset et al. 2009)。
# The Bak-Tang-Wiesenfeld (1987) sandpile on an L x L grid with open edges. Drop one grain on a random
# site; any site holding 4 or more grains topples, sending one grain to each neighbour (grains that
# leave the edge are lost). Avalanche size s = number of topplings caused by one grain.
# Nobody tunes the pile: driven slowly, it organises itself into a state with avalanches of all sizes.
import math, random
def sandpile(L=48, warm=10000, n=40000, seed=1):
rnd = random.Random(seed)
z = [[0] * L for _ in range(L)]
sizes = []
for k in range(warm + n):
i, j = rnd.randrange(L), rnd.randrange(L)
z[i][j] += 1
s, stack = 0, [(i, j)] if z[i][j] >= 4 else []
while stack:
a, b = stack.pop()
if z[a][b] < 4: continue
z[a][b] -= 4; s += 1
if z[a][b] >= 4: stack.append((a, b))
for c, d in ((a + 1, b), (a - 1, b), (a, b + 1), (a, b - 1)):
if 0 <= c < L and 0 <= d < L:
z[c][d] += 1
if z[c][d] == 4: stack.append((c, d))
if k >= warm: sizes.append(s)
return sizes, sum(map(sum, z)) / L ** 2
def mle_tau(xs, xmin):
"""Discrete power-law exponent, Clauset, Shalizi & Newman (2009) approximation, for x >= xmin.
The approximation is accurate only for xmin of about 6 or more; mle_exact below is exact."""
t = [x for x in xs if x >= xmin]
return 1 + len(t) / sum(math.log(x / (xmin - 0.5)) for x in t), len(t)
def hurwitz(tau, xmin, K=2000):
"""zeta(tau, xmin) = sum over k >= xmin of k^-tau: direct sum, then an Euler-Maclaurin tail."""
s = sum(k ** -tau for k in range(xmin, K))
return s + K ** (1 - tau) / (tau - 1) + 0.5 * K ** -tau + tau * K ** (-tau - 1) / 12
def mle_exact(xs, xmin):
"""Exact discrete MLE: maximise -tau * sum(ln x) - n ln zeta(tau, xmin) by golden-section search."""
t = [x for x in xs if x >= xmin]
n, sl = len(t), sum(math.log(x) for x in t)
f = lambda a: -a * sl - n * math.log(hurwitz(a, xmin))
lo, hi, g = 1.01, 3.0, (math.sqrt(5) - 1) / 2
for _ in range(60):
a, b = hi - g * (hi - lo), lo + g * (hi - lo)
if f(a) > f(b): hi = b
else: lo = a
return (lo + hi) / 2
if __name__ == "__main__":
sizes, dens = sandpile()
pos = [s for s in sizes if s > 0]
print(f"{len(sizes)} grains after warm-up: mean height {dens:.2f} grains per site; "
f"{len(sizes) - len(pos)} caused no toppling; largest avalanche {max(sizes)} topplings")
print("avalanche size share of avalanches (logarithmic bins)")
lo = 1
while lo <= max(pos):
hi = lo * 4
c = sum(1 for s in pos if lo <= s < hi)
print(f" {lo:6d}-{hi - 1:<6d} {c / len(pos):8.4f} " + "#" * round(60 * c / len(pos)))
lo = hi
for xmin in (2, 5, 10):
tau, m = mle_tau(pos, xmin)
print(f"maximum-likelihood exponent for s >= {xmin:2d}: exact {mle_exact(pos, xmin):.2f},"
f" approximation {tau:.2f} ({m} avalanches)")
top = sorted(sizes)[-len(sizes) // 100:]
print(f"the largest 1 % of avalanches account for {sum(top) / sum(sizes):.0%} of all topplings")
# 40000 grains after warm-up: mean height 2.09 grains per site; 22560 caused no toppling; largest avalanche 8063 topplings
# avalanche size share of avalanches (logarithmic bins)
# 1-3 0.2450 ###############
# 4-15 0.2014 ############
# 16-63 0.1899 ###########
# 64-255 0.1677 ##########
# 256-1023 0.1456 #########
# 1024-4095 0.0490 ###
# 4096-16383 0.0014
# maximum-likelihood exponent for s >= 2: exact 1.30, approximation 1.30 (15161 avalanches)
# maximum-likelihood exponent for s >= 5: exact 1.36, approximation 1.36 (12510 avalanches)
# maximum-likelihood exponent for s >= 10: exact 1.41, approximation 1.41 (10771 avalanches)
# the largest 1 % of avalanches account for 28% of all topplings
40,000 粒沙中，一半以上没有引起任何倒塌；最大的一次雪崩引起了八千多次倒塌。最大的 1% 的雪崩，贡献了全部倒塌次数的约 28%。这就是「小原因、大后果」的统计特征：绝大多数扰动无关紧要，极少数扰动造成了大部分破坏，而且事先无法从那一粒沙本身看出它会引起多大的雪崩。
但图中也能看到两个问题。第一，在有限的网格中，雪崩规模有一个由网格大小决定的上限，大雪崩一端急剧截断，分布并不是一条直线。第二，估计出的指数随着拟合的下限 从 2 提高到 10，从 1.30 升到 1.41。程序同时给出精确的离散最大似然估计（用赫尔维茨 ζ 函数归一化）与克劳塞特等人的近似公式；那个近似只在下限约为 6 以上时才准确，在这组数据中两者恰好一致。文献中二维沙堆模型雪崩规模指数的报告值不尽一致，它的标度性质至今仍有争议。「是不是幂律」「指数是多少」并不是看一眼双对数图就能回答的问题。
怎样估计幂律指数
一种很常见的做法是：画出直方图，取对数，用最小二乘法拟合一条直线，斜率就是指数。克劳塞特、沙利兹与纽曼系统地说明了这种做法为什么不可靠 (Clauset et al. 2009)：双对数直方图中，尾部的箱子只有很少几个数据点甚至为零，它们的噪声很大，却在最小二乘拟合中和数据充足的箱子权重相同；而空箱子只能被丢掉，这本身就造成偏差。
推导
从已知指数 的帕累托分布中抽取 300 个样本，每个 1,000 个数，用三种方法估计指数。最后，对一个不是幂律的对数正态样本，也照样「拟合」一个幂律。
# Estimating a power-law exponent. Samples are drawn from a known Pareto law, density
# p(x) ~ x^(-alpha) for x >= 1 with alpha = 2.5, then the exponent is estimated three ways:
# MLE: alpha = 1 + n / sum(ln x) (Clauset, Shalizi & Newman 2009)
# LS-linear: least squares on log(count) vs log(x) with equal-width bins (a common shortcut)
# LS-log: least squares on a histogram with logarithmic bins, counts divided by bin width
# 300 repeated samples of n = 1,000 show the bias and spread of each method.
import math, random, statistics as st
def slope(xs, ys):
mx, my = st.fmean(xs), st.fmean(ys)
return sum((a - mx) * (b - my) for a, b in zip(xs, ys)) / sum((a - mx) ** 2 for a in xs)
def ls_linear(x, width=1.0):
counts = {}
for v in x: counts[int((v - 1) // width)] = counts.get(int((v - 1) // width), 0) + 1
pts = [(math.log(1 + (k + 0.5) * width), math.log(c)) for k, c in counts.items()] # empty bins dropped
return -slope(*zip(*pts))
def ls_log(x, ratio=1.5):
edges = [ratio ** k for k in range(60)]
pts = []
for a, b in zip(edges, edges[1:]):
c = sum(1 for v in x if a <= v < b)
if c: pts.append((math.log(math.sqrt(a * b)), math.log(c / (b - a))))
return -slope(*zip(*pts))
def mle(x, xmin=1.0):
return 1 + len(x) / sum(math.log(v / xmin) for v in x)
rnd = random.Random(7)
ALPHA, N, R = 2.5, 1000, 300
res = {"MLE": [], "LS, equal-width bins": [], "LS, logarithmic bins": []}
for _ in range(R):
x = [(1 - rnd.random()) ** (-1 / (ALPHA - 1)) for _ in range(N)]
res["MLE"].append(mle(x)); res["LS, equal-width bins"].append(ls_linear(x)); res["LS, logarithmic bins"].append(ls_log(x))
print(f"true alpha = {ALPHA}; {R} samples of n = {N}")
for k, v in res.items():
print(f" {k:22s} mean {st.fmean(v):.2f} sd {st.stdev(v):.2f}")
print(f" MLE standard error formula (alpha - 1)/sqrt(n) = {(ALPHA - 1) / math.sqrt(N):.3f}")
# A log-normal sample also looks straight on a log-log plot over a limited range
y = [math.exp(rnd.gauss(0, 1.6)) for _ in range(20000)]
tail = [v for v in y if v >= 5]
print(f"log-normal sample (not a power law): tail above x = 5 holds {len(tail)} points;"
f" fitting a power law anyway gives alpha = {mle(tail, 5):.2f} (log-bin LS {ls_log(tail):.2f})")
# true alpha = 2.5; 300 samples of n = 1000
# MLE mean 2.50 sd 0.05
# LS, equal-width bins mean 1.51 sd 0.27
# LS, logarithmic bins mean 2.35 sd 0.15
# MLE standard error formula (alpha - 1)/sqrt(n) = 0.047
# log-normal sample (not a power law): tail above x = 5 holds 3138 points; fitting a power law anyway gives alpha = 2.19 (log-bin LS 2.19)最大似然估计的平均值正好是 2.50，标准差约 0.05，与理论上的标准误一致。等宽分箱的最小二乘把指数低估到约 1.5，而且每次结果相差很大；对数分箱的最小二乘好得多，但仍有偏差。最后一行更值得警惕：对数正态分布不是幂律，但只看它的尾部，任何方法都能给出一个看似合理的「指数」。克劳塞特等人因此建议：先用最大似然法估计指数和下限，再用拟合优度检验与似然比检验，和对数正态、带截断的幂律等替代分布比较 (Clauset et al. 2009)。
市场是一个沙堆吗
证据
金融收益率的重尾。 股票等资产的收益率分布比正态分布有更厚的尾部；在许多研究中，收益率尾部的指数在 2 到 5 之间 (Cont 2001)。加贝等人报告，股票价格变动的尾部指数约为 3（所谓「三次方定律」），并提出了一个与自组织临界性无关的解释：大机构投资者的规模本身服从幂律分布，它们的大额交易造成了价格的大幅变动 (Gabaix et al. 2003)。
所以，即使市场收益率有重尾，这也不能证明市场处于「自组织临界状态」：重尾可以来自许多不同的机制。沙堆模型里没有价格、没有交易者、没有预期，也没有中央银行。它是理解「缓慢驱动、阈值触发、连锁传递」这类系统的一个思想工具，用于市场时只是类比。
大缓和悖论
从 1980 年代中期到 2007 年，美国等经济体的产出与通胀波动明显下降，这一时期通常被称为「大缓和」。随后发生的是几十年来最严重的金融危机。这引出了一个问题：长期平静本身，会不会累积更大的风险？
一种简单的机制来自风险管理规则本身。许多金融机构按照风险价值（VaR）之类的规则确定杠杆：测得的波动越低，允许的杠杆越高。阿德里安与申铉松发现，美国券商的杠杆是顺周期的：资产价格上涨、测得的风险下降时，杠杆上升 (Adrian and Shin 2010)。
推导
一个按「99% 单日风险价值不超过权益的 4%」确定杠杆的投资者。比较不同的测得波动率下允许的杠杆，以及同样一次单日 10% 的下跌造成的权益损失。最后一列是按这条规则自己所假设的正态模型，单日下跌 10% 的概率。
# The "volatility paradox" in its simplest form: a leveraged investor sizes positions with a
# value-at-risk rule. Measured volatility sigma (annualised) is low in calm years, so the rule allows
# more leverage - and the same sudden shock then does more damage.
# Rule: 99 % one-day VaR must not exceed 4 % of equity -> L = 0.04 / (2.326 * sigma * sqrt(1/252)),
# capped at 10x. Shock: a one-day asset fall of 10 % before positions can be cut. Numbers are illustrative.
import math
Z99, H, CAP, SHOCK = 2.326, 1 / 252, 10.0, 0.10
def leverage(sigma):
return min(CAP, 0.04 / (Z99 * sigma * math.sqrt(H)))
print("measured vol leverage allowed equity lost in a 10 % fall probability of such a fall *")
for sigma in (0.06, 0.08, 0.12, 0.16, 0.20, 0.30):
L = leverage(sigma)
loss = min(1.0, L * SHOCK)
# * under the normal model the rule itself assumes, for one day:
z = SHOCK / (sigma * math.sqrt(H))
p = 0.5 * math.erfc(z / math.sqrt(2))
print(f" {sigma:4.0%} {L:5.2f}x {loss:5.0%}{' (wiped out)' if loss >= 1 else '':12s} {p:.1e}")
print("* the normal model calls a 10 % one-day fall almost impossible in calm years; real tails are fatter (Chapter 4).")
# measured vol leverage allowed equity lost in a 10 % fall probability of such a fall *
# 6% 4.55x 45% 1.5e-154
# 8% 3.41x 34% 6.3e-88
# 12% 2.27x 23% 3.0e-40
# 16% 1.71x 17% 1.7e-23
# 20% 1.36x 14% 1.0e-15
# 30% 0.91x 9% 6.1e-08
# * the normal model calls a 10 % one-day fall almost impossible in calm years; real tails are fatter (Chapter 4).测得的年化波动率为 6% 时，规则允许约 4.5 倍杠杆，一次 10% 的下跌就会损失近一半权益；测得波动率为 30% 时，杠杆不到 1 倍，同样的下跌只损失约 9%。而在规则所依据的正态模型中，平静时期单日下跌 10% 的概率小到没有意义，这正是模型给人的虚假安全感。第四章会说明真实收益率的尾部比正态分布厚得多。
证据
低波动之后的危机。 布伦纳迈尔与萨尼科夫的宏观金融模型中出现了「波动率悖论」：外生风险越低，金融部门越倾向于提高杠杆，内生的系统性风险反而上升 (Brunnermeier and Sannikov 2014)。丹涅尔松、巴伦苏埃拉与泽尔用数十个国家、约两个世纪的股市数据发现，波动率长期低于其趋势时，随后发生银行危机的概率上升；高波动本身则不是危机的好预测指标 (Danielsson et al. 2018)。
假说 · 尚未验证
压低小波动积累大风险假说（尚未确立）。在一个以测得波动率决定杠杆的金融体系中，持续压低日常波动（无论来自政策还是来自市场结构）会提高杠杆与资产价格对小冲击的敏感度，使随后的回撤更大。它可以用跨国面板数据检验：控制信贷增长之后，「低于趋势的波动率持续时间」是否仍能预测随后的最大回撤或危机。如果控制信贷增长与杠杆之后这一关系消失，或者在样本外不成立，这一假说就被削弱。
本章小结
沙堆模型显示，一个缓慢驱动、带有阈值的系统可以自发地产生各种规模的雪崩，少数大雪崩造成大部分破坏。估计幂律指数要用最大似然法并与替代分布比较，双对数直方图的最小二乘法会给出有偏的结果。金融收益率的重尾有多种可能的解释，沙堆只是类比。「长期平静累积更大风险」有模型机制与历史证据支持，但仍应作为假说。下一章讨论平静被打破时，保证金与抛售怎样互相放大。
参考文献
Fan, Y. (2026). Credit Cycles, Self-Organised Criticality and the Statistical Physics of Market Phase Transitions: From Minsky's Financial Instability Hypothesis and Sandpile Power-Law Avalanches to Implied-Volatility Skew, Read Beside Classical Chinese Thought on Rise and Decline [信贷周期、自组织临界性与市场相变的统计物理] (Version 1.0.0). 5D Library. https://doi.org/10.5281/zenodo.23754711 · Licensed under Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0).