几何布朗运动
如果没有几何布朗运动（Geometric Brownian Motion）就不会有布莱克-斯科尔斯-莫顿（BSM）欧式期权定价公式。
因此，金融数学（或者金融工程学或者量化金融学）的基础知识点之一就是要求学生熟练掌握几何布朗运动的所有性质。
在BSM欧式期权定价公式中，标的物–也就是股票价格是用几何布朗运动来描述的，这也就是为何我们说如果没有几何布朗运动就不会有布莱克-斯科尔斯-莫顿（BSM）欧式期权定价公式。
几何布朗运动较好的刻画了股票的以下特性：
几何布朗运动刻画的股票价格都是正数，不会出现负数的情形，这和现实中的股票价格的属性一致
几何布朗运动的走势有时会出现“疯涨”或者“疯跌”情形，这种特性和现实中的股票价格的特性相符合
将股票价格设定为几何布朗运动，并将股票价格作为衍生品的标的物，进而研究衍生品的定价，如此处理，较容易计算出衍生品的价格
当然，用几何布朗运动来刻画股票价格也是有明显的缺陷，比如：
几何布朗运动假设波动率（volatility）是一个常数，而在现实世界中，股票价格的波动率并非常数。【这也就是为何后来2003年诺贝尔经济学奖得主恩格尔教授一直倡导的ARCH模型的原因，ARCH家族模型一直将波动率作为一个动态变化过程来研究，具体可见：现代金融学的开端要从诺奖得主恩格尔说起】
股票价格市场有跳跃（jump）现象，而几何布朗运动无法刻画。
今天我们来演示一下几何布朗运动的数值表达式和计算机模拟实现，我们给出了Matlab、R和Python的代码实现。
一般而言我们有两种方法，第一种方法是利用的积分结果，即
（公式-1）。
。但我们不能直接用，我们都要做一些简单处理。
，因为我们知道，于是利用迭代方法可以得到的结果。
下面说说第二种方法，我们不妨假设于是 ，其中下面我们首先演示一下Malab代码的实现
clear
ti
c10; %Initial Stock price
10;
10e3;
N_sim=10e2;
0.0:dt:T];
0.05;
0.25;
1);
1)=S0*ones(N_sim,1);
subplot(1,2,1)
for j=1:N_sim
for i=1:N_time
1,1);
1)=S(j,i)*exp((mu-0.5*(vol^2))*dt+…
sqrt(dt)*z);
end
plot(t,S(j,?);
mu_s=S0*exp(mu*t);
exp(mu*t).*sqrt(exp((vol^2)*t)-1);
plot(t,mu_s+vol_s,’.k’);
xlabel(‘time’)
title(‘Simulated geometric Brownian Montion’)
100;
1,2,2)
hold on
histfit(S(:,end),nBins, ‘lognormal’);
set(gca, ‘xdir’, ‘reverse’)
270)
ylabel(‘frquencies’)
figure
1:end));
std(S(:,1:end));
plot(t,sample_mean_S,’.k’);
hold on
plot(t,mu_s);
plot(t,sample_mean_S+sample_std_S,’+r’);
plot(t,mu_s-vol_s,’.-k’);
grid on
xlabel(‘Time’)
ylabel(‘S_{t}’)
legend(‘Simluated sample mean of S_{t}^{i}’,…
‘Simulated sample +std’,’Simulated sample -std’,…
std’, ‘Theoretical -std’)
title(‘Simulated sample moments vs theoretical momnets’)
【图1几何布朗运动的模拟轨迹图与其对应的直方图分布】
图1模拟了1000次几何布朗运动从零时刻到T=10
，在左图中，我们还分别绘制了几何布朗运动的理论均值和正负一个理论标准差，几何布朗运动的一个理论标准差为。右图绘制了模拟1000次路径后，样本的直方图，因为我们知道
满足对数正态分布，使用了Matlab自带的“histfit”命令，让其将样本分布拟合对数正态分布，我们得到了一条非常完美的拟合曲线，从可视化角度我们认可验证，确实满足对数正态分布。
【图2几何布朗运动的样本和理论矩结果比较】
下面我们讨论第二种方法，我们用R语言同时实现了第一种和第二种方法，并比较了结果。可以发现，结果完全一致。下面我们给出一种方法的R语言的代码
mu<-0.3
sigma<-0.9
T0 <- 1
set.seed (123)
nn <-252
simulations.total <- 400
maturity <- T0
dt <- maturity/(nn+1)
timeline <- seq(0,maturity, dt)
f <- matrix(r0,(nn+2),simulations.total)
for(i in 2:(nn+2)){
f[i,j]<- f[i -1,j]*exp((mu-0.5*sigma^2)*dt + sigma*rnorm(1)*sqrt(dt))
}
plot(timeline,f[,1], ylim=range(f,thetaQ), type=”l”, col=”mediumorchid2″)
for(j in 2:simulations.total){
}
title(main=”Simulation of GMB”, col.main=”red”, font.main=4)
下面我们给出第二种方法的代码
mu<-0.3
sigma<-0.9
T0 <- 1
set.seed (123)
nn <-252
simulations.total <- 400
maturity <- T0
dt <- maturity/(nn+1)
timeline <- seq(0,maturity, dt)
f <- matrix(r0,(nn+2),simulations.total)
for(i in 2:(nn+2)){
f[i,j]<- f[i -1,j]*(1+mu*dt + sigma*rnorm(1)*sqrt(dt))
}
plot(timeline,f[,1], ylim=range(f,thetaQ), type=”l”, col=”mediumorchid2″)
for(j in 2:simulations.total){
}
title(main=”Simulation of GMB”, col.main=”red”, font.main=4)
【图3 两种模拟方法的显示结果】
最后我们还给出Python代码的实现
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
mu=1
n=50
dt=0.1
x0=100
x=pd.DataFrame()
np.random.seed(1)
for sigma in np.arange(0.8,2,0.2):
step=np.exp((mu-sigma**2/2)*dt)*np.exp(sigma*np.random.normal(0,np.sqrt(dt),(1,n)))
temp=pd.DataFrame(x0*step.cumprod())
x=pd.concat([x,temp],axis=1)
x.columns=np.arange(0.8,2,0.2)
plt.plot(x)
plt.legend(x.columns)
plt.xlabel(‘t’)
plt.ylabel(‘X’)
plt.title(‘Realizations of Geometric Brownian Motion with different variancesn mu=1’)
plt.show()
