1.#tushare导入
2.import pandas as pd
3.import datetime
4.import tushare as ts
5.import numpy as np
6.from math import log,sqrt,exp
7.from scipy import stats
8.import plotly
9.import plotly.graph_objects as go
10.import plotly.express as px
11.pro = ts.pro_api(token = '')#输入自己的tushare token
12.plotly.offline.init_notebook_mode(connected=True)
13.#数据提取
14.def extra_data(date): # 提取数据
15. # 提取沪深300ETF合约基础信息
16. df_basic = pro.opt_basic(exchange='SSE', fields='ts_code,name,call_put,exercise_price,list_date,delist_date')
17. df_basic = df_basic.loc[df_basic['name'].str.contains('300ETF')]
18. df_basic = df_basic[(df_basic.list_date<=date)&(df_basic.delist_date>date)] #提取当天市场上交易期权合约
19. df_basic = df_basic.drop(['name','list_date'],axis=1)
20. df_basic['date'] = date
21. # 提取日线行情数据
22. df_cal = pro.trade_cal(exchange='SSE', cal_date=date, fields = 'cal_date,is_open,pretrade_date')
23. if df_cal.iloc[0, 1] == 0:
24. date = df_cal.iloc[0, 2] # 判断当天是否为交易日，若否则选择前一个交易日
25. opt_list = df_basic['ts_code'].tolist() # 获取300ETF期权合约列表
26. df_daily = pro.opt_daily(trade_date=date,exchange = 'SSE',fields='ts_code,trade_date,close')
27. df_daily = df_daily[df_daily['ts_code'].isin(opt_list)]
28. df_300etf = pro.fund_daily(ts_code='510300.SH', trade_date = date,fields = 'close')
29. s = df_300etf.iloc[0, 0]
30. # 提取无风险利率数据（用一周shibor利率表示）
31. df_shibor = pro.shibor(date = date,fields = '2w')
32. rf = df_shibor.iloc[0,0]/100
33.
34. # 数据合并
35. df = pd.merge(df_basic,df_daily,how='left',on=['ts_code'])
36. df['s'] = s
37. df['r'] = rf
38. df = df.rename(columns={'exercise_price':'k', 'close':'c'})
39.
40.#隐含波动率求解模块：二分法
41.def bsm_value(s,k,t,r,sigma,q,option_type):
42. d1 = ( log( s/k ) + ( r- q+ 0.5*sigma**2 )*t )/( sigma*sqrt(t) )
43. d2 = ( log( s/k ) + ( r- q - 0.5*sigma**2 )*t )/( sigma*sqrt(t) )
44. if option_type.lower() == 'c':
45. value = s*exp( -q*t )*stats.norm.cdf( d1) - k*exp( -r*t )*stats.norm.cdf( d2)
46. else:
47. value = k * exp(-r * t) * stats.norm.cdf(-d2) - s*exp( -q*t )*stats.norm.cdf(-d1)
48. return value
49.
50.def black_value(F,k,t,r,sigma,option_type):
51. d1 = ( log( F/k ) + ( 0.5*sigma**2 )*t )/( sigma*sqrt(t) )
52. d2 = ( log( F/k ) + ( -0.5*sigma**2 )*t )/( sigma*sqrt(t) )
53. if option_type.lower() == 'c':
54. value = exp( -r*t )*(F*stats.norm.cdf( d1) - k*stats.norm.cdf( d2))
55. else:
56. value = exp(-r * t)*(k*stats.norm.cdf(-d2) - F* stats.norm.cdf(-d1))
57. return value
58.#二分法求隐含波动率
59.def bsm_imp_vol(s,k,t,r,c,q,option_type):
60. c_est = 0 # 期权价格估计值
61. top = 1 #波动率上限
62. floor = 0 #波动率下限
63. sigma = ( floor + top )/2 #波动率初始值
64. count = 0 # 计数器
65. while abs( c - c_est ) > 0.00000001:
66. c_est = bsm_value(s,k,t,r,sigma,q,option_type)
67. #根据价格判断波动率是被低估还是高估，并对波动率做修正
68. count += 1
69. if count > 1000: # 时间价值为0的期权是算不出隐含波动率的，因此迭代到一定次数就不再迭代了
70. sigma = 0
71. break
72.
73. if c - c_est > 0: #f(x)>0
74. floor = sigma
75. sigma = ( sigma + top )/2
76. else:
77. top = sigma
78. sigma = ( sigma + floor )/2
79. return sigma
80.def black_imp_vol(F,k,t,r,c,option_type):
81. c_est = 0 # 期权价格估计值
82. top = 1 #波动率上限
83. floor = 0 #波动率下限
84. sigma = ( floor + top )/2 #波动率初始值
85. count = 0 # 计数器
86. while abs( c - c_est ) > 0.00000001:
87. c_est = black_value(F,k,t,r,sigma,option_type)
88. #根据价格判断波动率是被低估还是高估，并对波动率做修正
89. count += 1
90. if count > 1000: # 时间价值为0的期权是算不出隐含波动率的，因此迭代到一定次数就不再迭代了
91. sigma = 0
92. break
93. if c - c_est > 0: #f(x)>0
94. floor = sigma
95. sigma = ( sigma + top )/2
96. else:
97. top = sigma
98. sigma = ( sigma + floor )/2
99. return sigma
100.def cal_iv(df): # 计算主程序
101. option_list = df.ts_code.tolist()
102.
103. #df['F']=0
104. #for i in range(len(df)):
105. #df['F'].iloc[i]=df['s'].iloc[i]*exp(df['r'].iloc[i]*df['r'].iloc[i])
106. df = df.set_index('ts_code')
107. alist = []
108.
109. for option in option_list:
110. s = df.loc[option,'s']
111. #F = df.loc[option,'F']
112. k = df.loc[option,'k']
113. t = df.loc[option,'t']
114. r = df.loc[option,'r']
115. c = df.loc[option,'c']
116. q = df.loc[option,'q']
117. option_type = df.loc[option,'call_put']
118. sigma = bsm_imp_vol(s,k,t,r,c,q,option_type)
119. #sigma = black_imp_vol(F,k,t,r,c,option_type)
120. alist.append(sigma)
121. df['iv'] = alist
122. return df
123. #print(df)
124.#隐含波动率求解方法：牛顿迭代法
125.def bsm_value(s,k,t,r,sigma,option_type):
126. d1 = ( log( s/k ) + ( r + 0.5*sigma**2 )*t )/( sigma*sqrt(t) )
127. d2 = ( log( s/k ) + ( r - 0.5*sigma**2 )*t )/( sigma*sqrt(t) )
128. if option_type.lower() == 'c':
129. value = (s*stats.norm.cdf( d1) - k*exp( -r*t )*stats.norm.cdf( d2))
130. else:
131. value = k * exp(-r * t) * stats.norm.cdf(-d2) - s* stats.norm.cdf(-d1)
132. return value
133.
134.def bsm_vega(s,k,t,r,sigma):
135. d1 = log( s/k ) + ( r + 0.5*sigma**2 )*t/( sigma*sqrt(t) )
136. vega = s*stats.norm.cdf(d1,0.,1.)*sqrt(t)
137. #print('vega',vega)
138. return vega
139.
140.#牛顿迭代法求隐含波动率，迭代次数设为100
141.def bsm_call_imp_vol_newton(s, k, t, r, c, sigma_est,option_type, it = 10000):
142. for i in range(it):
143. sigma_est =sigma_est-((bsm_value(s,k,t,r,sigma_est,option_type) - c)/
144. bsm_vega(s, k, t, r, sigma_est))
145. return sigma_est
146.
147.def cal_iv(df): # 计算主程序
148. option_list = df.ts_code.tolist()
149. df_shibor6 = pro.shibor(date = date,fields = '6m')
150. rf6 = df_shibor6.iloc[0,0]/100
151. df.loc[df['t']<=(186/365),'r']=rf6
152. df_shibor5 = pro.shibor(date = date,fields = '3m')
153. rf5 = df_shibor5.iloc[0,0]/100
154. df.loc[df['t']<=(93/365),'r']=rf5
155. df_shibor4 = pro.shibor(date =date,fields = '1m')
156. rf4 = df_shibor4.iloc[0,0]/100
157. df.loc[df['t']<=(31/365),'r']=rf4
158. df = df.set_index('ts_code')
159. alist = []
160. for option in option_list:
161. s = df.loc[option,'s']
162. k = df.loc[option,'k']
163. t = df.loc[option,'t']
164. r = df.loc[option,'r']
165. c = df.loc[option,'c']
166. sigma_est=0.5
167. option_type = df.loc[option,'call_put']
168. sigma = bsm_call_imp_vol_newton(s, k, t, r, c, sigma_est, option_type,it = 1000)
169. alist.append(sigma)
170. df['iv'] = alist
171. return df
172.def data_pivot(df): # 数据透视
173. df_c=df[(df['call_put']=='C')&(df['k']>df['s'])]
174. df_p=df[(df['call_put']=='P')&(df['k']<=df['s'])]
175. df=df_c.append(df_p,sort=False)
176. df = df.reset_index()
177. df = df.drop(['ts_code','trade_date','c','s','r','call_put'],axis=1)
178. df['t'] = df['t']*240
179. df['t'] = df['t'].astype(int)
180. df = df.pivot_table(index=["k"],columns=["t"],values=["iv"])
181. df.columns = df.columns.droplevel(0)
182. df.index.name = None
183. df = df.reset_index()
184. df = df.rename(columns={'index':'k'})
185. return df
186.
187.def fitting(df): # 多项式拟合
188. col_list = df.columns
189. for i in range(df.shape[1]-1):
190. x_col = col_list[0]
191. y_col = col_list[i+1]
192. df1 = df.dropna(subset=[y_col])
193. x = df1.iloc[:,0]
194. y = df1.iloc[:,i+1]
195. degree = 2
196. weights = np.polyfit(x, y, degree)
197. model = np.poly1d(weights)
198. predict = np.poly1d(model)
199. x_given_list = df[pd.isnull(df[y_col]) == True][x_col].tolist()
200. # 所有空值对应的k组成列表
201. for x_given in x_given_list:
202. y_predict = predict(x_given)
203. df.loc[df[x_col]==x_given, y_col] = y_predict
204. return df
205.
206.def im_surface(df): # 波动率曲面作图
207. df = fitting(df)
208. df = df.set_index('k')
209. y = np.array(df.index)
210. x = np.array(df.columns)
211. fig = go.Figure(data=[go.Surface(z=df.values, x=x, y=y)])
212. fig.update_layout(scene = dict(
213. xaxis_title='剩余期限',
214. yaxis_title='执行价格',
215. zaxis_title='隐含波动率'),
216. width=1400,
217. margin=dict(r=20, b=10, l=10, t=10))
218. plotly.offline.plot(fig)
219.
220.def smile_plot(df): # 波动率微笑作图
221. df = df.set_index('k')
222. df = df.stack().reset_index()
223. df.columns = ['k', 'days', 'iv']
224. fig = px.line(df, x="k", y="iv", color="days",line_shape="spline")
225. plotly.offline.plot(fig)
226.
227.def main():
228. date = '20210401' #
229. df = extra_data(date) # 提取数据
230. df = data_clear(df,date) # 数据清洗
231. df = cal_iv(df) # 计算隐含波动率
232. df = data_pivot(df) # 数据透视表
233. smile_plot(df) # 波动率微笑
234. im_surface(df)# 波动率曲面
235.
236.if __name__ == '__main__':
237. main()
238. return df
(Python)绘制波动率曲面
于 2023-04-28 21:28:05 首次发布
8563
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
