######------------------step1: 获取数据------------------#用requests库下载CSMAR网站的沪深300指数文件import urllib
import requests
import zipfile
import pandas as pd
url = 'http://119.147.213.32/ssis/csmardata/20170616/kquv1emc.zip'
urllib.request.urlretrieve(url, "kquv1emc.zip")
#获取文件下载后在本地的存放地址import os
address="%s" % os.getcwd()
file_address=address+'/IDX_Idxtrd.txt'#解压文件
zfile = zipfile.ZipFile('kquv1emc.zip','r')
try: #使用try-except可以跳过这个文件读写中的警示信息,直接得到我们想要的DataFramefor filename in zfile.namelist():
data = zfile.read(filename)[90::] #90的长度为第一行的长度,不读第一行(因为原文件第一行是变量行)with open(filename,'w') as f:
f.write(data.decode('utf-16'))
except:
import pandas as pd
import numpy as np
raw0=pd.read_table(file_address,names=['date','value','lowest','highest','return'])
raw0.head(3)
######------------------step2: 数据清理------------------
raw=raw0
raw.isnull().any()# 缺失值的查找:没有发现缺失值
raw.describe() #检查最大值最小中值有无异常## 1)变量设置:收益率R (使用对数价差RT=ln(Pt/Pt-1))import math
lt1=list(raw.value.values)
lt2=[math.log(lt1[i+1]/lt1[i]) for i in range(len(lt1)-1)]
lt3=[0]
[lt3.append(i) for i in lt2]
#说明:可以观察到我们用对数价差得到的日回报RT与原始数据return几乎一致,但RT相对更为保守,因此我们采用RT
raw['RT']=lt3
raw['return']=raw['return']/100## 2) 变量设置:滞后收益率lag_R
lt4=[0]
[lt4.append(i) for i in lt3[0:-1]]
raw['lag_RT']=lt4
## 3) 变量设置:指数日波动率sigma#说明:由于得不到每日指数的高频交易数据,因此我们只能采用上下限区间的极差来近似估计波动率sigma
raw['square_sigma']=(raw.highest/raw.lowest-1)*(raw.highest/raw.lowest-1)
## 4) 变量设置:滞后日波动率lag_sigma
lt5=[0]
lt6=list(raw.square_sigma.values)
[lt5.append(i) for i in lt6[0:-1]]
raw['square_lag_sigma']=lt5
raw.head()
## 5)变量设置:虚拟变量D的设置# 说明:用变量D表示股指期货的引入与否,1代表引入,0代表未引入
lt7=[]
for i in raw.date:
if str(i)>='2010-04-16':
lt7.append(1)
else:
lt7.append(0)
raw['D']=lt7
## 变量的标准中心化处理#原则:(观测值-均值)/标准差
raw['nom_RT']=(raw.RT-raw.RT.mean())/raw.RT.std()
raw['nom_lag_RT']=(raw.lag_RT-raw.lag_RT.mean())/raw.lag_RT.std()
raw['nom_square_sigma']=(raw.square_sigma-raw.square_sigma.mean())/raw.square_sigma.std()
raw['nom_square_lag_sigma']=(raw.square_lag_sigma-raw.square_lag_sigma.mean())/raw.square_lag_sigma.std()
raw.head()
######------------------step3: 数据探查 & 描述统计------------------import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.tsa.stattools as ts
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.arima_model import ARMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
## 描述统计:观察各指标#raw.describe() ## 描述统计:画出股指引入前后市场波动率对比图# 对主数据表raw进行拆分,before表为股指期货引入前的表,after表为股指期货引之后的表
before=raw[0:800]
before.set_index('date',inplace=True)
after=raw[800::]
after.set_index('date',inplace=True)
# 画出波动率对比图
Fig=plt.figure(figsize=(20,5),facecolor='white') # ————> 对应原文 Figure1 :波动率对比图
Ax1=Fig.add_subplot(121)
before['square_lag_sigma'].plot(color='gray',alpha=0.7,title='Before the introduction of stock index futures')
Ax2=Fig.add_subplot(122)
plt.ylim(0.000,0.016)
after['square_lag_sigma'].plot(color='cornflowerblue',alpha=0.7,title='After the introduction of stock index futures')
Fig.show()
## 数据探查:获得股指收益率的时间序列time_df
raw2=raw #使用备份的原始数据,命名为raw2
raw2.set_index('date',inplace=True) #将raw2改为用日期索引
time_df=raw2.iloc[:,4:5]
time_df.head()
## 数据探查:股指收益率的时间序列图
f=plt.figure(facecolor='white')
plt.title('Figure 1')
ax1=f.add_subplot(111)
time_df.plot(color='cornflowerblue', ax=ax1)
plt.show()
## 数据探查:自相关和偏相关图检验
raw2=raw #使用备份的原始数据,命名为raw2
time_df=raw2.iloc[:,4:5]
time_df.head()
f=plt.figure(facecolor='white')
plt.title("Figure1: test") # ————> 对应原文 Figure 2:自相关检验图
ax1=f.add_subplot(211)
plot_acf(time_df, lags=31, ax=ax1)
ax2=f.add_subplot(212)
plot_pacf(time_df, lags=31, ax=ax2)
plt.show()
## 数据探查:ADF检验
result_adf=ts.adfuller(time_df['RT'].values)
print("ADF检验结果:",result_adf)
## 数据探查:ARIMA(1,1)模型# 依据:根据figure2的结果,确定滞后阶数为1阶,因此采用ARIMA(1,1)模型来预测
raw3=raw #使用备份的原始数据,命名为raw3
time_df=raw3.iloc[:,4:5]
time_df.head()
from statsmodels.tsa.arima_model import ARMA
model=ARMA(time_df.values, order=(1, 1))
result_arma=model.fit(disp=-1, method='css') #进行ARMA回归,结果保存在result_arma中# 生成time_df2表,为下面的画图作准备
predict_ts=result_arma.predict()
predict_ts=list(predict_ts)
lt8=[0,0]
[lt8.append(i) for i in predict_ts[0:-1]]
time_df2=time_df
time_df2['pre']=lt8
# 画预测拟合效果对比图
Fig=plt.figure(figsize=(20,10),facecolor='white') # ————> 对应原文 Figure 3:预测拟合效果对比图
Ax1=Fig.add_subplot(121)
time_df2['RT'].plot(color='gray',alpha=0.7,title='Orignal')
Ax2=Fig.add_subplot(122)
plt.ylim(-0.08,0.10)
time_df2['pre'].plot(color='cornflowerblue',alpha=0.7,title='Predict')
Fig.show() #观察本图的拟合效果,发现与实际情况基本一致,说明拟合效果良好######------------------step 4: 实证分析(模型回归)------------------## 对GARCH-M的1)式进行回归import statsmodels.api as sm
x=raw.iloc[:,5:9]
x.drop(['square_sigma'],axis=1,inplace=True)
X=x.values
X=sm.add_constant(X)
Y=raw['RT'].values
est=sm.OLS(Y,X)
est=est.fit()
est.summary() # ————> 对此处产生的回归结果进行截图,对应文章中的表1## 说明:对于文章中GARCH-M模型的2)式,由于通过了前面的平稳性检验,因此可以认为残差服从正态随机游走过程## 对GARCH-M的3)式进行回归# 得到GARCH-M的1)式回归方程的残差序列resid、残差平方序列resid_square和滞后残差平方序列lag_resid_square
original_y=raw['RT'].values
predict_y=est.predict()
resid=original_y-predict_y #残差序列resid
resid_square=(original_y-predict_y)*(original_y-predict_y) #残差平方序列resid_square
resid_square2=list(resid_square)
lt9=[0]
[lt9.append(i) for i in resid_square2[0:-1]]
raw['lag_resid_square']=lt9 #滞后残差平方序列lag_resid_square# 开始回归
x1=raw.iloc[:,7::]
x1.drop(['nom_RT'],axis=1,inplace=True)
x1.drop(['nom_lag_RT'],axis=1,inplace=True)
x1.drop(['nom_square_sigma'],axis=1,inplace=True)
x1.drop(['nom_square_lag_sigma'],axis=1,inplace=True)
y1=raw.iloc[:,6:7]
X1=x1.values
X1=sm.add_constant(X1)
Y1=y1.values
import statsmodels.api as sm
est1=sm.OLS(Y1,X1)
est1=est1.fit()
est1.summary() # ————> 对此处产生的回归结果进行截图,对应文章中的表2## 对GJR-GARCH模型进行回归# 设置GJR-GARCH回归的其他变量:D2、cross
D2_lst=[] #建立第二个虚拟变量D2的list(规则:残差大于等于0,D2=1,代表好的冲击;残差小于0,D2=0,代表坏的冲击)for i in resid:
if i>=0:
D2_lst.append(1)
else:
D2_lst.append(0)
raw['D2']=D2_lst
raw['cross']=raw['D2'].values*resid_square #设置交互变量cross:及残差平方和与虚拟变量D2的乘积# 开始回归
x2=raw.iloc[:,12::]
x2.drop(['D2'],axis=1,inplace=True)
y2=raw.iloc[:,11:12]
X2=x2.values
X2=sm.add_constant(X2)
Y2=y2.values
est2=sm.OLS(Y2,X2)
est2=est2.fit()
est2.summary() # ————> 对此处产生的回归结果进行截图,对应文章中的表3
OLS Regression Results
Dep. Variable:
y
R-squared:
0.095
Model:
OLS
Adj. R-squared:
0.094
Method:
Least Squares
F-statistic:
68.19
Date:
Fri, 26 May 2017
Prob (F-statistic):
6.43e-42
Time:
17:55:34
Log-Likelihood:
5366.7
No. Observations:
1944
AIC:
-1.073e+04
Df Residuals:
1940
BIC:
-1.070e+04
Df Model:
3
Covariance Type:
nonrobust
coef
std err
t
P>|t|
[95.0% Conf. Int.]
const
-0.0058
0.001
-7.940
0.000
-0.007 -0.004
x1
0.1724
0.022
7.928
0.000
0.130 0.215
x2
0.3722
0.029
12.731
0.000
0.315 0.430
x3
0.0025
0.001
3.310
0.001
0.001 0.004
Omnibus:
358.175
Durbin-Watson:
2.232
Prob(Omnibus):
0.000
Jarque-Bera (JB):
1687.798
Skew:
-0.798
Prob(JB):
0.00
Kurtosis:
7.277
Cond. No.
101.
import pandas as pd
raw0=pd.read_table(IDX_Idxtrd.txt,names=['date','value','lowest','highest','return'])