陕 西 科 技 大 学
数理金融 课程设计任务书
学院 专业
题目： 南山铝业——股票价格及期权分析
课程设计从 2012 年 5 月 29 日起到 2012 年 6 月 5 日
1、课程设计的内容和要求（包括原始数据、技术要求、工作要求等）：
①绘制一个股价和时间的散点图（或拟合成曲线亦可）；②估计该股票价格
的二叉树模型中的增长因子u 和下降因子d （分别取∆t =1、3、5天）；并利
用该参数去求出15期（分别取∆t =1、3、5天）股票价格的二叉树模型；③如果用几何布朗运动模型来描述该股票价格，根据所给数据估计未知参
数μ和σ（以年为时间单位计算， 即∆t =1/365）；并用模型估计该股票在你的历史数据截止日前三个月的股票交易日的价格。并绘制该股票价格和
时间的散点图， 和①中的散点图对比（截止日前三个月）；④用几何布朗运动模型估计该在未来两个月的股票价格，并画出散点图。并分析其趋势
及相应统计特性；⑤取X 为历史数据中的平均值，r=0.05，利用以上数据
及相应估计参数分别计算历史数据截止日为起点的3个月后到期的欧式看
涨期权和看跌期权的价格。
根据所给的股票历史收盘价数据.---南山铝业
2、对课程设计成果的要求〔包括图表、实物等硬件要求〕：
3
指导教师：日期：
教研室主任：日期：
1 绘制股价和时间的图形
股票名称:南山铝业；
原始数据: 南山铝业1001个收盘价（2008年3月10日—2012年5月7日的历史行情数据）。
1.1 绘图原理
以时间为横坐标，以1001个股票实际收盘价为纵坐标，利用Matlab 编程实现画图。在画图过程中，用区间[1,1001]来代替时间作为横坐标。
1.2 绘制图形
图1-1 股价与时间图形
2.3 统计特征
股票价格平均值：mu=9.8821；方差：var= 7.2846；
2 二叉树模型对股票估值
2.1 计算原理
在风险中性的条件下，证券的预期收益率等于无风险利率r （r =10%），因此若该时段初证券价格为S ，则在小时间间隔∆t 段末的证券价格期望值为Se r ∆t 。参数p （价格上升的概率）、u 和d 的值必须满足这个要求，即：
Se r ∆t =pSu +(1-p ) Sd e r ∆t =pu +(1-p ) d （2-1） 二叉树模型也假设证券价格遵循几何布朗运动，在一个很小的时间段∆t 内证券价格变化的方差是S 2σ2∆t 。根据方差的定义，变量X 的方差等于X 2的期望值与X 期望值平方之差，因此：
S 2σ2∆t =pS 2u 2+(1-p ) S 2d 2-S 2[pu +(1-p ) d ]2
σ2∆t =pu 2+(1-p ) d 2-[pu +(1-p ) d ]2 （2-2）
从式（3-1）、（3-2）和u =1/d 可以求得，当∆t 很小时：
增长因子：u =e ；
下降因子：d =e -e r ∆t -d ； 价格上升的概率：p = ； u -d 价格下降的概率：q =1-p ；
2.2计算结果
（1）∆t =1时 u 1=1.1517 d 1=0.86825 p 1=0.46571 q 1=0.53429
（2）∆t =3时 u 2=1.2772 d 2=0.78295 p 2=0.44079 q 2=0.55921
（3）∆t =5时 u 3=1.3715 d 3=0.72914 p 3=0.42381 q 3=0.57619
2.3 二叉树模型的原理
二叉树模型首先把期权的有效期分为很多很小的时间间隔∆t ，并假设在每一个时间间隔∆t 内证券价格只有两种运动的可能：从开始的S 上升到原先的u 倍，即到达Su ；下降到原先的d 倍，即Sd 。其中，u >1，d
S
图2-1 ∆t 时间内资产价格的变动
相应地，期权价值也会有所不同，分别为f u 和f d 。在较大的时间间隔内，这种二值运
动的假设当然不符合实际，但是当时间间隔非常小的时候，比如在每个瞬间，资产价格只有这两个运动方向的假设是可以接受的。因此，二叉树模型实际上是在用大量离散的小幅度二值运动来模拟连续的资产价格运动。
2.4 二叉树模型估值结果
（1）∆t =1时股票价格15期的估值
表2-2 ∆t =1时股票价格15期的估值
（2）∆t =3时股票价格15期的估值
表2-3 ∆t =3时股票价格15期的估值
（3）∆t =5时股票价格15期的估值
表2-4 ∆t =5时股票价格15期的估值
2.5 二叉树模型估值结果分析
从表中可以看出, 从下往上, 股票的价格升高； 从上往下, 股票的价格降低, 这是由二叉树这种方法本身所决定的, 所以最后的股价要选择一个合适的价格, 一般在二叉树最后一层的中间取值较为合适。
3 几何布朗运动模型对股票估值
3.1 几何布朗运动模型中的参数
时间间隔∆t =1/365; 股票均值mu =7.2454; 标准差σ=0.5226;
3.2 几何布朗运动模型原理
布朗运动假设是现代资本市场理论的核心假设。现代资本市场理论认为证券期货价格具有随机性特征。这里的所谓随机性，是指数据的无记忆性，即过去数据不构成对未来数据的预测基础。同时不会出现惊人相似的反复。用
S (y)(0≤y ≤+∞) 表示y 时刻某证券的价格，若对任何非负实数y ，t 有：
（1）随机变量S (y+t)/S (y)独立于y 时刻及此前的所有价格；
（2）ln (S (y+t)/S (y))是均值为u ，t ，方差为t σ2的正态随机变量；
则称价格集服从漂移参数为u ，波动参数为σ的几何布朗运动。
如果证券价格遵循几何布朗运动，那么一旦u , σ的值确定了，影响未来价格概率分布的只是现在的价格，而与历史价格无关。涉及未来时刻t 以后的价格与当前价格比值的所有概率都与当前价格无关。比如一种证券在一个月后增长一倍的概率与该证券现在的价格是10还是25是没有关系的。
若随机变量为以y 为参数的对数正态随机变量，则
σ22E (Y ) =e μ+
（0）若已知证券的初始价格为S ，时刻t 价格的期望值仅依赖于几何布朗运动的漂移
参数和波动参数，即对于S (t ) 我们有
E [S (t )]=e t (μ+σ2
2) S (0)
用Δ表示一个小的时间增量，并假定，在每个Δ时间单位内，证券的价格或者以概率p 增长u 倍，或者以概率1-p 下跌d 倍，其中
布朗运动。
3.3 几何布朗运动模型估值结果
u =
e d =e -p =1(12当Δ取得越来越小时，价格的变化就越来越频繁，相应的价格集就近似为一个几何
（1）历史数据截止日前三个月的股票交易日的估值结果
6.[1**********]987 6.[1**********]105 6.[1**********]879 6.[1**********]046
6.[1**********]214 6.[1**********]716 6.[1**********]273 6.[**************].[1**********]943 5.[1**********]165 5.[1**********]165 6.[**************].[1**********]775 5.[1**********]332 6.[1**********]273 6.[**************].[1**********]426 6.[1**********]982 6.[1**********]426 6.[**************].[1**********]100 6.[1**********]416 6.[1**********]918 6.[**************].[1**********]081 6.[1**********]308 6.[1**********]751 6.[**************].[1**********]918 6.[1**********]140 6.[1**********]308 6.[**************].[1**********]746 6.[1**********]914 6.[1**********]470 6.[**************].[1**********]027 7.[1**********]579 7.[1**********]746 7.[**************].[1**********]968 7.[1**********]800 7.[1**********]130 7.[**************].[1**********]963 7.[1**********]342 7.[1**********]175 7.[**************].[1**********]121 7.[1**********]840 7.[1**********]781 7.[**************].[1**********]510 7.[1**********]229 7.[1**********]672 7.[**************].[1**********]451 7.[1**********]460 7.[1**********]125 7.[**************].[1**********]963 7.[1**********]012 7.[1**********]850 7.[**************].[1**********]687 7.[1**********]465 6.[1**********]194 6.[**************].[1**********]588 6.[1**********]194 6.[1**********]194 6.[**************].[1**********]805 6.[1**********]918 7.[1**********]524 7.[**************].[1**********]189 6.[1**********]249 7.[1**********]687 7.[**************].[1**********]239 7.[1**********]071 7.[1**********]741 7.[**************].[1**********]677 7.[1**********]736 7.[1**********]564 7.[**************].[1**********]510 7.[1**********]786
均值：mu=6.9453; 方差：var=0.2510;
（2）90天的股票交易日的真实数据
7.92 7.85 7.89 7.9 7.69 7.82 7.56 7.517.63 7.6 7.49 7.46 7.24 7.37 7.35 7.31
7.12 7.2 7.14 7.19 7.19 7 7 7.19 7.44 7.46 7.44 7.61 7.76 7.53 7.75 7.68 7.62
7.98 7.98 8 7.96 7.85 8.02 8.1 7.97 7.86 7.87 7.91 7.88 7.53 7.48 7.5 7.38 7.35
7.34 7.33 7.36 7.22 7.27 7.26 7.3 7.33 6.99 7.11 7.14 7.12 6.96 7.15 7.11 7.27
7.24 7.12 7.21 6.55 6.64 6.85 6.81 6.85 6.6 6.34 6.21 6.25 6.33 6.24 6.24 6.22
6.38 6.346.38 6.47 6.5 6.53 6.37 6.63
均值mu=7.2454; 方差var=0.2731;
(3) 估值结果与实际值的对比图
图3-1估值结果与实际值的对比图
3.4 几何布朗运动模型估值结果分析及与实际图形对比分析
从上图可以看出估值结果与实际股价的变化趋势大致相同，只是相对于真实值会较高或较低一些
4 用几何布朗运动模型估计未来股价
4.1模型参数
时间间隔∆t =1/365; 股票均值μ=9.8821; 标准差σ=2.6990;
4.2估值结果
(1)估值结果数据
8.[1**********]408 8.[1**********]851 7.[1**********]711
8.[1**********]327 11.[1**********]88 12.[1**********]92
7.[1**********]512 11.[1**********]67 10.[1**********]66 12.[1**********]22 10.[1**********]66 10.[1**********]72 14.[1**********]89 13.[1**********]64 18.[1**********]09 14.[1**********]97 18.[1**********]53 17.[1**********]73 23.[1**********]07 22.[1**********]16 14.[1**********]61 19.[1**********]31 11.[1**********]96 21.[1**********]80 23.[1**********]05 16.[1**********]47 16.[1**********]48 25.[1**********]69 17.[1**********]83 29.[1**********]61 25.[1**********]24 19.[1**********]31 22.[1**********]28 18.[1**********]30 21.[1**********]39 36.[1**********]07 41.[1**********]95 14.[1**********]64 27.[1**********]80 47.[1**********]47 24.[1**********]45 37.[1**********]98 42.[1**********]56 22.[1**********]27 30.[1**********]11
Mu=19.7221; var=92.8650
(2)估值结果与时间的趋势图
8.[**************].[**************]0.[1**********]628.[**************]7.[**************].[**************].[**************].[**************].[**************].[**************].[**************].[**************].[**************].[**************].[1**********]03
图4-1 几何布朗运动模型股价估值趋势图
4.3结果分析
由上图可以看出几何布朗运动模型预测的股价整体呈上升趋势，价格在随着时间上下波动。
5 计算股票欧式期权的价格
5.1蒙特卡洛模拟法原理
蒙特卡洛模拟方法亦称随机模拟方法，其基本思想是，求解科学、工程技术和经济金融等方面的问题。首先建立一个概率模型后随机过程，使其参数等于问题的解；然后通过对模型或过程的观察计算所求的统计特征，最后给出所求问题的近似值，解的精度可用估计值的标准差表示。
应用此方法求解可以分为两类:确定性问题和随机性问题。 解题步骤如下：
（1）根据提出的问题构造一个简单、适用的概率模型或随机模型，使问题的解对应于该模型中随机变量的某些特征（如概率、均值和方差等），所构造的模型在主要特征参量方面要与实际问题或系统相一致
（2）根据模型中各个随机变量的分布，在计算机上产生随机数，实现一次模拟过程所需的足够数量的随机数。通常先产生均匀分布的随机数，然后生成服从某一分布的随机数，方可进行随机模拟试验。
（3）根据概率模型的特点和随机变量的分布特性，设计和选取合适的抽样方法，并对每个随机变量进行抽样（包括直接抽样、分层抽样、相关抽样、重要抽样等）。
（4）按照所建立的模型进行仿真试验、计算，求出问题的随机解。 （5）统计分析模拟试验结果，给出问题的概率解以及解的精度估计。 5.2欧式看涨期权计算公式
c =max{0，S (0)e
(r -σ2/2)T +-X }
5.3 欧式看涨期权计算结果
(1) 下表分别为一般技术(MC)、对偶变量技术（MCAV ）、控制变量技术（MCCV ）的均值、方差的计算结果
解析解：BLprice=3.5467
表5-1 蒙特卡洛模拟法计算结果表
(2) 将蒙特卡洛模拟法计算结果用图形表示
图5-1 B-S公式和三种MC 技术方法算出的期权价格波动趋势图
图5-2 三种MC 法的置信区间变化图
图5-3 控制变量法的期权价格波动图
5.4欧式看跌期权计算公式
p =max{0，X -S (0)e (r -σ
2
/2)T +}
5.5欧式看跌期权计算结果
(1) 下表分别为一般技术(MC)、对偶变量技术（MCAV ）、控制变量技术（MCCV ）的均值、方差的计算结果
解析解：Put=5.3861; Call=3.5467;
表5-2 蒙特卡洛模拟法计算结果表
(2) 将蒙特卡洛模拟法计算结果用图形表示
图5-4 B-S公式和三种MC 技术方法算出的期权价格波动趋势图
图5-5 三种MC 法的置信区间变化图
图5-6 控制变量法的期权价格波动图
5.6 蒙特卡罗方法
蒙特卡罗方法的实质是模拟标的资产价格的随机运动，预测期权的平均回报，并由此得到期权价格的一个概率解。
蒙特卡罗模拟的主要优点包括：
（1）在大多数情况下，人们可以很直接地应用蒙特卡罗模拟方法，而无需对期权定价模型有深刻的理解，所用的数学知识也很基本；为了获得更精确的答案，只需要进
行更多的模拟；无需太多工作就可以转换模型。以上这些优点使得蒙特卡罗方法成为一
个相当广泛和强大的期权定价技术。
（2）蒙特卡罗模拟的适用情形相当广泛，其中包括： 期权的回报仅仅取决于标的变量的最终价值的情况；
期权的回报依赖于标的变量所遵循的路径，即路径依赖的情形；
期权的回报取决于多个标的变量的情况，尤其当随机变量的数量增加时，蒙特卡罗模拟的运算时间近似为线性增长而不象其他方法那样以指数增长，因此该方法对依赖三种以上风险资产的多变量期权模型很有竞争力。
因此，蒙特卡罗模拟可以适用于复杂随机过程和复杂终值的计算，在运算过程中蒙特卡罗模拟还能给出估计值的标准误差，这也是该方法的优点之一。
蒙特卡罗模拟的缺点主要是：
（1）只能为欧式期权定价，难以处理提前执行的情形。尝试使用蒙特卡罗模拟技巧来为美式期权定价，成为近年来这个领域的发展方向之一。
（2）为了达到一定的精确度，一般需要大量的模拟运算。尤其在处理三个以下的变量时，蒙特卡罗模拟相对于其他方法来说偏慢。
6 总结及体会
通过本次课程设计，加深了对于数理金融的理解，使得平日学习的理论课程得以在实践中应用。在设计的过程中，加深了对二叉树模型和欧式看涨，看跌期权的理解，特别是几何布朗运动的实验，学习到了一种预测股票未来价格走向的方法。在课程设计中，如何去发现问题、解决问题以及利用网络资源解决问题是最重要的。实验中应用不当就出错而且错误很费神，理论与实际差距的确很大，要想作出实际实用的东西来还需多多练习，不断改进，充分发挥所学知识，另外也应加强视野的开拓。设计过程中的解决问题的方法，让我明白了如何学习会更有效。如何学习才不会耽误太多的时间。设这次课程设计增加了自己动手的能力，也学会了解决问题的一般方法：向老师、同学请教，借助书籍网络等等。
此次课程设计，我也发现自己还有很多需要改进的地方。首先，在问题分析上还需要进一步加强，提高自己分析问题的能力，要学会挖掘题中的深层含义；其次，在编写
Matlab 程序时，需要有清晰思路，在编写Matlab 程序之前一定要学会写每一个程序的
流程；最后，要更加勤奋地学习理论知识，并且将理论与实践相结合，使自己能更加熟练地运用理论知识来解决实际问题。
通过本次课程设计, 我明白了一个道理:无论做什么事情, 都必需养成严谨, 认真, 善思的工作作风遇到问题最好的办法就是请教别人，因为每个人掌握的情况都不一样，一个人不可能做到处处都懂。
参 考 文 献
[1] 张元萍. 数理金融[M].北京：中国金融出版社，2004.
[2] 周博谢, 张宪海编著.Matlab 科学计算[M].北京：机械工业出版社，2010.
附录：
程序一
sp=[1001个收盘价] ; x=1001:-1:1; mu=mean(sp); var=var(sp); sigma=sqrt(var); plot(x,sp);
程序二
clear;clc; S0=7.92; r=0.1; sigma=2.6990;
dt=[1/365,3/365,5/365]; w=zeros(16,15,3); for n=1:3
u(n)=exp(sigma*sqrt(dt(n))); d(n)=1/u(n);
p(n)=(exp(r*dt(n))-d(n))/(u(n)-d(n)); q(n)=1-p(n); end for i=1:15 for j=i+1:-1:1 for n=1:3
w(i-j+2,i,n)=S0*u(n)^(j-1)*d(n)^(i-(j-1)); end end end
w1(:,:)=w(:,:,1); w2(:,:)=w(:,:,2); w3(:,:)=w(:,:,3);
程序三
clear;clc; d=[近90天收盘价]; x=90:-1:1; plot(x,d); hold on for i=1:90 l(i)=d(91-i); end l(i);
T=1/365;
mu=9.8821;
sigma=2.6990;
S=zeros(1,90);
nuT=(mu-0.5*sigma^2)*T;
siT=sigma*sqrt(T);
w=randn(1,1);
for i=1:90
S(i)=l(i)*exp(nuT+siT*w);
end
E=mean(S);
D=var(S);
x=1:90;
plot(x,S,'-.' );
程序四
clear;clc;
T=1/365;
mu= 9.8821;
sigma=2.6990;
S0=7.92;
S=zeros(10,60);
for i=1:60
nuT(i)=(mu-0.5*sigma^2)*T*i;
siT(i)=sigma*sqrt(T*i);
w=randn(10,1);
S(:,i)=S0*exp(nuT(i)+siT(i)*w);
[Price(i), VarPrice(i), CI(i,:)] = normfit(S(:,i));
end
E=mean(Price);
D=var(Price);
x=1:60;
plot(x,Price);
程序五
看涨期权蒙特卡洛模拟定价程序：
clear all ;
close all ;
clc;
NbSimu = 20000 : 20000: 200000;
NbPoints = numel(NbSimu);
ResAV = zeros(NbPoints,1);
ResMC = zeros(NbPoints,1);
ResMCCV = zeros(NbPoints,1);
CIAV = zeros(NbPoints,2);
CIMC = zeros(NbPoints,2);
CIMCCV = zeros(NbPoints,2);
VarPriceAV = zeros(NbPoints,1);
VarPriceMC = zeros(NbPoints,1);
VarPriceCV = zeros(NbPoints,1);
BLprice = blsprice(7.92,9.8821,0.05,0.25,2.6990,0) ;
for i = 1 :numel(NbSimu)
disp(['Computing Iteration N? num2str(i)' ' num2str(NbSimu(i))' 'Simulated Paths' ]);
[ResAV(i),VarPriceAV(i),CIAV(i,:)] =
BlsMCAV(7.92,9.8821,0.05,0.25,2.6990,NbSimu(i));
[ResMC(i),VarPriceMC(i),CIMC(i,:)] =
BlsMC(7.92,9.8821,0.05,0.25,2.6990,NbSimu(i));
[ResMCCV(i),VarPriceCV(i),CIMCCV(i,:)] =
BlsMCCV(7.92,9.8821,0.05,0.25,2.6990,NbSimu(i),1000);
end ;
h= figure; plot(NbSimu,[repmat(BLprice,NbPoints,1) ResAV ResMC ResMCCV]); legend({'Black Scholes' 'Antithetic' , 'Simple Monte Carlo', 'Monte Carlo with Control Variable'})
set(h,'WindowStyle' , 'Docked' );
grid on ;
xlabel('Number of simuation');
Ylabel('Vanilla option Price');
xlim([NbSimu(1) NbSimu(end)]);
h= figure;
plot(NbSimu,[diff(CIAV,1,2) diff(CIMC,1,2) diff(CIMCCV,1,2) ]);
legend({'Antithetic' , 'Simple Monte Carlo', 'Monte Carlo with Control Variable' , });
grid on ;xlim([NbSimu(1) NbSimu(end)]);
title('Comparing different Monte Carlo methods')
set(h,'WindowStyle' , 'Docked' );
xlabel('Number of simuation');
Ylabel('absolute confidence interval');
h= figure;
FillBetween(NbSimu',CIMCCV(:,1),CIMCCV(:,2),'g' ,0.3);
hold on ;
plot(NbSimu,ResMCCV,'LineWidth' ,1.5, 'Color' , 'r' );
plot(NbSimu,repmat(BLprice,NbPoints,1),'b' );
xlim([NbSimu(1) NbSimu(end)])
set(h,'WindowStyle' , 'Docked' );
title('Monte Carlo with Control Variable');
xlabel('Number of simulation');
Ylabel('Vanilla option Price');
grid on ;
看跌期权蒙特卡洛模拟定价程序：
clear all ;
close all ;
clc;
NbSimu =20000:20000:200000;
NbPoints = numel(NbSimu);
ResAV = zeros(NbPoints,1);
ResMC = zeros(NbPoints,1);
ResMCCV = zeros(NbPoints,1);
CIAV = zeros(NbPoints,2);
CIMC = zeros(NbPoints,2);
CIMCCV = zeros(NbPoints,2);
VarPriceAV = zeros(NbPoints,1);
VarPriceMC = zeros(NbPoints,1);
VarPriceCV = zeros(NbPoints,1);
[call,put] = blsprice(7.92,9.8821,0.05,0.25,2.6990,0);
for i = 1 :numel(NbSimu)
disp(['Computing Iteration N? num2str(i)' ' num2str(NbSimu(i))' 'Simulated Paths' ]);
[ResAV(i),VarPriceAV(i),CIAV(i,:)] =
BlsMCAV(7.92,9.8821,0.05,0.25,2.6990,NbSimu(i));
[ResMC(i),VarPriceMC(i),CIMC(i,:)] =
BlsMC(7.92,9.8821,0.05,0.25,2.6990,NbSimu(i));
[ResMCCV(i),VarPriceCV(i),CIMCCV(i,:)] =
BlsMCCV(7.92,9.8821,0.05,0.25,2.6990,NbSimu(i),1000);
end ;
h= figure;
plot(NbSimu,[repmat(put,NbPoints,1) ResAV ResMC ResMCCV]);
legend({'Black Scholes' 'Antithetic' , 'Simple Monte Carlo', 'Monte Carlo with Control Variable'})
set(h,'WindowStyle' , 'Docked' );
grid on ;
xlabel('Number of simuation');
Ylabel('Vanilla option Price');
xlim([NbSimu(1) NbSimu(end)]);
h= figure;
plot(NbSimu,[diff(CIAV,1,2) diff(CIMC,1,2) diff(CIMCCV,1,2) ]);
legend({'Antithetic' , 'Simple Monte Carlo', 'Monte Carlo with Control
Variable' , });
grid on ;xlim([NbSimu(1) NbSimu(end)]);
title('Comparing different Monte Carlo methods') set(h,'WindowStyle' , 'Docked' );
xlabel('Number of simuation');
Ylabel('absolute confidence interval');
h= figure;
FillBetween(NbSimu',CIMCCV(:,1),CIMCCV(:,2),'g' ,0.3); hold on ;
plot(NbSimu,ResMCCV,'LineWidth' ,1.5, 'Color' , 'r' ); plot(NbSimu,repmat(put,NbPoints,1),'b' );
xlim([NbSimu(1) NbSimu(end)])
set(h,'WindowStyle' , 'Docked' );
title('Monte Carlo with Control Variable');
xlabel('Number of simulation');
Ylabel('Vanilla option Price');
grid on ;
陕 西 科 技 大 学
数理金融 课程设计任务书
学院 专业
题目： 南山铝业——股票价格及期权分析
课程设计从 2012 年 5 月 29 日起到 2012 年 6 月 5 日
1、课程设计的内容和要求（包括原始数据、技术要求、工作要求等）：
①绘制一个股价和时间的散点图（或拟合成曲线亦可）；②估计该股票价格
的二叉树模型中的增长因子u 和下降因子d （分别取∆t =1、3、5天）；并利
用该参数去求出15期（分别取∆t =1、3、5天）股票价格的二叉树模型；③如果用几何布朗运动模型来描述该股票价格，根据所给数据估计未知参
数μ和σ（以年为时间单位计算， 即∆t =1/365）；并用模型估计该股票在你的历史数据截止日前三个月的股票交易日的价格。并绘制该股票价格和
时间的散点图， 和①中的散点图对比（截止日前三个月）；④用几何布朗运动模型估计该在未来两个月的股票价格，并画出散点图。并分析其趋势
及相应统计特性；⑤取X 为历史数据中的平均值，r=0.05，利用以上数据
及相应估计参数分别计算历史数据截止日为起点的3个月后到期的欧式看
涨期权和看跌期权的价格。
根据所给的股票历史收盘价数据.---南山铝业
2、对课程设计成果的要求〔包括图表、实物等硬件要求〕：
3
指导教师：日期：
教研室主任：日期：
1 绘制股价和时间的图形
股票名称:南山铝业；
原始数据: 南山铝业1001个收盘价（2008年3月10日—2012年5月7日的历史行情数据）。
1.1 绘图原理
以时间为横坐标，以1001个股票实际收盘价为纵坐标，利用Matlab 编程实现画图。在画图过程中，用区间[1,1001]来代替时间作为横坐标。
1.2 绘制图形
图1-1 股价与时间图形
2.3 统计特征
股票价格平均值：mu=9.8821；方差：var= 7.2846；
2 二叉树模型对股票估值
2.1 计算原理
在风险中性的条件下，证券的预期收益率等于无风险利率r （r =10%），因此若该时段初证券价格为S ，则在小时间间隔∆t 段末的证券价格期望值为Se r ∆t 。参数p （价格上升的概率）、u 和d 的值必须满足这个要求，即：
Se r ∆t =pSu +(1-p ) Sd e r ∆t =pu +(1-p ) d （2-1） 二叉树模型也假设证券价格遵循几何布朗运动，在一个很小的时间段∆t 内证券价格变化的方差是S 2σ2∆t 。根据方差的定义，变量X 的方差等于X 2的期望值与X 期望值平方之差，因此：
S 2σ2∆t =pS 2u 2+(1-p ) S 2d 2-S 2[pu +(1-p ) d ]2
σ2∆t =pu 2+(1-p ) d 2-[pu +(1-p ) d ]2 （2-2）
从式（3-1）、（3-2）和u =1/d 可以求得，当∆t 很小时：
增长因子：u =e ；
下降因子：d =e -e r ∆t -d ； 价格上升的概率：p = ； u -d 价格下降的概率：q =1-p ；
2.2计算结果
（1）∆t =1时 u 1=1.1517 d 1=0.86825 p 1=0.46571 q 1=0.53429
（2）∆t =3时 u 2=1.2772 d 2=0.78295 p 2=0.44079 q 2=0.55921
（3）∆t =5时 u 3=1.3715 d 3=0.72914 p 3=0.42381 q 3=0.57619
2.3 二叉树模型的原理
二叉树模型首先把期权的有效期分为很多很小的时间间隔∆t ，并假设在每一个时间间隔∆t 内证券价格只有两种运动的可能：从开始的S 上升到原先的u 倍，即到达Su ；下降到原先的d 倍，即Sd 。其中，u >1，d
S
图2-1 ∆t 时间内资产价格的变动
相应地，期权价值也会有所不同，分别为f u 和f d 。在较大的时间间隔内，这种二值运
动的假设当然不符合实际，但是当时间间隔非常小的时候，比如在每个瞬间，资产价格只有这两个运动方向的假设是可以接受的。因此，二叉树模型实际上是在用大量离散的小幅度二值运动来模拟连续的资产价格运动。
2.4 二叉树模型估值结果
（1）∆t =1时股票价格15期的估值
表2-2 ∆t =1时股票价格15期的估值
（2）∆t =3时股票价格15期的估值
表2-3 ∆t =3时股票价格15期的估值
（3）∆t =5时股票价格15期的估值
表2-4 ∆t =5时股票价格15期的估值
2.5 二叉树模型估值结果分析
从表中可以看出, 从下往上, 股票的价格升高； 从上往下, 股票的价格降低, 这是由二叉树这种方法本身所决定的, 所以最后的股价要选择一个合适的价格, 一般在二叉树最后一层的中间取值较为合适。
3 几何布朗运动模型对股票估值
3.1 几何布朗运动模型中的参数
时间间隔∆t =1/365; 股票均值mu =7.2454; 标准差σ=0.5226;
3.2 几何布朗运动模型原理
布朗运动假设是现代资本市场理论的核心假设。现代资本市场理论认为证券期货价格具有随机性特征。这里的所谓随机性，是指数据的无记忆性，即过去数据不构成对未来数据的预测基础。同时不会出现惊人相似的反复。用
S (y)(0≤y ≤+∞) 表示y 时刻某证券的价格，若对任何非负实数y ，t 有：
（1）随机变量S (y+t)/S (y)独立于y 时刻及此前的所有价格；
（2）ln (S (y+t)/S (y))是均值为u ，t ，方差为t σ2的正态随机变量；
则称价格集服从漂移参数为u ，波动参数为σ的几何布朗运动。
如果证券价格遵循几何布朗运动，那么一旦u , σ的值确定了，影响未来价格概率分布的只是现在的价格，而与历史价格无关。涉及未来时刻t 以后的价格与当前价格比值的所有概率都与当前价格无关。比如一种证券在一个月后增长一倍的概率与该证券现在的价格是10还是25是没有关系的。
若随机变量为以y 为参数的对数正态随机变量，则
σ22E (Y ) =e μ+
（0）若已知证券的初始价格为S ，时刻t 价格的期望值仅依赖于几何布朗运动的漂移
参数和波动参数，即对于S (t ) 我们有
E [S (t )]=e t (μ+σ2
2) S (0)
用Δ表示一个小的时间增量，并假定，在每个Δ时间单位内，证券的价格或者以概率p 增长u 倍，或者以概率1-p 下跌d 倍，其中
布朗运动。
3.3 几何布朗运动模型估值结果
u =
e d =e -p =1(12当Δ取得越来越小时，价格的变化就越来越频繁，相应的价格集就近似为一个几何
（1）历史数据截止日前三个月的股票交易日的估值结果
6.[1**********]987 6.[1**********]105 6.[1**********]879 6.[1**********]046
6.[1**********]214 6.[1**********]716 6.[1**********]273 6.[**************].[1**********]943 5.[1**********]165 5.[1**********]165 6.[**************].[1**********]775 5.[1**********]332 6.[1**********]273 6.[**************].[1**********]426 6.[1**********]982 6.[1**********]426 6.[**************].[1**********]100 6.[1**********]416 6.[1**********]918 6.[**************].[1**********]081 6.[1**********]308 6.[1**********]751 6.[**************].[1**********]918 6.[1**********]140 6.[1**********]308 6.[**************].[1**********]746 6.[1**********]914 6.[1**********]470 6.[**************].[1**********]027 7.[1**********]579 7.[1**********]746 7.[**************].[1**********]968 7.[1**********]800 7.[1**********]130 7.[**************].[1**********]963 7.[1**********]342 7.[1**********]175 7.[**************].[1**********]121 7.[1**********]840 7.[1**********]781 7.[**************].[1**********]510 7.[1**********]229 7.[1**********]672 7.[**************].[1**********]451 7.[1**********]460 7.[1**********]125 7.[**************].[1**********]963 7.[1**********]012 7.[1**********]850 7.[**************].[1**********]687 7.[1**********]465 6.[1**********]194 6.[**************].[1**********]588 6.[1**********]194 6.[1**********]194 6.[**************].[1**********]805 6.[1**********]918 7.[1**********]524 7.[**************].[1**********]189 6.[1**********]249 7.[1**********]687 7.[**************].[1**********]239 7.[1**********]071 7.[1**********]741 7.[**************].[1**********]677 7.[1**********]736 7.[1**********]564 7.[**************].[1**********]510 7.[1**********]786
均值：mu=6.9453; 方差：var=0.2510;
（2）90天的股票交易日的真实数据
7.92 7.85 7.89 7.9 7.69 7.82 7.56 7.517.63 7.6 7.49 7.46 7.24 7.37 7.35 7.31
7.12 7.2 7.14 7.19 7.19 7 7 7.19 7.44 7.46 7.44 7.61 7.76 7.53 7.75 7.68 7.62
7.98 7.98 8 7.96 7.85 8.02 8.1 7.97 7.86 7.87 7.91 7.88 7.53 7.48 7.5 7.38 7.35
7.34 7.33 7.36 7.22 7.27 7.26 7.3 7.33 6.99 7.11 7.14 7.12 6.96 7.15 7.11 7.27
7.24 7.12 7.21 6.55 6.64 6.85 6.81 6.85 6.6 6.34 6.21 6.25 6.33 6.24 6.24 6.22
6.38 6.346.38 6.47 6.5 6.53 6.37 6.63
均值mu=7.2454; 方差var=0.2731;
(3) 估值结果与实际值的对比图
图3-1估值结果与实际值的对比图
3.4 几何布朗运动模型估值结果分析及与实际图形对比分析
从上图可以看出估值结果与实际股价的变化趋势大致相同，只是相对于真实值会较高或较低一些
4 用几何布朗运动模型估计未来股价
4.1模型参数
时间间隔∆t =1/365; 股票均值μ=9.8821; 标准差σ=2.6990;
4.2估值结果
(1)估值结果数据
8.[1**********]408 8.[1**********]851 7.[1**********]711
8.[1**********]327 11.[1**********]88 12.[1**********]92
7.[1**********]512 11.[1**********]67 10.[1**********]66 12.[1**********]22 10.[1**********]66 10.[1**********]72 14.[1**********]89 13.[1**********]64 18.[1**********]09 14.[1**********]97 18.[1**********]53 17.[1**********]73 23.[1**********]07 22.[1**********]16 14.[1**********]61 19.[1**********]31 11.[1**********]96 21.[1**********]80 23.[1**********]05 16.[1**********]47 16.[1**********]48 25.[1**********]69 17.[1**********]83 29.[1**********]61 25.[1**********]24 19.[1**********]31 22.[1**********]28 18.[1**********]30 21.[1**********]39 36.[1**********]07 41.[1**********]95 14.[1**********]64 27.[1**********]80 47.[1**********]47 24.[1**********]45 37.[1**********]98 42.[1**********]56 22.[1**********]27 30.[1**********]11
Mu=19.7221; var=92.8650
(2)估值结果与时间的趋势图
8.[**************].[**************]0.[1**********]628.[**************]7.[**************].[**************].[**************].[**************].[**************].[**************].[**************].[**************].[**************].[**************].[1**********]03
图4-1 几何布朗运动模型股价估值趋势图
4.3结果分析
由上图可以看出几何布朗运动模型预测的股价整体呈上升趋势，价格在随着时间上下波动。
5 计算股票欧式期权的价格
5.1蒙特卡洛模拟法原理
蒙特卡洛模拟方法亦称随机模拟方法，其基本思想是，求解科学、工程技术和经济金融等方面的问题。首先建立一个概率模型后随机过程，使其参数等于问题的解；然后通过对模型或过程的观察计算所求的统计特征，最后给出所求问题的近似值，解的精度可用估计值的标准差表示。
应用此方法求解可以分为两类:确定性问题和随机性问题。 解题步骤如下：
（1）根据提出的问题构造一个简单、适用的概率模型或随机模型，使问题的解对应于该模型中随机变量的某些特征（如概率、均值和方差等），所构造的模型在主要特征参量方面要与实际问题或系统相一致
（2）根据模型中各个随机变量的分布，在计算机上产生随机数，实现一次模拟过程所需的足够数量的随机数。通常先产生均匀分布的随机数，然后生成服从某一分布的随机数，方可进行随机模拟试验。
（3）根据概率模型的特点和随机变量的分布特性，设计和选取合适的抽样方法，并对每个随机变量进行抽样（包括直接抽样、分层抽样、相关抽样、重要抽样等）。
（4）按照所建立的模型进行仿真试验、计算，求出问题的随机解。 （5）统计分析模拟试验结果，给出问题的概率解以及解的精度估计。 5.2欧式看涨期权计算公式
c =max{0，S (0)e
(r -σ2/2)T +-X }
5.3 欧式看涨期权计算结果
(1) 下表分别为一般技术(MC)、对偶变量技术（MCAV ）、控制变量技术（MCCV ）的均值、方差的计算结果
解析解：BLprice=3.5467
表5-1 蒙特卡洛模拟法计算结果表
(2) 将蒙特卡洛模拟法计算结果用图形表示
图5-1 B-S公式和三种MC 技术方法算出的期权价格波动趋势图
图5-2 三种MC 法的置信区间变化图
图5-3 控制变量法的期权价格波动图
5.4欧式看跌期权计算公式
p =max{0，X -S (0)e (r -σ
2
/2)T +}
5.5欧式看跌期权计算结果
(1) 下表分别为一般技术(MC)、对偶变量技术（MCAV ）、控制变量技术（MCCV ）的均值、方差的计算结果
解析解：Put=5.3861; Call=3.5467;
表5-2 蒙特卡洛模拟法计算结果表
(2) 将蒙特卡洛模拟法计算结果用图形表示
图5-4 B-S公式和三种MC 技术方法算出的期权价格波动趋势图
图5-5 三种MC 法的置信区间变化图
图5-6 控制变量法的期权价格波动图
5.6 蒙特卡罗方法
蒙特卡罗方法的实质是模拟标的资产价格的随机运动，预测期权的平均回报，并由此得到期权价格的一个概率解。
蒙特卡罗模拟的主要优点包括：
（1）在大多数情况下，人们可以很直接地应用蒙特卡罗模拟方法，而无需对期权定价模型有深刻的理解，所用的数学知识也很基本；为了获得更精确的答案，只需要进
行更多的模拟；无需太多工作就可以转换模型。以上这些优点使得蒙特卡罗方法成为一
个相当广泛和强大的期权定价技术。
（2）蒙特卡罗模拟的适用情形相当广泛，其中包括： 期权的回报仅仅取决于标的变量的最终价值的情况；
期权的回报依赖于标的变量所遵循的路径，即路径依赖的情形；
期权的回报取决于多个标的变量的情况，尤其当随机变量的数量增加时，蒙特卡罗模拟的运算时间近似为线性增长而不象其他方法那样以指数增长，因此该方法对依赖三种以上风险资产的多变量期权模型很有竞争力。
因此，蒙特卡罗模拟可以适用于复杂随机过程和复杂终值的计算，在运算过程中蒙特卡罗模拟还能给出估计值的标准误差，这也是该方法的优点之一。
蒙特卡罗模拟的缺点主要是：
（1）只能为欧式期权定价，难以处理提前执行的情形。尝试使用蒙特卡罗模拟技巧来为美式期权定价，成为近年来这个领域的发展方向之一。
（2）为了达到一定的精确度，一般需要大量的模拟运算。尤其在处理三个以下的变量时，蒙特卡罗模拟相对于其他方法来说偏慢。
6 总结及体会
通过本次课程设计，加深了对于数理金融的理解，使得平日学习的理论课程得以在实践中应用。在设计的过程中，加深了对二叉树模型和欧式看涨，看跌期权的理解，特别是几何布朗运动的实验，学习到了一种预测股票未来价格走向的方法。在课程设计中，如何去发现问题、解决问题以及利用网络资源解决问题是最重要的。实验中应用不当就出错而且错误很费神，理论与实际差距的确很大，要想作出实际实用的东西来还需多多练习，不断改进，充分发挥所学知识，另外也应加强视野的开拓。设计过程中的解决问题的方法，让我明白了如何学习会更有效。如何学习才不会耽误太多的时间。设这次课程设计增加了自己动手的能力，也学会了解决问题的一般方法：向老师、同学请教，借助书籍网络等等。
此次课程设计，我也发现自己还有很多需要改进的地方。首先，在问题分析上还需要进一步加强，提高自己分析问题的能力，要学会挖掘题中的深层含义；其次，在编写
Matlab 程序时，需要有清晰思路，在编写Matlab 程序之前一定要学会写每一个程序的
流程；最后，要更加勤奋地学习理论知识，并且将理论与实践相结合，使自己能更加熟练地运用理论知识来解决实际问题。
通过本次课程设计, 我明白了一个道理:无论做什么事情, 都必需养成严谨, 认真, 善思的工作作风遇到问题最好的办法就是请教别人，因为每个人掌握的情况都不一样，一个人不可能做到处处都懂。
参 考 文 献
[1] 张元萍. 数理金融[M].北京：中国金融出版社，2004.
[2] 周博谢, 张宪海编著.Matlab 科学计算[M].北京：机械工业出版社，2010.
附录：
程序一
sp=[1001个收盘价] ; x=1001:-1:1; mu=mean(sp); var=var(sp); sigma=sqrt(var); plot(x,sp);
程序二
clear;clc; S0=7.92; r=0.1; sigma=2.6990;
dt=[1/365,3/365,5/365]; w=zeros(16,15,3); for n=1:3
u(n)=exp(sigma*sqrt(dt(n))); d(n)=1/u(n);
p(n)=(exp(r*dt(n))-d(n))/(u(n)-d(n)); q(n)=1-p(n); end for i=1:15 for j=i+1:-1:1 for n=1:3
w(i-j+2,i,n)=S0*u(n)^(j-1)*d(n)^(i-(j-1)); end end end
w1(:,:)=w(:,:,1); w2(:,:)=w(:,:,2); w3(:,:)=w(:,:,3);
程序三
clear;clc; d=[近90天收盘价]; x=90:-1:1; plot(x,d); hold on for i=1:90 l(i)=d(91-i); end l(i);
T=1/365;
mu=9.8821;
sigma=2.6990;
S=zeros(1,90);
nuT=(mu-0.5*sigma^2)*T;
siT=sigma*sqrt(T);
w=randn(1,1);
for i=1:90
S(i)=l(i)*exp(nuT+siT*w);
end
E=mean(S);
D=var(S);
x=1:90;
plot(x,S,'-.' );
程序四
clear;clc;
T=1/365;
mu= 9.8821;
sigma=2.6990;
S0=7.92;
S=zeros(10,60);
for i=1:60
nuT(i)=(mu-0.5*sigma^2)*T*i;
siT(i)=sigma*sqrt(T*i);
w=randn(10,1);
S(:,i)=S0*exp(nuT(i)+siT(i)*w);
[Price(i), VarPrice(i), CI(i,:)] = normfit(S(:,i));
end
E=mean(Price);
D=var(Price);
x=1:60;
plot(x,Price);
程序五
看涨期权蒙特卡洛模拟定价程序：
clear all ;
close all ;
clc;
NbSimu = 20000 : 20000: 200000;
NbPoints = numel(NbSimu);
ResAV = zeros(NbPoints,1);
ResMC = zeros(NbPoints,1);
ResMCCV = zeros(NbPoints,1);
CIAV = zeros(NbPoints,2);
CIMC = zeros(NbPoints,2);
CIMCCV = zeros(NbPoints,2);
VarPriceAV = zeros(NbPoints,1);
VarPriceMC = zeros(NbPoints,1);
VarPriceCV = zeros(NbPoints,1);
BLprice = blsprice(7.92,9.8821,0.05,0.25,2.6990,0) ;
for i = 1 :numel(NbSimu)
disp(['Computing Iteration N? num2str(i)' ' num2str(NbSimu(i))' 'Simulated Paths' ]);
[ResAV(i),VarPriceAV(i),CIAV(i,:)] =
BlsMCAV(7.92,9.8821,0.05,0.25,2.6990,NbSimu(i));
[ResMC(i),VarPriceMC(i),CIMC(i,:)] =
BlsMC(7.92,9.8821,0.05,0.25,2.6990,NbSimu(i));
[ResMCCV(i),VarPriceCV(i),CIMCCV(i,:)] =
BlsMCCV(7.92,9.8821,0.05,0.25,2.6990,NbSimu(i),1000);
end ;
h= figure; plot(NbSimu,[repmat(BLprice,NbPoints,1) ResAV ResMC ResMCCV]); legend({'Black Scholes' 'Antithetic' , 'Simple Monte Carlo', 'Monte Carlo with Control Variable'})
set(h,'WindowStyle' , 'Docked' );
grid on ;
xlabel('Number of simuation');
Ylabel('Vanilla option Price');
xlim([NbSimu(1) NbSimu(end)]);
h= figure;
plot(NbSimu,[diff(CIAV,1,2) diff(CIMC,1,2) diff(CIMCCV,1,2) ]);
legend({'Antithetic' , 'Simple Monte Carlo', 'Monte Carlo with Control Variable' , });
grid on ;xlim([NbSimu(1) NbSimu(end)]);
title('Comparing different Monte Carlo methods')
set(h,'WindowStyle' , 'Docked' );
xlabel('Number of simuation');
Ylabel('absolute confidence interval');
h= figure;
FillBetween(NbSimu',CIMCCV(:,1),CIMCCV(:,2),'g' ,0.3);
hold on ;
plot(NbSimu,ResMCCV,'LineWidth' ,1.5, 'Color' , 'r' );
plot(NbSimu,repmat(BLprice,NbPoints,1),'b' );
xlim([NbSimu(1) NbSimu(end)])
set(h,'WindowStyle' , 'Docked' );
title('Monte Carlo with Control Variable');
xlabel('Number of simulation');
Ylabel('Vanilla option Price');
grid on ;
看跌期权蒙特卡洛模拟定价程序：
clear all ;
close all ;
clc;
NbSimu =20000:20000:200000;
NbPoints = numel(NbSimu);
ResAV = zeros(NbPoints,1);
ResMC = zeros(NbPoints,1);
ResMCCV = zeros(NbPoints,1);
CIAV = zeros(NbPoints,2);
CIMC = zeros(NbPoints,2);
CIMCCV = zeros(NbPoints,2);
VarPriceAV = zeros(NbPoints,1);
VarPriceMC = zeros(NbPoints,1);
VarPriceCV = zeros(NbPoints,1);
[call,put] = blsprice(7.92,9.8821,0.05,0.25,2.6990,0);
for i = 1 :numel(NbSimu)
disp(['Computing Iteration N? num2str(i)' ' num2str(NbSimu(i))' 'Simulated Paths' ]);
[ResAV(i),VarPriceAV(i),CIAV(i,:)] =
BlsMCAV(7.92,9.8821,0.05,0.25,2.6990,NbSimu(i));
[ResMC(i),VarPriceMC(i),CIMC(i,:)] =
BlsMC(7.92,9.8821,0.05,0.25,2.6990,NbSimu(i));
[ResMCCV(i),VarPriceCV(i),CIMCCV(i,:)] =
BlsMCCV(7.92,9.8821,0.05,0.25,2.6990,NbSimu(i),1000);
end ;
h= figure;
plot(NbSimu,[repmat(put,NbPoints,1) ResAV ResMC ResMCCV]);
legend({'Black Scholes' 'Antithetic' , 'Simple Monte Carlo', 'Monte Carlo with Control Variable'})
set(h,'WindowStyle' , 'Docked' );
grid on ;
xlabel('Number of simuation');
Ylabel('Vanilla option Price');
xlim([NbSimu(1) NbSimu(end)]);
h= figure;
plot(NbSimu,[diff(CIAV,1,2) diff(CIMC,1,2) diff(CIMCCV,1,2) ]);
legend({'Antithetic' , 'Simple Monte Carlo', 'Monte Carlo with Control
Variable' , });
grid on ;xlim([NbSimu(1) NbSimu(end)]);
title('Comparing different Monte Carlo methods') set(h,'WindowStyle' , 'Docked' );
xlabel('Number of simuation');
Ylabel('absolute confidence interval');
h= figure;
FillBetween(NbSimu',CIMCCV(:,1),CIMCCV(:,2),'g' ,0.3); hold on ;
plot(NbSimu,ResMCCV,'LineWidth' ,1.5, 'Color' , 'r' ); plot(NbSimu,repmat(put,NbPoints,1),'b' );
xlim([NbSimu(1) NbSimu(end)])
set(h,'WindowStyle' , 'Docked' );
title('Monte Carlo with Control Variable');
xlabel('Number of simulation');
Ylabel('Vanilla option Price');
grid on ;