[深度学习项目] - 时间序列预测 (3)

本文链接：https://blog.csdn.net/upr_rom/article/details/127053234

文章目录

移动平均法
ETS
时间序列的平稳性定义
Arima算法
Homewrok
- ETS 预测
- ARIMA

移动平均法

Moving Average smoothing(MA)
原理：用过去N个时刻的观测值的平均值作为对未来时刻的预测
$y_{t+h}=\frac{1}{N}\sum_{i=0}^{N-1} y_{t-i}$

使用过去N个时刻的均值来直接作为预测结果，不考虑外部变量的影响。可以用来判断预测的难易程度。

import pandas as pd

store_sales = pd.read_csv("store_sales.csv", parse_dates=["week"])
stock_price = pd.read_csv("stock_price.csv", parse_dates=["date"])

print(stock_price.head(3))

T = 10 # 10 days for stock price

date_end = pd.to_datetime("2022-01-01")
date_begin = date_end - pd.Timedelta(T-1, unit="D")  # 根据日期,往前推算 pd.Timedelta

FB_price = stock_price[(stock_price["company"]=="FB")&
                       (stock_price["date"]>=date_begin)&
                       (stock_price["date"]<=date_end)]["price"].mean()
print(FB_price)

当我们想要预测未来h天的结果时，我们会将该过去N个时刻的均值重复h次作为未来的预测结果。

ETS

简单ETS算法 (SETS)

Exponential Smoothing，指数平滑算法，简称ETS。

简单平滑预测可以看作 朴素预测法 和 平均预测法 的一种中和。

朴素预测法： $\hat{y}_{T+h|T}=y_T$ ，我们认为对未来的所有预测都等于该序列的最后一个观测值。
平均预测法： $\hat{y}_{T+h|T}=\frac{1}{T}\sum_{t=1}^T y_t$ ，我们认为对未来的所有预测与之前所有的观测值都有关，而且关联程度相等。【移动平均和平均预测类似，只不过是选用了一定的时间区间。】

而ETS的思想就是认为比较近的数据对于未来提供的信息更大一些，但是过去的数据也对于未来的数据提供一定的帮助。
因此，ETS：
$\hat{y}_{T+h|T}=\alpha y_T+\alpha(1-\alpha)y_{T-1}+\alpha(1-\alpha)^2 y_{T-2}+\alpha (1-\alpha)^3 y_{T-3}+\dots$
其中， $0\leq \alpha \leq 1$ 是平滑参数。

也可以写作:
Forecast equation: $\hat{y}_{T+h|T}=\ell_t$
Smoothing equation: $\ell_t=\alpha y_t +(1-\alpha) \ell_{t-1}$

带趋势的ETS算法

(霍尔特的线性趋势法)
我们将预测函数分别成两部分 水平项 和 趋势项。
Forecast equation: $\hat{y}_{t+h|t}=\ell_t+hb_t$
其中 $\ell_t$ 表示水平项，用于描述该时间序列的水平的估计值， $b_t$ 用于描述该时间序列在t时刻的趋势（斜率）的估计， $\alpha$ 是水平 $\leq \alpha \leq 1$ 的平滑参数， $\beta^*$ 表示趋势 $0\leq \beta^* \leq1$ 的平滑参数。
我们相当于将时间序列看作当前水平，和上升(下降)趋势【取决于 $b_t$ 的正负】。

水平和趋势的更新也源于数据的迭代：
Level equation: $\ell_t=\alpha y_t +(1-\alpha) (\ell_{t-1}+b_{t-1})$

水平项实际上取决于时刻t的真值和 t-1时刻在t时刻的预测值。

Trend equation: $b_t = \beta^* (\ell_t -\ell_{t-1}) + (1-\beta^*)b_{t-1}$
趋势项实际上取决于当前预测值与上一时刻预测值的差值（当前的上升趋势）和之前的趋势。

import pandas as pd
from statsmodels.tsa.exponential_smoothing.ets import ETSModel
import matplotlib.pyplot as plt

stock_price = pd.read_csv("stock_price.csv", parse_dates=["date"])


date_end = pd.to_datetime("2022-01-01")
data_train = stock_price[ (stock_price["company"]=="FB")&
                          (stock_price["date"]<=date_end)]

