期权定价问题的有限元Richardson外推法

N = 1e6 时的图形

期权定价问题的有限元Richardson外推法

德国布隆斯威克人。德国的数学家、物理学家和天文学家。高斯幼年时就显示出非凡的数学才能，得到Carl Wil-helm Ferdinand大公的赏识。在大公的支持下，1795—1798年在哥廷根(Gottingen)大学学习，1799年因证明代数学的基本定理而获得哈勒(Halle)大学的博士学位。从1807年到1855年逝世，一直担任哥廷根天文台台长兼大学教授。1796年用直尺圆规作出了正十七边形一自欧几里得以来两千年间几何作图的一个难题。接着又证明了数论中的欧勒猜想—二次互反律。 据说从此后他决心放弃古典文学而献身于数学。1801年用自己的行星轨道计称法和最小二乘法算出了意大利天文学家皮亚齐 (1746—1826) 发现的谷神星轨道; 谷神星的轨道计算使他一举名震世界。同年，出版经典著作《算术研究》，任职期间，高斯致力于数论、代数、几何、分析、复变函数、统计数学等多方面的研究、取得了一系列的成果。高斯定理、高斯公式、高斯函数等以他命名的多种发现至今仍在许多数学、科学部门中闪烁着光辉。高斯还涉足了大地测量工作。 为了进行长距离测量，发明了“目光反射器”，并在理论上创造了“大地问题解法”，导致他开创了曲面微分几何的理论。并由他的学生黎曼发展为黎曼几何。与德国的物理学家韦伯(wilhelm Eduavd Weber，1804—1891)一道建立了电磁学中的高斯单位制; 1833年还和韦伯一起发明了电磁电极。高斯的治学态度十分严谨。他的格言是 “瑰丽的大厦建成后，应拆除杂乱无章的脚手架。” 因此他发表的每篇著作都是经过仔细推敲、无懈可击的精品。因此，发表的论文比研究工作要少得多，但研究项目可在日记和书信中见到。全集包括日记、书信共计12卷。他是19世纪前半世纪最伟大的数学家。 [4]

蒙特·卡罗（Monte Carlo）法求定积分

MATLAB工作区 图示蒙特-卡罗方法

从上述结果我们可以看到，所得积分结果为 期权定价问题的有限元Richardson外推法 I=2.0075 ，可见其与真值 2 还是十分接近的。如果我们将 N 取得更大一些，所得结果可以更加逼近真值，比如我们取 N=10^6 将得到

N = 期权定价问题的有限元Richardson外推法 1e6 时的图形

我们容易知道这个定积分的值为 4 。但是采用蒙特·卡罗方法将得到

图示

此时的随机数数量已经取为 N=10^5 了，所得结果与真值仍然存在很大的差距，检查代码并结合上述图形发现计算过程应该没有问题，我也不知道问题出在哪里，希望广大知友能过指出错误之处，谢谢！

Romberg 算法

，不在 " width="" height="" />區間上的積分可以做相應的變量代換轉化到該區間上去。

我們先使用兩點梯形公式計算初始迭代值=>(f(0)+f(1))>" width="" height="" />，將區間不斷做二等分，得到一系列復化梯形公式 >=>T_>+>>\sum _^>f\!\left(>>\right),k=1,2,\cdots .>" width="" height="" /> 這些都是後續迭代的初始值，上式給出了初始值的通式。下面我們開始第一次迭代，令 >=>-T_>>>,k=0,1,2,\cdots .>" width="" height="" /> 上式的>>" 期权定价问题的有限元Richardson外推法 width="" height="" />是 Simpson 公式的復化自適應求積公式，進一步迭代 >=>-S_>>>,k=0,1,2,\cdots .>" width="" height="" /> 上式的>>" width="" height="" />是 Cotes 公式的復化自適應求積公式，進一步迭代 >=>-C_>>>,k=0,1,2,\cdots .>" width="" height="" /> 上式的>>" width="" height="" />是 Romberg 公式的復化自適應求積公式。

