#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Mar 6 17:45:07 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
data=pd.read_csv('/Users/vicky/Downloads/国债数据.csv')
xtrain=np.array(data.ix[:,1:])
#看哪些列包含nan
np.where(np.isnan(xtrain))[1]
#去除包含nan的列
xtrain=np.delete(xtrain,np.where(np.isnan(xtrain))[1],1)
#-------------PCA降维--------
(n,p)=np.shape(xtrain)
Xsta=(xtrain-np.tile(np.mean(xtrain,axis=0),(n,1)))/np.tile(np.std(xtrain,axis=0),(n,1));#标准化样本数据 X
U,D,V=la.svd(Xsta)#SVD分解
D=np.asmatrix(D).T
lam=np.multiply(D,D)/(n-1) #特征根
W=V.T/np.tile(np.sum(abs(V.T),axis=0),(p,1)) #标准化的特征向量
f=lam/sum(lam)#方差贡献率=特征值/所有特征值总和
#检验是否可以主成分降维
F=np.cumsum(lam)/sum(lam)#累计方差贡献率=前i个特征值总和/所有特征值总和
Index=np.where(F[0,:]>0.95)#取累计方差贡献达到95%的q个主成分
q=np.array(Index)[1,0]+1 #主成分个数
Xpc=np.dot(xtrain,W) #降维后的xpc矩阵＝标准化后的特征向量＊原x矩阵
pc=Xpc[:,0:q]#主成分矩阵＝降维后x矩阵的前q列
#主成分方差贡献率图像
plt.figure(figsize=(12,6))
plt.plot(f, 'k', linewidth=2)
plt.xlabel('n_components', fontsize=16)
plt.ylabel('explained_variance_', fontsize=16)
plt.title('PCA Plot',fontsize=14,)
plt.show()
#-----------时间序列--------
df=pd.DataFrame(pc)
df.columns=['level','slope','curvation']
#df.index = pd.to_datetime(data['trading_date']) #将字符串转换成时间
df.index = data['trading_date']
#原始数据 时序图，看数据的大概分布
df.plot(figsize=(12,6),title= 'Monthly Yield Curve of level',fontsize=8)
plt.show()
#有数据异常
#--------删去异常数据
df.loc[df['curvation']>0.0127].index #'2013-06-07', '2013-06-08', '2013-06-09', '2013-06-19', '2013-06-20','2013-06-21'
df.loc[df['level']<-0.04].index # '2013-06-24', '2013-06-25'
df.loc[df['slope']>0.0045].index #'2013-06-14', '2013-06-26'
df=df.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)
#删去异常数据后的时序图
df.index = pd.to_datetime(df.index) #将index转换成datetime格式
df.plot(figsize=(12,6),title= 'Monthly Yield Curve',fontsize=8)
plt.show()
#生成pd.Series对象
ts_level=pd.Series(df['level'].values, index=df.index)
ts_slope=pd.Series(df['slope'].values, index=df.index)
ts_curvation=pd.Series(df['curvation'].values, index=df.index)
#-------确定差分阶数
#----level,一阶
#平稳性检测，未通过
adfuller(ts_level) #p值>0.05
#返回值依次为：adf, pvalue p值， usedlag, nobs, critical values临界值 , icbest, regresults, resstore
#单位检测统计量对应的p值显著大于0.05，不能拒绝H0,认为该序列不是平稳序列
#adf_test的返回值：
#Test statistic：代表检验统计量
#p-value：代表p值检验的概率
#Lags used：使用的滞后k，autolag=AIC时会自动选择滞后
#Number of Observations Used：样本数量
#Critical Value(5%) : 显著性水平为5%的临界值。
#(1)假设是存在单位根，即不平稳；
#(2)显著性水平，1%：严格拒绝原假设；5%：拒绝原假设，10%类推。
#(3)看P值和显著性水平a的大小，p值越小，小于显著性水平的话，就拒绝原假设，认为序列是平稳的；大于的话，不能拒绝，认为是不平稳的
#(4)看检验统计量和临界值，检验统计量小于临界值的话，就拒绝原假设，认为序列是平稳的；大于的话，不能拒绝，认为是不平稳的
#自相关性图，未通过
plot_acf(ts_level)
plt.show()
#白噪声检验，通过
acorr_ljungbox(ts_level, lags= 1)#差分序列的白噪声检验结果，返回统计量和 p 值
# 差分序列的白噪声检验结果： (array([200.4273461]), array([1.68490686e-45])) p值为第二项， 远小于 0.05
#p值小于0.05，拒绝纯随机的原假设，为非随机序列，可建模
#一阶level_d
ts_level_d = ts_level.diff().dropna()
ts_level_d.columns = ['diff']
#序列图
ts_level_d.plot(title= 'Difference Curve of Level')
plt.show()
#平稳性检测，通过
adfuller(ts_level_d)
#自相关性图，通过
plot_acf(ts_level_d)
plt.show()
#白噪声检验，通过
acorr_ljungbox(ts_level_d, lags= 1)#p值<0.05
#----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)
#自相关性图，通过
plot_acf(ts_slope_d)
plt.show()
#白噪声检验，通过
acorr_ljungbox(ts_slope_d, lags= 1)#p值<0.05
#----curvation，二阶
#平稳性检测，通过
adfuller(ts_curvation) #p值<0.05
#自相关性图，未通过
plot_acf(ts_curvation,lags=100)
plt.show()
#白噪声检验，通过
acorr_ljungbox(ts_curvation, lags=1)#p值<0.05
#一阶curvation_d
ts_curvation_d = ts_curvation.diff().dropna()
ts_curvation_d.columns = ['diff']
#序列图
ts_curvation_d.plot(title= 'Difference Curve of Curvation')
plt.show()
#平稳性检测，通过
adfuller(ts_curvation_d)
#自相关性图，通过
plot_acf(ts_curvation_d)
plt.show()
#白噪声检验，未通过
acorr_ljungbox(ts_curvation_d, lags= 1)#p值=0.013,，未通过随机性检验
#二阶curvation_d
ts_curvation_dd = ts_curvation_d.diff().dropna()
ts_curvation_dd.columns = ['diff']
#序列图
ts_curvation_dd.plot(title= 'Second Order Difference Curve of Curvation')
plt.show()
#平稳性检测，通过
adfuller(ts_curvation_dd)
#自相关性图，通过
plot_acf(ts_curvation_dd)
plt.show()
#白噪声检验，通过
acorr_ljungbox(ts_curvation_dd, lags= 1)#p值<0.01
#--------确定ARIMA的阶数
#----利用自相关图和偏自相关图
def threeplots(ts,titlename):
#时序图
ts.plot(figsize=(12,6),title= 'Difference Curve of '+titlename,fontsize=8)
plt.show()
#自相关图 和 偏自相关图
#不限制lags值时图看不清，限制lags值再画
plot_acf(ts,lags=100)
plt.title('Acf Plot of '+titlename)
plt.show()
plot_pacf(ts,lags=100)
plt.title('Pacf Plot of '+titlename)
plt.show()
return 0
threeplots(ts_level_d,'Level') #p=3,d=1,q=1
threeplots(ts_slope_d,'Slope') #p=2,d=1,q=2
threeplots(ts_curvation_dd,'Curvation') #p=2,d=1,q=5
#----借助AIC、BIC统计量自动确定
#在statsmodels包里还有更直接的函数：
order = st.arma_order_select_ic(ts_level_d,max_ar=5,max_ma=5,ic=['aic', 'bic', 'hqic'])
order.aic_min_order
#输出(1,0)
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)
order = st.arma_order_select_ic(ts_curvation_dd,max_ar=5,max_ma=5,ic=['aic', 'bic', 'hqic'])
order.aic_min_order
#输出(1,1)
#我们常用的是AIC准则，AIC鼓励数据拟合的优良性但是尽量避免出现过度拟合(Overfitting)的情况。所以优先考虑的模型应是AIC值最小的那一个模型。
#为了控制计算量，我们限制AR最大阶不超过5，MA最大阶不超过5。 但是这样带来的坏处是可能为局部最优。
#timeseries是待输入的时间序列，是pandas.Series类型，max_ar、max_ma是p、q值的最大备选值。
#order.bic_min_order返回以BIC准则确定的阶数，是一个tuple类型
##借助AIC统计量自动确定,可换成bic
def proper_model(data_ts, maxLag):
init_aic = float("inf")
init_p = 0
init_q = 0
init_properModel = None
for p in range(1,maxLag):
for q in range(1,maxLag):
model = ARMA(data_ts, order=(p, q))
try:
results_ARMA = model.fit(disp=-1, method='css')
except:
continue
aic = results_ARMA.aic
print('p:',p,'q:',q,'aic:',aic)
if aic < init_aic:
init_p = p
init_q = q
init_properModel = results_ARMA
init_aic = aic
return init_aic, init_p, init_q, init_properModel
proper_model(ts_level_d,5) #1,1
proper_model(ts_slope_d,5) #1,1
proper_model(ts_curvation_dd,10) #1,1
def findmodel(ts,maxLag):
#init_aic = float("inf")
init_p = 0
init_q = 0
init_properModel = None
init_wucha=float("inf")
for p in range(1,maxLag):
for q in range(1,maxLag):
model = ARIMA(ts, order=(p, 1, q)) #ARIMA(p，d，q)
try:
results_AR = model.fit(disp=-1) #disp为-1代表不输出收敛过程的信息，True代表输出
#print('p:',p,'q:',q,'AIC:',results_AR.aic)
except:
continue
predictions_diff = pd.Series(results_AR.fittedvalues, copy=True)
predictions_diff_cumsum = predictions_diff.cumsum() #从第一个开始累计想加
predictions = pd.Series(ts.ix[0], index=ts.index)
predictions = predictions.add(predictions_diff_cumsum,fill_value=0)
#偏差率，平均绝对误差
wucha=sum(abs((predictions-ts)/ts))/len(ts)
print('p:',p,'q:',q,'wucha:',wucha)
#aic = results_AR.aic
if wucha < init_wucha:
init_p = p
init_q = q
init_properModel = results_AR
init_wucha = wucha
return init_p, init_q, init_properModel,init_wucha
findmodel(ts_level,6) #3,2
findmodel(ts_slope,6) #4,4
findmodel(ts_curvation_d,5) #1,1
def findmodel(ts,maxLag):
#init_aic = float("inf")
init_p = 0
init_q = 0
init_properModel = None
init_wucha=float("inf")
for p in range(1,maxLag):
for q in range(1,maxLag):
model = ARIMA(ts, order=(p, 0, q)) #ARIMA(p，d，q)
try:
results_AR = model.fit(disp=-1) #disp为-1代表不输出收敛过程的信息，True代表输出
#print('p:',p,'q:',q,'AIC:',results_AR.aic)
except:
continue
predictions = pd.Series(results_AR.fittedvalues, copy=True)
#偏差率，平均绝对误差
wucha=sum(abs((predictions-ts)/ts))/len(ts)
print('p:',p,'q:',q,'wucha:',wucha)
#aic = results_AR.aic
if wucha < init_wucha:
init_p = p
init_q = q
init_properModel = results_AR
init_wucha = wucha
return init_p, init_q, init_properModel,init_wucha
#---------建模预测
#ARIMA(p,d,q)模型，RSS是残差平方和
#----建模
ts=ts_level
ts_d=ts_level_d
ts=ts_slope
ts_d=ts_slope_d
ts=ts_curvation
ts_d=ts_curvation_d
def arima_ts(ts,ts_d,p,d,q):
model = ARIMA(ts, order=(p, d, q)) #ARIMA（p，d，q）
results_AR = model.fit(disp=-1) #disp为-1代表不输出收敛过程的信息，True代表输出
# results_AR.summary2()
# print('AIC:',results_AR.aic)
plt.plot(ts_d)
plt.plot(results_AR.fittedvalues, color='red')#红色线代表差分的预测值
rss=sum(((results_AR.fittedvalues-ts_d).dropna())**2)
print('RSS:',rss)
plt.title('RSS: %.4f'% rss)#残差平方和
plt.show()
#计算模型历史预测值和偏差率
predictions_diff = pd.Series(results_AR.fittedvalues, copy=True)
predictions_diff_cumsum = predictions_diff.cumsum() #从第一个开始累计想加
predictions = pd.Series(ts.ix[0], index=ts.index)
predictions = predictions.add(predictions_diff_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((predictions-ts)/ts)
wucha=sum(res)/len(ts)
print('wucha',wucha)
#正负一致率
rate=signrate(ts,predictions)
print(rate)
#回测曲线与原始数据对比图
plt.figure()
plt.plot(ts)
plt.plot(predictions,color='red')
# plt.title('RMSE: %.4f'% np.sqrt(sum((predictions-ts)**2)/len(ts)))
plt.title('MAE: %.4f'% wucha)
plt.show()
return (predictions,wucha,results_AR,rate)
#差分正负个数一致率
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_level,wucha_level,results_AR_level,rate_level=arima_ts(ts_level,ts_level_d,3,1,1) #0.0753,0.634
pre_slope,wucha_slope,results_AR_slope,rate_slope=arima_ts(ts_slope,ts_slope_d,2,1,2) #5.06,0.604
#pre_curvation,wucha_curvation,results_AR_curvation=arima_ts(ts_curvation_d,ts_curvation_dd,2,1,5) #
#---level
def arima_ts_curvation(ts,ts_d,p,d,q):
model = sm.tsa.statespace.SARIMAX(ts_d,order=(p,0,q),seasonal_order=(0,0,0,0),
enforce_stationarity=False, enforce_invertibility=False).fit()
#We can use SARIMAX model as ARIMAX when seasonal_order is (0,0,0,0) .
#注意这里用的是ts_curvation_d算，因为SARIMAX包的fittedvalues直接给出原始拟合值而不是差分拟合值
#差分还原
predictions_d = pd.Series(model.fittedvalues, copy=True)
predictions_d_cumsum = predictions_d.cumsum() #从第一个开始累计想加
predictions = pd.Series(ts.ix[0], index=ts_curvation.index)
predictions = predictions.add(predictions_d_cumsum,fill_value=0)
return(predictions,model)
pre_level2,results_AR_level2=arima_ts_curvation(ts_level,ts_level_d,1,1,3)
plt.figure()
plt.plot(pre_level)
plt.plot(pre_level2)
pre_level-pre_level2
differ=results_AR_level.fittedvalues
differ2=results_AR_level2.fittedvalues
plt.figure()
plt.plot(differ)
plt.plot(differ2)
differ-differ2
model = sm.tsa.statespace.SARIMAX(ts_level,order=(1,1,3),seasonal_order=(0,0,0,0),
enforce_stationarity=False, enforce_invertibility=False).fit()
pre_level22=model.fittedvalues
plt.figure()
plt.plot(ts_level)
plt.plot(pre_level)
plt.plot(pre_level2)
plt.plot(pre_level22)
#---curvation
model = ARIMA(ts_curvation, order=(2,2,4)) #ARIMA（p，d，q）
results_AR_curvation = model.fit(disp=-1) #disp为-1代表不输出收敛过程的信息，True代表输出
#results_AR_curvation.summary2()
plt.plot(ts_curvation_dd)
plt.plot(results_AR_curvation.fittedvalues, color='red')#红色线代表差分的预测值
rss=sum(((results_AR_curvation.fittedvalues-ts_curvation_dd).dropna())**2)
print('RSS:',rss)
plt.title('RSS: %.4f'% rss)#残差平方和
plt.show()
#计算模型历史预测值和偏差率
predictions_dd = pd.Series(results_AR_curvation.fittedvalues, copy=True)
predictions_dd_cumsum = predictions_dd.cumsum() #cumsum从第一个开始累计加和
predictions_d = pd.Series(ts_curvation_d.ix[0], index=ts_curvation_d.index)
predictions_d = predictions_d.add(predictions_dd_cumsum,fill_value=0) #同序列值相加
plt.figure()
plt.plot(ts_curvation_d)
plt.plot(predictions_d,color='red')
plt.show()
predictions_d_cumsum = predictions_d.cumsum() #cumsum从第一个开始累计加和
predictions = pd.Series(ts_curvation.ix[0], index=ts_curvation.index)
predictions = predictions.add(predictions_d_cumsum,fill_value=0) #同序列值相加
#偏差率，平均绝对误差
res=abs((predictions-ts_curvation)/ts_curvation)
wucha=sum(res)/len(ts_curvation)
print('wucha',wucha)
#回测曲线与原始数据对比图
plt.figure()
plt.plot(ts_curvation)
plt.plot(predictions,color='red')
plt.title('MAE: %.4f'% wucha)
plt.show()
#预测结果还原
predict_ts = results_AR_curvation.predict()
# 一阶差分还原
diff_shift_ts = ts_curvation_d.shift(1) #序号不变，值往下移一个
diff_recover_1 = predict_ts.add(diff_shift_ts)
# 再次一阶差分还原
rol_shift_ts = ts_curvation_dd.shift(1)
diff_recover = diff_recover_1.add(rol_shift_ts)
plt.figure()
plt.plot(ts_curvation)
plt.plot(diff_recover,color='red')
plt.show()
#偏差率，平均绝对误差
res=abs((predictions-ts_curvation)/ts_curvation)
wucha=sum(res)/len(ts_curvation)
print('wucha',wucha)
#回测曲线与原始数据对比图
plt.figure()
plt.plot(ts_curvation_d)
plt.plot(predictions_d,color='red')
# plt.title('RMSE: %.4f'% np.sqrt(sum((predictions-ts)**2)/len(ts)))
plt.title('MAE: %.4f'% wucha)
plt.show()
#----------白噪声检验
results_AR=results_AR_level
results_AR=results_AR_slope
def jianyan(results_AR):
#残差的自相关图
resid = results_AR.resid#残差
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = plot_acf(resid.values.squeeze(), lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = plot_pacf(resid, lags=40, ax=ax2)
#做D-W检验
print('D-W检验结果：',sm.stats.durbin_watson(results_AR.resid.values))
#0≤DW≤４,其中 DW＝O＝＞ρ＝１即存在正自相关性; DW＝４＜＝＞ρ＝－１即存在负自相关性; DW＝２＜＝＞ρ＝０即不存在（一阶）自相关性
#所以，当DW值显著的接近于O或4时，则存在自相关性，而接近于2时，则不存在（一阶）自相关性
#检验结果是1.998，说明不存在自相关性
#观察是否符合正态分布
resid = results_AR.resid#残差
fig = plt.figure()
ax = fig.add_subplot(111)
fig = qqplot(resid, line='q', ax=ax, fit=True)
#Ljung-Box检验
r,q,p = sm.tsa.acf(resid.values.squeeze(), qstat=True)
data = np.c_[range(1,41), r[1:], q, p]
table = pd.DataFrame(data, columns=['lag', "AC", "Q", "Prob(>Q)"])
print(table.set_index('lag'))
#6期和12期的p值均大于0.05，所以残差序列不存在自相关性
return 0
jianyan(results_AR_level)
jianyan(results_AR_slope)
#----预测未来走势
# forecast方法会自动进行差分还原，当然仅限于支持的1阶和2阶差分
forecast_n = 1000#预测未来n天走势
forecast_AR = results_AR.forecast(forecast_n)[0] #返回预测值，标准差，置信区间
print(forecast_AR)
stderr_AR=sum(results_AR.forecast(forecast_n)[1])
print(stderr_AR)
#将预测的数据和原来的数据绘制在一起，为了实现这一目的，我们需要增加数据索引，使用开源库arrow:
import arrow
def get_date_range(start, limit, level='day',format='YYYY-MM-DD'):
start = arrow.get(start, format)
result=(list(map(lambda dt: dt.format(format) , arrow.Arrow.range(level, start,limit=limit))))
dateparse2 = lambda dates:pd.datetime.strptime(dates,'%Y-%m-%d')
return map(dateparse2, result)
#预测从训练数据最后一个数据的后一个日期开始
new_index = get_date_range('2019-02-25', forecast_n)
forecast_ARIMA= pd.Series(forecast_AR, copy=True, index=new_index)
print(forecast_ARIMA.head())
#绘图如下
plt.plot(ts,label='Original',color='blue')
plt.plot(forecast_ARIMA, label='Forcast',color='red')
plt.legend(loc='best')
plt.title('forecast')
#### 预测
predicts_ARIMA = results_AR.predict(start =pd.to_datetime('2012-01-10'),
end=pd.to_datetime('2019-02-25'), dynamic = True,typ='levels')
#print (len(predicts_ARIMA))
compare_ARIMA = pd.DataFrame()
compare_ARIMA['original'] = ts
compare_ARIMA['predict_ARIMA'] = predicts_ARIMA
compare_ARIMA
compare_ARIMA.plot()
#### 拟合优度检验
delta_ARIMA = results_AR.fittedvalues - ts
score_ARIMA = 1 - delta_ARIMA.var()/ts[1:].var()
print (score_ARIMA)
扫一扫
1116
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
