Python金融分析（三）：随机过程模拟

发表于 2019-03-15 分类于他山攻错阅读次数： Disqus：

前言

如今，随机数学已经成为金融中最重要的数学分支之一。现代金融的开端，可以追溯至20世纪七八十年代，当时研究的主要目标是寻找给定模型的封闭解。

近年来，金融研究的目标也发生了巨大的变化，不仅要求对单一的金融工具给予准确估值，还要求对不同衍生品组合的估值具有一致性。类似的，对于风险管理而言，也需要提出一致的风险衡量指标，例如在险价值（Value-at-risk, VaR）和信用风险评估调整（Credit valuation adjustments, CVA）。这些任务难以寻找精确的封闭解，转而寻求通过数值方法解决。因此，随机数学，尤其是蒙特卡洛模拟的作用越发重要。

接下来，我们将从Python的角度介绍随机数、随机模拟、估值、风险度量等内容。

随机数

在Python中，与随机数产生有关的模块为numpy.random。但是，需要注意的是，所谓的随机数，只是伪随机数（pseudo-random），此外还有一个类似的概念，叫做准随机数（quasi-random），两个兄弟长得很像，但是并不能混为一谈。

伪随机过程是一个看似随机但实际并非如此的过程。伪随机序列通常表现出统计上的随机性，尽管实际上是由完全确定的生成过程产生的。

准随机序列，是多位情况下均匀分布的点，具有低差异（low-discrepancy, LD）的特性。准随机数并不是为了看上去随机，也不一定能够通过随机性的假设检验，它的目的在于均与填充要采样的整个空间。

因此，伪随机数和准随机数都使用计算算法来生成随机序列，差异在于空间中的分布。 在金融中，准随机数和伪随机数都有应用，但是似乎还是伪随机数为王。

1	%matplotlib inline

import math
import numpy as np
import numpy.random as npr
import matplotlib as mpl
import matplotlib.pyplot as plt

plt.style.use('ggplot')
mpl.rcParams['figure.figsize']= [15, 9]

1
2
3

npr.seed(0)
np.set_printoptions(precision=4)
npr.rand(10)

array([0.5488, 0.7152, 0.6028, 0.5449, 0.4237, 0.6459, 0.4376, 0.8918,
       0.9637, 0.3834])

Numpy可以直接产生简单的随机数（如下表所示），此外num.random包还可以产生不同分布的随机数，例如beta、二项式、卡方、负二项式、帕累托等等……当然，应用最多的还是（标准）正态分布。

Function	Parameters	Returns/result
rand	d0, d1, …, dn
randn	d0, d1, …, dn	A sample (or samples) from the standard normal
randint	low[, high, size]	Random integers from low (inclusive) to high (exclusive)
random_integers	low[, high, size]	Random integers between low and high, inclusive
random_sample	[size]	Random floats in the half-open interval [0.0, 1.0)
random	[size]	Random floats in the half-open interval [0.0, 1.0)
ranf	[size]	Random floats in the half-open interval [0.0, 1.0)
sample	[size]	Random floats in the half-open interval [0.0, 1.0)
choice	a[, size, replace, p]	Random sample from a given 1D array
bytes	length	Random bytes

尽管关于正态分布在金融中的使用（甚至是滥用）的批评声音不绝于耳，但是正态分布还是不可或缺的工具。之所以如此，大概有两个原因：

许多金融模型或多或少用到正态分布或者对数正态分布的条件；
不依赖正态分布的模型可以通过正态分布离散化或者近似来模拟。

伪随机数矩校正

既然我们是用伪随机数来近似随机数，还是会造成统计矩的微小偏差，例如我们生成一组标准正态分布的样本，计算其样本均值和方差，可以发现，虽然随着样本数目的增加，统计矩的偏差在缓慢缩小，但是样本非常大，偏差仍然存在。（实际上，即使是从真实的分布中抽样，同样还会存在这个问题）。

