

用python实现股票价格几何布朗运动的Monte Carlo method

数据皮皮侠

2020-04-24

导读：前言在经典的期权定价理论中，有一个假设是股票价格在短期内服从几何布朗运动，以及为基础发展出了Black-Sc

前言

在经典的期权定价理论中，有一个假设是股票价格在短期内服从几何布朗运动，以及为基础发展出了Black-Scholes等期权定价公式，但在现实中股票价格是否真能如假设一样服从几何布朗运用？我国股市，特别是股市中的个股是否是符合这一规律？为此我使用python构建 Monte Carlo method，对股票价格进行模拟，与真实对比。

理论基础

影响股票价格变化的因素主要有股票价格随时间上涨的趋势和股票价格的平均波动率。前者对股票价格增长的贡献取决于时间的长短；后者只取决于布朗运动造成的随机波动。

如果用s表示股票价格，表示股票预期收益率，表示波动率，且均为常数，t代表时间，z为标准布朗运动，则有：

好了下面让我们进入代码部分

具体代码

第一步导入一些需要的库

# -*- coding: utf-8 -*-import numpy as npimport pandas as pdimport matplotlib.pyplot as pltfrom time import time

主要使用到的是numpy，pandas,matplotlib.pyplot和计算时间的time

第二步：我们要将收益率对数化并且计算μ和σ：

data=pd.read_csv(r'C:\Users\10547\Desktop\Midterm\data\002436.SZ.CSV',encoding='gbk')S_true=pd.DataFrame(data,columns=['日期','前收盘价(元)'])In_S_true=list()for i in range(len(S_true.iloc[:,1])-1):    In_S_true.append(np.log(S_true.iloc[i+1,1]/S_true.iloc[i,1]))mu=np.mean(In_S_true)sigma=np.sqrt(np.var(In_S_true))

这里面我使用的是来自Wind数据库的股票数据，其默认格式为”股票代码.交易所代码.CSV”，在导入时有一点需要特别说明，下载下来的数据里有中文列名，在导入时需要加上encoding='gbk'，否则在后续代码中可能会报错。在代码中S_true为真实股价，In_S_true为股价的对数收益率，我们可以根据上述公式求出需要的参数。

上面的代码里还涉及到一些Dataframe里数据索引的操作，与R里面的不太一样，诸位读者可以自己多尝试尝试。

第三步：进行Monte Carlo 模拟路径：

#几何布朗运动的蒙特卡洛模拟M=len(In_S_true) #时间步数I=2000 #模拟次数np.random.seed(19990801)#Tip：想要出来的模拟图不一样改变这个seed括号里面的数值就行start=time()S_model=np.zeros((M+1,I))S_model[0]=S_true.iloc[0,1]for j in range(0,I-1):    for i in range(1,len(In_S_true)):

S_model[i,j]=S_model[i-1,j]*np.exp((mu-0.5*sigma**2)+sigma*np.random.standard_normal(1))#S_model中每一列都是一种模拟情况下的股票价格变动#S_model=S_true.iloc[0,1]*np.exp(np.cumsum((mu-0.5*sigma**2)+sigma*np.random.standard_normal((M,I)),axis=0))end=time()print('total time is %.6f second'%(end-start))

这里一开始使用了for循环进行模拟，不过for循环是在太弟弟了，大家可以尝试一下我在end=time()上一行的语句，把#去掉，同时把S_model=np.zeros((M+1,l))开始到刚刚那行前的全部注释，你就可以使用全向量化的方法，效率会高更多。start=time()和end=time()就是我为了直观比较程序运行速度而写的代码，每个人电脑性能可能不同，在我这边for循环运行了68秒，而全向量化的只用了0.35秒。这一行的具体效果是这样的：生成一个随机变量的数组，M行，I列，同时计算出每一条路径，每一个时间点的指数水平的增量， np.cumsum(axis=0)，在列的方向上进行累加得到每一个时间步数上的指数水平。

第四步是绘图

可以先看看2000条可能路径长什么样

plt.figure(figsize=(12,7))plt.grid(True)plt.xlabel('Time')plt.ylabel('Stock Price')for i in range(I):plt.plot(S_model[:,i])

然后将2000条路径平均，与真实股价进行比较，注意两次绘图在Ipython console 里进行，放在同一个脚本里可能会图形绘到同一张里面去。

S_modelavg=list()
for i in range(len(In_S_true)):
    S_modelavg.append(np.average(S_model[i,:]))
plt.figure(figsize=(12,7))
plt.grid(True)
plt.xlabel('Time')
plt.ylabel('Stock Price')
plt.plot(S_modelavg)
plt.plot(S_true.iloc[:,1])