DOI: 10.12677/AAM.2020.98158, PDF, HTML, XML,  被引量    科研立项经费支持
作者: 贾翔宇：唐山学院基础教学部，河北 唐山
关键词: 跳–扩散模型；蒙特卡罗法；加权极大似然法；Jump-Diffusion Models； Monte Carlo Method； Weighted Maximum Likelihood Method

文章引用：贾翔宇. 跳–扩散模型参数校准的极大似然法[J]. 应用数学进展, 2020, 9(8): 1338-1345. https://doi.org/10.12677/AAM.2020.98158

1. 引言

随着计量经济理论的发展，各种参数估计方法不断涌现出来，极大似然法是求估计的一种方法，这一方法首先由德国数学家Gauss于1821年提出，但未得到足够的重视。英国的统计学家Fisher在1922年的论文“On the mathematical foundations of theoretical statistics, reprinted in Contributions to Mathematical Statistics” [1] 中再次提出了极大似然法的思想，并且首先探讨了这种方法的一些性质，如一致性，渐进正态性，渐进有效性和不变性等。极大似然估计这一名称也是Fisher命名的，这是一种目前仍然得到广泛应用的方法。

人们常用极大似然法来估计跳–扩散模型，这也是统计学里很经典的方法。但在实现操作过程中，极大似然法会有很大难度，因为只有少数扩散过程的似然函数具有封闭形式，如Black-Scholes-Merton (1973) [2]、Vasicek (1977) [3]、CIR (1985) [4] 等模型，而只有具有封闭形式的似然函数才可用来估计，这就引发了扩散模型的转移密度函数的研究。求解转移密度函数的方法一般分为两类，一类是直接将扩散模型离散化并假设模型的参差符合正态分布，在这种假设条件下，模型的参数估计可以直接通过求解非线性计量经济模型参数的方法求得，如Euler法、Nowman法等。另一类是从扩散模型满足的偏微分方程入手，通过数值方法求解偏微分方程来获得转移密度函数的近似解，如SMLE法、Crank-Nicoslon法等。

本文中，我们讨论了跳–扩散模型参数校准问题的极大似然法。首先模拟股票价格路径的生成，然后对生成的路径建立对数似然函数，并构造权重函数，利用加权极大似然法校准模型中的波动率，最后数值模拟试验表明了该算法的可行性。

2. 模型的建立

在风险中性测度下，假设股票价格 $S_t$ 满足如下随机微分方程

$\frac{\text{d}{S}_t}{S_t}=\left(\mu -\lambda \kappa \right)\text{d}t+\sigma \text{d}{W}_t+\left(\eta -1\right)\text{d}{Q}_t,$ (1)

其中 $\mu$ 是漂移率， $\sigma$ 是股票的波动率， $W_t$ 是标准布朗运动， $\text{d}{Q}_t$ 是一个强度为 $\lambda$ 的泊松过程， $\eta -1$ 是股票发生跳跃时的相对跳跃高度， $\kappa =E\left[\eta -1\right]$，其中 $E[\cdot ]$ 表示期望算子。

在Merton跳–扩散模型下， $\eta$ 服从对数正态分布，若用 $g\left(\eta \right)$ 表示 $\eta$ 的概率密度函数，则有

$g\left(\eta \right)=\frac{1}{\sqrt{2\pi }{\sigma }_J\eta }\exp \left(-\frac{{\left(\log \eta -\nu \right)}^2}{2{\sigma }_J^2}\right),$ (2)

其中 $\kappa =E\left[\eta -1\right]=\exp \left(\nu +\frac{{\sigma }_J^2}{2}\right)-1$， $\nu$ 和 ${\sigma }_J$ 表示跳跃的均值和方差。

而在Kou跳–扩散模型下， $g\left(\eta \right)$ 满足对数双指数分布

$g\left(\eta \right)=p{\alpha }_1{\text{e}}^{-{\alpha }_1\eta }H\left(\eta \right)+q{\alpha }_2{\text{e}}^{{\alpha }_2\eta }H\left(-\eta \right),$ (3)