print('%15s %15s' % ('Mean', 'Std. Deviation'))
print(31 * '-')
for i in range(1, 31, 2):
    npr.seed(100)
    sn = npr.standard_normal(i ** 2 * 10000)
    print('%15.12f %15.12f' % (sn.mean(), sn.std()))
i ** 2 * 10000

           Mean  Std. Deviation
-------------------------------
 0.001150944833  1.006296354600
 0.002841204001  0.995987967146
 0.001998082016  0.997701714233
 0.001322322067  0.997771186968
 0.000592711311  0.998388962646
-0.000339730751  0.998399891450
-0.000228109010  0.998657429396
 0.000295768719  0.998877333340
 0.000257107789  0.999284894532
-0.000357870642  0.999456401088
-0.000528443742  0.999617831131
-0.000300171536  0.999445228838
-0.000162924037  0.999516059328
 0.000135778889  0.999611052522
 0.000182006048  0.999619405229

8410000

不过，我们可以通过一些简单的方差减少技术改善正态分布的前两阶矩的偏差：

第一种方法是采用 对偶变量，也就是先生成所需数量一半的样本，另一半则通过取已有样本的相反数获得，这样可以保证平均数为0，但是不能减小方差。
第二种方法是采用 矩匹配(moment matching) 方法，即首先生成所有样本，然后再对样本进行标准化（减去样本均值，除以样本标准差）。

def gen_sn(M, I, method=None):
    ''' Function to generate random numbers for simulation.
    Parameters
    ==========
    M: int
     number of time intervals for discretization
    I: int
    method: available variance reduction techniques-'anti_paths' and 'mo_match'
    '''
    if method is 'anti_paths':
        sn = npr.standard_normal((M + 1, int(I / 2)))
        sn = np.concatenate((sn, -sn), axis=1)
    else:
        sn = npr.standard_normal((M + 1, I))
        if method is 'mo_match':
            sn = (sn - sn.mean()) / sn.std()
    return sn

sn = gen_sn(0, 10000,'mo_match')
print(sn.mean(), sn.std())

2.966515921798418e-17 0.9999999999999999

仿真模拟（Simulation）

蒙特卡洛模拟（Monte Carlo simulation, MCS）大概是金融中最重要的数值方法之一。之所以能够被广泛使用，主要得益于它是估计数学表达式（例如积分），尤其是在金融衍生品估值方面最为灵活有效的方法。当然，这种灵活性也不是毫无代价的：蒙特卡洛模拟常常为了估计仅仅一个值就要进行成千上万次的复杂运算，导致其计算复杂性较高。使用NumPy进行矢量化是一种在Python中实现蒙特卡罗模拟算法的自然，简洁且高效的方法。但是，使用NumPy矢量化通常会带来更大的内存占用。后面会介绍如何在节约内存的情况下提高运算速度。

随机变量

布莱克-舒尔斯模型（英语：Black-Scholes Model），简称BS模型，又称布莱克-舒尔斯-墨顿模型（Black–Scholes–Merton model），是一种为期权或权证等金融衍生工具定价的数学模型，由美国经济学家迈伦·舒尔斯与费雪·布莱克首先提出，并由罗伯特·C·墨顿修改模型于有派发股利时亦可使用而更完善。（来源：https://en.wikipedia.org/wiki/Black–Scholes_model ）

我们以期权定价中的称布莱克-舒尔斯-墨顿（BSM）模型（Black–Scholes–Merton model）为例：

S_T = S_0\exp((r-\frac{1}{2}\sigma^2)T+\sigma \sqrt{T}\zeta)

其中：

$S_T$ – T时期的指数水平；
$S_0$ – 现在的指数水平；
$r$ – 无风险利率；
$\sigma$ – S的内在波动性；
$\zeta$ – 正态分布随机变量。

对于BSM模型我们用Python进行模拟，其中只有 $\zeta$ 是随机变量，模拟结果如下：

S0 = 100
r = 0.05
sigma = 0.25
T = 2.0
# simulation loops
I = 10000
ST = S0*np.exp((r - 0.5*sigma**2) * T +
               sigma * math.sqrt(T) * npr.standard_normal(I))

