蒙特卡洛估计

用随机抽样的平均来回答

金融中的许多问题没有整洁的公式：组合一年内亏损 20% 的概率，收益取决于整条路径的期权价格。蒙特卡洛用蛮力回答它们。把所求量写成平均值 E[f(X)]，模拟 X 的大量抽样，对每个抽样应用 f 再取平均。由大数定律，样本均值收敛到真实期望。概率只是「是/否」指示变量的期望，所以「事件发生的模拟占比」就是它的估计。

import math
import numpy as np

# Estimate pi: the quarter circle x^2 + y^2 <= 1 covers pi/4 of the unit square
rng = np.random.default_rng(42)
x, y = rng.random(100000), rng.random(100000)
inside = x ** 2 + y ** 2 <= 1                        # True/False for every point
print(round(4 * inside.mean(), 4))                    # 3.1378   close to 3.1416, with 100,000 points

模拟价格与概率

要模拟一只股票，先抽取标准正态冲击，再代入上一节的对数正态公式：全部 n 个终值价格由一个数组表达式给出。之后关于价格分布的任何问题都变成对数组求平均。收于 80 以下的概率是 np.mean(prices < 80)，期望价格是 prices.mean()，亏损的分位数是 np.quantile(...)。把生成器 rng 作为参数传进函数而不要在内部创建，这样调用者可以控制可重复性。

def simulate_terminal_prices(rng, s0, mu, sigma, t, n):
    z = rng.standard_normal(n)
    return s0 * np.exp((mu - 0.5 * sigma ** 2) * t + sigma * math.sqrt(t) * z)

prices = simulate_terminal_prices(np.random.default_rng(1), 100, 0.05, 0.2, 1, 200000)
print(round(float(np.mean(prices < 80)), 4))          # 0.1027   P(price below 80 after a year), about 10%
print(round(float(np.mean(prices > 130)), 4))         # 0.1216   P(price above 130)

为期权定价

欧式看涨期权在到期时支付 max(S - K, 0)。在风险中性世界里，所有资产都以无风险利率 rate 增长，所以公平价格是折现后的期望收益：用 mu = rate 模拟终值价格，计算每个的收益，取平均，再乘以 exp(-rate * t)。这个简单期权有闭式公式（布莱克-斯科尔斯，下面这组数字是 10.4506），模拟结果与之接近。蒙特卡洛的意义在于：同样的几行代码，对完全没有公式的收益依然适用。

def mc_call_price(rng, s0, strike, rate, sigma, t, n):
    prices = simulate_terminal_prices(rng, s0, rate, sigma, t, n)    # risk-neutral: the price grows at the risk-free rate
    payoff = np.maximum(prices - strike, 0)                           # max(S - K, 0) for every simulated path
    return math.exp(-rate * t) * payoff.mean()                        # discounted average payoff

price = mc_call_price(np.random.default_rng(1), 100, 100, 0.05, 0.2, 1, 200000)
print(round(price, 4))        # 10.3993   the exact Black-Scholes price is 10.4506: the difference is Monte Carlo error

# More points, less error (1 / sqrt(n)): the estimate of pi with the same seed
for n in (100, 10000, 1000000):
    g = np.random.default_rng(42)
    x, y = g.random(n), g.random(n)
    print(n, round(4 * float((x * x + y * y <= 1).mean()), 4))
# 100 3.36
# 10000 3.1488
# 1000000 3.1424

• 蒙特卡洛：把所求量写成平均 E[f(X)]，模拟大量抽样，应用 f，取平均。
• 概率是比较结果的平均：np.mean(prices < 80)；把生成器作为参数传入，使运行可重复。
• 期权价格 = exp(-rate * t) * mean(max(S - K, 0))，价格按 mu = rate 模拟；同时报告误差。

下一节：蒙特卡洛误差与收敛，估计有多精确，以及如何以更少的代价获得更多精度。

实验：评分内容

• 用 estimate_pi 向正方形投点来求 pi25 分
• 用 estimate_probability 统计事件发生的频率25 分
• 用 simulate_terminal_prices 模拟股票可能的终点25 分
• 用 mc_call_price 通过模拟为期权定价25 分