import pandas as pd import numpy as np import matplotlib.pyplot as plt from statsmodels.tsa.stattools import adfuller # 导入ADF检验函数 from statsmodels.tsa.seasonal import seasonal_decompose # 导入季节性分解函数，将数列分解为趋势、季节性和残差三部分 from statsmodels.stats.diagnostic import acorr_ljungbox # 导入白噪声检验函数 from statsmodels.graphics.tsaplots import plot_pacf, plot_acf # 导入自相关和偏自相关的绘图函数 from matplotlib.ticker import MaxNLocator # 导入自动查找到最佳的最大刻度函数 from statsmodels.tsa.arima_model import ARIMA # 导入ARIMA模型 from statsmodels.tsa.statespace.sarimax import SARIMAX from sklearn import svm # 单位根检验(test_stationarity，ADF检验)，用于检验序列是否是平稳的 # 季节性分解函数(seasonal_decompose)，通过分解后的趋势、季节性确认序列是否是平稳的 # 白噪声检验函数 # 自相关性和偏自相关性(acf_pacf)，通过截尾或拖尾的lag值，初步确认p,q。也可以用来检验序列是否平稳 # ADF检验：这是一种检查数据稳定性的统计测试。无效假设：时间序列是不稳定的。 # 测试结果由测试统计量和一些置信区间的临界值组成。 # 如果“测试统计量”少于“临界值”，我们可以拒绝无效假设，并认为序列是稳定的。 # 当p-value<0.05，且Test Statistic显著小于Critical Value (5%)时，数列稳定 # 主要看p-value，显著小于的判断不精确 def test_stationarity(timeseries, window=12): rolmean = timeseries.rolling(window=window, center=False).mean() rolstd = timeseries.rolling(window=window, center=False).std() # 旧版方法，即将被移除 # rolmean = pd.rolling_mean(timeseries, window=window) # rolstd = pd.rolling_std(timeseries, window=window) # 设置原始图，移动平均图和标准差图的式样 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() print('ADF检验结果：') dftest = adfuller(timeseries, autolag='AIC') # 使用减小AIC的办法估算ADF测试所需的滞后数 # 将ADF测试结果、显著性概率、所用的滞后数和所用的观测数打印出来 dfoutput = pd.Series(dftest[0:4], index=['Test Statistic', 'p-value', 'Num Lags Used', 'Num Observations Used']) for key, value in dftest[4].items(): dfoutput['Critical Value (%s)' % key] = value print(dfoutput) # 按照滑动均值将维经过指数转化的时间序列分为趋势(trend), 季节性(seasonality)和残差(residual)三部分 def decompose_plot(series, title=''): decomposition = seasonal_decompose(series) # 季节性分解函数 trend = decomposition.trend # 分解出的趋势，包含NaN值，原因不明 seasonal = decomposition.seasonal # 分解出的季节性 residual = decomposition.resid # 分解出的残差，包含NaN值，原因不明 fig = decomposition.plot() fig.set_size_inches(12, 6) fig.suptitle(title) fig.tight_layout() fig2 = acf_pacf(residual, title='Residuals', figsize=(12, 6)) # 分析残差的自相关，偏自相关 # test_stationarity(residual.dropna()) # Dropna后才能做稳定性检验 # 原数据的残差一般是不稳定的，经过差分后的数据残差可能是平稳的 # 定义一个画自相关，偏自相关的函数 # series 输入的时间序列 # lags 自相关和偏自相关函数的滞后取值范围 def acf_pacf(series, lags=40, title=None, figsize=(12, 6)): # 求自相关函数 fig = plt.figure(figsize=figsize) # figure指设置图形的特征。figsize为设置图形大小，单位为inch ax1 = fig.add_subplot(211) # 子图2行1列的第一张,大于10后用逗号分隔，如(3,4,10) ax1.set_xlabel('lags') # 横坐标为滞后值 ax1.set_ylabel('AutoCorrelation') # 纵坐标为自相关系数 ax1.xaxis.set_major_locator(MaxNLocator(integer=True)) # 设置主坐标轴为自动设置最佳的整数型坐标标签 plot_acf(series.dropna(), lags=lags, ax=ax1) # 没有title参数，需要删除title # 求偏自相关函数 ax2 = fig.add_subplot(212) ax2.set_xlabel('lags') ax2.set_ylabel('Partial AutoCorrelation') ax2.xaxis.set_major_locator(MaxNLocator(integer=True)) plot_pacf(series.dropna(), lags=lags, ax=ax2) # 没有title参数，需要删除title plt.tight_layout() # 设置为紧凑布局 plt.show() # #生成一个伪随机白噪声用于测试acorr_ljungbox的可靠性 # from random import gauss # from random import seed # from pandas import Series # # seed random number generator # # seed(1) #指定生成伪随机数的种子，指定后每次生成的随机数相同 # # create white noise series # whitenoise = Series([gauss(0.0, 1.0) for i in range(1000)])#创建一个高斯分布的白噪声 # acf_pacf(whitenoise, lags=40, title='White Noise Series') # print(u'白噪声检验结果为：', acorr_ljungbox(whitenoise, lags=1))#检验结果：平稳度，p-value。p-value>0.05为白噪声 # 骑自行车人数预测 df = pd.read_csv('D:\\portland-oregon-average-monthly-.csv') print(df.head(), '\nindex type:\n', type(df.index)) df['Month'] = pd.to_datetime(df['Month'], format='%Y-%m') # 索引并resample为月 indexed_df = df.set_index('Month') ts = indexed_df['riders'] ts = ts.resample('M').sum() # ts.plot(title='Monthly Num. of Ridership') # print(ts.head(), '\nindex type:\n', type(ts.index)) # #原始数据分析 # test_stationarity(ts) # decompose_plot(ts, title='ts decompose') # print(u'白噪声检验结果为：', acorr_ljungbox(ts, lags=1))# 是否为白噪声测试,p-value<0.05非白噪声。平稳且非白噪声函数，可用于时间序列建模 # acf_pacf(ts, lags=20, title='ts ACF&PACF')# 自相关和偏自相关初步确认p,q # 季节差分12个月 ts_sdiff = ts.diff(12) test_stationarity(ts_sdiff.dropna()) decompose_plot(ts_sdiff.dropna(), title='ts_sdiff decompose') print(u'白噪声检验结果为：', acorr_ljungbox(ts_sdiff.dropna(), lags=1)) # 是否为白噪声测试,p-value<0.05非白噪声。平稳且非白噪声函数，可用于时间序列建模 acf_pacf(ts_sdiff.dropna(), lags=20, title='ts_sdiff ACF&PACF') # 自相关和偏自相关初步确认p,q # 趋势差分1个月 ts_diff_of_sdiff = ts_sdiff.diff(1) test_stationarity(ts_diff_of_sdiff.dropna()) decompose_plot(ts_diff_of_sdiff.dropna(), title='ts_diff_of_sdiff decompose') print(u'白噪声检验结果为：', acorr_ljungbox(ts_diff_of_sdiff.dropna(), lags=1)) # 是否为白噪声测试,p-value<0.05非白噪声。平稳且非白噪声函数，可用于时间序列建模 acf_pacf(ts_diff_of_sdiff.dropna(), lags=20, title='ts_diff_of_sdiff ACF&PACF') # 自相关和偏自相关初步确认p,q # 得出ACF&PACF后，开始计算使用的p,q,P,Q # PACF(lag=k)=0，则p=k-1。如果PACF=0时，ACF仍然显著>0，则再增加一些p # ACF(lag=k)=0,则q=k-1。如果ACF=0时，PACF仍然显著>0，则再增加一些q # P: 如果季节差分后序列的ACF(lag=季节周期)显著为正, 考虑增加P # Q: 如果季节差分后序列的ACF(lag=季节周期)显著为负, 考虑增加Q # 绘出拟合后的ARIMA的拟合残差的ACF和PACF, 按上述规则调整p, q, P, Q # 通常P, Q最多为1 # 这里考虑p=0, q=0，季节周期=12, 季节差分后的ACF(12)显著为负, 可以考虑P=0, Q=1 pdq = (5, 1, 0) PDQ = (0, 1, 1, 12) # # SARIMA的参数: order = p, d, q, seasonal_order=P, D, Q, season的周期=12 # model = SARIMAX(ts, order=pdq,seasonal_order=PDQ).fit()# 已拟合的模型 # print(model.summary()) # ts.plot() # model.fittedvalues.plot(color='red') # residuals = pd.DataFrame(model.resid)# 拟合的数据和原始数据的残差 # acf_pacf(residuals, title='Residuals', figsize=(10,4)) # plt.title('RSS: %.4f'% sum((model.fittedvalues-ts)**2))# 计算拟残差平方和(RSS) # residuals.plot(kind='kde')# 残差的核密度估计 # print('\n\nResiduals Summary:\n', residuals.describe()) # plt.show() test_point = 30 # 测试点数 train_size = int(len(ts) - test_point) # 测试集大小 train, test = ts[:train_size], ts[train_size:] # 切分数据集 model_train = SARIMAX(train, order=pdq, seasonal_order=PDQ).fit() # 用于测试的已拟合的模型 predict_train = model_train.forecast(test_point) fitted_train = model_train.fittedvalues.append(predict_train) # 拟合曲线和预测曲线 model_run = SARIMAX(ts, order=pdq, seasonal_order=PDQ).fit() # 用于测试的已拟合的模型 predict_run = model_run.forecast(test_point) fitted_run = model_run.fittedvalues.append(predict_run) # 拟合曲线和预测曲线 # 通过支持向量机，将预测有偏差的残差再进行预测 residual = predict_train - test # for n in np.arange(3,15): #n从3到15测试最佳值 n = 13 residual_list = list() for i in np.arange(n - 1, len(residual)): # n-1之前的值无法计算 for j in np.arange(-(n - 1), 1): residual_list.append(residual[i + j]) # df_residual=pd.DataFrame(np.repeat(np.array([np.nan]),n).reshape(1,n)) #组成n+1列的数据框，用前n列来预测n+1列的残差 df_residual = pd.DataFrame() df_residual = df_residual.append(pd.DataFrame(np.array(residual_list).reshape(len(residual) - (n - 1), n))) df_residual.index = residual[n - 1:].index matrix_residual = df_residual.as_matrix() # 训练SVM data_train = matrix_residual[:int(0.5 * len(matrix_residual)), :] data_test = matrix_residual[int(0.5 * len(matrix_residual)):, :] x_train = data_train[:, 0:n - 1] y_train = data_train[:, n - 1] # .astype(int) # 必须是离散的int型才能预测 x_test = data_test[:, 0:n - 1] y_test = data_test[:, n - 1] # .astype(int) # 必须是离散的int型才能预测 model_svm = svm.SVR(epsilon=0.2) # (C=1.0, kernel='rbf', gamma='auto') model_svm.fit(x_train, y_train) y_pred = model_svm.predict(x_test[0]) y_pred = model_svm.predict(x_test[1]) y_pred = model_svm.predict(x_test[2]) y_pred = model_svm.predict(x_test[3]) residual_pred = y_test - y_pred print(y_pred) print('n=', n, '残差预测误差均值=', np.mean(residual_pred), '残差预测误差标准差=', np.std(residual_pred)) # # SVM最终训练 # x_final = matrix_residual[:, 0:n-1] # y_final = matrix_residual[:, n-1].astype(int) #必须是离散的int型才能预测 # x_test = data_test.as_matrix()[:, 0:n-1] # y_test = data_test.as_matrix()[:, n-1].astype(int) #必须是离散的int型才能预测 # model_svm = svm.SVC() # (C=1.0, kernel='rbf', gamma='auto') # model_svm.fit(x_final, y_final) # matrix_residual_new=model_svm.predict(x_test) # print (n) # print(matrix_residual_new) print('end') predict_train[-len(y_pred):]=predict_train[-len(y_pred):]-y_pred fitted_train_svm = model_train.fittedvalues.append(predict_train) # 拟合曲线和预测曲线 plt.figure(figsize=(12, 6)) plt.title('Riders') plt.plot(ts, 'o', label='observed') plt.plot(fitted_train, 'g', label='forecast', color='green') plt.plot(fitted_train_svm, 'g', label='forecast', color='pink') plt.plot(fitted_run, 'g', label='forecast', color='red') plt.legend(loc='upper left') plt.show() print('end')
07-24
1002
1002
10-13
2万+
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