plt.hist(ST, bins=50)
plt.xlabel('index level')
plt.ylabel('frequency')
plt.show()

output_17_0

BSM模型有5个重要假设，其中一条就是，金融资产价格服从对数正态分布，即金融资产的对数收益率服从正态分布。

根据BSM模型可知，BSM模型背后的原理是几何布朗运动，所以我们的模拟实际上就是通过正态分布得到一个几何布朗运动的截面。而且，我们知道S_T服从对数正态分布（从上面的表达式可以看出来），那么我们可以直接用对数正态分布进行模拟：

ST_2 = S0 * npr.lognormal((r - 0.5 * sigma ** 2) * T,  # the mean
                          sigma * math.sqrt(T),  # the standard deviation
                          size=I)   

plt.hist(ST, bins=50)
plt.hist(ST_2,bins=50, alpha=0.5)
plt.xlabel('index level')
plt.ylabel('frequency')
plt.legend(['by normal', 'by lognormal'])
plt.show()

output_19_0

比较两个分布的直方图，可以发现，两个计算结果非常相近，我们可以使用KS检验判断两组数据是否属于不同的分布，但是简单期间，我们先看一下两者的统计量：

import scipy.stats as scs

def print_statistics(a1, a2):
    ''' Prints selected statistics.

    Parameters
    ==========
    a1, a2: ndarray objects
     results objects from simulation
    '''
    sta1 = scs.describe(a1)  
    sta2 = scs.describe(a2)  
    print('%14s %14s %14s' %
     ('statistic', 'data set 1', 'data set 2'))
    print(45 * "-")
    print('%14s %14.3f %14.3f' % ('size', sta1[0], sta2[0]))
    print('%14s %14.3f %14.3f' % ('min', sta1[1][0], sta2[1][0]))
    print('%14s %14.3f %14.3f' % ('max', sta1[1][1], sta2[1][1]))
    print('%14s %14.3f %14.3f' % ('mean', sta1[2], sta2[2]))
    print('%14s %14.3f %14.3f' % ('std', np.sqrt(sta1[3]),
                               np.sqrt(sta2[3])))
    print('%14s %14.3f %14.3f' % ('skew', sta1[4], sta2[4]))
    print('%14s %14.3f %14.3f' % ('kurtosis', sta1[5], sta2[5]))

1	print_statistics(ST, ST_2)

     statistic     data set 1     data set 2
---------------------------------------------
          size      10000.000      10000.000
           min         26.810         25.893
           max        391.321        401.891
          mean        110.759        110.422
           std         41.112         40.731
          skew          1.231          1.174
      kurtosis          2.657          2.512

当然两者还是存在一定的差异，可以认为是采样误差引起的。此外，在实际应用中，还有一类误差是离散化误差，即对连续随机过程进行离散化采样引入的误差。

随机过程

刚上大学的时候，看过著名“影片”——《电子系的复兴》，很多台词都非常经典，其中一句“随机过程随机过，量子力学量力学”，讲的就是随机过程。我先后在本科和研究生一共学了两遍随机过程，遗憾的是还是处于随机理解的状态。

简单地说，随机过程就是随机变量随时间变化的过程，也就是在样本产生中引入了先后的概念，那么不同时刻产生的样本也就有可能存在某种关系。如果随机过程满足的关系是—— 当前时刻的样本取值只与上一时刻的样本值有关，而与更早的样本值无关 ——马尔科夫性（Markov property）（此时为1阶）也称为无记忆性（memoryless），那么称为马尔科夫过程。马氏过程在金融分析中得到广泛应用。

随机微分方程

随机微分方程（stochastic differential equation, SDE）是一个微分方程，其中一个或多个项是随机过程，从而产生一个也是随机过程的解。SDE用于模拟各种现象，例如不稳定的股票价格或受热波动影响的物理系统。通常，SDE包含表示随布朗运动或维纳过程的导数计算的随机白噪声的变量。但是，其他类型的随机行为也是可能的，例如跳转过程。

