所有由邓一硕发布的文章

关于邓一硕

专注于R语言在金融投资分析和计量经济学中的应用;《R in a Nut Shell》、《R Graphics Cookbook》、《Introductory Statistics with R》等书的译者。

蒙特卡洛方法与定积分计算

本文讲述一下蒙特卡洛模拟方法与定积分计算,首先从一个题目开始:设$0\leq f(x) \leq 1$,用蒙特卡洛模拟法求定积分$J=\int_{0}^{1}f(x)dx$的值。

随机投点法

设$(X,Y)$服从正方形 $\{0\leq x \leq 1,0\leq y\leq 1\}$上的均匀分布,则可知 $X,Y$分别服从[0,1]上的均匀分布,且$X,Y$相互独立。记事件$A=\{Y\leq f(X)\}$,则$A$的概率为

$
P(A)=P(Y\leq f(X))=\int_{0}^{1}\int_{0}^{f(x)}dydx=\int_{0}^{1}f(x)dx=J
$

即定积分$J$的值 就是事件$A$出现的频率。同时,由伯努利大数定律,我们可以用重复试验中$A$出现的频率作为 $p$的估计值。即将$(X,Y)$看成是正方形$\{0\leq x \leq 1,0\leq y \leq 1\}$内的随机投点,用随机点落在区域${y\leq f(x)}$中的频率作为定积分的近似值。这种方法就叫随机投点法,具体做法如下:

图1 随机投点法示意图

1、首先产生服从 $[0,1]$上的均匀分布的$2n$个随机数( $n $为随机投点个数,可以取很大,如 $n=10^4$ )并将其配对。

2、对这$n$对数据 $(x_i,y_i),i=1,2,…,n$ ,记录满足不等式$y_i\leq f(x_i)$的个数,这就是事件 $A$ 发生的频数$\mu_n$,由此可得事件$A$发生的频率$\frac{\mu_n}{n}$,则$J\approx \frac{\mu_n} {n}$。

举一实例,譬如要计算$\int_{0}^{1}e^{-x^2/2}/\sqrt{2\pi}dx$,模拟次数$n=10^4$时,R代码如下:

n=10^4;
 x=runif(n);
 y=runif(n);
 f=function(x)
 {
 exp(-x^2/2)/sqrt(2*pi)
 }
 mu_n=sum(y<f(x));
 J=mu_n/n;
 J

模拟次数 $n=10^5$ 时,令$n=10^5$,其余不变。

定积分$\int_{0}^{1}e^{-x^2/2}/\sqrt{2\pi}dx$的精确值和模拟值如下:

表1

精确值 $10^3$ $10^4$ $10^5$ $10^6$ $10^7$
0.3413447 0.342 0.344 0.34187 0.341539 0.341302
注:精确值用integrate(f,0,1)求得

扩展

如果你很细心,你会发现这个方法目前只适用于积分区间$[0,1]$ ,且积分函数 $f(x)$ 在区间$[0,1]$上的取值也位于 $[0,1]$ 内的情况。那么,对于一般区间 $[a,b]$ 上的定积分$J’=\int_{a}^{b}g(x)dx$ 呢?一个很明显的思路,如果我们可以将 $J’$ 与$J$ 建立代数关系就可以了。

首先,做线性变换,令 $y=(x-a)/(b-a)$ ,此时,

$x=(b-a)y+a$, $J’=(b-a)\int_{0}^{1}g[(b-a)y+a]dy$。

进一步如果在区间$[a,b]$上有$c\leq g(x) \leq d $ ,令

$f(y)=\frac{1}{d-c}{g(x)-c}=\frac{1}{d-c}{g[a+(b-a)y]-c} $,

则$0\leq f(y) \leq 1$。此时,可以得到$J’=\int_{a}^{b}g(x)dx=S_0J+c(b-a) $。

其中,$S_0=(b-a)(d-c) $, $J=\int_{0}^{1}f(y)dy $。

这说明,用随机投点法计算定积分方法具有普遍意义。

举一个实例,求定积分 $J’=\int_{2}^{5}e^{-x^2/2}/\sqrt{2\pi}dx $。

显然$a=2,b=5$,$c=min\{g(x)|2\leq x \leq 5\},d=max\{g(x)|2\leq x \leq 5\}$,由于$g(x)=e^{-x^2/2}/\sqrt{2\pi}$在 $[2,5]$上时单调减函数,所以$c=g(5),d=g(2)$,$S_0=(b-a)(d-c)$。R中代码为

a=2;
 b=5;
 g=function(x)
 {
 exp(-x^2/2)/sqrt(2*pi);
 }
 f=function(y)
 {
 (g(a+(b-a)*y)-c)/(d-c);
 }
 c=g(5);d=g(2);s_0=(b-a)*(d-c);
 n=10^4;
 x=runif(n);y=runif(n);
 mu_n=sum(y<=f(x));
 J=mu_n/n;
 J_0=s_0*J+c*(b-a);