其中 ${\alpha }_1>1$， ${\alpha }_2,p,q>0$， $p+q=1$， $H\left(\eta \right)$ 是赫维赛德阶跃函数。

随机微分方程(1)的解为

$S\left(t\right)=S\left(0\right){\text{e}}^{\left(a-\frac{1}{2}{\sigma }^2\right)t+\sigma W\left(t\right)}{\displaystyle \underset{j=1}{\overset{Q\left(t\right)}{\prod }}{\eta }_j},$ (4)

其中令 $a=\mu -\lambda \kappa$。

3. 模拟样本路径的蒙特卡洛方法

固定时间 $0=t_0<t_1<\cdots <t_n$，并令 $\Delta t_i=t_{i+1}-t_i,i=0,1,\cdots ,n-1$。Glasserman (2004) [5] 给出了模拟带跳的几何布朗扩散过程的方法。

假设Q是一个泊松过程， ${\eta }_1,{\eta }_2,\cdots$ 独立同分布， $Q,W$ 和 $\left\{{\eta }_1,{\eta }_2,\cdots \right\}$ 相互独立。为了模拟 $S_t$ 在时刻 $t_1,\cdots ,t_n$ 的值，我们推广(4)式为

$S\left(t_{i+1}\right)=S\left(t_i\right){\text{e}}^{\left(a-\frac{1}{2}{\sigma }^2\right)\Delta t_i+\sigma \left[W\left(t_{i+1}\right)-W\left(t_i\right)\right]}{\displaystyle \underset{j=Q\left(t_i\right)+1}{\overset{Q\left(t_{i+1}\right)}{\prod }}{\eta }_j},$ (5)

若 $Q\left(t_{i+1}\right)=Q\left(t_i\right)$，则乘积项 ${\displaystyle \underset{j=Q\left(t_i\right)+1}{\overset{Q\left(t_{i+1}\right)}{\prod }}{\eta }_j}=1$，即退化为几何布朗运动。令 $x\left(t\right)=\log S\left(t\right)$，则有

$x_{i+1}=x_i+\left(a-\frac{1}{2}{\sigma }^2\right)\Delta t_i+\sigma \left(W_{i+1}-W_i\right)+{\displaystyle \underset{j=Q_i+1}{\overset{Q_{i+1}}{\sum }}\log {\eta }_j},$ (6)

其中 $x_i=x\left(t_i\right)$， $W_i=W\left(t_i\right)$， $Q_i=Q\left(t_i\right)$。

Merton跳–扩散模型，一般模拟(6)式从时刻 $t_i$ 到 $t_{i+1}$ 的方法有以下步骤：

1) 生成服从标准正态分布的随机变量 $Z\left(Z~N\left(0,1\right)\right)$ ；

2) 生成服从参数为 $\lambda \Delta t_i$ 的泊松分布的随机变量Q；如果 $Q=0$，则令 $M=0$，并跳到第4步；

3) 按照 $\eta$ 的分布生成 $\log {\eta }_1,\log {\eta }_2,\cdots ,\log {\eta }_Q$，令 $M=\log {\eta }_1+\log {\eta }_2+\cdots +\log {\eta }_Q$ ；

4) 令

$x_{i+1}=x_i+\left(a-\frac{1}{2}{\sigma }^2\right)\Delta t_i+\sigma \sqrt{\Delta t_i}Z+M.$ (7)

这种方法充分依赖于泊松过程的两个性质：一是增量 $Q_{i+1}-Q_i$ 服从均值为 $\lambda \Delta t$ 的泊松过程；二是具有独立增量，即 $Q_{i+1}-Q_i$ 与 $Q_i$ 独立。作进一步的假设，设 ${\eta }_j$ 服从对数正态分布 $LN\left(\nu ,{\sigma }_J^2\right)$，那么 $\log {\eta }_j~N\left(\nu ,{\sigma }_J^2\right)$，则有

${\displaystyle \underset{j=1}{\overset{n}{\sum }}\log {\eta }_j}~N\left(n\nu ,n{\sigma }_J^2\right)=n\nu +{\sigma }_J\sqrt{n}N\left(0,1\right).$ (8)