随机微分方程是微分方程的扩展。一般微分方程的对象为可导函数，并以其建立等式。然而，随机过程函数本身的导数不可定义，所以一般解微分方程的概念不适用于随机微分方程。随机微分方程多用于对一些多样化现象进行建模，比如不停变动的股票价格，部分物理现象如热扰动等。

随机微分方程的概念最早以布朗运动的形式，由阿尔伯特·爱因斯坦在《热的分子运动论所要求的静液体中悬浮粒子的运动》论文中提出。这项研究随后由保罗·朗之万继续。此后伊藤清和鲁斯兰斯特拉托诺维奇完善了随机微分方程的数学基础，使得这门领域更加的科学严谨。

一般而言，随机微分方程的解是一随机过程函数，但解方程需要先定义随机过程函数的微分。随机微分方程特别是随机偏微分方程的数值解是一个相对而言的年轻领域。几乎所有用于求解常微分方程的算法对于SDE都非常有效，具有非常差的数值收敛性。常用的方法包括Euler-Maruyama方法（我们后面会用到），Milstein方法和Runge-Kutta方法等

几何布朗运动

几何布朗运动（Geometric Brownian motion, GBM），也叫做指数布朗运动，是连续时间情况下的随机过程，其中随机变量的对数遵循布朗运动，也称为维纳过程。几何布朗运动在金融中的应用之一，就是在前面我们提到的BSM模型中模仿股票价格。

对于满足以下随机微分方程（stochastic differential equation, SDE）的随机过程 $S_t$ ，认为其遵循几何布朗运动：

d{S_t} = \mu S_t dt + \delta S_t d{Z_t} \tag{1}

其中， $Z_t$ 是一个维纳过程（布朗运动）， $\mu$ （百分比drift）和 $\delta$ （百分比volatility）为常量。

公式 $(1)$ 可以通过欧拉-丸山法精确地离散化：

S_t = S_{t-\Delta t}\exp ((r-\frac{1}{2}\sigma^2)\Delta t+\sigma \sqrt{\Delta t}\zeta_t) \tag{2}

其中， $\Delta t$ 是离散化间隔， $\zeta_t$ 是服从正态分布的随机变量。

我们前面采用静态仿真的方法生成了随机变量的样本，那么在这里我们引入时间顺序的概念，按照随机过程的概念进行动态模拟：

# simulation loops
I = 10000
# simulation interval
M = 50
dt = T / M
S = np.zeros((M+1, I))
S[0] = S0
for t in range(1, M+1):
    S[t] = S[t-1] * np.exp((r - 0.5 * sigma ** 2) * dt +
                           sigma * math.sqrt(dt) * npr.standard_normal(I))
ST_3 = S[-1]
plt.hist(ST, bins=50, alpha=0.5)
plt.hist(ST_2,bins=50, alpha=0.5)
plt.hist(ST_3,bins=50, alpha=0.5)
plt.xlabel('index level')
plt.ylabel('frequency')
plt.legend(['by normal', 'by lognormal', 'by GBM'])
print_statistics(ST_3, ST_2)

     statistic     data set 1     data set 2
---------------------------------------------
          size      10000.000      10000.000
           min         22.313         25.893
           max        352.963        401.891
          mean        110.057        110.422
           std         39.967         40.731
          skew          1.078          1.174
      kurtosis          1.972          2.512

output_30_1

使用动态模拟的好处是，我们可以直接根据离散化的过程画出一条路径（每一次蒙特卡洛对应一条路径），而且还可以对美式期权或者百慕大期权进行定价（因为他们的行权和收益是与路径有关的），也就是得到了一个随时间推移的全景图：

1
2
3

plt.plot(S[:, :10], lw=1.5)
plt.xlabel('time')
plt.ylabel('index level')

output_32_1

平方根扩散（Square-root diffusion）

均值回归过程是金融数学中另外一类重要的随机过程，可以用于模拟诸如短期利率或波动率过程。其中一种流行且广泛使用的模型是由Cox，Ingersoll和Ross（1985）提出的平方根扩散模型(Cox–Ingersoll–Ross model, CIR model) ，其随机微分方程如下：

