傅里叶变换期权定价的数值方法

0.通过傅里叶变换计算BSM 定价的call 期权

Carr/Madan payoff 傅立叶变换法

underlying asset X :

dX=\mu S dt+\sigma S dw\\ ln(X_T) \sim N(x_0+(\mu-\frac{\sigma^2}{2}）T,\sigma^2T) ln(X_T)=x_T 的特征方程为

\Phi_{x_T}(\nu)=exp(i(x_0+(\mu -\frac{\sigma^2}{2})T)\nu -\frac{\sigma^2 \nu ^2}{2}T)

期权价格的 C_T

\Psi _T(v)=\frac{C}{(\alpha+iv)(\alpha+iv+1)}\Phi_{x_T}(v-(\alpha+1)i)\\ C_T(k)=\frac{e^{-\alpha k}}{\pi} \int e^{-ivk} \Psi _T(v) dv\\C =e^{-rT} 如果使用特定的特征方程，就可以得到对应的随机过程下的定价, 这里可以发现Carr/Madan 方法的第一个好处，知道underlying assest的 随机过程或者特征方程，就可以获得期权价格的计算方程。对于BSM 定价假设的随机过程的特征函数： \Phi_{x_T}(\nu)=exp(i(x_0+(\mu -\frac{\sigma^2}{2})T)\nu -\frac{\sigma^2 \nu ^2}{2}T)

1. 数值积分计算

简单的数值逼近积分（Simposon's rule）

\int_a^b f(x)\approx\frac{b-a}{6} (f(a)+4f(\frac{a+b}{2})+f(b))\\ =\frac{\Delta x}{3}[f(a)+4(f1+f3+f5+...)+2(f2+f4+f6 +...)+f(b)]

我们设定 \Delta b=\frac{B}{N}=\eta, v_j=j\eta， for\ \ j=0...N.

\begin{align*} C_T(k)&=\frac{e^{-\alpha k}}{\pi} \int_0^\infty e^{-ivk} \Psi _T(v) dv\\ &\approx \frac{e^{-\alpha k}}{\pi} \int_0^B e^{-ivk} \Psi _T(v) dv\\ &= \frac{e^{-\alpha k}}{\pi} \sum^N_{j=0} e^{-i v_jk}\Psi _T(v_j) *w_j \end{align*}

根据simposon rule

w_j=\frac{\eta}{3}[1,4,2,4,2,4,2,4.....1]\approx \frac{\eta}{3}(3+(-1)^j-\delta_j) \\ \text{for } \begin{equation} \delta_j= \begin{cases} 1, & \text{if}\ j=1 \\ 0, & \text{otherwise} \end{cases} \end{equation} 这是一个单纯的数值积分方法。。 B代表着 对于连续的积分区间 替换成有限的积分区间B= N*\eta 。 \eta 代表采样频率。 由于积分term 内有 e^{iv_jk} ，整个积分项会指数的衰减。所以更大 的B带来的精确度提升是可以忽略不计的。N 越大(采样点越多)， \eta 越小（积分间距），数值解就精度更高。不考虑效率可以十分奢侈的设定 N 极大， \eta 极小。

下面是数值积分方法的算法demo。这里我们使用了 \alpha=1.5,N=400000,B=100 (有关\alpha 的取值会在FFT详细讨论)

def BSM_characteristic_function(v, x0, T, r, sigma):
    cf_value = np.exp(((x0 / T + r - 0.5 * sigma ** 2) * 1j * v
                - 0.5 * sigma ** 2 * v ** 2) * T)
    return cf_value
def BSM_call_characteristic_function(v,alpha, x0, T, r, sigma):
    res=np.exp(-r*T)/((alpha+1j*v)*(alpha+1j*v+1))\
        *BSM_characteristic_function((v-(alpha+1)*1j), x0, T, r, sigma)
    return res
    
def SimpsonW(N,eta):
    delt = np.zeros(N, dtype=np.float)
    delt[0] = 1
    j = np.arange(1, N + 1, 1)
    SimpsonW = eta*(3 + (-1) ** j - delt) / 3
    return SimpsonW
    

def Simposon_numerical_integrate(S0, K, T, r, sigma):
    k = np.log(K)
    x0 = np.log(S0)
    N=400000
    B=100
    eta=B/N
    W=SimpsonW(N,eta)
    
    alpha=1.5
    sumx=0
    for j in range(N):
        v_j=j*eta
        temp=np.exp(-1j*v_j*k)*\
            BSM_call_characteristic_function(v_j,alpha, x0, T, r, sigma)*\
            W[j]            
        sumx+=temp.real

        
    return sumx*np.exp(-alpha*k)/np.pi

以下面的call期权为例子，做一个简单数值解和解析解的比较

S0 = 100.0  # index level
K = 90.0  # option strike
T = 1.0  # maturity date
r = 0.05  # risk-less short rate
sigma = 0.3  # volatility

print BSM_call_value(S0, K, T, r, sigma)
print Simposon_numerical_integrate(S0, K, T, r, sigma)

19.6974420868
19.6974420868

缺点： 慢！

2.FFT 使用快速傅里叶变换加速

我们会把FFT 当作黑箱子使用，然后分析在定价上如何实现，以及使用FFT 定价的优点和缺点。以及参数的问题。

FFT 的黑箱子

我们不讨论FFT 的实现原理，只描述其作用，输入，输出。

如果出现形如

w(m)=\sum^N_{j=0}e^{-\frac{2\pi}{N}j m}x(j) \ \ \text{for } m=0..N

我们可以使用FFT 使得求w(m)的函数值从 O(n^2) 复杂度降低为 O(n\log \ n) ，注意的是我们获得的是整个w函数的数值解，而不是特定m 的函数值。这也意味着我们可以一次求得一组只有k 不同的期权定价。

FFT 的使用

重新回到我们之前的定价公式。

\begin{align*} C_T(k)&=\frac{e^{-\alpha k}}{\pi} \int_0^\infty e^{-ivk} \Psi _T(v) dv\\ &\approx\frac{e^{-\alpha k}}{\pi} \sum^N_{j=0} e^{-i v_jk}\Psi _T(v_j) *w_j \end{align*}

把 C_T(k) 视为k 和期权的函数关系，我们要做的第一步是对k 进行离散化。采样N 个点，采样间隔为 \lambda . 这里采用一个相对合理的下界，我们设定下界 \beta=ln(X_0)-\frac{\lambda N}{2} , 这样可以使得at the money 的期权在定义域的中间 。即 , 积分的采样和之前相似 for j=0...N

\begin{align*} C_T(k_m)&\approx\frac{e^{-\alpha k_m}}{\pi} \sum^N_{j=0} e^{-i v_jk_m}\Psi _T(v_j) *w_j\\ &=\frac{e^{-\alpha k_m}}{\pi} \sum^N_{j=0} e^{-i(j\eta) (\beta+m\lambda)}\Psi _T(v_j) *w_j\\ &=\frac{e^{-\alpha k_m}}{\pi} \sum^N_{j=0} e^{-i(\lambda \eta )(jm)}e^{-i (\beta v_j)}\Psi _T(v_j) *w_j\\ \end{align*}

令x(j)=e^{-i (\beta v_j)}\Psi _T(v_j) *w_j, \lambda \eta=\frac{2\pi}{N} ， 我们就可以使用FFT 了。

代码demo

def BSM_call_value_FFT(S0, K, T, r, sigma):
    k = np.log(K)
    x0 = np.log(S0)
    N =2**10
    alpha=1.5
    
    eta=0.15
    lambda_ = 2 * np.pi / (N *eta)
    beta=x0-lambda_*N/2
    km=np.asarray([beta+i*lambda_ for i in range(N)])
    W=SimpsonW(N,eta)
    v=np.asarray([i*eta for i in range(N)])
    Psi=np.asarray([BSM_call_characteristic_function(vj,alpha, x0, T, r, sigma)  for vj in v])
    FFTFunc=Psi*np.exp(-1j*beta*v)*W
    y=fft(FFTFunc).real
    cT=np.exp(-alpha*km)*y/np.pi
    
    return np.exp(km),cT

使用之前的 期权做定价比较

print BSM_call_value(S0, K, T, r, sigma)
print Simposon_numerical_integrate(S0, K, T, r, sigma)
k,c=BSM_call_value_FFT(S0, K, T, r, sigma)
print np.interp(K, k, c)

19.6974420868
19.6974420868
19.7192567818

我们使用相同参数的FFT 对一组k 不同的期权定价，并和BSM 得到的期权价格比较

结果不是特别令人满意。下面我们讨论误差产生的原因。

插值误差

上图可以很好解释差值误差，不再采样点上k 我们只能通过插值来求得。在采样点相对准确。

一个细节是 我们一直使用的 k 是 \log K ,我们的 k_m 序列虽然是均匀的采样。我们日常更关心的真实的执行价格K，在转换成 K=e^k 后，采样会膨胀的厉害，随着k变大，不再采样点上k的插值误差会越来越大。这导致我们一开始满怀期待的一次求得一组只有K不同的期权定价 变得不实际。能不能通过增加采样点，提高采样频率避免差值误差呢，也不是能的，原因在参数选取部分讨论。

参数的选取

alpha

\alpha 存在的意义是为了是积分收敛。 \alpha 最终通过 \frac{e^{-\alpha k_m}}{\pi} 消除影响。不影响最终结果的计算才对。但实际上并不是这样。 重新回顾期权的价格公式。

\begin{align*} C_T(k)&=\frac{e^{-\alpha k}}{\pi} \int_0^\infty e^{-ivk} \Psi _T(v) dv\\ &\approx \frac{e^{-\alpha k}}{\pi} \int_0^B e^{-ivk} \Psi _T(v) dv\\ &= \frac{e^{-\alpha k}}{\pi} \sum^N_{j=0} e^{-i v_jk}\Psi _T(v_j) *w_j \end{align*}

\alpha 的不同取值会影响数值积分的准确率。我们先plot 出不同alpha 情况下

一个直观的结论： alpha 越大，我们要积分的项震荡的就越厉害，要求我们提供更多的采样点，才能很好的处理数值积分。而且收敛到0 的速度也更慢，也就需要更大的B（ \eta N ).即我们需要更大的N 和更小的 \eta 。在Simposon数值积分方法中我们不考虑效率的情况下 我们也正是这样做的。但是在FFT 中 更大的N 和更小 \eta 的会导致另一个矛盾，K 的采样距离膨胀。我们接下来讨论。