则可以替换第3步为

3’) 生成 $Z_2~N\left(0,1\right)$，设定 $M=\nu Q+{\sigma }_J\sqrt{Q}{Z}_2$。

参考Ait-Sahalia (2002) [6]，我们模拟Merton跳–扩散模型的路径生成时，选取参数为 $a=0.2,\sigma =0.3,\lambda =1/3,\nu =0,{\sigma }_J=0.2,x_0=0$，我们对该过程进行一次样本路径的模拟，并令 $n=1000,\Delta t=1/52$，见图1。

而对于Kou跳–扩散模型， $\log {\eta }_j$ 满足伽马分布(实际上是指数分布)，且 $\log {\eta }_j$ 大于零的概率为p，小于零的概率为q。在这种情况下，当泊松随机变量 $Q=n$ 时， $\log {\eta }_j$ 取正值的个数满足参数为n和p的二项分布。因此第3步应改为

3a’) 生成 $K~b\left(Q,p\right)$，其中 $b\left(Q,p\right)$ 表示以n和p为参数的二项分布；

3b’) 生成 $R_1~Ga\left(K{\alpha }_1,{\alpha }_2\right)$， $R_2~Ga\left(\left(Q-K\right){\alpha }_1,{\alpha }_2\right)$，令 $M=R_1-R_2$。

在第3b’步中，若 $K=0$，则令 $R_1=0$ ；若 $K=Q$，则令 $R_2=0$。在第3a’步中，当 $Q=n$ 时，二项分布K的密度函数为

$P\left(K=k\right)=\frac{n!}{k!\left(n-k\right)!}{p}^k{q}^{n-k},k=0,1,\cdots ,n$ (9)

4. 跳–扩散模型参数校准的加权极大似然法

Press (1967) [7] 已经证明，对于Merton跳–扩散模型， $x_t$ 的转移概率密度函数是一个独立的泊松过程和一个正态概率密度函数的合成，即为

$f_M\left(x_{i+1}|x_i,\sigma \right)={\displaystyle \underset{j=0}{\overset{\infty }{\sum }}\frac{{\left(\lambda \Delta t\right)}^j{\text{e}}^{-\lambda \Delta t}}{j!}}\cdot \frac{1}{\sqrt{2\pi \left({\sigma }^2\Delta t+j{\sigma }_J^2\right)}}\cdot \exp \left\{-\frac{{\left(x_{i+1}-x_i-a\Delta t+\frac{1}{2}{\sigma }^2\Delta t-j\nu \right)}^2}{2\left({\sigma }^2\Delta t+j{\sigma }_J^2\right)}\right\},$ (10)

其中j表示发生跳跃的次数。对于 $\left\{x_i\right\},i=1,\cdots ,n$，可得对数似然函数为

$L_{Merton}\left(\sigma \right)=\log f_M\left(x_0,\sigma \right)+{\displaystyle \underset{i=0}{\overset{n-1}{\sum }}\log f_M}\left(x_{i+1}|x_i,\sigma \right).$ (11)

而对于Kou跳–扩散模型， $x_t$ 的转移概率密度函数为

$\begin{array}{c}{f}_K\left(x_{i+1}|x_i,\sigma \right)={\displaystyle \underset{j=0}{\overset{L}{\sum }}\frac{L!}{j!\left(L-j\right)!}{p}^j{q}^{L-j}\frac{1}{\sqrt{2\pi \left({\sigma }^2\Delta t-j\frac{\left(Q-K\right){\alpha }_1}{{\alpha }_2^2}+j\frac{K{\alpha }_1}{{\alpha }_2^2}\right)}}}\\ \cdot \exp \left\{-\frac{{\left(x_{i+1}-x_i-a\Delta t+\frac{1}{2}{\sigma }^2\Delta t+j\frac{\left(Q-K\right){\alpha }_1}{{\alpha }_2}-j\frac{K{\alpha }_1}{{\alpha }_2}\right)}^2}{2\left({\sigma }^2\Delta t-j\frac{\left(Q-K\right){\alpha }_1}{{\alpha }_2^2}+j\frac{K{\alpha }_1}{{\alpha }_2^2}\right)}\right\},\end{array}$ (12)

