slope.py-20190423

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Mar 27 17:06:32 2019

@author: vicky
"""

import numpy as np
import pandas as pd
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 statsmodels.api as sm
from statsmodels.graphics.api import qqplot
from statsmodels.stats.diagnostic import acorr_ljungbox
import arch 
from arch import arch_model

#----检验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
#用一阶差分


#------------1. arma(1,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=4

#粗略选择均值模型为AR(1)模型，波动率模型选择ARCH(4)模型
def arch_ts_slope(ts,ts_d,n_ar,n_arch):
    am_ARCH = arch.arch_model(ts_d, mean='AR', lags=n_ar, vol='ARCH', p=n_arch) 
    res_ARCH = am_ARCH.fit()
    #res_ARCH.summary()
    #res_ARCH.params
    
    #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
    #resid=pre_ARCH.residual_variance ＃残差
    #var=pre_ARCH.variance ＃方差
    pre =pd.Series(pre_mean['h.1'], index=ts_d.index)
    
    #偏差率，平均绝对误差
    res=(abs((ts_d-pre)/pre)).dropna()
    wucha=sum(res)/len(res.dropna())
    print('wucha',wucha)
    
    #回测曲线与原始数据对比图
    plt.figure()
    plt.plot(ts_d)
    plt.plot(pre, color='red')#红色线代表差分的预测值
    rss=sum(((pre-ts_d).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.ix[0], index=ts.index) 
    pres = pres.add(pres_d_cumsum,fill_value=0)
    
    #偏差率，平均绝对误差
    res=abs((pres-ts)/ts)
    wucha_arch=sum(res)/len(res)
    print('wucha',wucha_arch)
    #6.43
    
    #回测曲线与原始数据对比图
    plt.figure()
    plt.plot(ts)
    plt.plot(pres,color='red')
    plt.title('MAE: %.4f'% wucha)
    plt.show()
    
    r_arch=signrate(ts,pres)
    print(r_arch)
    
    return(pres,res_ARCH,wucha_arch,r_arch)

#差分正负个数一致率
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)

pre_arch_slope,results_ARCH_slope,wucha_arch_slope,rate_arch_slope=arch_ts_slope(ts_slope,ts_slope_d,1,4)  
#误差6.427，一致率0.979



#------------2. arma(1,1)-garch(1,0,4)模型------------
#对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) #q=4
plot_pacf(at2,lags = 50) #p=1

#粗略选择均值模型为ARMA(1,1)模型，波动率模型选择GARCH(4,0,1)模型
def garch_ts_slope(ts,ts_d,p_arima,d_arima,q_arima,p_garch,o_garch,q_garch):
    model2=arch_model(ts_d, mean='HAR',lags=[p_arima,d_arima,q_arima],vol='GARCH',
                      p=p_garch, o=o_garch, q=q_garch)
    res_garch=model2.fit()
    #print(res_garch.summary())
    #print(res_garch.params)
    
    #标准偏差和条件波动率
    res_garch.plot()
    #残差图显示方差没了集群效应，波动率（条件方差图）显示条件方差的起伏与波动率的大小一致
    res_garch.hedgehog_plot()
    #可以看出，方差的还原还不错。
    
    pre_garch = res_garch.forecast(horizon=1,start=1,method='simulation')
    pre_garch_mean=pre_garch.mean
    #resid=pre_garch.residual_variance
    #var=pre_garch.variance
    
    pre_g_d=pd.Series(pre_garch_mean['h.1'], index=ts.index)
    
    #偏差率，平均绝对误差
    res=(abs((ts_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_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.ix[0], index=ts.index) 
    pres_g = pres_g.add(pres_d_cumsum2,fill_value=0)

    #偏差率，平均绝对误差
    res=abs((pres_g-ts)/ts)
    wucha_garch=sum(res)/len(res)
    print('wucha',wucha_garch)
    #7.72
    
    #回测曲线与原始数据对比图
    plt.figure()
    plt.plot(ts)
    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()
    
    #差分正负个数一致率
    r_garch=signrate(ts,pres_g)
    print('rate',r_garch)
    #0.979

    return(pres_g,res_garch,wucha_garch,r_garch)

pre_garch_slope,result_garch_slope,wucha_garch_slope,rate_garch_slope=garch_ts_slope(ts_slope,ts_slope_d,1,1,1,1,0,4)
#误差6.427,一致率0.979


#-------------3. arima(2,1,2)----------
#----确定arima系数，函数来源level.py
threeplots(ts_slope_d,'Slope') #p=1,d=1,q=1

order = st.arma_order_select_ic(ts_slope_d,max_ar=5,max_ma=5,ic=['aic', 'bic', 'hqic'])
order.aic_min_order
#aic输出(3,3)

proper_model(ts_slope_d,5) #1,1

findmodel(ts_slope,6) #4,4

#----计算arima(2，1，2）误差率
pre_slope,wucha_slope,results_AR_slope,rate_slope=arima_ts_level(ts_slope,ts_slope_d,1,1,1) 
#误差5.0052,差分一致率0.601

#白噪声检验
jianyan(results_AR_slope)#通过检验

#------------三种模型对比：ar-arch,arma-garch,arima
#对比arima，arch和garch
plt.figure()
plt.plot(ts_slope,label = 'Slope')
plt.plot(pre_slope,color='orange',label = 'ARIMA')
plt.plot(pre_arch_slope,color='red',label = 'AR-ARCH')
plt.plot(pre_garch_slope,color='green',label = 'ARMA-GARCH') #曲线和ar-arch是一样的
plt.legend(loc='upper right')

#综合误差和正负一致率情况，选ar(1)-arch(4)