# 将日期变成 训练集的index
data_train = data_train[ ["date","price"] ].set_index("date")

# 构建模型
model = ETSModel( data_train["price"] )
fit = model.fit()

# 获取数据
data_train["fit_price"] = fit.fittedvalues
data_plot = data_train[ (data_train.index>="2021-09-18")  ]

plt.plot( data_plot["price"], ".-", label="actual" )
plt.plot( data_plot["fit_price"], ",-", label="fit" )
plt.legend()
plt.xticks(rotation=90)
plt.show()

# 预测未来h天
h = 10
pred_price = fit.get_prediction(start=data_train.shape[0]+1,
                                end=data_train.shape[0]+h)

print(pred_price.summary_frame()[ ["mean"] ])

>>>
>            mean
2423  336.763047
2424  336.763047
2425  336.763047
2426  336.763047
2427  336.763047
2428  336.763047
2429  336.763047
2430  336.763047
2431  336.763047
2432  336.763047

在这里插入图片描述

带季节性的ETS算法

（预测效果不佳）

时间序列的平稳性定义

平稳性定义( stationary )：如果截取某个时间序列中任意一段数值(起码大于5/10个时间点)，它的均值，方差不随时间发生变化，那么这个时间序列是平稳的。

平稳时间序列的好处是可以用移动平均的方式做预测。

可以采用自相关系数来估计序列的平稳性。如果自相关系数没有迅速衰减到0，那么序列是不平稳的。

使用差分的方式，获取平稳的时间序列

$y_t^1=y_t-y_{t-1}$

########### difference
FB_price = stock_price[ stock_price["company"]=="FB" ].sort_values("date").reset_index()
FB_price["price_last"] = FB_price["price"].shift(1)
FB_price["price_diff"] = FB_price["price"] - FB_price["price_last"]

plt.plot( FB_price["date"], FB_price["price_diff"] )
plt.show()

plot_acf(FB_price["price_diff"])
plt.show()

FB difference
通过差分将一直上升的 FB股价 转换成了 均值近似恒定 (虽然方差变化仍不小) 的时间序列。

FB_difference_ACF
也可以看到 ACF 自相关系数非常小。

随机散步模型(random walk model)

如果时序的差分结果是稳定的，那么结果就是带有固定噪音，并且均值恒定的白噪声，我们将白噪声可以记作 $c+\epsilon_t$ 【也叫残差，表示随机变量影响所造成的偏差。】其中c表示的是噪声的均值， $\epsilon_t$ 表示方差。
如果一阶差分是稳定的，那么：
$y^1_t=y_t-y_{t-1}=c+\epsilon_t$
那么，我们就可以得到：
$y_t=y_{t-1}+c+\epsilon_t$

也就是当前时刻的值 只取决于过去一个时刻的值，和一个上升/下降的趋势 ( c ).

季节差分

当前观测值和过去相同周期的观测值之间的差值是季节差分：
$y^1_t=y_t-y_{t-m}$
如果季节差分数据是白噪声，那么我们就可以将模型看作
$y_t=y_{t-m}+\epsilon_t+c$

有时，数据整体上有向上或者向下的趋势，同时也有季节性变化，我们可以既进行季节性差分，也进行一阶差分。如果有明显的季节性因素 一般我们先进行季节性差分。如果差分后结果稳定，那么就不再进行差分，否则再进行差分。【另外，我们在差分之前可以对数据进行处理，Log指数化处理等，让数据的变化趋势呈现线性】

多阶差分

为了获取平稳的时间序列，我们通常会用差分的方式
一阶差分：
$y^1_t=y_t-y_{t-1}=y_t-By_t=(1-B)y_t$
其中B表示 “backshift”，表示回退
$By_t$ 表示由 $y_t$ 回退一个时刻，所以表示 $y_{t-1}$
二阶差分：
$y^2_t=y^1_t-y^1_{t-1}=y_t-y_{t-1}-y_{t-1}+y_{t-2}=y_t-2By_t+B^2y_t=(1-B)^2 y_t$
d阶差分：
$y^d_t=(1-B)^d y_t$