定积分 $J’=\int_{2}^{5}e^{-x^2/2}/\sqrt{2\pi}dx$ 的精确值和模拟值如下:

表2

真实值 $10^3$ $10^4$ $10^5$ $10^6$ $10^7$
0.02274985 0.02332792 0.02311736 0.02262659 0.02284152 0.02278524
注:精确值用integrate(g,2,5)求得)

平均值法

蒙特卡洛模拟法计算定积分时还有另一种方法,叫平均值法。这个原理也很简单:设随机变量 $X$ 服从$[0,1]$上的均匀分布,则$Y=f(X)$的数学期望为

$E(f(X))=\int_{0}^{1}f(x)dx=J $

所以估计$J$的值就是估计$f(X)$的数学期望值。由辛钦大数定律,可以用$f(X)$的观察值的均值取估计$f(X)$的数学期望。具体做法:

先用计算机产生$n$个服从$[0,1]$上均匀分布的随机数:$x_i,i=1,2,…,n$。

对每一个$x_i$ ,计算$f(x_i)$。

计算$J\approx \frac{1}{n}\sum_{i=1}^{n}f(x_i)$。

譬如,计算 $J=\int_{0}^{1}e^{-x^2/2}/\sqrt{2\pi}dx $,R中的代码为

n=10^4;
 x=runif(n);
 f=function(x)
 {
 exp(-x^2/2)/sqrt(2*pi)
 }
 J=mean(f(x));

其精确值和模拟值如下:
表3

真实值 $10^3$ $10^4$ $10^5$ $10^6$ $10^7$
0.3413447 0.3405831 0.3410739 0.3414443 0.3414066 0.3413366

平均值法与随机投点法类似可以进行扩展,这里不再赘述。

结论

用蒙特卡洛模拟法计算定积分具有普遍意义。蒙特卡洛模拟方法为我们提供了一个看待世界的新思路,即一个不具随机性的事件可以通过一定的方法用随机事件来模拟或逼近。

风险评价的VaR方法简介

风险管理作为商业银行维持其正常经营的重要手段,19世纪初就已在世界范围内得到共识。西方发达国家早已建立起一套成熟的风险管理体系,其运作的依 据通常都是某种统计或者数学模型。虽然中国的风险管理才刚刚起步,然而,近年来在风险管理方面也进行了如火如荼的探究。由于中国过去的风险管理理念是建立 在定性和主观经验的基础之上,因此现阶段对定量的数学模型的引进及探讨十分得宠。VaR作为风险管理的最重要的模型之一,尤其需要学者和风险管 理者的重视,本文对其进行简单介绍。

一、VaR模型的由来

1993年7月G30国成员曾发表了一个关于金融衍生工具的报告,首次建议用“风险价值系统”(Value at Risk System,简称VaRS)来评估金融风险。1999年的新巴塞尔协议征求意见稿中, 新巴塞尔委员会又极力提倡商业银行用VaR 模型度量其所面临的信用风险, 在2004年发布的新巴塞尔协议中委员会把风险管理的对象扩大到市场风险、信用风险和操作风险的总和,并进一步主张用VaR模型对商 业银行面临的风险进行综合管理。此外,委员会也鼓励商业银行在满足监管和审计要求的前提下, 可以自己建立以VaR 为基础的内部模型。此后,VaR 模型作为一个很好的风险管理工具开始正式在新巴塞尔协议中获得应用和推广,并逐步奠定了其在风险管理领域的元老地位。

二、VaR模型的含义

VaR(Value at Risk),通译为“风险价值”。其内涵是在正常情况下,一定时期内$\Delta t$内,一定的置信水平$1-\alpha$某种资产组合面临的最大损失。统计学表达为:
$Prob(\Delta p\leq VaR)=1-\alpha$
其中$\Delta p$是指在一定的时期$\Delta t$内某种资产组合市场价值的变化,$1-\alpha$为给定的概率。即在一定的持有期$\Delta t$内,给定的置信水平$1-\alpha$下,该资产组合的最大损失不会超过VaR。用VaR进行风险衡 量时,首先要确定持有期和置信水平,巴塞尔委员会规定的持有期标准为10天,置信水平为99%,但各个商业银行可以确定自己的标准。如 J.P.Morgan 公司在1994年的年报中,规定的持有期为1天,置信水平为95%,VaR值为1500万美元。其含义即为J.P.Morgan公 司在一天内,其所持有的风险头寸的损失小于1500万的概率为95%,超过1500 万的概率为5%。当然,使用者也可以根据自己的喜好来选择持有期和置信水平,一般而言,置信水平直接反应使用者的风险偏好水平。