\eta ，\lambda 和N

我们简单的解释这些参数的作用

越小 ，数值积分越精细，能获得更好的local accuracy

越小，我们的k 采样就越密集，到时候插值计算就越准确

N 我们的采样点个数，越多自然越好 N\eta=B, N \lambda 则是k domin 的长度。

为了应用FFT 的时候我们获得一个附加条件： \lambda \eta=\frac{2\pi}{N} ， \lambda 和 \eta 是反比关系。我们如果要求 \eta 越小， \lambda 会变大，我们就更容易出现插值误差，k 的domin 会放大，采样稀疏，这也不符合实际。反之亦然， \lambda 越小，积分的准确率就会下降。 我们在讨论FFT 时尽量避开了时域频域的描述。这是傅里叶变换本身的一个特性决定的： bandlimit 和 timelimit 不可能同时成立，一个信号不能在时空域和频域上同时过于集中。

3总结

通过对FFT 定价敏感性分析，采样参数和alpha对定价十分敏感。鱼和熊掌不可兼得。相比直接的数值积分，FFT 会丢失更多的准确性。我们使用FFT 目的是更快计算出期权价格。实际上给了我们一个渠道，和bloom filter 一样，准确率换取效率，而这个兑换比率是通过参数设定可控的。

alpha 越大，我们的数值解就越好。但是代价时需要 更大的N 和更小的 \eta ，但这会导致k 采样膨胀，经验上取值在1~1.5 是合适的。

我们为了使用FFT， 被迫追加了新条件 \lambda \eta=\frac{2\pi}{N} ，能不能使得 \lambda 和 \eta 独立，不丢失准确性，并且也能获得 O(n\log n) 效率呢。能。Fractional FFT。。paper 正在读。。

参考书目：

1.Computational Methods in Finance

2.lecture Note for EE261 the fourier transform and its applications