季节差分和一阶差分同时作用：
$1-B)(1-B^m)y_t=(1-B-B^m+B^{m+1})y_t=y_t-y_{t-1}-y_{t-m}+y_{t-m-1}$

Arima算法

Autoregressive Integrated Moving Average Model [差分整合移动平均自回归模型]

Autoregression models

我们可以使用自回归模型来对模型进行预测， $A R (p)$ 表示 p阶的自回归模型。

$y_t =c+\phi_1 y_{t-1}+\phi_2 y_{t-2} +\dots +\phi_p y_{t-p} +\epsilon_t$

p阶数表示了 $y_t$ 跟多远的自身变量有关。

在AR(1)中，如果 $\phi_1=0$ ，那么模型就是白噪声。 $\phi_1=1,c=0$ ，那么 $y_t$ 就是随机漫步算法。 $\phi_1<0$ ，那么 $y_t$ 就在均值c上震荡。
AR(1)中，平稳性的条件: $-1<\phi_1<1$
AR(2)中，平稳性的条件: $-1<\phi_2<1,\phi_1+\phi_2<1,\phi_2-\phi_1<1$

Moving Average models

该Moving average models和之前的 MA算法完全不同。
该moving average model使用过去预测的误差来进行回归预测。MA(q) 表示q阶的模型。
$y_t=c+\epsilon_t+\theta_1 \epsilon_{t-1}+ \theta_2 \epsilon_{t-2} + \dots + \theta_q \epsilon_{t-q}$
每一个 $y_t$ 都可以看作过去预测误差的移动平均，

我们可以将任意的 AR(q) 看作 $MA(\infty)$ 模型：
$y_t=\phi_1 y_{t-1}+\epsilon_t =\phi_1(\phi_1 y_{t-2}+\epsilon_{t-1})+\epsilon_t=\phi_1^2 y_{t-2}+\phi_1 \epsilon_{t-1}+\epsilon_t=\phi_1^3 y_{t-3}+\phi_1^2 \epsilon_{t-2}+\phi_1 \epsilon_{t-1}+\epsilon_t=\epsilon_t +\phi_1 \epsilon_{t-1}+\phi_1^2 \epsilon_{t-2}+\phi_1^3 \epsilon_{t-3}+\dots$

Non-seasonal ARIMA models

$y'_t=c+\phi_1 y'_{t-1}+ \phi_2 y'_{t-2}+\dots + \phi_p y'_{t-p}+\theta_1 \epsilon_{t-1}+\dots+\theta_{q} \epsilon_{t-q}+\epsilon_t$
$y'_t$ 已经差分了d次，我们将上述ARIMA 记作 ARIMA(p,d,q)。

ARIMA 结合了差分，自回归和 moving average 三部分。

不同的模型	对应的ARIMA参数
White noise	ARIMA(0,0,0) $y_t=\epsilon_t$
Random Walk	ARIMA(0,1,0) with no constant
Random Walk with drift	ARIMA(0,1,0) with a constant
Autoregression	ARIMA(p,0,0)
Moving average	ARIMA(0,0,q)

我们使用 backshift 对上式进行修改
$(1-\phi_1 B-\dots-\phi_p B^p) (1-B)^dy_t =c +(1+\theta_1 B+\dots+\theta_q B^q)\epsilon_t$
AR ( p ) · d differences = c + MA(q)

meaning of different parameters

d=0 说明了 $y_t$ 本身就是平稳的，会围绕c上下波动。
d=1 说明了 $y_t'$ 是稳定的，当c=0时，说明 $y_t'$ 会趋于0，说明 $y_t$ 会趋于稳定值；如果 $c\neq0$ ， $y_t'$ 趋于固定值， $y_t$ 会呈现直线。
d=2 说明了 $y_t''$ 是稳定的，当c=0时，说明 $y_t''$ 会趋于0，然后 $y_t'$ 会趋于稳定值， $y_t$ 会呈现直线；如果 $c\neq 0$ ，那么 $y_t$ 会呈现二次曲线。

determine p,q