三、VaR的主要性质

VaR的 主要性质:
(1)变换不变性:$a\in R$,满足$VaR(X+a)=VaR(X)+a$
(2)正齐次性:$h<0$,满足$VaR(hX)=hVaR(X)$,保证资产的风 险与其持有的头寸成正比。
(3)协单调可加性:$VaR(X_1+X_2)=VaR(X_1)+VaR(X_2)$ 。
(4)不满足次可加性和凸性。不满足前者意味着资产组合的风险不一定小于各资产风险之和。这一点是不合理的,因为它意味着,一个金融机构不能通过计算其分支机构 的VaR来 推导整个机构的VaR。 不满足凸性意味着以VaR为 目标函数的规划问题一般不是凸规划,其局部最优解不一定是全局最优解。所以在基于VaR对资产组合进行优化时,可能存在多个局部极值。因此,也 就无法求得最佳资产组合。这也是VaR用于资产组合风险研究时候的主要障碍。
(5)满足一阶随机占优。
(6)VaR 关于概率水平$1-\alpha$不是连续的。
Artzner(1999)指出好的风险度量方法需满足一致性,满足一致性意味着其能同时满足次可加性、正齐次性、单调性和变换不变性四个性质。显然, 不具有这个特性。

三、VaR模型的优点

VaR模 型作为现代商业银行风险管理最重要的方法之一,存在不少优点。首先,VaR 使用规范的数理统计技术和现代工程方法来度量银行风险,与以往靠定性和主观经验的风险度量技术相比更具客观性;其次,它使用单一指标对风险进行衡量,具有 直观性,即使没有专业背景的投资者和管理者,也能通过这一指标评价风险的大小;再次,它不仅可以衡量单一的金融资产的风险,还能衡量投资组合的风险;而 且,它对风险的衡量具有前瞻性,是对未来风险的衡量,不像以往对风险的衡量都是在事后进行;最后,VaR 把对未来预期损失的规模和发生的可能性结合起来,是管理者不仅能了解损失的规模,还能了解在这一规模上损失的概率。并且通过对不同的置信区间的选择可以得 到不同的最大损失规模,便于管理者了解在不同可能程度上的风险大小。

四、VaR模型的缺点

VaR模 型当然也有一些缺点,首当其冲的是,VaR模型是对正常的市场环境中的金融资产的风险的衡量,但一旦金融环境出现动荡,即当极端情况发生时,VaR模型所代表 的风险大小就失去了参考价值。然而,很多时候恰恰是这些极端值,决定了为了完善对所有市场状况下的风险衡量。因此,为了完善对所有市场状况下的风险衡量, 通常在VaR模 型的基础上可以引入Stress Test即压力试验和极值分析两种方法进行辅助。压力试验主要是在违背模型假设的极端市场情景下,对资产组合收益的不利影响进行评价;而对VaR进行辅助的 极值分析的方法有两种,分别是BMM和POT,其中最常用的一种模型是POT模型。它是在当风险规模超过某一最大值的情况下进行的建模。它直接处理风险概 率分布的尾部,事先并不对数据的分布做任何假设,在利用设定参数建立的模型的基础上,对极端情况下的风险规模和概率进行衡量。
第二,VaR模型是在收益分布为正态分布的情况下的衡量。但事实表明,资产的收益的尾部比正态分布的尾部更厚,通常成为厚尾性,且其与正态分布的对称性也并不一致。当这 种情况出现的时候,VaR模 型就不会产生一致性度量的结果。所谓一致性风险度量即是风险衡量得出的度量值的大小与风险的实际大小具有一致性。对风险大的金融资产衡量得出的风险值大于 对风险小的金融资产衡量得出的风险值,具有相同风险的金融资产具有相同的风险度量值,具有不同风险的金融资产具有不同的风险度量值。针对VaR模型度量的 不一致性可以引进ES(Expected shortfall)或者CVaR(conditional value at risk)模型对VaR模型进行修 正,两者是同一概念的不同说法。CVaR称为条件风险价值。它是当资产组合的损失超过某个给定的VaR的情况下,资产组合的损失的期望。用数学公式表达为:
$CVaR_{\alpha}=E(-X|-X\leq VaR_{\alpha}(x))$
其中X表示资产的损益。
根据CVaR的定义可以证明CVaR满足以下性质:
(1)是一致连续的。
(2)满足次可加性。
$\forall X,Y$,满足$\rho(X+Y)\leq\rho(X)+\rho(Y)$。
(3)满足二阶随机占优。
(4)满足单调性。
$\forall X,Y$满足$\rho(X)\leq \rho(Y)$,如果 $ X\leq Y $。