时间序列模型

时间序列预测模型

• 以上海客源预测为例

定义函数

• 自相关图 | 偏自相关图
• 稳定性检验（单位根）
• 白噪声检验

# 自相关和偏自相关图，默认阶数为31
def draw_acf_pacf(ts, subtitle, lags=31):
    print('自相关图和偏自相关图，maxlags={0}'.format(lags))
    f = plt.figure(facecolor='white', figsize=(18,4))
    ax1 = f.add_subplot(121)
    plot_acf(ts, lags=lags, ax=ax1, title='ACF\n{}'.format(subtitle))
    ax2 = f.add_subplot(122)
    plot_pacf(ts, lags=lags, ax=ax2, title='PACF\n{}'.format(subtitle))
    plt.show()

from statsmodels.tsa.stattools import adfuller
# 平稳性检验(单位根检验)
def test_stationarity(ts):
    print('Results of Dickey-Fuller Test:')
    dftest = adfuller(ts, autolag='AIC')
    dfoutput = pd.Series(dftest[0:4],index=['统计量值','P值','阶数','观测值数'])
    for key,value in dftest[4].items():
        dfoutput['Critical Value({})'.format(key)] = value
    print(dfoutput)

# 白噪声检验（ljung-box test）
def randomness(ts, lags=31):
    rdtest = acorr_ljungbox(ts, lags=lags)
    rddata = np.c_[range(1,lags+1),rdtest[1:][0]]
    rdoutput = pd.DataFrame(rddata, columns=['lags','p-value'])
    return rdoutput.set_index('lags')

Box-Jenkins建模

平稳性

图示：明显的周期性，周期大概是一年，即365天；杭州还有明显的趋势性

相关图和自相关图

图示：看不出明显的指数下降和截尾，加时序图，肯定就是非平稳序列了，但还是做下单位跟检验

单位根检验

北京P值小于0.05，但看统计值，虽小于5%的置信水平，但距离99%的置信区间很近。上海则99%的置信度下接受H0，杭州更明显，5%的置信水平下接受H0。所以3个就作为非平稳序列吧

差分

图示：看起来平稳很多

# 先进行步长差分，去掉周期性(365)
df_bj1 = df_bj.北京.diff(365)
df_sh1 = df_sh.上海.diff(365)
df_hz1 = df_hz.杭州.diff(365)
df_bj1[0:364]

单位根检验：还是差挺大的

# 进行阶次差分
df_bj2 = np.diff(df_bj1, n=1)
df_sh2 = np.diff(df_sh1, n=1)
df_hz2 = np.diff(df_hz1, n=1)

非常好，p值非常小，是妥妥的平稳序列了

AIC判断模型

best_aic = np.inf#设初始值
best_order = None
best_mdl = None

pq_rng = range(2)
d_rng = range(1)#1阶差分后已经是平稳序列了，所以d没有必要了，但这里还是写出来，便于下面写ARIMA模型的order参数

for p in pq_rng:
    for d in d_rng:
        for q in pq_rng:
            try:
                tmp_mdl = smt.ARIMA(df_hz2,order=(p,d,q)).fit()
                tmp_aic = tmp_mdl.aic
                print('aic: {:6.5f}| order: {}'.format(tmp_aic,(p,d,q)))
                if tmp_aic < best_aic:
                    best_aic = tmp_aic
                    best_order = (p,d,q)
                    best_mdl = tmp_mdl
            except: 
                continue

print('\naic: {:6.5f}| order: {}'.format(best_aic,best_order))

SARIMAX 季节性的ARIMA模型

df_bjarima_final = smt.SARIMAX(df_bj,order=(1,1,1),seasonal_order=(0,1,0,365)).fit()

Holt-Winters建模

北京

df_bjtrain_HW = ExponentialSmoothing(np.asarray(df_bjtrain['北京']), seasonal="add", seasonal_periods =365).fit()

# 测试集一共438期
pre_bj = df_bjtrain_HW.forecast(438)
# pre

# 均方根误差
np.sqrt(np.mean(np.square(abs(df_bjtest.北京-pre_bj))))

df_bjtest['HW'] = pre_bj
df_bjtest

走势和变化还可以

上海

可以看到，21年旺季，市场是超预期的，但10月末开始走向平稳，预测也更准确。