d{x_t}=\kappa(\theta-x_t)dt+\delta\sqrt{x_t}d{Z_t} \tag{3}

其中，

$x_t$ – t时期的指数水平；
$\kappa$ – 均值回归银子；
$\theta$ – 该过程的长期均值；
$\sigma$ – 波动性常数；
$Z_t$ – 标准布朗运动。

众所周知， $x_t$ 服从卡方分布。然而，如前所述，许多金融模型可以通过使用正态分布（即所谓的欧拉-丸山法）来离散化和近似。但是，虽然欧拉法对于几何布朗运动是精确的，它对于大多数其他随机过程是有偏差的。即便如此，而且即使存在可用的精确解（稍后给出平方根扩散的精确解），还是可能会使用欧拉法。例如，下面就是一种欧拉法的离散化（此外还有其他形式）：

\begin{aligned} \tilde{x}_t &= \tilde{x}_S + \kappa(\delta - \tilde{x}_S^+) \Delta t + \delta\sqrt{\tilde{x}_S^+}\sqrt{\Delta t}d{Z_t}\\ x_t &= \tilde{x}_t^+ \end{aligned} \tag{4}

其中，

$s=t-\Delta t$
$x^+\equiv\max(x,0)$

平方根扩散具有使用方便、符合实际（其值严格为证）的特点，但是由于欧拉法本身无法排除负值，所以在模拟的随机过程中加以限制。那么Python的仿真程序如下：

x0 = 0.15  # 初始利率
kappa = 3.0  # 回归系数
theta = 0.02  # 长期均值
sigma = 0.1  # 波动性
I = 10000
M = 100
dt = T / M

def srd_euler():
    xh = np.zeros((M + 1, I))
    x = np.zeros_like(xh)
    xh[0] = x0
    x[0] = x0
    for t in range(1, M + 1):
        xh[t] = (xh[t - 1] +
                 kappa * (theta - np.maximum(xh[t - 1], 0)) * dt +
                 sigma * np.sqrt(np.maximum(xh[t - 1], 0)) *
                 math.sqrt(dt) * npr.standard_normal(I))  
    x = np.maximum(xh, 0)
    return x

x1 = srd_euler()

plt.hist(x1[-1], bins=50)
plt.xlabel('value')
plt.ylabel('frequency')
plt.show()

output_35_0

既然叫均值回归过程，我们通过其随时间推移的动态路径观察是否存在，可以看到这个趋势是非常明显的：

plt.plot(x1[:, :10], lw=1.5)
plt.xlabel('time')
plt.ylabel('index level')
plt.show()

output_37_0

如果我们使用平方根扩散随机微分方程离散化的精确解：

x_t = \frac{\delta^2(1-e^{-\kappa\Delta t})}{4\kappa}\chi'^2_d \left(\frac{4\kappa e^{-\kappa\Delta t}}{\delta^2(1-e^{-\kappa\Delta t})}x_s \right)

比较两个的计算结果和运算时间，我们可以发现，采用欧拉法近似与精确解的结果十分接近，但是速度可以提高接近一倍。

def srd_exact():
    x = np.zeros((M + 1, I))
    x[0] = x0
    for t in range(1, M + 1):
        df = 4 * theta * kappa / sigma ** 2  
        c = (sigma ** 2 * (1 - np.exp(-kappa * dt))) / (4 * kappa)  
        nc = np.exp(-kappa * dt) / c * x[t - 1]  
        # Exact discretization scheme, making use of npr.noncentral_chisquare().
        x[t] = c * npr.noncentral_chisquare(df, nc, size=I)
    return x

x2 = srd_exact()

I = 100000
%time x1 = srd_euler()
%time x2 = srd_exact()

print_statistics(x1[-1], x2[-1])

Wall time: 662 ms
Wall time: 1.02 s
     statistic     data set 1     data set 2
