Python实现蒙特卡洛模拟的期权估值

Posted on | Edited on | In | Comments: | Views:

期权定价原理

常用的期权定价模型是BSM模型,适用于欧式期权：

本项目用蒙特卡洛模拟方法既可以对欧式期权定价，也可以对美式期权定价。蒙特卡洛方法的理论基础是概率论与数理统计，其实质是通过模拟标的资产价格路径预测期权的平均回报并得到期权价格估计值。

实现过程

生成随机数

随机数的生成是随机分析的开始，蒙特卡洛随机模拟的精度很大程度上依赖于所选的随机数。这里通过numpy.random提供的函数实现，但由于函数产生的是伪随机数，而且我们抽样的数量是变化的，因此得到的随机数集合通常某些统计量并不是精确的满足我们的要求。例如，要生成均值为0方差为1的标准正态分布随机数，但实际生成的数据是有一定误差的。这里使用2种方差降低技术来更为精确的匹配标准正态分布的前两阶矩。
1、对偶方法（antithetic variates），先抽取一半数量的随机数，然后加上相同数字的相反数。这种方式很容易实现一阶矩的匹配，因为这样抽样是能够保证均值为0的，但对于二阶矩标准差没有影响。
2、矩匹配（moment matching），将每个随机数减去所有随机数的均值，再除以标准差。可以准确的拟合一阶和二阶矩。

1
2
3
4
5
6
7
8
9
10
11
12
 
import numpy as np
def sn_random_numbers(shape,antithetic = True,moment_matching = True,fixed_seed=False):
    if fixed_seed:
        np.random.seed(1000)
    if antithetic:
        ran = np.random.standard_normal((shape[0],shape[1], round(shape[2]/2)))
        ran = np.concatenate((ran,-ran),axis=2)
    else:
        ran = np.random.standard_normal(shape)
    if moment_matching:
        ran = ran-np.mean(ran)
        ran = ran/np.std(ran)

生成模拟路径

将时间间隔[0, T]分为等距的、长度为d𝑡的子时段，M是时段个数，I是生成的路径数，根据公式计算每条路径的在t时的价格，存放在paths数组中，模拟的结果可以通过时段数M和路径数I来进行控制。 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
 
def generate_paths(self,fixed_seed = False,day_count=365.):
       if self.time_grid is None:
           self.generate_time_grid()
       M = len(self.time_grid)
       I = self.paths
       paths = np.zeros((M,I))
       paths[0]=self.initial_value
       rand = sn_random_numbers((1,M,I),fixed_seed=fixed_seed)
       short_rate = self.discount_curve.short_rate
       for t in range(1,len(self.time_grid)):
           ran = rand[t]
           dt = (self.time_grid[t]-self.time_grid[t-1]).days/day_count
           paths[t] = paths[t-1]*np.exp((short_rate-0.5*self.volatility**2)*dt+self.volatility*np.sqrt(dt)*ran)
       self.instrument_values = paths

当给定输入条件:
无风险利率r=5%,定价日为2017-1-1，到期日2017-12-31，波动率volatility=20%，标的资产初始价格S=36，模拟频率按月，路径数取10000条时，看到模拟结果：

对比波动率分别取20%和50%时，资产价格的走势： 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
 

men = market_environment('men',dt.datetime(2017,1,1))
men.add_constant('initial_value',36.)
men.add_constant('volatility',0.2)
men.add_constant('final_date',dt.datetime(2017,12,31))
men.add_constant('currency','EUR')
men.add_constant('frequency','M')
men.add_constant('paths',10000)
csr = constant_short_rate('csr',0.05)
men.add_curve('discount_curve',csr)

gbm = geometric_brownian_motion('gbm',men)
gbm.generate_time_grid()
path_1=gbm.get_instrument_values()
gbm.update(volatility=0.5)
path_2 = gbm.get_instrument_values()
plt.figure(figsize=(8,4))
p1 = plt.plot(gbm.time_grid,path_1[:,:10],'b')
p2 = plt.plot(gbm.time_grid,path_2[:,:10],'r-.')
plt.grid(True)
l1 = plt.legend([p1[0],p2[0]],['low_volatility','high_volatility'],loc=2)
plt.gca().add_artist(l1)
plt.savefig('simulation_volatility.jpg')
plt.show()

