大三的一个大作业，时间有点久远，这里把代码块和解释都尽量补齐，毕竟是耗费两个月的作品还是有感情。我用的是 jupyter notebook， 先下一个 Anaconda 就可以，以后要用的时候就终端输入 jupyter notebook 就能自动跳转页面
首先 task 是
dS/S = r*dt + sigma*dW
dB/B = r*dt
a.1) Simulate S with MC reduction methods
a.2) Simulate European option price with MC and BS model
b. Conduct performance analysis by changing M and N
c. Estimate the probability functions of strike price
d) implement static delta-neutral strategy in MC
e) implement static delta-gamma-neutral strategy in MC
随机过程的公式就是上面的 dS/S，股价随机波动，当然具体的sigma，r 等参数后面是自己设定的
这里用人话解释一下是什么意思（中外合办大学，英文教学两行泪）
a1 首先用蒙特卡洛的 reduction 方法模拟出股票也就是 Stock 的价格变化
a2 然后 MC 和 BS 的方法算出欧洲期权的实时价格，也就是一个分离和连续时间的区别，这里用一下公式就行了
b 通过交换 M 和 N（M 是有几种股价 jump 的方向，N 是进行多少个周期也就是 steps），来监测股价波动的表现
c 画出 strike price（K）的概率图
d 用 monte-carlo 方法实施静态 delta 中性对冲策略
e 用 monte-carlo 方法实施静态 delta-gamma 中性对冲策略
首先调包 numpy, pandas, matplotlib 基本操作，防止后面出现一些奇怪的报错
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm
然后是一些函数定义，之后要用到就直接拿来用了
def EulerMilsteinMCStock(scheme, parameters):
np.random.seed(1000)
# time setup
T = parameters['setup']['T'] # total time/maturity
numSteps = parameters['setup']['numSteps'] # number of steps
numPaths = parameters['setup']['numPaths'] # number of simulated paths
dt = parameters['setup']['dt']
# model parameters
S_0 = parameters['model']['S0'] # initial value
sigma = parameters['model']['sigma'] # initial value
rf = parameters['model']['rf'] # initial value
# simulation
S = np.zeros((numSteps + 1, numPaths),dtype=float)
S[0,:] = np.log(S_0)
################ simluations for asset price S ########
for i in range(numPaths):
for t_step in range(1, numSteps+1):
# the 2 stochastic drivers for variance V and asset price S and correlated
Zs = np.random.normal(0, 1, 1)
if scheme == 'Euler':
S[t_step,i] = S[t_step-1,i] + (rf-sigma**2/2)*dt + sigma*np.sqrt(dt)*Zs
elif scheme == 'Milstein':
# Euler and Milstein have the same discretization scheme for Log(S) due to dsigma(t,X)/dX =0
S[t_step,i] = S[t_step-1,i] + (rf-sigma**2/2)*dt + sigma*np.sqrt(dt)*Zs
return np.exp(S)
这个 Euler 和 Milstein 是两种 MC simulation 的方法，各有利弊和简便之分，这个我们一般用 Euler 就可以了。np.random.seed() 是随机生成种子数，了解蒙特卡洛原理的就知道这个本身就是要很多次实验下产生的结果更加准确，你想要一亿两亿都可以。
然后是time的设定，总时间，总步数（N），还有一次jump的路径数量（M），dt 就是 T/N
model 的参数设定，Stock 初始价格，波动参数，无风险利率，都是后面自己设定的，不要太夸张给个2%就行。
然后就是算法原理的代码实践。输出股价和时间的对应。
第二个函数定义，BS model，这个公式是量化金融常出现的。用来测算期权价格应用很广泛，定义之后就可以直接用了。
def Black_ScholesPrice(parameters):
S0 = parameters['model']['S0']
sigma = parameters['model']['sigma']
K = parameters['asset']['K']
rf = parameters['model']['rf']
T = parameters['setup']['T']
option_type = parameters['asset']['optype']
DF = np.exp(-rf*T)
d1 = (np.log(S0/K)+(rf+0.5*sigma**2)*T)/(sigma*np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
if option_type == 1: # call
option_price = S0 * norm.cdf(d1) - K * DF * norm.cdf(d2)
elif option_type == -1: # put
option_price = K * DF * norm.cdf(-d2) - S0 * norm.cdf(-d1)
return option_price
至此，我们就可以做 a1 a2 小题了，因为公式都定义好了，之后用 main 环境就可以直接运行。定义一下参数们，用一下之前定义好的函数。这个等会再一起和对冲码出来。
接下来就是 b 的要求，变化 M 和 N 的数值，看一下股价的表现，在不同方法的模拟下，这一块由于我们的电脑实在带不动大的 M 和 N，你想 M=10000 就意思是一个价格在某一时刻有 10000 种上升或下降的不同价格，配上一个大 N，我的破 Mac 根本带不动，于是我们用小的数据慢慢的一个个试，从 M=200，N=200 开始，一直到 M=5000, N=200, 我们发现其实 N 太大了会影响作图的清晰度所以我们就保持在200 就刚刚好，M 上了千之后其实目测的效果都是差不多的，大概在 50-200，的价格波动（初始价格 S0 =100），头铁的话可以试试建个数组让 M 和 N 分别变化看 performance，我们的小电脑一下就卡死了。形而上学，看看就好：
M = np.arange(2000,50000,2000)
Option_Price = np.zeros(len(M))
BS_line = np.empty(len(M))
for i in range(len(M)):
BS_line[i] = BS_price
#print(BS_line)
for i in range(0,len(M)):
M_new=M[i]
parameters = {'model':{'S0':S0, 'sigma':sigma, 'rf':rf},
'asset':{'K':K, 'optype': Optype},
'setup':{'T':T, 'numSteps':N, 'dt': T/N, 'numPaths':M_new}
}
Sim_S = EulerMilsteinMCStock('Euler', parameters)
S_T = Sim_S[-1,:]
Option_Price[i] = MC_option_price(S_T, parameters)
#print(Option_Price)
plt.figure(figsize=(12,8))
plt.plot(M,Option_Price)
plt.plot(M,BS_line)
plt.ylabel('Monte Carlo Option Price')
plt.xlabel('Number of Paths')
plt.show()
N = np.arange(500,44500,2000)
Option_Price_1 = np.zeros(len(N))
#print(N)
BS_line = np.empty(len(N))
for i in range(len(N)):
BS_line[i] = BS_price
#print(BS_line)
for i in range(0,len(N)):
N_new=N[i]
parameters = {'model':{'S0':S0, 'sigma':sigma, 'rf':rf},
'asset':{'K':K, 'optype': Optype},
'setup':{'T':T, 'numSteps':N_new, 'dt': T/N_new, 'numPaths':M}
}
Sim_S = EulerMilsteinMCStock('Euler', parameters)
S_T = Sim_S[-1,:]
Option_Price_1[i] = MC_option_price(S_T, parameters)
#print(Option_Price_1)
plt.figure(figsize=(12,8))
plt.plot(N,Option_Price_1)
plt.plot(N,BS_line)
plt.ylabel('Monte Carlo Option Price')
plt.xlabel('Number of time steps')
plt.show()
这块其实不是那么重要，因为大家都跑不出来，你有理说一个你觉得好的 M 和 N 就可以了，毕竟咱们也没那么好的电脑是不是。
其实到这 b 就解决了，出来一个像 K 线一样的图。
股价波动图 N=300
这里我们用了一个很小的 N 就是因为发现大了之后图真看不清了，还不如小点清晰度高，反正之前都验证过了大点小点没什么特别大的关系，有关系我们也看不出来，也要对应大交易量才能看出来。
那 c 的话就是求一个 St=K 的概率, call option 的话当你 St>K 那肯定就要成交了，所以我们在这个地方需要用到 call option 的公式，经过一系列的数学推导（提示，option price 的公式中有用到期望值，期望值可以用积分代替，这样就可以转化成某个区域的 probability），代码如下：
#estimate cumulative/density function
K_set = np.arange(150,50,-0.5, float)
num_K = len(K_set)
MC_price_Ks = np.zeros(num_K)
for i in range(num_K):
parameters['asset']['K'] = K_set[i]
MC_price_Ks[i] = MC_option_price(S_T, parameters)
# convert <K_set, MC_price_Ks> into dataframe
MC_K_vs_price = pd.DataFrame(np.transpose([K_set,MC_price_Ks]),
index =list(range(0, num_K)),
columns=['K', 'Price'])
#print(MC_K_vs_price)
# estimate culumative probs
HigherThanK_prob = ProbHigherThanK(K_set, MC_price_Ks, parameters)
# estimate density probs
den_probs = Density_prob(K_set, MC_price_Ks, parameters)
# save all data
All_info = pd.DataFrame(np.transpose([K_set,MC_price_Ks, HigherThanK_prob, den_probs]),
index =list(range(0, num_K)),
columns=['K', 'Price', 'Prob(ST>K)', 'Prob(ST=K)'])
# print(All_info)
All_info.to_excel(r"all_MC_info.xlsx")
# plot
plt.figure(figsize=(60,45))
plt.subplot(1,3,1)
plt.plot(K_set, den_probs, 'b')
plt.grid(True)
plt.xlabel('K')
plt.ylabel('Prob(ST=K)')
plt.subplot(1,3,2)
plt.plot(K_set, HigherThanK_prob, 'b')
plt.grid(True)
plt.xlabel('K')
plt.ylabel('Prob(ST>K)')
plt.subplot(1,3,3)
plt.plot(K_set, 1-HigherThanK_prob, 'b')
plt.grid(True)
plt.xlabel('K')
plt.ylabel('Prob(ST<K)')
用的是subplot（1，3）所以我们就出来一行三个图
St=K 概率图
从这个图中我们其实可以看出啊，K=80 左右的时候，St 最有可能相等，S0=100，当你 K 定一个 80 左右的价钱时候是很有可能在股价这么波折波折中赚钱的，厉害吧。
接下来就是要用量化对冲的两个策略，所谓什么叫delta对冲和gamma对冲呢，就是（期权价格/股价）的一阶导和二阶导，是为了用来对冲的股市变动带来资产变化风险的（具体原理可以百度，数学原理还蛮复杂的，不过大白话意思就是导的越多肯定对冲风险越精确），那d和e就是要看一下两种对冲策略下股价的表现以及资产变化怎么样。
首先要定义函数 Black-Scholes 的在量化对冲上的使用，这些都是有数学公式做支撑的，算法就是翻译一下，在这里就不糊眼睛了。
def Black_Scholes_Delta(parameters):
# based on the lecture, delta = N(d1)
S0 = parameters['model']['S0']
sigma = parameters['model']['sigma']
rf = parameters['model']['rf']
K = parameters['option']['K']
T = parameters['option']['T']
option_type = parameters['option']['optype']
d1 = (np.log(S0/K)+(rf+0.5*sigma**2)*T)/(sigma*np.sqrt(T))
if option_type == 1: # call
bs_delta = norm.cdf(d1)
else:
bs_delta = norm.cdf(d1) - 1
return bs_delta
def Black_Scholes_Gamma(parameters):
# gamma = phi(d1)/S_t/(sigma*sqrt(T))
S0 = parameters['model']['S0']
sigma = parameters['model']['sigma']
rf = parameters['model']['rf']
K = parameters['option']['K']
T = parameters['option']['T']
d1 = (np.log(S0/K)+(rf+0.5*sigma**2)*T)/(sigma*np.sqrt(T))
bs_gamma = norm.pdf(d1)/S0/(sigma*np.sqrt(T))
return bs_gamma
那再把主环境设置一下，参数自行设置好。方便电脑快速运转，参数都调了小的数。把 a 的两小问也放进去了，反正就是函数调用一下的事情。
if __name__=="__main__": # main function for d,e
#==================================================
# Environment Setup
#==================================================
# underlying info
S0 = 100
sigma = 0.40
rf = 0.025
# option info
K = 110
T = 3 # maturity
Optype = 1 # 1: call -1: put
# discrete setup
N = 20 # steps
dt = T/N
M = 10
# info structure
parameters = {'model':{'S0':S0, 'sigma':sigma, 'rf':rf},
'asset':{'K':K, 'optype':Optype},
'option':{'K':0, 'optype':0, 'T':0}, # save all option info
'setup':{'T':T, 'numSteps':N, 'dt': T/N, 'numPaths':M},
'strategy':{'n_c':0, 'n_s':0, 'n_b':0, 'k':0}
}
#==================================================
# simulation of stock price paths
#==================================================
Sim_S = EulerMilsteinMCStock('Euler', parameters)
Stock_DF = pd.DataFrame(Sim_S)
Stock_DF.to_excel(r"Stock_paths.xlsx")
# plot simulated stock prices
plt.figure(figsize=(12, 4))
plt.plot(Sim_S)
plt.ylabel('Simulated stock price paths')
plt.xlabel('time step')
plt.tight_layout()
plt.show()
#==================================================
# Task 1 -- delta-neutral strategy (static)
# Key ideas: 1) the strategy is established at t=0
# 2) hold it until maturity T wihout rebalancing
# 3) replication errors are recored during [0,T]
#==================================================
# setup
#option_1 at t =0
parameters_1 = parameters
parameters_1['option']['K'] = K
parameters_1['option']['optype'] = Optype # call
parameters_1['option']['T'] = T
BS_price_1 = Black_ScholesPrice(parameters_1)
delta_1 = Black_Scholes_Delta(parameters_1)
gamma_1 = Black_Scholes_Gamma(parameters_1)
parameters_1['strategy']['n_c'] = -1 # -1: a short position in Option 1
parameters_1['strategy']['n_s'] = delta_1
parameters_1['strategy']['n_b'] = -(parameters_1['strategy']['n_c']*BS_price_1
+ parameters_1['strategy']['n_s']* parameters_1['model']['S0']
)
parameters_1['strategy']['k'] = 0
print('parameter_1=\n')
print(parameters_1)
# perform delta-neutral hedging
PI_DN = DeltaGammaNeutralHedging(Sim_S, parameters_1)
PI_DN_DF = pd.DataFrame(PI_DN)
PI_DN_DF.to_excel(r"PI_DN_paths.xlsx")
#==================================================
# Task 2 -- delta-gamma-neutral strategy (static)
# Key ideas: 1) the strategy is established at t=0
# 2) hold it until maturity T wihout rebalancing
# 3) replication errors are recored during [0,T]
#==================================================
# setup
#option_2 at t =0
parameters_2 = parameters
parameters_2['option']['K'] = K
parameters_2['option']['optype'] = Optype # call
parameters_2['option']['T'] = T + 1 # with longer maturity
BS_price_2 = Black_ScholesPrice(parameters_2)
delta_2 = Black_Scholes_Delta(parameters_2)
gamma_2 = Black_Scholes_Gamma(parameters_2)
parameters_2['strategy']['n_c'] = -1 # -1: a short position in Option 1
parameters_2['strategy']['k'] = -parameters_2['strategy']['n_c']*gamma_1/gamma_2
parameters_2['strategy']['n_s'] = -parameters_2['strategy']['n_c']*delta_1 - parameters_2['strategy']['k']* delta_2
parameters_2['strategy']['n_b'] = -(parameters_2['strategy']['n_c']*BS_price_1
+ parameters_2['strategy']['n_s']* parameters_2['model']['S0']
+ parameters_2['strategy']['k']* BS_price_2
)
print('parameter_2=\n')
print(parameters_2)
print('option_1=%f vs option_2 = %f'%(BS_price_1, BS_price_2))
print('delta_1=%f vs delta_2=%f'%(delta_1, delta_2))
print('gamma_1=%f vs gamma_2=%f'%(gamma_1, gamma_2))
#perform delta-gamma-neutral hedging
PI_DGN = DeltaGammaNeutralHedging(Sim_S,parameters_2)
PI_DGN_DF = pd.DataFrame(PI_DGN)
PI_DGN_DF.to_excel(r"PI_DGN_paths.xlsx")
#==================================================
# Plot results
#==================================================
plt.figure(figsize=(12,8))
plt.subplot(2,1,1)
plt.plot(PI_DN)
plt.grid(True)
#plt.xlabel('Paths(M)')
plt.ylabel('Replication Errors (in cash units)')
plt.title('Delta-Neutral Hedging')
plt.subplot(2,1,2)
plt.plot(PI_DGN)
plt.grid(True)
plt.xlabel('Paths(M)')
plt.ylabel('Replication Errors (in cash units)')
plt.title('Delta-Gamma-Neutral Hedging')
plt.figure(figsize=(12,4))
data_DN = np.transpose([Sim_S[-1,:],PI_DN[-1,:]]) # convert data into the one with two columns
#print(data_DN[np.lexsort(data_DN[:,::-1].T)]) # sort data based on the first column and change the data in 2nd column accordingly
new_data_DN = data_DN[np.lexsort(data_DN[:,::-1].T)] # save the sorted data to a new variable
plt.plot(new_data_DN[:,0],new_data_DN[:,1],"b",label='Delta-Neutral') # plot data in the way we want
data_DGN = np.transpose([Sim_S[-1,:],PI_DGN[-1,:]]) # convert data into the one with two columns
#print(data_DGN[np.lexsort(data_DGN[:,::-1].T)]) # sort data based on the first column and change the data in 2nd column accordingly
new_data_DGN = data_DGN[np.lexsort(data_DGN[:,::-1].T)] # save the sorted data to a new variable
plt.plot(new_data_DGN[:,0],new_data_DGN[:,1],"r-*",label='Delta-Gamma-Neutral') # plot data in the way we want
plt.legend(loc=0)
plt.grid(True)
plt.xlabel('Stock Price (S_T)')
plt.ylabel('Replication Errors (in cash units)')
plt.title('(Static) Hedging Performance at Maturity(T)')
plt.figure(figsize=(12,4))
plt.plot(np.mean(PI_DN,axis=1),"b", label='Delta-Neutral')
plt.plot(np.mean(PI_DGN,axis=1),'r-*',label='Delta-Gamma-Neutral')
plt.legend(loc=0)
plt.grid(True)
plt.xlabel('Paths(M)')
plt.ylabel('Averaged Replication Errors (in cash units)')
plt.title('(Static) Hedging Performance over tim [0,T]')
股价波动图
那其实这个能看出来这个股价是会有不停的波动的，如果我们不去管他，在这个市场上，我们的资产价值就一直在波动，那肯定是控制不好风险的，那看看经过对冲之后表现如何。delta 只用一个 option，gamma 要用两个 option。
随步数变化的股价波动值
从这个图其实我们已经能看出来了，单单是delta hedging无法让风险被很好的对冲掉，一开始是有些线重合在一起了好像挺稳定的，但也一直在往上走，到过了一段时间的时候差不多又毫无章法了，反观delta-gamma的对冲效果就明显好很多，error一直保持在 6/每现金单位，就算之后失去控制也没有那么夸张。
对冲表现
这张图就更明显了，落实到平均上看一下，不管股价怎么动，delta-gamma的对冲效果始终能保证误差在一个小范围内可控，但是delta波动较大，下面的地方，红色的双option对冲策略一条直线，谁更平稳显而易见。
以上是之前大作业的回顾，有些地方如果有错误欢迎指正～毕竟我也不是很专业，这门课是金融工程，第一次用 python 代码老师给了很多，第一次学我们以应用和理解为主。讲真为什么 delta-gamma hedging 这么好还要用那个 delta hedging 呢，那肯定是是因为好的更贵更麻烦呀，这么复杂的过程，落实到真正的股市上，那计算量是翻倍的涨啊，得考虑到经济效益的问题。