---------------------------------------------
          size     100000.000     100000.000
           min          0.003          0.004
           max          0.062          0.053
          mean          0.020          0.020
           std          0.006          0.006
          skew          0.580          0.576
      kurtosis          0.575          0.479

随机波动（Stochastic volatility）

BS模型假设波动性恒定，但是在实际情况中波动性既不是常数，也不是确定性的——它是随机的。因此，在20世纪90年代早期，随机波动率模型的引入，在金融建模方面取得了重大进展，其中最受欢迎的模型之一是Heston(1993)提出的，模型如下：

\begin{aligned} dS_t &= rS_td + \sqrt{\nu_t}S_tdZ^1_t\\ d\nu_t &= \kappa_\nu(\theta_\nu - \nu_t)dt + \delta_\nu\sqrt{\nu_t}dZ^2_t\\ dZ^1_tdZ^2_t &= \rho \end{aligned} \tag{5}

根据前面的几何布朗运动和平方根扩散可以推断 $(5)$ 中变量和参数的含义：参数 $\rho$ 表示两个标准布朗运动之间的瞬时相关性。这使我们能够解释国外金融市场杠杆效应的现象：波动性在压力时期（市场下跌）上升，在牛市（市场上涨）时下降（国内市场似乎表现为反杠杆效应）。

考虑以下模型的参数化。为了解释两个随机过程之间的相关性，需要确定相关矩阵的Cholesky分解：

S0 = 100
r = 0.05
v0 = 0.1 # 初始波动性
kappa = 3.0
theta = 0.25
sigma = 0.1
rho = 0.6 # 两个布朗运动的相关性
T = 2.0

corr_mat = np.zeros((2,2))
corr_mat[0,:] = [1.0, rho]
corr_mat[1,:] = [rho, 1.0]
cho_mat = np.linalg.cholesky(corr_mat)
cho_mat

array([[1. , 0. ],
       [0.6, 0.8]])

在随机过程模拟开始之前，生成两个过程的整组随机数：一个作为指数过程，一个作为波动率过程。首先模拟波动率过程，通过平方根扩散采用欧拉法建模，通过Cholesky矩阵引入相关性：

M = 50
I = 10000
dt = T / M
ran_num = npr.standard_normal((2, M+1, I))

v = np.zeros_like(ran_num[0])
vh = np.zeros_like(v)
S = np.zeros_like(ran_num[0])

v[0] = v0
vh[0] = v0
S[0] = S0

for t in range(1, M+1):
    ran = np.dot(cho_mat, ran_num[:, t, :])
    vh[t] = (vh[t - 1] +
             kappa * (theta - np.maximum(vh[t - 1], 0)) * dt +
             sigma * np.sqrt(np.maximum(vh[t - 1], 0)) *
             math.sqrt(dt) * ran[1])
    v[t] = np.maximum(vh[t], 0)
    # 用欧拉法模拟几何布朗运动
    S[t] = S[t-1] * np.exp((r - 0.5 * v[t]) * dt +
                           np.sqrt(v[t]) * ran[0] * np.sqrt(dt))


fig, (ax1, ax2) = plt.subplots(1,2, figsize=(15,6))
# 期末水平
ax1.hist(S[-1], bins=50)
ax1.set_xlabel('index level')
ax1.set_ylabel('frequency')
ax2.hist(v[-1], bins=50)
ax2.set_xlabel('volatility')
plt.show()

output_44_0

这说明了使用欧拉法进行平方根扩散的另一个优点：由于只需要标准的正态分布随机数，因此容易且一致地考虑相关性。而使用混合方法（即指数过程用欧拉法，波动率过程使用非中心卡方的精确方法）没有简单的方法来实现相同的方法。

1	print_statistics(S[-1], v[-1])

     statistic     data set 1     data set 2
---------------------------------------------
          size      10000.000      10000.000
           min          5.789          0.173
           max       1120.679          0.330
          mean        116.469          0.250
           std         90.146          0.021
          skew          2.802          0.125
      kurtosis         14.669          0.057