欧式期权定价

有了标的资产的价格模拟路径，期权定价就容易了，以欧式看涨期权为例，计算到期日的每条路径的payoff（即max（S-X,0)），再对所有payoff取平均，折现到起始日期即得到欧式期权价格。

1
2
3
4
5
6
 
payoff_func='np.maximum(maturity_value-strike,0)'
    def present_value(self,accuracy=6,fixed_seed=False):
        cash_flow = self.generate_payoff(fixed_seed=fixed_seed)
        discount_factor=self.discount_curve.get_discount_factors((self.pricing_date,self.maturity))[0,1]
        result=discount_factor * np.sum(cash_flow)/len(cash_flow)
        return round(result,accuracy)

取执行价格X=40时，计算期权价格：

1
2
3
4
5
6
7
8
9
 
gbm = geometric_brownian_motion('gbm',men)
me_call=market_environment('me_call',men.pricing_date)
me_call.add_constant('strike',40.)
me_call.add_constant('maturity',dt.datetime(2017,12,31))
me_call.add_constant('currency','EUR')
payoff_func='np.maximum(maturity_value-strike,0)'
eur_call=valuation_mcs_european('eur_cal',underlying=gbm,mar_env=me_call,payoff_func=payoff_func)
call_value=eur_call.present_value()
print(call_value)

得到欧式看涨期权价格为：2.184223

美式期权定价

欧式期权只有到期日1个行权日，美式期权可以在到期前任意时刻行权，定价相对更复杂。需要在每个时间点比较行权的价值和不行权的价值，选择二者的高值，而行权价值是期权的内在价值（以美式看跌期权为例）：max（X-S，0），不行权的价值是剩余时间期权产生的价值的折现，也称为持续价值，这里用最小二乘法对期权价值和标的资产价格做二项式回归，计算期权的持续价值。 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
 
def present_value(self,accuracy=6,fixed_seed=False,bf=5,full=False):
       instrument_values,inner_values,time_index_start,time_index_end =self.generate_payoff(fixed_seed=fixed_seed)
       time_list=self.underlying.time_grid[time_index_start:time_index_end+1]

       discount_factors=self.discount_curve.get_discount_factors(time_list,dtobjects=True)
       V=inner_values[-1]
       for t in range(len(time_list)-2,0,-1):
           df=discount_factors[t,1]/discount_factors[t+1,1]   #向前折一期的折现因子
           rg=np.polyfit(instrument_values[t],V*df,bf)
           C=np.polyval(rg,instrument_values[t])
           #比较持续价值和内在价值，取高值
           V=np.where(inner_values[t]>C,inner_values[t],V*df)
       df=discount_factors[0,1]/discount_factors[1,1]
       result=np.sum(V*df)/len(V)
       return round(result,accuracy)

假设标的资产和执行价格都与欧式期权相同，按周模拟，路径数取50000条，计算期权价格： 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
 
men = market_environment('men',dt.datetime(2017,1,1))
men.add_constant('initial_value',36.)
men.add_constant('volatility',0.2)
men.add_constant('final_date',dt.datetime(2017,12,31))
men.add_constant('currency','EUR')
men.add_constant('frequency','W')
men.add_constant('paths',50000)

csr = constant_short_rate('csr',0.05)
men.add_curve('discount_curve',csr)

gbm = geometric_brownian_motion('gbm',men)
payoff_func='np.maximum(strike-instrument_values,0)'

me_am_put=market_environment('me_am_put',men.pricing_date)
me_am_put.add_constant('strike',40.)
me_am_put.add_constant('maturity',dt.datetime(2017,12,31))
me_am_put.add_constant('currency','EUR')

am_put=valuation_mcs_american('am_put',underlying=gbm,mar_env=me_am_put,payoff_func=payoff_func)
temp=am_put.present_value(fixed_seed=True,bf=5)
print(temp)

得到美式看跌期权价格为:4.577837。
相同条件下，美式期权价值一定高于欧式期权。