其中 $L=0,1,\cdots$。对于 $\left\{x_i\right\},i=1,\cdots ,n$，可得对数似然函数为

$L_{Kou}\left(\sigma \right)=\log f_K\left(x_0,\sigma \right)+{\displaystyle \underset{i=0}{\overset{n-1}{\sum }}\log f_K}\left(x_{i+1}|x_i,\sigma \right).$ (13)

由于当观察次数足够多时，似然函数中第一项 $\log f\left(x_0,\sigma \right)$ 的权重很小，因此在实际的估计过程中往往省略。求似然函数 $L\left(\sigma \right)$ 的最大化，即可得到参数的估计值。对于j的取值，Ball和Torous (1985) [8] 发现：在 $j=10$ 处截断可以得到精确的极大似然估计量。

由于在实际情况中，我们获取的时间序列数据都包含着观测误差和市场噪音，因此我们假设股票的市场观测值 $S_i^M$ 与模型所得价格之间存在一定的高斯误差，即

$S\left(t_i,\sigma \right)=S_i^M+u_i,i=0,1,\dots ,n$ (14)

其中 $S\left(t_i,\sigma \right)$ 是在参数 $\sigma$ 已知的条件下，通过模拟路径生成的值， $u_i~N\left(0,{\varphi }_i\right)$， ${\varphi }_i$ 是已知的。

因此，为了使得估计得到的参数能更好地拟合市场的动态过程，构造如下权重函数

$g\left(x_{i+1}|x_i,\sigma \right)=\frac{1}{2\pi {\varphi }_{i+1}}\cdot \exp \left\{-\frac{{\left(S_{i+1}-S_{i+1}^M\right)}^2}{2{\varphi }_{i+1}}\right\},i=0,1,\cdots ,n-1$ (15)

综合(11)和(15)式，并省略常数项，我们得到Merton跳–扩散模型的加权对数似然函数为

$\begin{array}{l}{F}_{Merton}\left(\sigma \right)={\displaystyle \underset{i=0}{\overset{n-1}{\sum }}\log f_M\left(x_{i+1}|x_i,\sigma \right)\cdot g\left(x_{i+1}|x_i,\sigma \right)}\\ ={\displaystyle \underset{i=0}{\overset{n-1}{\sum }}\log \left\{{\displaystyle \underset{j=0}{\overset{\infty }{\sum }}\frac{{\left(\lambda \Delta t\right)}^j{\text{e}}^{-\lambda \Delta t}}{j!}\cdot \frac{1}{\sqrt{2\pi \left({\sigma }^2\Delta t+j{\sigma }_J^2\right)}}\cdot \exp \left\{-\frac{{\left(x_{i+1}-x_i-a\Delta t+\frac{1}{2}{\sigma }^2\Delta t-j\nu \right)}^2}{2\left({\sigma }^2\Delta t+j{\sigma }_J^2\right)}\right\}}\right\}}\\ \cdot \frac{1}{2\pi {\varphi }_{i+1}}\cdot \exp \left\{-\frac{{\left(S_{i+1}-S_{i+1}^M\right)}^2}{2{\varphi }_{i+1}}\right\}.\end{array}$ (16)

综合(13)和(15)式，并省略常数项，我们得到Kou跳–扩散模型的加权对数似然函数为