查看两个过程的模拟路径,我们看到在随机波动率的情况下，期末指数的最大值远远高于固定波动率（其他条件不变）情况下的几何布朗运动。

fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True,
                               figsize=(15, 12))
ax1.plot(S[:, :10], lw=1.5)
ax1.set_ylabel('index level')
ax2.plot(v[:, :10], lw=1.5)
ax2.set_xlabel('time')
ax2.set_ylabel('volatility')
plt.show()

output_48_0

跳跃扩散（Jump diffusion）

扩散过程，是指具有连续样本轨道（trajectory）的马尔科夫过程，在时间、状态上均连续。相对应地，跳跃过程是指样本轨道存在跳跃点的随机过程，例如计数过程，通常建模为简单或复合泊松过程。那么，跳跃扩散过程则是同时具有跳跃和离散的随机过程，在凝聚态物理、模式识别、期权定价等方面都有重要应用。

除了前面提到的随机波动率和杠杆效应外，还有一个经验现象时资产价格或者波动性的跳跃。在期权定价中，跳跃扩散模型是混合模型的一种形式，混合了跳跃过程和扩散过程。 Robert C. Merton（也就是BSM模型中的M）引入了跳跃扩散模型作为跳跃模型的扩展。由于它们的计算易处理性，基本仿射跳跃扩散的特殊情况在一些信用风险和短期模型中很受欢迎。

在风险中性的情况下，Merton跳跃扩散过程的随机微分方程如下：

dS_t = (r-r_J)S_tdt + \delta S_tdZ_t + J_tS_tdN_t \tag{6}

其中，

$S_t$ – t时期的资产价格
$r$ – 无风险利率
$r_J \equiv \lambda e^{\mu_j+\sigma^2-1}$ – 为保持风险中性对跳跃的漂移校正
$\delta$ – 波动率常数
$Z_t$ – 标准布朗运动
$J_t$ – t时刻的跳跃（服从lognormal）
$N_t$ – 抵达率为 $\lambda$ 的泊松过程

其欧拉法数值解为：

S_t = S_{t-\Delta t}\left(e^{(r-r_J-\Delta^2/2) \sigma t + \sigma \sqrt{\Delta t}\zeta^1_t} + \left(e^{\mu_J+\delta\zeta^2_t}-1 \right)y_t\right)

其中， $\zeta^n_t$ 服从标准正态分布， $y_t$ 服从参数为 $\lambda$ 的泊松分布。

S0 = 100
r = 0.05
sigma = 0.2
lamb = 0.75 # 跳跃到达率
mu = -0.6 # 平均跳跃大小
delta = 0.25  # 跳跃波动性
# 漂移校正
rj = lamb * (math.exp(mu + 0.5 * delta ** 2) - 1)

T = 1.0
M = 50
I = 10000
dt = T / M

S = np.zeros((M+1, I))
S[0] = S0
sn1 = npr.standard_normal((M+1, I))
sn2 = npr.standard_normal((M+1, I))
poi = npr.poisson(lamb * dt, (M+1, I))

for t in range(1, M+1, 1):
    S[t] = S[t - 1] * (np.exp((r - rj - 0.5 * sigma ** 2) * dt +
                              sigma * math.sqrt(dt) * sn1[t])+
                       (np.exp(mu + delta * sn2[t]) - 1) * poi[t])
    S[t] = np.maximum(S[t], 0)

plt.hist(S[-1], bins=50)
plt.xlabel('value')
plt.ylabel('frequency')
plt.show()

output_51_0

plt.plot(S[:, :10], lw=1.5)
plt.xlabel('time')
plt.ylabel('index level')
plt.title('Dynamically simulated jump diffusion process paths')
plt.show()

output_52_0

小结

本文简单介绍了金融分析中常用到的随机数，尤其是正态分布及校正技巧；然后介绍了金融衍生品定价中用到的几种随机过程（几何布朗、平方根扩散、随机波动率、跳跃扩散）的离散化方法和蒙特卡洛模拟方法。

前言

随机数