5D Library › Credit Cycles, Self-Organised Criticality and the Statistical Physics of Market Phase Transitions · 中文
Chapter 2 Sandpiles and the Great Moderation: Self-Organised Criticality
An avalanche set off by one grain
In 1987 Bak, Tang and Wiesenfeld proposed an extremely simple model (Bak et al. 1987). Each cell of a square grid holds some grains of sand; at each step one grain is added to a random cell; any cell holding 4 grains “topples”, passing one grain to each of its four neighbours, and grains that fall off the edge are lost. A neighbour that receives grains may itself reach 4 and topple in turn, producing an “avalanche”.
What is striking about the model is that, with nobody adjusting any parameter, the pile evolves by itself into a special state in which adding one grain may cause nothing or may cause an avalanche spreading across the whole grid, and the distribution of avalanche sizes approximates a power law:
Bak called the phenomenon self-organised criticality and in later work suggested that it might explain the frequency of “large events” in earthquakes, extinctions and even economic fluctuations (Bak 1996).
Calculation
A 48 × 48 grid. First add 10,000 grains to let the pile settle into its steady state, then add 40,000 more and record how many topplings each one causes. Finally estimate the power-law exponent by maximum likelihood (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
Of the 40,000 grains, more than half caused no toppling at all; the largest avalanche caused more than eight thousand topplings. The largest 1 % of avalanches accounted for about 28 % of all topplings. That is the statistical signature of “small causes, large effects”: most disturbances do not matter, a very few cause most of the damage, and nothing about a grain beforehand tells you how large an avalanche it will set off.
But the figure also shows two problems. First, on a finite grid avalanche size has an upper limit set by the size of the grid, so the large end is cut off sharply and the distribution is not a straight line. Second, the estimated exponent rises from 1.30 to 1.41 as the lower limit of the fit, , rises from 2 to 10. The program gives both the exact discrete maximum-likelihood estimate (normalised with the Hurwitz zeta function) and the approximation of Clauset and colleagues; that approximation is accurate only for lower limits of about 6 or more, and for these data the two happen to agree. Values reported in the literature for the avalanche-size exponent of the two-dimensional sandpile are not consistent, and its scaling properties are still debated. “Is it a power law?” and “what is the exponent?” are not questions that a glance at a log–log plot can answer.
How to estimate a power-law exponent
A very common approach is to plot a histogram, take logarithms, fit a straight line by least squares and read the exponent off the slope. Clauset, Shalizi and Newman explained systematically why this is unreliable (Clauset et al. 2009): in a log–log histogram the bins in the tail hold very few points or none, so they are very noisy, yet in a least-squares fit they carry the same weight as well-filled bins; and empty bins can only be dropped, which itself causes bias.
Calculation
Draw 300 samples of 1,000 numbers each from a Pareto distribution with known exponent , and estimate the exponent three ways. Finally, “fit” a power law in the same way to a log-normal sample, which is not a power law.
# 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)The maximum-likelihood estimate averages exactly 2.50, with a standard deviation of about 0.05, matching the theoretical standard error. Least squares on equal-width bins underestimates the exponent at about 1.5, and the result varies greatly from sample to sample; least squares on logarithmic bins is much better but still biased. The last line deserves more caution: a log-normal distribution is not a power law, yet looking only at its tail, any method yields a plausible-looking “exponent”. Clauset and colleagues therefore recommend estimating the exponent and the lower cut-off by maximum likelihood, then using goodness-of-fit tests and likelihood-ratio tests against alternatives such as the log-normal or a power law with a cut-off (Clauset et al. 2009).
Is the market a sandpile?
Evidence
Fat tails in financial returns. The return distributions of stocks and other assets have fatter tails than a normal distribution; in many studies the tail exponent of returns lies between 2 and 5 (Cont 2001). Gabaix and colleagues report a tail exponent of about 3 for stock price changes (the “cubic law”) and propose an explanation unrelated to self-organised criticality: the sizes of large institutional investors themselves follow a power law, and their large trades cause large price moves (Gabaix et al. 2003).
So even if market returns have fat tails, that does not prove the market is in a “self-organised critical state”: fat tails can come from many different mechanisms. The sandpile has no prices, no traders, no expectations and no central bank. It is a thinking tool for systems that are “slowly driven, triggered at thresholds and connected in chains”, and applied to markets it is only an analogy.
The paradox of the Great Moderation
From the mid-1980s to 2007 the volatility of output and inflation in the United States and other economies fell markedly, a period usually called the “Great Moderation”. It was followed by the worst financial crisis in decades. That raises a question: can long calm itself build up larger risk?
One simple mechanism lies in risk-management rules themselves. Many financial institutions set leverage by rules such as value at risk (VaR): the lower the measured volatility, the higher the leverage allowed. Adrian and Shin found that the leverage of US broker-dealers is procyclical: it rises when asset prices rise and measured risk falls (Adrian and Shin 2010).
Calculation
An investor sets leverage by the rule “99 % one-day value at risk must not exceed 4 % of equity”. Compare the leverage allowed at different measured volatilities, and the loss of equity caused by the same one-day fall of 10 %. The last column is the probability of a one-day fall of 10 % under the normal model that the rule itself assumes.
# 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).At a measured annual volatility of 6 % the rule allows about 4.5 times leverage, and a single fall of 10 % costs nearly half the equity; at a measured volatility of 30 % leverage is below 1 and the same fall costs only about 9 %. Under the normal model the rule relies on, the probability of a one-day fall of 10 % in calm times is too small to mean anything, which is exactly the false sense of security the model gives. Chapter 4 shows that the tails of real returns are much fatter than normal.
Evidence
Crises after low volatility. In Brunnermeier and Sannikov’s macro-financial model a “volatility paradox” appears: the lower the exogenous risk, the more the financial sector raises leverage, and endogenous systemic risk rises (Brunnermeier and Sannikov 2014). Using stock-market data for dozens of countries over about two centuries, Danielsson, Valenzuela and Zer found that when volatility stays below its trend for a long time, the probability of a subsequent banking crisis rises, while high volatility itself is not a good predictor of crises (Danielsson et al. 2018).
Hypothesis · not established
Suppressed-volatility hypothesis (not yet established). In a financial system where leverage is set by measured volatility, persistently suppressing day-to-day volatility (whether by policy or by market structure) raises leverage and the sensitivity of asset prices to small shocks, making later drawdowns larger. It can be tested with cross-country panel data: after controlling for credit growth, does “the length of time volatility has stayed below trend” still predict later maximum drawdowns or crises? If the relation disappears once credit growth and leverage are controlled for, or fails out of sample, the hypothesis is weakened.
Summary
The sandpile model shows that a slowly driven system with thresholds can produce avalanches of every size by itself, with a few large avalanches causing most of the damage. Power-law exponents should be estimated by maximum likelihood and compared with alternative distributions; least squares on a log–log histogram gives biased results. The fat tails of financial returns have several possible explanations, and the sandpile is only an analogy. “Long calm builds up larger risk” has a model mechanism and historical evidence behind it, but should still be treated as a hypothesis. The next chapter looks at what happens when calm breaks: how margins and fire sales amplify each other.
References
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).