$\begin{array}{l}{F}_{Kou}\left(\sigma \right)={\displaystyle \underset{i=0}{\overset{n-1}{\sum }}\log f_K\left(x_{i+1}|x_i,\sigma \right)\cdot g\left(x_{i+1}|x_i,\sigma \right)}\\ ={\displaystyle \underset{i=0}{\overset{n-1}{\sum }}\log \left\{{\displaystyle \underset{j=0}{\overset{L}{\sum }}\frac{L!}{j!\left(L-j\right)!}{p}^j{q}^{L-j}\frac{1}{\sqrt{2\pi \left({\sigma }^2\Delta t-j\frac{\left(Q-K\right){\alpha }_1}{{\alpha }_2^2}+j\frac{K{\alpha }_1}{{\alpha }_2^2}\right)}}}\right\}}\\ \cdot \exp \left\{-\frac{{\left(x_{i+1}-x_i-a\Delta t+\frac{1}{2}{\sigma }^2\Delta t+j\frac{\left(Q-K\right){\alpha }_1}{{\alpha }_2}-j\frac{K{\alpha }_1}{{\alpha }_2}\right)}^2}{2\left({\sigma }^2\Delta t-j\frac{\left(Q-K\right){\alpha }_1}{{\alpha }_2^2}+j\frac{K{\alpha }_1}{{\alpha }_2^2}\right)}\right\}\\ \cdot \frac{1}{2\pi {\varphi }_{i+1}}\cdot \exp \left\{-\frac{{\left(S_{i+1}-S_{i+1}^M\right)}^2}{2{\varphi }_{i+1}}\right\}.\end{array}$ (17)

通过求解 $F\left(\sigma \right)$ 的最大化( $F\left(\sigma \right)$ 表示 $F_{Merton}\left(\sigma \right)$ 和 $F_{Kou}\left(\sigma \right)$ )，便可以得到参数 $\sigma$ 的估计，即

${\sigma }^{*}=\arg \max F\left(\sigma \right).$ (18)

从 $F\left(\sigma \right)$ 的表达式可以看出，加入权重函数(15)的目的是对于能够较好拟合市场期权价格的点，赋予的权重大一些，而对于那些与市场期权价格观测值差距较大的点，则赋予的权重要小一些。

5. 数值实验

下面给出算法的具体过程。

首先由第三节的算法步骤生成一条样本路径，即生成 $x_i,i=1,2,\cdots ,n$ 的值，同时算得 $S_i$ 和 $S_i^M$ 的值，然后通过线性搜索的方法求解(18)式，并得到参数 $\hat{\sigma }$ 的估计值。由于生成 $x_i$ 和 $S_i^M$ 时都有随机项的存在，因此我们引入蒙特卡洛方法进行H次模拟，最终得到 $\sigma$ 的平均估计值 $\bar{\sigma }$，并定义均方根误差为

$\text{RMSE}\left(\sigma \right)=\sqrt{{\displaystyle \underset{i=1}{\overset{H}{\sum }}{\left({\hat{\sigma }}_i-{\sigma }_0\right)}^2/H}},$

其中 ${\hat{\sigma }}_i$ 为校准所得到的数值解。

对于Merton跳–扩散模型，我们选取参数 $a=0.2,\lambda =\frac{1}{3},\nu =0,{\sigma }_J=0.2,x_0=0$，选取等距的离散时间点，令 $\Delta t=0.01$，并假设参数 $\sigma$ 的真实值为 ${\sigma }_0=0.3$。

从表1可以看出，当 $\varphi ={10}^{-3}S$ 时，波动率估计值 $\hat{\sigma }$ 的平均值 $\bar{\sigma }$ 在 $n=10$ 时得到的波动率估计值的均值有一定的误差，并且均方根误差(RMSE)较大，而在 $n=50,100,200$ 三种情况下都能得到较好的估计，其中n表示模拟生成路径上点的个数。在这四种情况下，波动率估计值的均方根误差(RMSE)会随着n的增大而减小，并且减小到一定程度后，就只会缓慢减小。这表明随着路径上点的个数的增加，我们估计得到的波动率的值更稳定趋于真实值，但增加到一定程度，效果就会减慢。

表2表明，当 $n=100$ 时，在 $\varphi ={10}^{-2}S,{10}^{-3}S,{10}^{-4}S$ 三种情况下波动率估计值 $\hat{\sigma }$ 的平均值 $\bar{\sigma }$ 与均方根误差的数值结果。实验结果表明，三种情况下得到的均值和均方根误差相差不大，都能得到较好的估计，且均方根误差都较小。

对于Kou跳–扩散模型，我们选取参数 $a=0.2,\lambda =\frac{1}{3},{\alpha }_1=3.0456,{\alpha }_2=3.0775,p=0.3445,x_0=0$，选取等距的离散时间点，令 $\Delta t=0.01$，并假设参数 $\sigma$ 的真实值为 ${\sigma }_0=0.3$。