可以使用 ACF 和 PACF 来确定p和q。
PACF 可以在判断 $y_t$ 与 $y_{t-k}$ 之间相关性时排除掉因为 $y_{t-1},y_{t-2},\dots,y_{t-k+1}$ 造成的影响。
有时候我们可以使用 auto arima来获取参数，而不是自己直接定义 p,d,q。但是auto arima并不能探索到整个范围内，结果可能不是最优解。

Homewrok

分别使用 ETS算法，ARIMA算法预测所有公司的2022-01-01及以后的股价，每次往后预测1天，并且计算预测MAPE。

ETS 预测

[‘FB’ ‘TWTR’ ‘ETSY’ ‘SNAP’ ‘PINS’]

import pandas as pd
from statsmodels.tsa.exponential_smoothing.ets import ETSModel
from sklearn.metrics import mean_absolute_percentage_error as mape
#####################  ETS
stock_price = pd.read_csv("stock_price.csv", parse_dates=["date"])
COMPANY_NAME = "FB"

date_end = pd.to_datetime("2022-01-01")
pre_data_train = stock_price[ stock_price["company"]==COMPANY_NAME].sort_values("date")
all_predict_date = pre_data_train[ pre_data_train["date"]>=date_end  ]["date"]
# print(predict_date)
predict_result = []
true_result = []
mape_errors = []


for i,predict_date in all_predict_date.items():
    data_train = pre_data_train[pre_data_train["date"]<predict_date][ ["date","price"] ].set_index("date")
    model = ETSModel( data_train["price"] )
    fit = model.fit()
    price = fit.get_prediction(start=data_train.shape[0]+1,end=data_train.shape[0]+1)
    predict_result.append(price.predicted_mean)
    true_result.append(pre_data_train[ pre_data_train["date"]==predict_date ]["price"])
    mape_errors.append(  mape(true_result[-1],predict_result[-1]) )

plt.subplot(2,1,1)
plt.plot(all_predict_date,predict_result,".-",label="fit")
plt.plot(all_predict_date,true_result,".-",label="actual")
plt.xticks(rotation=90)
plt.ylabel(COMPANY_NAME+" stock price predict")
plt.xlabel("date")
plt.legend()

plt.subplot(2,1,2)
plt.plot(all_predict_date,mape_errors,".-")
plt.ylabel("MAPE")
plt.xlabel("date")
plt.xticks(rotation=90)

plt.show()

在这里插入图片描述

ARIMA

from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

########### difference
COMPANY = "TWTR"
date_end = pd.to_datetime("2022-01-01")
data_price = stock_price[ (stock_price["company"]==COMPANY)&
                          (stock_price["date"]<date_end)  ].sort_values("date").reset_index()
data_price["price_last"] = data_price["price"].shift(1)
data_price["price_diff"] = data_price["price"] - data_price["price_last"]

######## first difference
plt.figure(1)
plt.ylabel("first difference of "+COMPANY)
plt.plot( data_price["date"], data_price["price_diff"] )

###### acf & pacf
# 确定 d, q
plot_acf(data_price["price_diff"])


# 决定 p
plot_pacf(data_price["price"])
plt.show()

first diff
acf
pacf

from statsmodels.tsa.arima.model import ARIMA
predict_result = []
true_result = []
mape_errors = []

for i,predict_date in all_predict_date.items():
    data_train = pre_data_train[pre_data_train["date"]<predict_date]["price"].values
    model = ARIMA(data_train, order=(1, 1, 0)) # p, d, q
    fit = model.fit()

    price = fit.forecast(1)[0]
    predict_result.append(price)
    true_result.append(pre_data_train[ pre_data_train["date"]==predict_date ]["price"])
    mape_errors.append(  mape([true_result[-1]],[predict_result[-1]]) )


plt.subplot(2,1,1)
plt.plot(all_predict_date,predict_result,".-",label="fit")
plt.plot(all_predict_date,true_result,".-",label="actual")
plt.xticks(rotation=90)
plt.ylabel(COMPANY+" stock price predict")
plt.xlabel("date")
plt.legend()

plt.subplot(2,1,2)
plt.plot(all_predict_date,mape_errors,".-")
plt.ylabel("MAPE")
plt.xlabel("date")
plt.xticks(rotation=90)

plt.show()

在这里插入图片描述