yieldcurve.py-20190417

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Feb 25 08:51:13 2019

@author: vicky
"""

#导入数值计算库
import numpy as np
from numpy import linalg as la
from sklearn import linear_model
import pandas as pd
import math 
from sklearn.decomposition import PCA

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列

#from sklearn.decomposition import PCA
#x= np.mat(xtrain)
#pca=PCA(n_components=3)
#pca.fit(x)
#pc2=pca.fit_transform(x)

#------------PCA降维--------
from sklearn.decomposition import PCA
pca = PCA(n_components=0.95)
pca.fit(xtrain)
print(pca.explained_variance_ratio_) #方差贡献率
#print(pca.explained_variance_)#方差
print(pca.n_components_)#主成分个数
pc2=pca.transform(xtrain)
#p=np.array(df4)
#x=np.linalg.inv(df3)
#w=np.dot(np.array(df4),np.linalg.inv(np.array(df3)))

import matplotlib.pyplot as plt
plt.figure(figsize=(12,6),title= 'PCA Plot',fontsize=14)
plt.plot(pca.explained_variance_, 'k', linewidth=2)
plt.xlabel('n_components', fontsize=16)
plt.ylabel('explained_variance_', fontsize=16)
plt.show()


#-----------时间序列--------
df=pd.DataFrame(pc)
df.columns=['level','slope','curvation']
df.index = pd.to_datetime(data['trading_date'])  #将字符串转换成时间

#--------原始数据 时序图，看数据的大概分布
df.plot(figsize=(12,6),title= 'Monthly Yield Curve',fontsize=14)
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)

#原始数据 自相关性图
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
plot_acf(ts_level)
plt.show()
 
#原始数据 平稳性检测
from statsmodels.tsa.stattools import adfuller
adfuller(ts_level)
#返回值依次为：adf, pvalue p值， usedlag, nobs, critical values临界值 , icbest, regresults, resstore 
#adf分别大于3中不同检验水平的3个临界值，单位检测统计量对应的p值显著大于0.05，说明序列可以判定为非平稳序列



#进行ADF检验
#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)看检验统计量和临界值，检验统计量小于临界值的话，就拒绝原假设，认为序列是平稳的；大于的话，不能拒绝，认为是不平稳的

##滚动统计
#def rolling_statistics(timeseries):
#    #Determing rolling statistics
#    rolmean = pd.rolling_mean(timeseries, window=12)
#    rolstd = pd.rolling_std(timeseries, window=12)
# 
#    #Plot rolling statistics:
#    orig = plt.plot(timeseries, color='blue',label='Original')
#    mean = plt.plot(rolmean, color='red', label='Rolling Mean')
#    std = plt.plot(rolstd, color='black', label = 'Rolling Std')
#    plt.legend(loc='best')
#    plt.title('Rolling Mean & Standard Deviation')
#    plt.show(block=False)
# 
###ADF检验
#from statsmodels.tsa.stattools import adfuller
#def adf_test(timeseries):
#    rolling_statistics(timeseries)#绘图
#    print ('Results of Augment Dickey-Fuller Test:')
#    dftest = adfuller(timeseries, autolag='AIC')
#    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
#    for key,value in dftest[4].items():
#        dfoutput['Critical Value (%s)'%key] = value   #增加后面的显著性水平的临界值
#    print (dfoutput)
#    
#adf_test(ts_level)   #从结果中可以看到p值为0.1097>0.1，不能拒绝H0,认为该序列不是平稳序列


##自相关图ACF和偏相关图PACF
#import statsmodels.api as sm
#def acf_pacf_plot(ts_log_diff):
#    sm.graphics.tsa.plot_acf(ts_log_diff,lags=500) #ARIMA,q
#    sm.graphics.tsa.plot_pacf(ts_log_diff,lags=30) #ARIMA,p
#    
#acf_pacf_plot(ts_level)   #查看数据的自相关图和偏自相关图


#--------确定ARIMA的阶数
#----利用自相关图和偏自相关图
#一阶差分 时序图
ts_d = ts_level.diff().dropna()
ts_d.columns = ['diff']
ts_d.plot(figsize=(12, 6),title= 'Difference Curve', fontsize=14) 
plt.show()

#一阶差分 自相关图 和 偏自相关图
#plot_acf(ts_d)   
#plt.show()
#plot_pacf(ts_d)  
#plt.show()

#不限制lags值时图看不清，限制lags值再画
plot_acf(ts_d,lags=100) 
plt.show()
plot_pacf(ts_d,lags=100) 
plt.show()
#p=2,q=2

#一阶差分 平稳性检验 ADF检验结果
adfuller(ts_d)
#返回值依次为：adf, pvalue p值， usedlag, nobs, critical values临界值 , icbest, regresults, resstore 
#adf分别大于3中不同检验水平的3个临界值，单位检测统计量对应的p值显著小于0.05，说明序列可以判定为平稳序列


#----借助AIC、BIC统计量自动确定
#在statsmodels包里还有更直接的函数：
import statsmodels.tsa.stattools as st
order = st.arma_order_select_ic(ts_d,max_ar=5,max_ma=5,ic=['aic', 'bic', 'hqic'])
order.bic_min_order
#输出(1,0)

#我们常用的是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
from statsmodels.tsa.arima_model import ARMA
def proper_model(data_ts, maxLag): 
    init_aic = float("inf")
    init_p = 0
    init_q = 0
    init_properModel = None
    for p in np.arange(maxLag):
        for q in np.arange(maxLag):
            model = ARMA(data_ts, order=(p, q))
            try:
                results_ARMA = model.fit(disp=-1, method='css')
            except:
                continue
            aic = results_ARMA.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_d,5)


#---------建模预测
#ARIMA模型，d=1,p=2,q=2
#RSS是残差平方和

#----建模
from statsmodels.tsa.arima_model import ARIMA
model = ARIMA(ts_level, order=(2, 1, 2))  #ARIMA（p，d，q）
results_AR = model.fit(disp=-1)  #disp为-1代表不输出收敛过程的信息，True代表输出
results_AR.summary2()
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)#残差平方和


#----计算模型误差
#ARIMA拟合的其实是一阶差分ts_log_diff，predictions_ARIMA_diff[i]是第i个月与i-1个月的ts_log的差值。
#由于差分化有一阶滞后，所以第一个月的数据是空的，
predictions_diff = pd.Series(results_AR.fittedvalues, copy=True)
print(predictions_diff.head())
#累加现有的diff，得到每个值与第一个月的差分（同log底的情况下）。
#即predictions_ARIMA_diff_cumsum[i] 是第i个月与第1个月的ts_log的差值。
predictions_diff_cumsum = predictions_diff.cumsum()
#先ts_log_diff => ts_log=>ts_log => ts 
#先以ts_log的第一个值作为基数，复制给所有值，然后每个时刻的值累加与第一个月对应的差值(这样就解决了，第一个月diff数据为空的问题了)
#然后得到了predictions_ARIMA_log => predictions_ARIMA
predictions = pd.Series(ts_level.ix[0], index=ts_level.index)
predictions = predictions.add(predictions_diff_cumsum,fill_value=0)

plt.figure()
plt.plot(ts_level)
plt.plot(predictions)
plt.title('RMSE: %.4f'% np.sqrt(sum((predictions-ts_level)**2)/len(ts_level)))


#----预测未来走势
# forecast方法会自动进行差分还原，当然仅限于支持的1阶和2阶差分
forecast_n = 1000#预测未来30个天走势
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)

#预测从2017-12-03开始，也就是我们训练数据最后一个数据的后一个日期
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_level,label='Original',color='blue')
plt.plot(forecast_ARIMA, label='Forcast',color='red')
plt.legend(loc='best')
plt.title('forecast')