参考链接：常用7种时间序列预测模型
运用ARIMA进行时间序列建模的基本步骤：
- 1）加载数据：构建模型的第一步当然是加载数据集。
- 2）预处理：根据数据集定义预处理步骤。包括创建时间戳、日期/时间列转换为d类型、序列单变量化等。
- 3）序列平稳化：为了满足假设，应确保序列平稳。这包括检查序列的平稳性和执行所需的转换。
- 4）确定d值：为了使序列平稳，执行差分操作的次数将确定为d值。
- 5）创建ACF和PACF图：这是ARIMA实现中最重要的一步。用ACF PACF图来确定ARIMA模型的输入参数。
- 6）确定p值和q值：从上一步的ACF和PACF图中读取p和q的值。
- 7）拟合ARIMA模型：利用我们从前面步骤中计算出来的数据和参数值，拟合ARIMA模型。
- 8）在验证集上进行预测：预测未来的值。
- 9）计算RMSE：通过检查RMSE值来检查模型的性能，用验证集上的预测值和实际值检查RMSE值。
ARMA模型公式：
信息准则定阶
- AIC(Akaike Information Criterion)
L是数据的似然函数，k=1表示模型考虑常数c，k=0表示不考虑。最后一个1表示算上误差项，所以其实第二项就是2乘以参数个数。
- AICc（修正过的AIC）
- BIC(Bayesian Information Criterion)
注意事项：
- 信息准则越小，说明参数的选择越好，一般使用AICc或者BIC。
- 差分d，不要使用信息准则来判断，因为差分会改变了似然函数使用的数据，使得信息准则的比较失去意义，所以通常用别的方法先选择出合适的d。
- 信息准则的好处是可以在用模型给出预测之前，就对模型的超参做一个量化评估，这对批量预测的场景尤其有用，因为批量预测往往需要在程序执行过程中自动定阶。
数据平稳性检验
# 测试序列平稳性
from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries):
#Determing rolling statistics
rolmean = timeseries.rolling(window=12).mean()
rolstd = timeseries.rolling(window=12).std()
#Plot rolling statistics:
fig = plt.figure(figsize=(12, 8))
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()
#Perform Dickey-Fuller test:
#迪基-福勒检验还可以扩展为增广迪基-福勒检验（Augmented Dickey-Fuller test），简称ADF检验，可检验模型是否存在单位根（unit root）
print('Results of 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)
要将非平稳时间序列转为平稳序列，有如下几种方法：
- Deflation by CPI
- Logarithmic（取对数）
- First Difference（一阶差分）
- Seasonal Difference（季节差分）
- Seasonal Adjustment
这里会尝试取对数、一阶查分、季节差分三种方法，先进行一阶差分，去除增长趋势后检测稳定性：
可以看到图形上看上去变稳定了，但p-value的并没有小于0.05。再来看看12阶查分（即季节查分），看看是否稳定：
从图形上，比一阶差分更不稳定（虽然季节指标已经出来了），我们再来将一阶查分和季节查分合并起来，再来看下：
可以看到，在一阶差分后再加上季节差分，数据终于稳定下了。P-value小于了0.05，并且Test Statistic的值远远小于Critical Value (5%)的值。
对原始数据取对数有两个用处：第一个是将指数增长转为线性增长；第二可以平稳序列的方差
取对数后，从图形上看有些稳定了，但p-value还是没有接近0.05，没能达稳定性要求。在此基础上我们再来做下一阶差分，去除增长趋势：
还是没能满足需求，再来试试在对数基础上进行一阶差分+季节差分的效果：
可以看到稳定了满足需求了，但是稳定性相比不取对数的情况下，并没有提升而是下降了。
#原始序列 test_stationarity(df['riders']) #一阶差分 df['first_difference'] = df['riders'].diff(1) test_stationarity(df['first_difference'].dropna(inplace=False)) #季节差分 df['seasonal_difference'] = df['riders'].diff(12) test_stationarity(df['seasonal_difference'].dropna(inplace=False)) #一阶差分+季节差分 df['seasonal_first_difference'] = df['first_difference'].diff(12) test_stationarity(df['seasonal_first_difference'].dropna(inplace=False)) #对数 df['riders_log']= np.log(df['riders']) test_stationarity(df['riders_log']) #对数+一阶差分 df['log_first_difference'] = df['riders_log'].diff(1) test_stationarity(df['log_first_difference'].dropna(inplace=False)) #对数+季节差分 df['log_seasonal_difference'] = df['riders_log'].diff(12) test_stationarity(df['log_seasonal_difference'].dropna(inplace=False)) #对数+一阶差分+姐姐查分 df['log_seasonal_first_difference'] =df['log_first_difference'].diff(12) test_stationarity(df['log_seasonal_first_difference'].dropna(inplace=False))
如何确定差分阶数和常数项：
1.假如时间序列在很高的lag（10以上）上有正的ACF，则需要很高阶的差分
2.假如lag-1的ACF是0或负数或者非常小，则不需要很高的差分，假如lag-1的ACF是-0.5或更小，序列可能overdifferenced。
3.最优的差分阶数一般在最优阶数下标准差最小（但并不是总数如此）
4.模型不查分意味着原先序列是平稳的；1阶差分意味着原先序列有个固定的平均趋势；二阶差分意味着原先序列有个随时间变化的趋势
5.模型没有差分一般都有常数项；有两阶差分一般没常数项；假如1阶差分模型非零的平均趋势，则有常数项
如何确定AR和MA的阶数：
1.假如PACF显示截尾或者lag-1的ACF是正的（此时序列仍然有点underdifferenced），则需要考虑AR项；PACF的截尾项表明AR的阶数
2.假如ACF显示截尾或者lag-1的ACF是负的（此时序列有点overdifferenced），则需要加MA项，ACF的截尾项表明AR的阶数
3.AR和MA可以相互抵消对方的影响，所以假如用AR-MA模型去拟合数据，仍然需要考虑加些AR或MA项。尤其在原先模型需要超过10次的迭代去converge。
假如在AR部分有个单位根（AR系数和大约为1），此时应该减少一项AR增加一次差分
假如在MA部分有个单位根（MA系数和大约为1），此时应该减少一项AR减少一次差分
假如长期预测出现不稳定，则可能AR、MA系数有单位根
如何确定季节性部分：
1.假如序列有显著是季节性模式，则需要用季节性差分，但是不要使用两次季节性差分或总过超过两次差分（季节性和非季节性）
2.假如差分序列在滞后s处有正的ACF，s是季节的长度，则需要加入SAR项；假如差分序列在滞后s处有负的ACF，则需要加入SMA项，如果使用了季节性差异，则后一种情况很可能发生，如果数据具有稳定和合乎逻辑的季节性模式。如果没有使用季节性差异，前者很可能会发生，只有当季节性模式不稳定时才适用。应该尽量避免在同一模型中使用多于一个或两个季节性参数（SAR + SMA），因为这可能导致过度拟合数据和/或估算中的问题。
模型构建
from statsmodels.tsa.arima_model import ARIMA # 1,1,2 ARIMA Model model = ARIMA(df.value, order=(1,1,2)) model_fit = model.fit(disp=0) print(model_fit.summary())
中间的表格列出了训练得到的模型各项和对应的系数，如果系数很小，且‘P>|z|’ 列下的P-Value值远大于0.05，则该项应该去掉，比如左图中的ma部分的第二项，系数是-0.0010,P-Value值是0.998,那么可以重建模型为ARIMA(1,1,1)，从右图可以看到，修改阶数后的模型的AIC等信息准则都有所降低
通常会检查模型拟合的残差序列，即训练数据原本的序列减去训练数据上的拟合序列后的序列。该序列越符合随机误差分布(均值为0的正态分布)，说明模型拟合的越好，否则，说明还有一些因素模型未能考虑。
# Plot residual errors residuals = pd.DataFrame(model_fit.resid) fig, ax = plt.subplots(1,2) residuals.plot(title="Residuals", ax=ax[0]) residuals.plot(kind='kde', title='Density', ax=ax[1]) plt.show() # Actual vs Fitted model_fit.plot_predict(dynamic=False) plt.show()
from statsmodels.tsa.stattools import acf
# Create Training and Test
train = df.value[:85]
test = df.value[85:]
# Build Model
model = ARIMA(train, order=(3,2,1))
# model = ARIMA(train, order=(1, 1, 1)) #预测效果很差，选用上面的效果有提升
fitted = model.fit(disp=-1)
# Forecast
fc, se, conf = fitted.forecast(15, alpha=0.05) # 95% conf
# Make as pandas series
fc_series = pd.Series(fc, index=test.index)
lower_series = pd.Series(conf[:, 0], index=test.index)
upper_series = pd.Series(conf[:, 1], index=test.index)
# Plot
plt.figure(figsize=(12,5), dpi=100)
plt.plot(train, label='training')
plt.plot(test, label='actual')
plt.plot(fc_series, label='forecast')
plt.fill_between(lower_series.index, lower_series, upper_series,
color='k', alpha=.15)
plt.title('Forecast vs Actuals')
plt.legend(loc='upper left', fontsize=8)
plt.show()
AUTO-ARIMA模型 （python - pmdarima；r - forecast）
通过预测结果来推断模型阶数的好坏毕竟还是耗时耗力了些，一般可以通过计算AIC或BIC的方式来找出更好的阶数组合。pmdarima模块的auto_arima方法就可以让我们指定一个阶数上限和信息准则计算方法，从而找到信息准则最小的阶数组合。
from statsmodels.tsa.arima_model import ARIMA
import pmdarima as pm
df = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/wwwusage.csv', names=['value'], header=0)
model = pm.auto_arima(df.value, start_p=1, start_q=1,
information_criterion='aic',
test='adf', # use adftest to find optimal 'd'
max_p=3, max_q=3, # maximum p and q
m=1, # frequency of series
d=None, # let model determine 'd'
seasonal=False, # No Seasonality
start_P=0,
D=0,
trace=True,
error_action='ignore',
suppress_warnings=True,
stepwise=True)
print(model.summary())
# Forecast
n_periods = 24
fc, confint = model.predict(n_periods=n_periods, return_conf_int=True)
index_of_fc = np.arange(len(df.value), len(df.value)+n_periods)
# make series for plotting purpose
fc_series = pd.Series(fc, index=index_of_fc)
lower_series = pd.Series(confint[:, 0], index=index_of_fc)
upper_series = pd.Series(confint[:, 1], index=index_of_fc)
# Plot
plt.plot(df.value)
plt.plot(fc_series, color='darkgreen')
plt.fill_between(lower_series.index,
lower_series,
upper_series,
color='k', alpha=.15)
plt.title("Final Forecast of WWW Usage")
plt.show()
如何自动构建季节性ARIMA模型？
如果模型带有季节性，则除了p,d,q以外，模型还需要引入季节性部分：
与非季节性模型的区别在于，季节性模型都是以m为固定周期来做计算的，比如D就是季节性差分，是用当前值减去上一个季节周期的值，P和Q和非季节性的p,q的区别也是在于前者是以季节窗口为单位，而后者是连续时间的。
上节介绍的auto arima的代码中，seasonal参数设为了false，构建季节性模型的时候，把该参数置为True，然后对应的P，D，Q,m参数即可，代码如下：
# !pip3 install pyramid-arima
import pmdarima as pm
# Seasonal - fit stepwise auto-ARIMA
smodel = pm.auto_arima(data, start_p=1, start_q=1,
test='adf',
max_p=3, max_q=3, m=12,
start_P=0, seasonal=True,
d=None, D=1, trace=True,
error_action='ignore',
suppress_warnings=True,
stepwise=True)
smodel.summary()
注意这里的stepwise参数，默认值就是True，表示用stepwise algorithm来选择最佳的参数组合，会比计算所有的参数组合要快很多，而且几乎不会过拟合，当然也有可能忽略了最优的组合参数。所以如果你想让模型自动计算所有的参数组合，然后选择最优的，可以将stepwise设为False。
如何在预测中引入其它相关的变量？
在时间序列模型中，还可以引入其它相关的变量，这些变量称为exogenous variable(外生变量，或自变量)，比如对于季节性的预测，除了之前说的通过加入季节性参数组合以外，还可以通过ARIMA模型加外生变量来实现，那么这里要加的外生变量自然就是时间序列中的季节性序列了(通过时间序列分解得到)。需要注意的是，对于季节性来说，还是用季节性模型来拟合比较合适，这里用外生变量的方式只是为了方便演示外生变量的用法。因为对于引入了外生变量的时间序列模型来说，在预测未来的值的时候，也要对外生变量进行预测的，而用季节性做外生变量的方便演示之处在于，季节性每期都一样的，比如年季节性，所以直接复制到3年就可以作为未来3年的季节外生变量序列了。
%matplotlib inline
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import datetime
from dateutil.relativedelta import relativedelta
import seaborn as sns
import statsmodels.api as sm
from statsmodels.tsa.stattools import acf,pacf
from statsmodels.tsa.seasonal import seasonal_decompose
decomposition = seasonal_decompose(df['riders'], freq=12)
def load_data():
"""
航司乘客数时间序列数据集
该数据集包含了1949-1960年每个月国际航班的乘客总数。
"""
from datetime import datetime
date_parse = lambda x: datetime.strptime(x, '%Y-%m-%d')
data = pd.read_csv('https://www.analyticsvidhya.com/wp-content/uploads/2016/02/AirPassengers.csv', index_col='Month', parse_dates=['Month'], date_parser=date_parse)
# print(data)
# print(data.index)
ts = data['value']
# print(ts.head(10))
# plt.plot(ts)
# plt.show()
return ts,data
# 加载时间序列数据
_ts,_data = load_data()
# 时间序列分解
result_mul = seasonal_decompose(_ts[-36:], # 3 years
model='multiplicative',
freq=12,
extrapolate_trend='freq')
_seasonal_frame = result_mul.seasonal[-12:].to_frame()
_seasonal_frame['month'] = pd.to_datetime(_seasonal_frame.index).month
# seasonal_index = result_mul.seasonal[-12:].index
# seasonal_index['month'] = seasonal_index.month.values
print(_seasonal_frame)
_data['month'] = _data.index.month
print(_data)
_df = pd.merge(_data, _seasonal_frame, how='left', on='month')
_df.columns = ['value', 'month', 'seasonal_index']
print(_df)
print(_df.index)
_df.index = _data.index # reassign the index.
print(_df.index)
build_arima(_df,_seasonal_frame,_data)
# SARIMAX Model
sxmodel = pm.auto_arima(df[['value']],
exogenous=df[['seasonal_index']],
start_p=1, start_q=1,
test='adf',
max_p=3, max_q=3, m=12,
start_P=0, seasonal=False,
d=1, D=1, trace=True,
error_action='ignore',
suppress_warnings=True,
stepwise=True)
sxmodel.summary()
# Forecast
n_periods = 36
fitted, confint = sxmodel.predict(n_periods=n_periods,
exogenous=np.tile(seasonal_frame['y'].values, 3).reshape(-1, 1),
return_conf_int=True)
index_of_fc = pd.date_range(data.index[-1], periods = n_periods, freq='MS')
# make series for plotting purpose
fitted_series = pd.Series(fitted, index=index_of_fc)
lower_series = pd.Series(confint[:, 0], index=index_of_fc)
upper_series = pd.Series(confint[:, 1], index=index_of_fc)
# Plot
plt.plot(data['y'])
plt.plot(fitted_series, color='darkgreen')
plt.fill_between(lower_series.index,
lower_series,
upper_series,
color='k', alpha=.15)
plt.title("SARIMAX Forecast of a10 - Drug Sales")
plt.show()