這樣一直進行下去，精度不斷提高，一般意義下的遞推公式為 =T_-T_>-1>>,j,k=0,1,2,\cdots .>" width="" height="" /> 在上面的計算中，>=T_,S_>=T_,C_>=T_,R_>=T_.>" width="" height="" /> 上述計算的精度可以控制為 -I=O(h^).>" width="" 期权定价问题的有限元Richardson外推法 height="" /> 判斷迭代停止的條件可以是先驗的，即預先給出絕對誤差0>" width="" height="" 期权定价问题的有限元Richardson外推法 />後判斷停止的條件可以是 -T_| ^.>" width="" height="" /> 這時可以以>" width="" height="" />作為積分" width="" height="" />的數值解。

数值积分常用方法

(1) 式的几何意义即为：底为 b-a 而高为 f(\xi) 的矩形的面积恰等于所求曲边梯形的面积 I 。因此，要想求出 (1) 式左端积分，我们只需要知道三个值： a,b,f(\xi) 即可，这里 a,b 是显然的，问题在于 \xi 的具体位置一般是不清楚的，从而 f(\xi) 未知。我们暂且将 f(\xi) 称为区间 [a,b] 上的平均高度，我们的目标就是寻求一种求出平均高度 f(\xi) 的算法，这样我们的问题也就解决了。

首先我们想到的是用区间两端点的“高度” f(a) 与 f(b) 的算术平均值作为平均高度 f(\xi) 的近似值，从而导出下述求积公式

\int_a^b f(期权定价问题的有限元Richardson外推法 x)\mathrmx\approx\frac[f(a)+f(b)] \\ (2)

(2) 式即为我们我们熟知的梯形公式。写出MATLAB函数为

而若用区间中点 \frac 的“高度” f(\frac) 近似取代平均高度 f(\xi) ，则可导出中矩形公式（简称矩形公式）

总结：更一般地，我们可以在区间 [a,b] 上适当选取某些节点 x_k ，然后用 f(x_k) 的加权平均得到平均高度 f(\xi) 的近似值，这样构造出来的求积公式具有下列通式

其中 x_k 期权定价问题的有限元Richardson外推法 期权定价问题的有限元Richardson外推法 称为求积节点， A_k 称为求积系数，也称作伴随节点 x_k 的权。

这种数值积分方法通常称为机械求积，其优势在于将积分求值问题归结为被积函数值的计算，很适合在计算机上使用。

前面提到的梯形公式 (2) 和矩形公式 (3) 均具有一次代数精度。

求积公式的收敛性与稳定性

在求积公式 (4) 中，若 \lim_\sum\limits_^A_kf(x_k)=\int_a^bf(x)\mathrmx,\\ (5) 其中 h=\max\limits_<1\leqslant i\leqslant n>\\> ，则称求积公式 (4) 是收敛的。
对任给 \varepsilon>0 ，若 \exists \delta>0 ，只要 |f(x_k)-\tilde_k|\leqslant \delta(k=0,1,2,\cdots,n) 就有 |I_n(f)-I_n(\tilde)|=\left|\sum\limits_^A_k[f(x_k)-\tilde_k]\right|\leqslant \varepsilon,\\ (期权定价问题的有限元Richardson外推法 6) 成立，其中 \tilde_k 为计算 f(x_k) 期权定价问题的有限元Richardson外推法 时产生误差 \delta_k 后实际得到的值，即 f(x_k)=\tilde_k+\delta_k ，则称求积公式 (4) 是稳定的。

定理 \mathbf 期权定价问题的有限元Richardson外推法 若求积公式 (4) 中系数 A_k>0(k=0,1,\cdots,n) ，则此求积公式是稳定的。

牛顿-柯特斯(Newton-Cotes)公式

将积分区间 [a,b] 划分为 n 等份，此时步长为 h=\frac ，选取等距节点 x_k=a+kh 构造出的插值型求积公式

称为牛顿-柯特斯(Newton-Cotes)期权定价问题的有限元Richardson外推法 公式，其中 C_k^ 称为柯特斯系数。可以通过下式确定

