statsmodels是统计建模分析的核心工具包,其包括了几乎所有常见的各种回归模型、非参数模型和估计、时间序列分析和建模以及空间面板模型等
statsmodels统计建模分析
时间序列预测7种方法
1.朴素法
如果数据集在一段时间内都很稳定,我们想预测第二天的价格,可以取前面一天的价格,预测第二天的值。这种假设第一个预测点和上一个观察点相等的预测方法就叫朴素法。即
y
^
t
+
1
=
y
t
\hat{y}_{t+1}=y_t
y^t+1=yt
2.简单平均法
物品价格会随机上涨和下跌,平均价格会保持一致。我们经常会遇到一些数据集,虽然在一定时期内出现小幅变动,但每个时间段的平均值确实保持不变。这种情况下,我们可以预测出第二天的价格大致和过去天数的价格平均值一致。这种将预期值等同于之前所有观测点的平均值的预测方法就叫简单平均法。
y
^
x
+
1
=
1
x
∑
i
x
y
i
\hat{y}_{x+1}=\frac{1}{x}\sum_{i}^{x}y_i
y^x+1=x1i∑xyi
3.移动平均法
物品价格在一段时间内大幅上涨,但后来又趋于平稳。,比如价格或销售额某段时间大幅上升或下降。如果我们这时用之前的简单平均法,就得使用所有先前数据的平均值,但在这里使用之前的所有数据是说不通的,因为用开始阶段的价格值会大幅影响接下来日期的预测值。因此,我们只取最近几个时期的价格平均值。很明显这里的逻辑是只有最近的值最要紧。这种用某些窗口期计算平均值的预测方法就叫移动平均法。
计算移动平均值涉及到一个有时被称为滑动窗口的大小值P。使用简单的移动平均模型,我们可以根据之前数值的固定有限数p的平均值预测某个时序中的下一个值。这样,对于所有的 i >p:
y
^
l
=
1
p
(
y
i
−
1
+
y
i
−
2
+
.
.
.
+
y
i
−
p
)
\hat{y}_l = \frac{1}{p}(y_{i-1}+y_{i-2}+...+y_{i-p})
y^l=p1(yi−1+yi−2+...+yi−p)
在上文移动平均法可以看到,我们对“p”中的观察值赋予了同样的权重。但是我们可能遇到一些情况,比如“p”中每个观察值会以不同的方式影响预测结果。将过去观察值赋予不同权重的方法就叫做加权移动平均法。加权移动平均法其实还是一种移动平均法,只是“滑动窗口期”内的值被赋予不同的权重,通常来讲,最近时间点的值发挥的作用更大了。即:
y
^
l
=
1
m
(
w
1
∗
y
i
−
1
+
w
2
∗
y
i
−
2
+
.
.
.
+
w
m
y
i
−
m
)
\hat{y}_l = \frac{1}{m}(w_1*y_{i-1}+w_2*y_{i-2}+...+w_my_{i-m})
y^l=m1(w1∗yi−1+w2∗yi−2+...+wmyi−m)
这种方法并非选择一个窗口期的值,而是需要一列权重值(相加后为1)。例如,如果我们选择[0.40, 0.25, 0.20, 0.15]作为权值,我们会为最近的4个时间点分别赋给40%,25%,20%和15%的权重。
4.简单指数法
我们注意到简单平均法和加权移动平均法在选取时间点的思路上存在较大的差异。我们就需要在这两种方法之间取一个折中的方法,在将所有数据考虑在内的同时也能给数据赋予不同非权重。例如,相比更早时期内的观测值,它会给近期的观测值赋予更大的权重。按照这种原则工作的方法就叫做简单指数平滑法。它通过加权平均值计算出预测值,其中权重随着观测值从早期到晚期的变化呈指数级下降,最小的权重和最早的观测值相关:
y
^
T
+
1
=
α
y
T
+
α
(
1
−
α
)
y
T
−
1
+
α
(
1
−
α
)
2
y
T
t
−
2
+
.
.
.
\hat{y}_{T+1} = αy_T+α(1-α)y_{T-1}+α(1-α)^2y_{Tt-2}+...
y^T+1=αyT+α(1−α)yT−1+α(1−α)2yTt−2+...
其中0≤α≤1是平滑参数。对时间点T+1的单步预测值是时序y1,…,yT的所有观测值的加权平均数。权重下降的速率由参数α控制,预测值yx是α⋅yt与(1−α)⋅yx的和。
因此,它可以写成:
y
^
T
+
1
=
α
y
T
+
(
1
−
α
)
y
T
−
1
\hat{y}_{T+1} = αy_T + (1-α)y_{T-1}
y^T+1=αyT+(1−α)yT−1
所以本质上,我们是用两个权重α和1−α得到一个加权移动平均值,让表达式呈递进形式。
from statsmodels.tsa.api import SimpleExpSmoothing
y_hat_avg = test.copy()
fit = SimpleExpSmoothing(np.asarray(train['Count'])).fit(smoothing_level=0.6,optimized=False)
y_hat_avg['SES'] = fit.forecast(len(test))
// 模型中使用的α值为0.6,我们可以用测试集继续调整参数以生成一个更好的模型。
5.Holt 线性趋势法
假设y轴表示某个物品的价格,x轴表示时间(天)。
如果物品的价格是不断上涨的(见上图),我们上面的方法并没有考虑这种趋势,即我们在一段时间内观察到的价格的总体模式。
每个时序数据集可以分解为相应的几个部分:趋势(Trend),季节性(Seasonal)和残差(Residual)。任何呈现某种趋势的数据集都可以用霍尔特线性趋势法用于预测。
import statsmodels.api as sm
sm.tsa.seasonal_decompose(train['Count']).plot()
result = sm.tsa.stattools.adfuller(train['Count'])
plt.show()
该数据集呈上升趋势。因此我们可以用霍尔特线性趋势法预测未来价格。该算法包含三个方程:一个水平方程,一个趋势方程,一个方程将二者相加以得到预测值y:
我们在上面算法中预测的值称为水平(level)。正如简单指数平滑一样,这里的水平方程显示它是观测值和样本内单步预测值的加权平均数,趋势方程显示它是根据 ℓ(t)−ℓ(t−1) 和之前的预测趋势 b(t−1) 在时间t处的预测趋势的加权平均值。
我们将这两个方程相加,得出一个预测函数。我们也可以将两者相乘而不是相加得到一个乘法预测方程。当趋势呈线性增加和下降时,我们用相加得到的方程;当趋势呈指数级增加或下降时,我们用相乘得到的方程。实践操作显示,用相乘得到的方程,预测结果会更稳定,但用相加得到的方程,更容易理解。
from statsmodels.tsa.api import Holt
y_hat_avg = test.copy()
fit = Holt(np.asarray(train['Count'])).fit(smoothing_level=0.3, smoothing_slope=0.1)
y_hat_avg['Holt_linear'] = fit.forecast(len(test))
// 若调整参数,结果会更好
6.Holt-Winters季节性预测模型
在应用这种算法前,我们先介绍一个新术语。假如有家酒店坐落在半山腰上,夏季的时候生意很好,顾客很多,但每年其余时间顾客很少。因此,每年夏季的收入会远高于其它季节,而且每年都是这样,那么这种重复现象叫做“季节性”(Seasonality)。如果数据集在一定时间段内的固定区间内呈现相似的模式,那么该数据集就具有季节性。
Holt-Winters季节性预测模型:模型在预测时考虑到数据集的季节性,它是一种三次指数平滑预测,其背后的理念就是除了水平和趋势外,还将指数平滑应用到季节分量上。
Holt-Winters季节性预测模型由预测函数和三次平滑函数——一个是水平函数ℓt,一个是趋势函数bt,一个是季节分量 st,以及平滑参数α,β和γ ∈(0,1)
# 水平函数为季节性调整的观测值和时间点t处非季节预测之间的加权平均值。 # 趋势函数和霍尔特线性方法中的含义相同。
# 季节函数为当前季节指数和去年同一季节的季节性指数之间的加权平均值。
在本算法,我们同样可以用相加和相乘的方法。当季节性变化大致相同时,优先选择相加方法,而当季节变化的幅度与各时间段的水平成正比时,优先选择相乘的方法。
from statsmodels.tsa.api import ExponentialSmoothing
y_hat_avg = test.copy()
fit1 = ExponentialSmoothing(np.asarray(train['Count']), seasonal_periods=7, trend='add', seasonal='add', ).fit()
y_hat_avg['Holt_Winter'] = fit1.forecast(len(test))
# 搭建这个模型的时候用的是默认参数。你可以试着调整参数来优化模型。
7.ARIMA 自回归移动平均模型
指数平滑模型都是基于数据中的趋势和季节性的描述,而自回归移动平均模型的目标是描述数据中彼此之间的关系。ARIMA的一个优化版就是季节性ARIMA。它像Holt-Winters季节性预测模型一样,也把数据集的季节性考虑在内。
import statsmodels.api as sm
y_hat_avg = test.copy()
fit1 = sm.tsa.statespace.SARIMAX(train.Count, order=(2, 1, 4), seasonal_order=(0, 1, 1, 7)).fit()
y_hat_avg['SARIMA'] = fit1.predict(start="2013-11-1", end="2013-12-31", dynamic=True)
plt.figure(figsize=(16, 8))
plt.plot(train['Count'], label='Train')
plt.plot(test['Count'], label='Test')
plt.plot(y_hat_avg['SARIMA'], label='SARIMA')
plt.legend(loc='best')
plt.show()
ARIMA的开发流程为:
1.数据可视化,识别平稳性。
2.对非平稳的时间序列数据做差分,得到平稳序列。
3.建立合适的模型:
- 1.平稳化处理后,若PACF是截尾的,而ACF是拖尾的,则建立AR模型;
- 2.若PACF是拖尾的,而ACF是截尾的,则建立MA模型;
- 3.若PACF和ACF均是拖尾的,则序列适合ARMA模型;
4.模型的阶数在确定之后,对ARMA模型进行参数估计,比较常用是最小二乘法进行参数估计。
5.假设检验,判断(诊断)残差序列是否为白噪声序列。
6.利用已通过检验的模型进行预测。
项目代码基本步骤:
1.用pandas读取和加载数据,并同时做好数据预处理,处理空值以及时间点缺失。
2.使用可视化对数据进行探索,通过图表了解基本的统计和可能存在的季节性趋势。
3.选择模型,规定网格搜寻的参数范围,以最低AIC得分穷举最佳参数组合。
4.设置需要预测的时间段,然后用最优模型来做predict。