5.1.1 蒙特卡洛

什么是蒙特卡洛？

蒙特卡洛（Monte Carlo）是一种思想或一类方法，它通过随机采样来估计数学问题的数值解，本质上是通过大量采样来暴力解决问题，特别适用于复杂的积分、概率计算和优化问题等。

例子：如何利用蒙特卡洛方法估计圆周率$\pi$？

1. 问题描述：在一个边长为 2 的正方形中，内嵌一个半径为 1 的圆，我们希望估计圆的面积。

2. 生成随机样本：在正方形内随机生成$n$个点，记录落在圆内的点数$m$

3. 计算估计值：

   • 圆的面积与正方形面积的比值是$\frac{\pi r^2}{(2r)^2} = \frac{\pi}{4}$

   • 通过 $\frac{m}{n}$估计这个比值，$\pi \approx 4 \times \frac{m}{n}$

4. 增加样本数量：随着样本数量$n$增加，$\pi$的估计值会越来越接近真实值

在未知$\pi$的情况下，我们无法直接计算圆的面积，于是通过暴力采样和计算点数量的方式估计圆的面积，进而估计$\pi$的值。

蒙特卡洛如何帮助我们处理复杂后验分布？

在求解后验分布的时候，我们经常要在分布上做一些计算，比如求期望、方差、或其他积分。但是如果碰上复杂分布、高维变量，就会让计算变得特别困难。这时候我们就可以尝试绕过解析解，通过大量采样的方法直接得到统计量。

例如，我们想对一个连续随机变量$X$求期望，概率分布函数为：

$$p(x) = e^{-0.5 \cdot x^2} + \sin(x) + x^3$$

这是一个非标准概率密度函数，我们无法直接求积分，但是我们可以从$p(x)$中采样$x_1, x_2, x_3, ...x_i$，再利用大数定理（在本书中不做详细解释）近似期望：

$$\mathbb{E}[X] = \int x_i\cdot p(x_i)\,d(x)\approx \frac{1}{n} \sum_{i}^{n} x_i$$

同样，对离散连续变量$X$求期望，虽然不涉及积分困难，但当变量的维度很高的时候$(X \in \mathbb{R}^D )$，会出现维数灾难，求和计算难度随着维度增加而指数级增加。我们也可以通过从$P(x)$采样来近似期望：

$$\mathbb{E}[X] = \sum_{i}^{n} x_i\cdot P(x_i)\approx \frac{1}{n} \sum_{i}^{n}x_i$$

一句话总结：遇到棘手的概率分布形式，避开难求的解析解，通过大量采样进行暴力计算估计。

注意：这里我们已知概率函数$p(x)$, 只是这个$p(x)$可能比较复杂，直接积分比较难求。

蒙特卡洛的实际应用案例

接下来我们借着一个模拟股票走势的代码案例，展示蒙特卡洛的用法。

假定某个股票目前的股价是$S_t$, 预测时间$\Delta t$之后的股价$S_{t+1}$，金融工程中有公式：

$$S_{t+1} = S_t + \hat{\mu} S_t \Delta t + \sigma S_t \epsilon \sqrt{\Delta t}$$

假定初始股价$S_0=10$，我们可以模拟一下，一年之后 (244 个交易日 ) 的股价分布。

1. 首先，定义初始状态和模拟次数。

from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt

s0 = 10. #刚买入的股价，初始状态
mu = 0.15
sigma = 0.2 
nSimu = 10000 # 模拟多少次
nStep = 244 # 每一次模拟走多少步, 1年只有244个交易日
dt = 1/nStep # 由于公式中的单位是年，而一年中有244个交易日，所以dt是理想化的每次股价变化间隔时间。

模拟10000次，得到一年后的股价分布。

price = np.empty((nSimu, nStep))
e = norm()
for i in range (nSimu):
    s = s0
    price[i, 0] = s
    for j in range(nStep):
        et = e.rvs() # 
        s = s + mu*s*dt + sigma*s*et*np.sqrt(dt) #股价的计算公式
        price[i, j] = s

画出一年后股价的分布。

plt.hist(price[:,-1], bins=30, density=True)
plt.ylabel('pdf')
plt.xlabel('Price')