yieldcurve4.py-20190328

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Mar 12 17:13:38 2019

@author: vicky
"""

import numpy as np
from numpy import linalg as la
import pandas as pd
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import adfuller
import statsmodels.tsa.stattools as st
from statsmodels.tsa.arima_model import ARMA
from statsmodels.tsa.arima_model import ARIMA
import copy
import statsmodels.api as sm
from statsmodels.graphics.api import qqplot
import warnings
from statsmodels.stats.diagnostic import acorr_ljungbox
import arch 
from scipy import  stats
import statsmodels.api as sm  # 统计相关的库
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import arch  # 条件异方差模型相关的库
from arch import arch_model
import datetime as dt

#----检验slope用ar还是arima
#slope平稳性检测，不平稳，需一阶差分
adfuller(ts_slope) #p值<0.05
#自相关性图，未通过
plot_acf(ts_slope)
plt.show() 
#白噪声检验，通过
acorr_ljungbox(ts_slope, lags= 1)#p值<0.05

#一阶slope_d
ts_slope_d = ts_slope.diff().dropna()
ts_slope_d.columns = ['diff']
#序列图
ts_slope_d.plot(title= 'Difference Curve of Slope') 
plt.show()
#平稳性检测，通过
adfuller(ts_slope_d) #p值<0.05
#自相关性图，通过
plot_acf(ts_slope_d) #明显三角对称，非平稳序列 
plt.show() 
#白噪声检验，通过
acorr_ljungbox(ts_slope_d, lags= 1)#p值<0.05
#用一阶差分



#------------ar(1)-arch(4)模型
#对slope_d建立 AR(p) 模型，并判断阶数
fig = plt.figure(figsize=(15,6))
ax1 = fig.add_subplot(111)
fig = plot_pacf(ts_slope_d,lags = 50,ax=ax1) #p=1
#根据图，建立 AR(1)模型,即均值方程
model = ARMA(ts_slope_d, order=(1,0)).fit(disp=-1)  #ARMA（p，q）

#----ARCH 效应检验
#计算均值方程的残差：a_t = r_t - mu_t,并作出残差和残差平方图像：
at = ts_slope_d - model.fittedvalues
at2 = np.square(at)
plt.figure(figsize=(10,6))
plt.subplot(211)
plt.plot(at,label = 'at')
plt.legend()
plt.subplot(212)
plt.plot(at2,label='at^2')
plt.legend(loc=0)

#残差自相关下性检验：原假设H0:序列没有相关性，备择假设H1：序列具有相关性
#Ljung-Box检验
acf,q,p = sm.tsa.acf(at2,nlags=20,qstat=True)  ## 计算自相关系数及p-value
out = np.c_[range(1,21), acf[1:], q, p]
output=pd.DataFrame(out, columns=['lag', "AC", "Q", "P-value"])
output = output.set_index('lag')
output
#p<0.05，认为残差序列具有相关性，符合异方差的假设，因此具有 ARCH 效应。
#参数均显著，说明残差有长期自相关性，需要很大的q值。

#----ARCH模型建立
#ARCH模型的阶次，可以用{at^2}序列的偏自相关函数PACF来确定：
figure = plt.figure(figsize=(20,5))
ax1 = figure.add_subplot(111)
fig = plot_pacf(at2,lags = 50, ax=ax1) #q=1

#train_ARCH = ts_slope['2012-01-04':'2019-01-04']
#test_ARCH = ts_slope['2019-01-05':'2019-02-25']

#train_ARCH = ts_slope[:-10]
#test_ARCH = ts_slope[-10:]

#粗略选择均值模型为AR(1)模型，波动率模型选择ARCH(4)模型
train_ARCH = ts_slope_d

am_ARCH = arch.arch_model(train_ARCH, mean='AR', lags=1, vol='ARCH', p=4) 
res_ARCH = am_ARCH.fit()
res_ARCH.summary()
res_ARCH.params
#日收益率期望大约为0.005%， AdjR^2值＝0.969，模型拟合效果不错

#ARCH模型的拟合效果
res_ARCH.hedgehog_plot()
#可以看出，虽然具体值差距挺大，但是均值和方差的变化相
res_ARCH.hedgehog_plot(type='volatility')
res_ARCH.plot()


pre_ARCH = res_ARCH.forecast(horizon=1,start=1,method='simulation')
#pre_mean=pre_ARCH.mean.iloc[0:]
pre_mean=pre_ARCH.mean
print(pre_mean)
resid=pre_ARCH.residual_variance
print(resid)
var=pre_ARCH.variance
print(var)

pre =pd.Series(pre_mean['h.1'], index=train_ARCH.index)
#偏差率，平均绝对误差
res=(abs((train_ARCH-pre)/pre)).dropna()
wucha=sum(res)/len(res.dropna())
print('wucha',wucha)

#回测曲线与原始数据对比图
plt.figure()
plt.plot(train_ARCH)
plt.plot(pre, color='red')#红色线代表差分的预测值
rss=sum(((pre-train_ARCH).dropna())**2) 
print('RSS:',rss)
plt.title('RSS: %.4f'% rss)#残差平方和
plt.show()

#计算模型历史预测值和偏差率
pres_d = pd.Series(pre, copy=True)
pres_d_cumsum = pres_d.cumsum() #从第一个开始累计想加
pres = pd.Series(ts_slope.ix[0], index=ts_slope.index) 
pres = pres.add(pres_d_cumsum,fill_value=0)

#    predictions.index=data.trading_date
#    predictions=predictions.drop(labels=['2013-06-07', '2013-06-08', '2013-06-09','2013-06-13','2013-06-14', '2013-06-26',
#    '2013-06-19','2013-06-20','2013-06-21','2013-06-24', '2013-06-25','2013-06-27'],axis=0)
#    predictions.index=pd.to_datetime(predictions.index)

#偏差率，平均绝对误差
res=abs((pres-ts_slope)/ts_slope)
wucha_arch=sum(res)/len(res)
print('wucha',wucha_arch)
#6.43

#回测曲线与原始数据对比图
plt.figure()
plt.plot(ts_slope)
plt.plot(pres,color='red')
#    plt.title('RMSE: %.4f'% np.sqrt(sum((predictions-ts)**2)/len(ts)))
plt.title('MAE: %.4f'% wucha)
plt.show()





#------------arma-garch模型

#对slope_d建立 ARMA(p,q) 模型，并判断阶数
plot_acf(ts_slope_d,lags = 50) #p=1
plt.show()
plot_pacf(ts_slope_d,lags = 50) #q=1
plt.show()
#根据图，建立 arma(1,1)模型,即均值方程
model = ARMA(ts_slope_d, order=(1,1)).fit(disp=-1)  #ARMA（1，1）

#----GARCH效应检验
#计算arma均值方程的残差：a_t = r_t - mu_t,并作出残差和残差平方图像：
at = ts_slope_d - model.fittedvalues
at2 = np.square(at)
plt.figure()
plt.subplot(211)
plt.plot(at,label = 'at')
plt.legend()
plt.subplot(212)
plt.plot(at2,label='at^2')
plt.legend(loc=0)

#残差自相关下性检验：原假设H0:序列没有相关性，备择假设H1：序列具有相关性
#Ljung-Box检验
acf,q,p = sm.tsa.acf(at2,nlags=20,qstat=True)  ## 计算自相关系数及p-value
out = np.c_[range(1,21), acf[1:], q, p]
output=pd.DataFrame(out, columns=['lag', "AC", "Q", "P-value"])
output = output.set_index('lag')
output
#p<0.05，认为残差序列具有相关性，符合异方差的假设，因此具有 ARCH 效应。
#参数均显著，说明残差有长期自相关性，需要很大的q值。


#----GARCH模型建立
#ARCH模型的阶次，可以用{at^2}序列的偏自相关函数PACF来确定：
plot_acf(at2,lags = 50) #p=4
plot_pacf(at2,lags = 50) #q=1

#粗略选择均值模型为ARMA(1,1)模型，波动率模型选择GARCH(4,0,1)模型
#train = ts_slope_d

model2=arch_model(ts_slope_d, mean='HAR',lags=[1,1,1],vol='GARCH', p=4, o=0, q=1)
res_garch=model2.fit()
print(res_garch.summary())
print(res_garch.params)

#标准偏差和条件波动率
res_garch.plot()
#残差图显示方差没了集群效应，波动率（条件方差图）显示条件方差的起伏与波动率的大小一致
res_garch.hedgehog_plot()
#可以看出，方差的还原还不错。

#res_garch.conditional_volatility
#res_garch.hedgehog_plot(type='volatility')
#res_garch.plot()


pre_garch = res_garch.forecast(horizon=1,start=1,method='simulation')
#pre_mean=pre_ARCH.mean.iloc[0:]
pre_garch_mean=pre_garch.mean
print(pre_garch_mean)
resid=pre_garch.residual_variance
print(resid)
var=pre_garch.variance
print(var)

pre_g_d=pd.Series(pre_garch_mean['h.1'], index=ts_slope.index)

#偏差率，平均绝对误差
res=(abs((ts_slope_d-pre_g_d)/pre_g_d)).dropna()
wucha=sum(res)/len(res.dropna())
print('wucha',wucha)

#回测曲线与原始数据对比图
plt.figure()
plt.plot(ts_slope_d)
plt.plot(pre_g_d, color='red')#红色线代表差分的预测值
rss=sum(((pre_g_d-ts_slope_d).dropna())**2) 
print('RSS:',rss)
plt.title('RSS: %.4f'% rss)#残差平方和
plt.show()

#计算模型历史预测值和偏差率
pres_d2 = pd.Series(pre_g_d, copy=True)
pres_d_cumsum2 = pres_d2.cumsum() #从第一个开始累计想加
pres_g = pd.Series(ts_slope.ix[0], index=ts_slope.index) 
pres_g = pres_g.add(pres_d_cumsum2,fill_value=0)

#偏差率，平均绝对误差
res=abs((pres_g-ts_slope)/ts_slope)
wucha_garch=sum(res)/len(res)
print('wucha',wucha_garch)
#7.72

#回测曲线与原始数据对比图
plt.figure()
plt.plot(ts_slope)
plt.plot(pres_g,color='red')
#    plt.title('RMSE: %.4f'% np.sqrt(sum((predictions-ts)**2)/len(ts)))
plt.title('MAE: %.4f'% wucha)
plt.show()

#差分正负个数一致率
def signrate(a,b):
    a=np.sign(a.diff())
    b=np.sign(b.diff())
    c=(a-b).dropna()
    return (c[c==0].count())/len(c)

r_arch=signrate(ts_slope,pres)
print(r_arch)
#0.9791549295774647

r_garch=signrate(ts_slope,pres_g)
print(r_garch)
#0.979

r_arima=signrate(ts_slope,pre_slope)
print(r_arima)
#0.604

#对比arima，arch和garch
plt.figure()
plt.plot(ts_slope,label = 'Slope')
plt.plot(pre_slope,color='orange',label = 'ARIMA')
plt.plot(pres,color='red',label = 'AR-ARCH')
plt.plot(pres_g,color='green',label = 'ARMA-GARCH')
plt.legend(loc='upper right')