• 当 n=1 时， 期权定价问题的有限元Richardson外推法 C_0^=C_1^=\frac12, 这时的求积公式便是之前的梯形公式 (2) 。
• 当 n=2 时， C_0^=\frac14\int_0^2(t-1)(t-2)\mathrmt=\frac16,\\ C_1^=-\frac12\int_0^2t(t-2)\mathrmt=\frac46,\\ C_2^=\frac14\int_0^2t(t-1)\mathrmt=\frac16. \\ (9) 此时对应的求积公式便是辛普森(Simpson)公式： S=\frac\left[f(a)+4f\left(\frac\right)+f(b)\right]. \\ (10) 写成MATLAB函数为

• 当 n=4 时的牛顿-柯特斯公式则特别称为柯特斯公式，其形式为 C=\frac[7f(x_0)+32f(x_1)+12f(期权定价问题的有限元Richardson外推法 x_2)+32f(x_3)+7f(x_4)]. \\ (11) 这里 h=\frac,x_k=a+kh(k=0,1,2,3,4) 。

复合求积公式

由于牛顿-柯特斯公式在 n\geqslant 8 时不具有稳定性，故不可能再通过提高阶的方法来提高求积精度。复合求积法便是通过把积分区间分成若干个子区间（通常是等分），再在每个子区间上使用低阶求积公式的方法来提高精度的。

其实细心的同学可以发现，前面我所展示的MATLAB代码使用的便是对每个求积公式使用复合法完成的。常用的有复合梯形公式和复合辛普森公式。

龙贝格(Romberg)求积公式

其中 T_ 表示在 T_n 基础上步长 h 减半后的积分值。

定理 \mathbf 设 f(x)\in C^\infty[a,b] ，则有

T(h)=I+\alpha_1h^2+\alpha_2h^4+\cdots+\alpha_l h^+\cdots, \\ (12) 其中 T(h)=T_n ，系数 \alpha_l(l=1,2 ,\cdots) 与 h 无关。

从 (12) 式可以看出， T(h)\approx I 是 O(h^2) 阶，若用 \frac h2 代替 h ，有
T(\frac h2)=I+\alpha_1\frac4+\alpha_2\frac+\cdots+\alpha_l\left(\frac h2\right)^+\cdots, \\ (13)

再用 4 乘 (13) 式再减去 (12) 期权定价问题的有限元Richardson外推法 式后再除以 3 得

这里 \beta_1,\beta_2,\cdots 是与 h 无关的系数。

从 (14) 式可看出，用近似积分值 I ，其误差阶为 O(h^4) ，显然误差阶是提高了。类似这种将计算 I 的近似值的误差阶由 O(h^2) 提高到 O(h^4) 的方法称为外推算法，也称为理查德森(Richardson)外推算法。

龙贝格(Romberg)算法

此方法常称为理查德森外推加速方法。

设用 T_0^ 表示二分 k 次后求得的梯形值，且以 T_m^ 表示序列 \ 的 m 次加速值，则由递推公式 (15) 可得

上次则称为龙贝格求积算法，计算过程如下：

1. 取 k=0,h=b-a ，求 T_0^=\frac h2[f(a)+f(b)] 。令 1\rightarrow k ( k 记区间 [a,b] 的二分次数)。
2. 求梯形值 T_0\left(\frac\right) ，即按递推公式 (11) 计算 T_0^ 。
3. 求加速值。
4. 若 |T_k^-T_|<\varepsilon (预先给定的精度)，则终止计算，并且 T_k^\approx I ；否则令 k+1\rightarrow k 转 2. 继续计算。

自适应积分方法

复合求积方法通常适用于被积函数变化不太大的积分，如果在求积区间中被积函数变化很大，有的部分函数值变化剧烈，另一部分却变化平缓。这时如果统一将区间等分再用复合求积公式计算积分将会导致计算量很大，我们想实现的是在满足误差要求的前提下，对变化剧烈部分将区间细分，而平缓部分则可使用大步长，也即针对被积函数在区间上不同情形采用不同的步长，使得再满足精度前提下积分计算工作量尽可能小。其算法技巧是在不同区间上预测被积函数变化的剧烈程度确定相应步长。这就是自适应积分方法的基本思想。