从表3可以看出，当 $\varphi ={10}^{-3}S$ 时，波动率估计值 $\hat{\sigma }$ 的平均值 $\bar{\sigma }$ 在 $n=10$ 时，得到的估计会有一定的误差，并且均方根误差(RMSE)较大，而当 $n=50,100,200$ 时，能得到波动率 $\sigma$ 的较好估计，且其均方根误差(RMSE)会随着n的增大而减小，并且减小到一定程度后，就只会缓慢减小。这表明随着路径上点的个数的增加，我们估计得到的波动率的值更稳定趋于真实值，但增加到一定程度，效果就会减慢。

表4表明，当 $n=100$ 时，在 $\varphi ={10}^{-2}S,{10}^{-3}S,{10}^{-4}S$ 三种情况下波动率估计值 $\hat{\sigma }$ 的平均值 $\bar{\sigma }$ 与均方根误差的数值结果。实验结果表明，三种情况下得到的均值和均方根误差相差不大，都能得到较好的估计，并且每次的估计值基本都稳定在真值附近。

我们对Merton跳–扩散模型和Kou跳–扩散模型波动率校准进行比较。表5表示在其余参数相同的情况下，当 $n=10,50,100,200$ 时的计算结果。可以看出，当n较小时，两种模型得到的估计值都会有一定的偏差，并且Kou模型得到的结果要比Merton模型得到的更精确且稳定一些，而随着n的增加，两种模型都能得到比较精确且稳定的估计值，且n较大时，Merton模型得到的估计值更稳定一些。对于我们所举的数值算例，若路径上的点较少时，我们使用Kou跳–扩散模型进行波动率的校准，而当路径上的点较多时，我们采用Merton模型进行校准。

6. 结论

本文中，我们讨论了Merton跳–扩散模型和Kou跳–扩散模型参数校准问题的极大似然法。首先用蒙特卡洛罗法模拟股票价格路径的生成，然后找到生成的路径的转移概率密度函数，并建立对数似然函数。令跳–扩散模型的股票市场价格与模型价格之间的误差为高斯误差，并构造权重函数，通过线性搜索的方法求解加权极大似然函数的最大化，从而得到波动率的校准。最后给出数值模拟实验，实验结果表明了该算法的有效性和可行性。

基金项目

由唐山学院博创基金“金融中的反问题及其数值计算”资助。

参考文献

[1]	Fisher, R.A. (1950) On the Mathematical Foundations of Theoretical Statistics, Reprinted in Contributions to Mathematical Statistics. J. Wiley & Sons, New York.
[2]	Black, F. and Scholes, M. (1973) The Pricing of Options and Corporate Liabilities. Journal of Political Economy, 81, 637-659. https://doi.org/10.1086/260062
[3]	Vasicek, O. (1977) An Equilibrium Characterization of the Term Structure. Journal of Financial Economics, 5, 177-188. https://doi.org/10.1016/0304-405X(77)90016-2
[4]	Cox, J.C., Ingersoll, J.E. and Ross, S.A. (1985) A Theory of the Term Structure of Interest Rates. Econometrica, 53, 385-407. https://doi.org/10.2307/1911242
[5]	Glasserman, P. (2004) Monte Carlo Methods in Financial Engineering. Spring Verlag, New York. https://doi.org/10.1007/978-0-387-21617-1
[6]	Ait-Sahalia, Y. (2002) Maximum Likelihood Estimation of Discretely Sampled Diffusions: A Closed-Form Approximation Approach. Econometrica, 70, 223-262. https://doi.org/10.1111/1468-0262.00274
[7]	Press, S.J. (1967) A Compound Events Model for Security Prices. Journal of Business, 40, 317-335. https://doi.org/10.1086/294980
[8]	Ball, C.A. and Torous, W.N. (1985) On Jumps in Common Stock Prices and Their Impact on Call Option Pricing. Journal of Finance, 40, 155-173. https://doi.org/10.1111/j.1540-6261.1985.tb04942.x