基于ARMA差分还原的客流量时间序列预测 完整代码数据 毕设

# -*- coding:utf-8 -*- # import statsmodels.tsa.stattools as st from statsmodels.tsa.arima_model import ARMA import numpy as np import pandas as pd import matplotlib.pylab as plt from datetime import datetime from statsmodels.graphics.tsaplots import plot_acf, plot_pacf # acf pacf from statsmodels.tsa.stattools import adfuller # adf检验 # 示例数据 date = [10930, 10318, 10595, 10972, 7706, 6756, 9092, 10551, 9722, 10913, 11151, 8186, 6422, 6337, 11649, 11652, 10310, 12043, 7937, 6476, 9662, 9570, 9981, 9331, 9449, 6773, 6304, 9355, 10477, 10148, 10395, 11261, 8713, 7299, 10424, 10795, 11069, 11602, 11427, 9095, 7707, 10767, 12136, 12812, 12006, 12528, 10329, 7818, 11719, 11683, 12603, 11495, 13670, 11337, 10232, 13261, 13230, 15535, 16837, 19598, 14823, 11622, 19391, 18177, 19994, 14723, 15694, 13248, 9543, 12872, 13101, 15053, 12619, 13749, 10228, 9725, 14729, 12518, 14564, 15085, 14722, 11999, 9390, 13481, 14795, 15845, 15271, 14686, 11054, 10395] date_index = pd.date_range('2017-01-01', periods=90, Freq='D') date = pd.Series(date, index=date_index) timeSeries = date def draw_ts(timeSeries): f = plt.figure(facecolor='white') timeSeries.plot(color='blue') plt.show() def draw_acf_pacf(timeSeries, lags=31): '''acf pacf 默认为31阶差分 ''' # bestdiff = best_diff(timeSeries) # 差分运算 # ts_diff = timeSeries.diff(int(bestdiff)) # ts_diff.dropna(inplace=True) # 绘图 f = plt.figure(facecolor='white') ax1 = f.add_subplot(211) plot_acf(timeSeries, lags=lags, ax=ax1) ax2 = f.add_subplot(212) plot_pacf(timeSeries, lags=lags, ax=ax2) plt.show() def outliers(timeSeries, threshold=3): '''异常值处理 采用移动中位数方法''' ts_rolling = rolling_median(timeSeries, window=3, center=True) ts_dropna = ts_rolling.fillna(method='bfill').fillna(method='ffill') return ts_dropna def best_diff(timeSeries, maxdiff=8): '''判断最优差分''' p_set = {} for i in range(0, maxdiff): temp = pd.DataFrame(timeSeries).copy() # 每次循环前，重置 if i == 0: temp['diff'] = temp[temp.columns[0]] else: temp['diff'] = temp[temp.columns[0]].diff(i) temp = temp.drop(temp.iloc[:i].index) # 差分后，前几行的数据会变成nan，所以删掉 pvalue = adfuller(temp['diff'])[1] p_set[i] = pvalue p_df = pd.DataFrame.from_dict(p_set, orient="index") p_df.columns = ['p_value'] i = 0 while i < len(p_df): if p_df['p_value'][i] < 0.01: bestdiff = i break i += 1 return bestdiff def draw_trend(timeSeries, size=12): '''移动平均图''' f = plt.figure(facecolor='white') # 对size个数据进行移动平均 rol_mean = timeSeries.rolling(window=size).mean() # 对size个数进行加权移动平均 rol_weighted_mean = pd.ewma(timeSeries, span=size) # ggplot style import matplotlib matplotlib.style.use('ggplot') timeSeries.plot(color='blue', label='Original') rol_mean.plot(color='red', label='Rolling Mean') rol_weighted_mean.plot(color='black', label='Weighted Rolling Mean') plt.legend(loc='best') plt.title('Rolling Mean') plt.show() def testStationarity(timeSeries): # adf 检验 dftest = adfuller(timeSeries) # 语义描述 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 return dfoutput def diff(timeSeries): '''ADF检验时间序列不是平稳序列 进行差分运算 并绘制图形''' date = timeSeries fig = plt.figure(facecolor='white', figsize=(12, 8)) ax1 = fig.add_subplot(221) date_plot = date.plot() # diff_1 ax2 = fig.add_subplot(222) diff_1 = timeSeries.diff(1) date_diff_1_plot = diff_1.plot() # diff_2 ax3 = fig.add_subplot(223) diff_2 = timeSeries.diff(2) date_diff_2_plot = diff_2.plot() # diff_3 ax4 = fig.add_subplot(224) diff_3 = timeSeries.diff(3) date_diff_3_plot = diff_3.plot() plt.show() return diff_1, diff_2, diff_3 # 数据平稳后，需要对模型定阶，即确定p、q的阶数 def proper_model(data_ts, maxLag): '''对于个数不多的时序数据，我们可以通过观察自相关图和偏相关图来进行模型识别. 对于数量较多的时序数据，依据BIC准则识别模型的p, q值，通常认为BIC值越小的模型相对更优。 BIC准则，它综合考虑了残差大小和自变量的个数，残差越小BIC值越小，自变量个数越多BIC值越大''' import sys from statsmodels.tsa.arima_model import ARMA init_bic = sys.maxint 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 bic = results_ARMA.bic if bic < init_bic: init_p = p init_q = q init_properModel = results_ARMA init_bic = bic return init_bic, init_p, init_q, init_properModel def disassemble(data_ts): # 分解 from statsmodels.tsa.seasonal import seasonal_decompose decomposition = seasonal_decompose(data_ts, model='additive') # 趋势项 trend = decomposition.trend teend_plot = decomposition.trend.plot(color='black', label='Trend') # 季节因素项 seasonal = decomposition.seasonal seasonal_plot = decomposition.seasonal.plot(color='red', label='Seasonal') # 剩余项 residual = decomposition.resid residual_plot = decomposition.resid.plot(color='blue', label='Residual') plt.show() return trend, seasonal, residual if __name__ == "__main__": # 1.可视化时间序列 draw_ts(timeSeries) # 2.可视化 acf pacf draw_acf_pacf(timeSeries) # 3.异常值处理 判断最优差分 timeSeries = outliers(timeSeries) maxLag = best_diff(timeSeries) # 模型识别 # rol_mean = ts.rolling(window=12).mean() # rol_mean.dropna(inplace=True) # diff_1 # ts_diff_1 = rol_mean.diff(1) # ts_diff_1.dropna(inplace=True) # diff_2 # ts_diff_2 = ts_diff_1.diff(1) # ts_diff_2.dropna(inplace=True) # 2阶差分 adf # testStationarity(ts_diff_2) # 拟合ARMA order = st.arma_order_select_ic(timeSeries, max_ar=3, max_ma=3, ic=['aic', 'bic', 'hqic']) order.bic_min_order model = ARMA(ts_diff_2, order=(1, 1)) result_arma = model.fit(disp=-1, method='css') # 模型拟合完后，我们就可以对其进行预测了 predict_ts = result_arma.predict() # 预测的Y值 还原 '''由于ARMA拟合的是经过相关预处理后的数据，故其 预测值需要通过相关逆变换进行还原（根据相关的变换情况）''' # 一阶差分还原 diff_shift_ts = ts_diff_1.shift(1) diff_recover_1 = predict_ts.add(diff_shift_ts) # 再次一阶差分还原 rol_shift_ts = rol_mean.shift(1) diff_recover = diff_recover_1.add(rol_shift_ts) # 移动平均还原 rol_sum = ts_log.rolling(window=11).sum() rol_recover = diff_recover * 12 - rol_sum.shift(1) # 对数还原 log_recover = np.exp(rol_recover) log_recover.dropna(inplace=True) # 使用均方根误差（RMSE）来评估模型样本内拟合的好坏 ts = ts[log_recover.index] # 过滤没有预测的记录 plt.figure(facecolor='white') log_recover.plot(color='blue', label='Predict') ts.plot(color='red', label='Original') plt