一、实验目的
使学生了解蒙特卡罗方法在金融学中的应用。
二、实验内容
1.年金保险中的蒙特卡罗模拟；
2.基于EViews的期权定价蒙特卡罗模拟。
三、实验工具
EViews软件。
四、年金保险中的蒙特卡罗模拟
在这节里，我们用年金保险的例子来介绍蒙特卡罗模拟的方法。年金保险是寿险公司最常用的一种险种。选用Eviews软件，来完成任意一次蒙特卡洛模拟，并由此方法来计算年金保险案例中保险公司的净收入情况。
假定某保险公司推出一项纯粹保险的业务，该保险业务要求投保人保费支付期为20年（即从41岁开始支付保费，到61岁停止支付保费），每年支付保费200元，若投保人在保费支付期内死亡，则保险公司无需对投保人进行补偿。投保人支付保费20年后，若投保人仍然生存，则保险公司在其生存期内每年给他（她）支付500元，另外，假定被保险人41-60岁期间的死亡概率为3.5%，61岁以后死亡的概率为6％。金融市场的无风险贴现率为4.5％。用蒙特卡罗模拟法计算该保险公司在开展这项保险业务出现资不抵债的概率以及公司的期望收益。
我们的解题过程大致分为两步，首先解出一个客户时保险公司的期望收益，二是用蒙特卡罗方法模拟出保险公司的期望收益。为了反映出每一步模拟时的情况，我们先模拟一个客户保险的期望收益，其步骤如下：
4.1程序设计思路
在程序设计时，我们通过蒙特卡洛模拟的方法产生一系列随机序列，并且假设这些序列是呈现出正态分布的，进而可以得出这些随机序列的概率，当投保人的死亡概率大于其中的某一随机概率时，我们则认为该投保人已经死亡，相应地，投保人保费以及保险年金的现金流则开始中断。
4.2操作具体步骤如下
首先开启Eviews程序，在其命令框中输入：
Program insurance这时Eviews就会打开一个名为insurance的程序设计框，在该程序设计框中，输入以下程序，其中括号内的文字为程序解释文字）。
'this programm is about life insurance case （’为程序说明文字解释符号）
workfile insurance u 1 70 （定义非时序工作文件insurance，其样本容量为70）
series u=@rnorm （生成一个白噪声序列，定义为u）
series v=@dnorm(u) （求出白噪声序列u的相应概率值，定义为序列v）
series age （定义序列age，为投保人的年龄）
series rf （定义序列rf，为无风险收益率）
rf=0.05 （设该无风险收益率恒为5%）
series discount （定义序列discount，为折现率）
smpl 1 70 （把样本容量控制为1-70）
age(1)=41 序列age第一个序列值为41，即投保人的交年费的首年年龄）
discount(1)=1/(1+rf(1)) （定义序列discount的第一个序列值为1/(1+rf(1))）
smpl 2 70 （把样本容量控制为2-70）
age=age(-1)+1 （定义序列age的产生方式）
discount=discount(-1)/(1+rf(1)) （定义序列discount的产生方式）
series income （定义序列income，为保险公司的当年收支现金流）
series prob （定义序列prob，为投保人死亡的概率）
smpl 1 20 （把样本容量控制为1-20）
income=200 （序列income的前20个值为200，表示投保人所交的年金）
prob=0.035 （序列prob的前20个值为0.035，表示投保人前20年死亡的概率）
smpl 21 70 （把样本容量控制为21-70）
income=-500 （序列income的后50个值为-500，表示保险公司所支付的年金）
prob=0.06 （序列prob的后50个值为0.035，表示投保人后50年死亡的概率）
series pv （定义序列pv，表示保险公司的年金现值）
smpl 1 70 （把样本容量控制为1-70）
pv(1)=income(1) （定义序列pv的第一个值为序列income的第一个值，为200）
smpl 2 70 （把样本容量控制为2-70）
!N=70 （定义控制变量!N为70）
for !i=2 to !N （!i为控制变量，通过for…next语句使得计算循环!N次）
if pv(!i-1)=0 then （这是一个条件结构程序，如果pv的前一个值为0，那么
pv(!i)=0 pv的当前值也为0）
else （下面一个嵌套条件结构程序）
if v(!i)>prob(!i) then （如果随机概率大于投保人的死亡概率，那么 pv(!i)=income(!i)*discount(!i)+pv(!i-1) 投保人仍然活着，pv值为现金流的折现值）
else （条件结构if…then，else…endif）
pv(!i)=0 （如果投保人死亡，那么pv当前值为0）
endif （结束嵌套条件语句）
endif （结束条件语句）
next （与for…，next循环结构对应）
所有程序输入完毕以后，如图4.1所示，点击程序窗口的Save键保存程序，然后点击Run，运行程序。程序运行结果会产生一个insurance的工作文件（workfile），打开序列（series）pv，可以看到对其中的某个投保人的模拟的结果，pv序列显示的是保险公司的现金流现值，其结果如图4.2所示。由于蒙特卡洛模拟是的原理是基于随机过程的，因此，每一次程序运行的结果不会相同，而保险公司只需要增大样本运行的次数，根据模拟的大量本就可以得出该项保险产品的期望现金流状况。
图4.1 Eviews中程序输入界面
图4.2 Eviews程序运行结果界面
五、基于EViews 的蒙特卡罗模拟
5.1相关系数分布的蒙特卡罗模拟
例题1、两个相互独立O(l)变量相关系数分布的蒙特卡罗模拟
1、程序设计思路
首先生成两个相互独立的O(l)变量（样本容量为100)，然后计算这两个O(l)变量的相关系数。将这个过程重复5万次，每次得到的相关系数都存储在rxy里面，最后生成列rxy的直方图以观察和分析相关系数的分布特征。
2、操作具体步骤如下
首先编辑程序文件（program）。打开Eviews， 在主菜单中点击File，选择New / Program , 图5.1所示，然后输入如下程序：
图5.1 Eviews 的操作界面
“The distribution of r of two indePendent O(l) series” <程序的说明文字必须以符号“ ‘” 开头>
!N =50000 （定义数值变量!N，表示模拟50,000 次）
workfile simu5 u 1 !N （定义一个非时序工作文件，imus，容量为!N）
series rxy （定义序列rxy，储存每次计算得到的相关系数值）
for !i=1 !N （!i为控制变量，通过一个for ……next语句使计算循环N 次）
smpl 1 100 （把样本容量设置成100）
series u=＠nrnd （生成一个白噪声序列，定义为u）
series v=＠nrnd （生成一个白噪声序列，定义为v）
u(l)=0 （定义序列u 的初始值为0）
v(1)=0 （定义序列u 的初始值为0）
series x （定义序列 x）
series y （定义序列 y）
x(1 =0 （定义序列X 的初始值为0）
y(1)=O （定义序列y 的初始值为0）
smpl 2 100 （把样本容量设置成2-100）
x=x(-1)+u （生成一阶单积序列（随机游走序列）x）
y=y(-1)+v （生成一阶单积序列（随机游走序列）y）
smpl 1 100 （把样本容量恢复成l-100）
scalar r2=@cor(x,y) （计算序列X 和y 的相关系数，并将结果放到标量r2 中）
rxy(!i)=r2 （将相关系数值r2 放到序列rxy中）
next （与上面的for 相对应，表示至此完成！N 次循环计算）
smpll !N （工作文件simu5 的样本范围（sample）恢复到l-50,000）
rxy.hist （显示序列rxy的直方图以及有关统计量的值）
在输入完毕后，点击程序窗口的Save键保存程序，然后点击图5.3程序窗口中的Run键，在随后弹出的对话框中的执行方式（Execution mode）选择区选Quiet（安静）运行方式。点击OK 键，大约需要4分钟的时间就能完成5万次的模拟运行。最后显示出相关系数r的分布直方图（见图5.4）。
图5.2 程 序 界 面
图5.3 运行程序界面
两个独立平稳变量的相关系数应该服从均值为零的正态分布但是从图6来看，均值虽然接近零，但分布为倒U字型。JB=2595.172 的相应概率为0.0000，也说明相关系数值不服从正态分布。从这个模拟结果可以看出，非平稳变量之间存在虚假相关问题。
图5.4 相关系数r的分布直方图
5.2单位根检验DF统计量分布的蒙特卡罗模拟
Dickey和Fuller于1976年提出了用于单位根检验的DF统计量，可以利用Eviews程序对DF统计量的分布进行蒙特卡罗模拟。
- DF统计量的分布特征
先给出三个简单的自回归模型：
$$y_{t}=\beta y_{t-1}+u_{t}, y_{0} =0$$ （1）
$$y_{t}=\mu+\beta y_{t-1}+u_{t},y_{0}=0$$ （2）
$$y_{t}= \mu + \alpha t + \beta y_{t-1} + u_{ t} , y_{0} = 0$$ （3）
其中μ是位移项（也称漂移项）,αt是趋势项。ut满足弱条件，即允许ut存在暂时性自相关和异方差。显然，对于以上三个模型，当︱β︱<1时，yt是稳的；当︱β︱=1时，yt是非平稳的。以模型(1)式为例，若β＝0时，统计量t(β^)=(β^)/s(β^)～(t（t-1）)的极限分布为标准正态分布；若︱β︱<1时，统计t(β^)=(β^)-(β)/s(β^))渐进服从标准正态分布；而当︱β︱=1时，定义此条件下的统计量为：DF=t(β^)=(β^-1)/s(β^)。
2、程序设计思路
首先，在原假设条件下生成一个容量为100的I(l)序列yt，然后利用这个I(l)序列进行带漂移项的DF检验，提取相应的DF统计量值，并保存到给定的序列中。将这个过程重复1万次，画出该序列（DF值）的直方图，并显示其第1、5和10百分位数值。
3、操作具休步骤
首先生成一个新的程序文件(Program)。打开Eviews，在主菜单中点击File，选择New/ Progrom，然后输入下列程序语句。
‘The distribution of DF
!N=10000 ——(定义数值变量N ，表示模拟次数为10,000次)
workfile dfvalue u 1 !N ——(定义一个非时序工作文件，名称为dfralue，容量为10,000)
series df ——(定义一个序列，名为df )
equation eq01 ——(定义一个方程，名为eq01)
fo i＝1 to !N ——(!1为控制变量，通过一个for …… next语句使计算循环10,000次)
smpl 1 100 ——(把样本范围设置成1-100)
series u=@nrnd ——(生成一个白噪声序列，定义为u)
u(l)=0 ——(定义序列u的初始值为0)
series y ——(定义一个序列y)
y(1)=0 ——(定义序列y的初始值为0)
smpl 2 100 ——(把样本范围设置成2-100)
y=y(-l)+ u —(定义序列y等于其前一期值加上一个白噪声，即随机游走)
eq01.ls d(y) y(-l) c ——(估计带有漂移项的DF检验式，即)
smpl ＠all ——(将样本范围改为整个工作文件的范围)
df(!i)=eq01.@tstats(l) ——(提取回归方程eq0l中的第1个回归系数对应的t统计量
的值，即DF值，保存在系列df的相应位置)
next ——(与for对应,表示至此完成一次循环计算，共要完成1万次)
show df.hist ——(显示系列df的直方图以及有关统计量的值)
matrix(2,3) cvalue ——(定义2×3矩阵，名为cvalue)
cvalue.fill(b=r) 0.01 ,0.05 , 0.10 ——(按行往矩阵中填充3个数据)
!K =1 ——定义控制变量!k = l)
For ％ 1 0.01 0.05 0.10 ——(将这3个数据作为字符串进行循环)
cvalue(2,!k)=＠quantile(df,% 1) —(计算df分布的3个分位数值，放到矩
阵cvalue的第2行第!k列中)
!K=!k+l ——(!k的数值加1)
next ——(与for相对应，至此完成一次字符串的循环计算)
show cvalue ——(显示矩阵Cvalue，即DF统计量的第1,5,10 百分位数值)
在program文件窗口中点击run键，按照“Quiet ( fast ) no screen/status line updates”模式运行(具体界面见上例中的图示), 1分钟左右能完成(运行过程中按Fl键或Esc键可以随时终止运行)。DF统计量的分布直方图见图5.5。
从图5.6可知，对应于1％、5％和10％置信度的临界值分别为：-3.520502、-2.914110和-2.592803(每次模拟的结果会略有差异)。
图5.5 DF统计量的分布直方图
图5.6 检验水平的临界值
5.3 股票价格波动的蒙特卡罗模拟
在这一部分，我们用股票价格波动的例子来介绍蒙特卡罗模拟的方法在金融学中的应用。在本例中将尝试使用Eviews中的相关技术对股票价格波动进行模拟测算。Eviews自身既为我们提供了方便快捷的交互式窗口操作方式，也提供了强大的程序语句，借助于其程序语句，可以完成更多的复杂任务。在本次实验中，我们主要运用程序语言的方式对数据进行处理和分析。
现以上证综指从2000年1月4日到2009年6月17日每个交易日的收盘价（指数）为历史样本数据（见图5.7)，然后用蒙特卡罗模拟10 , 000个股票价格在单位日内可能出现的变化率分布。
1、计算原理
选择几何布朗运动模型作为民生银行股票价格变化的随机过程，它是股票价格行为中最常用的模型之一。几何布朗运动模型的离散形式为：
$$\Delta S / S=\mu\Delta t+\sigma\varepsilon\Delta t^{1/2}$$ （4）
图5.7上证综指历史交易收盘价（已隐藏部分数据）
其中：变量ΔS为短时间Δt后股票价格S的变化，ε为标准正态随机变量，σ为股票收益的标准差（波动率），μ为单位时间内股票的预期收益率，σ和μ假设为常数。
方程（4）式等号左边是短时间Δt后股票的收益比，右边的μΔt项是这一收益的期望值，σεΔt1/2项是收益的随机部分，随机部分的方差（即整个收益的方差）为了σ2Δt。方程（4）表明ΔS/S是均值为μΔt，标准差为σΔt的正态分布，即：
$$\Delta S/S\sim\varphi(\mu\Delta t,\sigma\Delta t^{1/2})$$ （5）
由于估计的是单位日内股票价格波动，所以当Δt=1时，股价变化的随机模型为：
$$\Delta S/S\sim\varphi(\mu,\sigma)$$ （6）
在计算收益率的时候，用对数收益率代替百分比收益率，因为对数收益率可以保持连续性。将对数收益率表示为rt=ln（pt/pt-1），所以，(2)式中的μ和σ分别为历史样本的每天股票的对数收益率的均值和标准差。
现在，利用从2000年1月4日到2009年6月17日每天股票的收盘价（一共2247个样本），计算可得μ=m，σ=s，从而可估算出上证综指收益率变化服从的分布。
$$\Delta S/S\sim\varphi(m,s)$$ （7）
其中，φ（m,s）表示均值为m、标准差为S的正态分布。
由以上所得上证综指收益率的均值m和标准差s计算股价变化率序列。从标准正态分布(即u～φ(0,l))中取样值，将其转换为φ(m,s)中的样本点，来模拟股票价格的一条轨迹，该操作的第1步是从标准正态分布（即u～φ(0,l)）中取样值u1,然后将其转换为φ(m,s)中样本点x1， 利用公式：
$$x_1=\mu+\sigma\varepsilon=\mu+\sigma u_1$$ （8）
$$\Delta S_{1}=S_{1}(\mu+\sigma u_{1})=S_{1}x_{1}$$ （9）
$$\Delta S_1 / S_1=\mu+\sigma u_1=x_1$$ （10）
再依次类推，可得上证综指2247个样本的股价变化率序列。
同理，我们可以利用所得股价收益率均值m和标准差s参数和公式对上证综指10 ,000个股票价格的变化进行蒙特卡罗模拟，从而得出其分布及其分布统计量。因此我们解题大致分为两步：一是先将Excel 数据导入Eviews软件中，计算出上证综指对数收益率的均值和标准差，以确定上证综指变化模拟所用的参数。二是用蒙特卡罗方模拟方法模拟出10,000个股票价格单位日内可能出现的变化率分布及其分布统计量。
2、计算步骤
首先导入已知数据。打开Eviews，并在其中新建一个工作文件，将Excel表内的数据导入到工作文件中，生成价格序列price。数据表中包括从2000年1月4日到2009年6月17日期间一共2247个交易日的收盘价。我们根据数据的特点，在Eviews中选择Undated 类型建立工作文件，并在样本值中键入总样本数2247。如图5.8所示。
图5.8 生成序列price
其次计算出上证综指的对数收益率的均值和标准差，作为蒙特卡罗模拟参数。在打开的工作文件中建立一个新的程序文件，新建一个程序是在主菜单中选择File / New / Program，屏幕会显示一个文本编辑窗口，也可以键入命令program 打开程序编辑窗口，后面还可以输入程序名称，如p1:
Program p1
再在文本编辑窗口键入如下程序：
series price
series m
series v
series s
series sp
series pl
genr pl=price()/price (-1)
“表示根据表达式生成一组新的序列pl”
series r
genr r=log(pl)
series u=@nrnd
“表示生成一个均值为O，方差为1的标准正态分布随机数序列，也即白噪声序列了
m=＠mean(r)
“表示生成r的均值序列m
v=＠var(r)
表示生成r的方差序列v
s=＠sqr(v)
“表示生成r 的平方根序列,
series x=m+s*u
series sp=price*x
按Save键保存并按Run运行该程序生成上证综指对数收益率序列r的均值和标准差序列，从运行结果可知：上证综指2247个历史样本每天股票的对数收益率的均值m =μ=0.0002829,标准差s=σ=0.01662878。如图5.9所示
图5.9 上证综指均值m和标准差s序列
3. 利用以上所得参数蒙特卡罗模拟10,000个股票价格的可能变化率。
利用从2000年1月4日到2009年6月17日每天股票的收盘价（一共2247个样本）可计算可得参数m =μ=0.0002829，s=σ=0.01662878，模拟股票价格变化公式为：
$$\Delta S_{t} = S_{t}( m + s u_{t}) = S_{t} x_{t}$$
$$\Delta S_{t} / S_{t} =(m+su_{t})=x_{t}$$
在参数和公式都具备的时候，现在在Eviews软件中对上证综指价格变化进行10,000次蒙特卡洛模拟。
我们在上面建立的工作文件中，新建另一个程序文件（见图5.9），并在程序文件的编辑窗口中键入如下程序内容：
series ssp
for !k=1 to !n
“!k为控制循环变量，通过一个for…next 语句使计算循环!n次”
series u=@nrnd
“表示生成一个白噪声序列u”
series x
x(1)=(0.000282998884054344)+(0.016628780888637)*u(1)
series price
price(1)=1406.37+u(1)
series sp
sp(1)=price(1)*x(1)
ssp(1)=sp(1)/price(1)
if !k>=2 then
x(!k)=(0.000282998884054344)+(0.016628780888637)*u(!k)
price(!k)=1406.37+u(!k)
sp(!k)=price(!k)*x(!k)
ssp(!k)=sp(!k)/price(!k)
endif
next
smpl 1 !n
ssp.hist
按Save键保存并按Run运行该程序生成民生股票价格序列Price、随机游走u列、股价增量sp序列和股价变化率ssp序列，如图5.10所示：
图5.10 蒙特卡罗模拟所得各序列
运行以上程序得出上证综指股票价格变化率ssp序列的分布及其分布统计量。如图13 所示。从图5.11中ssp序列分布的统计值可知：偏度（Skewness）接近0，峰度 (Kurtosis)约为3.072062且JB统计量=2.507517，所对应的概率约为0.285430，可推断10,000次蒙特卡罗模拟出的上证综指股价变化率ssp序列值是服从均值Mean约为0.0002829，标准差Std.Dev为0.01662878 的正态分布。
图5.11上证综指股价变化率ssp序列值分布及其统计值
发表回复