【时间序列】时间序列的平稳性

本专栏将会介绍时间序列、预测和不同的基本原理我们将在讨论的不同章节中使用。具体来说，本笔记本将讨论：

1. 时间序列
2. 时间序列预测
3. 随机过程
   我们使用Python及相关Python第三方库来分析时间序列

# Required Libraries
import numpy as np
import numpy.random as rng
import pandas as pd
import scipy.stats as st
import matplotlib.pyplot as plt
import statsmodels.graphics.tsaplots as tg
%matplotlib inline

plt.rcParams['figure.figsize'] = [15, 2]

General Introduction

我们大多数人都会听说过市场上的新嗡嗡声，即加密货币。我们中的许多人也会投资他们的硬币。但是，将资金投资于如此波动的货币安全吗？我们如何确保现在投资这些硬币肯定会在未来产生健康的利润？我们不能确定，但我们肯定可以根据之前的价格生成一个近似值。时序模型是预测它们的一种方法。

What is a Time Series?

引用两篇博客中对时间序列的定义：

1. 时间序列数据（time-series data) ，亦即单一变量按时间的先后次序产生的数据。
2. 时间序列是一系列索引值，其中每个值都是随机变量的结果。换句话说，时间序列是相应过程的一种实现。

时间序列的一个例子是在德国耶拿马克斯普朗克生物地球化学研究所的气象站记录的天气时间序列数据集。在这个数据集中，几年内每 14 分钟记录 10 个不同的量（如气温、大气压力、湿度、风向等）。原始数据可以追溯到 2003 年，但此示例仅限于 2009-2016 年的数据。此数据集非常适合学习使用数值时间序列。我们将使用此数据集来构建我们的模型。以下是数据中的实际变量：

• Date Time
• p (mbar) atmospheric pressure
• T (degC) temperature
• Tpot (K) potential temperature
• Tdew (degC) dew point temperature
• rh (%) relative humidity
• VPmax (mbar) saturation water vapor pressure
• VPact (mbar) actual water vapor pressure
• VPdef (mbar) water vapor pressure deficit
• sh (g/kg) specific humidity
• H2OC (mmol/mol) water vapor concentration
• rho $(g/m^3)$ air density
• wv (m/s) wind velocity
• max. wv (m/s) maximum wind velocity
• wd (deg) wind direction

Time Series Forecasting

预测是在给定所有可用信息（例如历史数据和可能影响数据的任何可能变量）的情况下预测时间序列的未来值。除了加密货币之外，还有多个重要领域使用时间序列预测，例如：预测销售，呼叫中心的呼叫量，太阳能活动，海洋潮汐，股票市场行为等等。假设一家酒店的经理想要预测他预计明年有多少访客来相应地调整酒店的库存，并对酒店的收入做出合理的猜测。根据前几年/月/日的数据，他可以使用时间序列预测并得到访问者的近似值。游客的预测价值将有助于酒店管理资源并相应地计划事情。
在我们预测时间序列之前需要考虑时间序列的几个点：

1. 数据类型（趋势、季节性等）。
2. 影响他的因素分析（多元分析、因果关系分析等）。
3. 数据量的大小。
   在接下来的章节中，我们将讨论如何使用它们的不同预测方法。每种方法都有自己的一组假设和局限性，但一般的想法是我们假设有一个基本的数学过程来定义我们的时间序列。了解这些数学过程有助于我们找到最佳方法或算法（以及相应的超参数），以更好地预测我们的时间序列。
   在下一节中我们讨论时间序列的第一种表述方式，随机过程。

stochastic process

随机过程是一系列索引随机变量，在时间序列分析的背景下，我们将使用、离散的时间作为索引：$t=0,1,2,\cdots$
{ X 0 , X 1 , X 2 , ⋯   } \{X_0,X_1,X_2,\cdots\} {X0​,X1​,X2​,⋯}表示一个随机过程，其中$X_0,X_1,X_2\cdots$表示随机变量。
在此上下文中，时间序列是一系列索引值，其中每个值都是随机变量的结果。换句话说，时间序列是相应过程的一种实现。
示例 1：伯努利试验。回想一下，伯努利试验有两种可能的结果或实现。让我们把它们当作.让我们生成几个值：{−1,+1}

bern_outcomes = [-1., 1.]
size = 100
# Generate 100 Bernoulli trials
flips = rng.choice(bern_outcomes, size=size)
series = pd.Series(flips)
fig, ax = plt.subplots(figsize=(15, 2))
series.plot(ax=ax) # here we implicitly index them; let the index stand for time
plt.show()

示例2：高斯白噪声。高斯白噪声过程生成一个时间序列，其中每个值都来自均值和单位方差为零的高斯分布：$x_t\sim Gaussian(0,1)$

size = 1000
wn = st.norm.rvs(size=size)
wn_series = pd.Series(wn)
fig, ax = plt.subplots(figsize=(15, 2))
wn_series.plot(ax=ax)
plt.show()

高斯白噪声具有独立且分布相同的特性，这导致了线性时间序列的一些关键特性，下面将介绍。
示例 3：一维布朗运动。一维布朗运动或维纳过程可以被认为是白噪声时间序列的累积和。序列在$t$时刻的值是所有先前值的总和，可以认为是随机步行者在时间$t$的位移：$X_t=X_{t1}+X_{t-2}+\cdots ={\sum }_{i=0}^t{X}_i$

size = 100000
bms = []
bms_msq = []
for i in range(100):
    wn = st.norm.rvs(size=size)
    bm = np.cumsum(wn)
    bms.append(bm)
    bms_msq.append([x**2 for x in bm])
bm_ave = np.mean(bms, axis=1)
bm_msq = np.mean(bms_msq, axis=1)

fig, ax = plt.subplots(figsize=(15, 2))
bms_series = pd.Series(bm_ave)
bms_msq_series = pd.Series(bm_msq)
bms_series.plot(label="Displacement", ax=ax)
bms_msq_series.plot(label="MSD", ax=ax)
plt.legend()
plt.show()

示例4：移动平均`给定一个白噪声过程见上，在$t$时刻的一个q阶的移动平均过程是当前值和过去几个值的加权和：
$X_t={ϵ}_t+a_1{ϵ}_{t-1}+\cdots +a_q{ϵ}_{t-q}$,其中 { ϵ t } ∼ W N ( 0 , σ 2 ) \{\epsilon_t\}\sim WN(0,\sigma^2) {ϵt​}∼WN(0,σ2)。我们写作 { X t } ∼ M A ( q ) \{X_t\}\sim MA(q) {Xt​}∼MA(q)，该移动平均过程是一个线性过程，假如我们有一个参数$a_1=0.5$和$a_2=0.25$的二阶移动平均过程：
$X_t={ϵ}_t+a_1{ϵ}_{t-1}+a_2{ϵ}_{t-2}$

def ma2(wn):
    x = []
    # First element is just epsilon_t
    x.append(wn[0])
    # For the second element, let's take just the first coefficient (admittedly a judgment call)
    x.append(0.5*wn[0] + wn[1])
    # For the rest, use the full equation
    for t in range(2, len(wn)):
        x.append(wn[t] + 0.5*wn[t - 1] + 0.25*wn[t - 2])
    return np.array(x)

size = 1000
wn = st.norm.rvs(size=size)
ma2_series = pd.Series(
    ma2(wn)
)
wn_series = pd.Series(wn)
fig, ax = plt.subplots(figsize=(15, 2))
wn_series.plot(label="WN(0,1)", ax=ax)
ma2_series.plot(label=r"MA(2)", ax=ax)
plt.legend(loc="best")
plt.show()

Stationarity

时序数据具有平稳性，意味着其统计特性（如均值、方差）在时间上是常数或不会随时间变化。
给定一个时间序列 { X t } \{X_t\} {Xt​}，在时间$t_1+\tau ,t_2+\tau ,\cdots ,t_n+\tau$具有相应的时间分布$P(X_{t_{1+\tau }},X_{t_{2+\tau }},\cdots ,X_{t_{n+\tau }})$。
对任何$\tau ,t_1\cdots ,t_n\text{和}n$来说，如果$P(X_{t_{1+\tau }},X_{t_{2+\tau }},\cdots ,X_{t_{n+\tau }})=P(X_{t1}+,X_{t2},\cdots ,X_{tn})$，那么时间序列$X_t$是平稳的。
本质上讲，平稳过程的统计量对时间平移是不变的。
如果时间序列 { X t } \{X_t\} {Xt​}满足：

1. 对任意$\tau$来说$\mu (X_t)=\mu (X_{t+\tau })=\mu$。
2. 对任意$t$来说$E[|X_t{|}^2]<\infty$。
3. $\Gamma (t_1,t_2)=\Gamma (t_1-t_2)$

上式中$\Gamma (t_1,t_2)$表示$t_1$和$t_2$的协方差定义为：$\Gamma (t_1,t_2)=cov(t_1,t_2)=E[(X_{t1}-\mu X_{t1})(X_{t2}-\mu X_{t2})]$。
平稳相关性意味着延迟参数$\tau$，可用于参数化自协方差，从而：
$\Gamma (\tau )=E[(X_t-\mu )(X_{t+\tau }-\mu )]={\gamma }_{t-s}$。
{ γ t } \{\gamma_t\} {γt​}为序列 { X t } \{X_t\} {Xt​}的自协方差函数
通过标准差归一化$t_1$和$t_2$时刻的自协方差，并产生自相关：
$$P(t_1,t_2)=\frac{cov(t_1,t_2)}{{\sigma }_{t_1}{\sigma }_{t_2}}=\frac{E[(X_{t_1}-\mu (X_{t_1}))(X_{t_2}-\mu (X_{t_2}))]}{{\sigma }_{t_1}{\sigma }_{t_2}}$$
对于弱平稳性时间序列来说：
$$P(\tau )=\frac{E[(X_t-\mu )(X_{t+\tau }-\mu )]}{{\sigma }^2}$$

Autocovariance of a stationary process

平稳时间序列的自协方差函数具有特定的形式，让$\tau =|t_1-t_2|$表示时间滞后lag，对于具有$n$个时间点的平稳时间序列 { X t } \{X_t\} {Xt​}来说，它的自协方差函数矩阵是对称和正定的：

1. 对称性 ${\gamma }_k={\gamma }_{-k},\forall k\in \mathbb{Z}$。
2. 正定性：
   $${\Gamma }_n=({\gamma }_{k-j}{)}_{k,j=1}^n=\left(\begin{array}{ccccc}{\gamma }_0& {\gamma }_1& {\gamma }_2& \cdots & {\gamma }_{n-1}\\ {\gamma }_1& {\gamma }_0& {\gamma }_1& \cdots & {\gamma }_{n-2}\\ {\gamma }_2& {\gamma }_1& {\gamma }_0& \cdots & {\gamma }_{n-3}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {\gamma }_{n-1}& {\gamma }_{n-2}& \cdots & {\gamma }_{n-3}& {\gamma }_0\end{array}\right)$$

Moving Average Process

移动平均模型（MA模型）描述的是当前时间点的数据与过去噪声的关系。严格定义上来讲：其模型的定义是基于白噪声序列的假设。白噪声是一种特殊的时间序列模型，每个时间点的数据都是独立且服从相同分布的，且具有常数的均值和方差。$q$阶的移动平均模型的函数表达式为：
$$Y_t={ϵ}_t+a_1{ϵ}_{t-1}+\cdots +a_q{ϵ}_{t-q}$$
这一公式被称为q阶的移动平均模型，写作$MA(q)$，其中：

1. $Y_t$是我们感兴趣的时间序列，在时间点$t$的观察值。
2. $\mu$是时间序列的均值或者期望值，这个值对所有的时间点都是相同的。
3. ${ϵ}_t,{ϵ}_{t-1},{ϵ}_{t-2},\dots ,{ϵ}_{t-q}$这些是所谓的白噪声项，每个时间点的值都是独立同分布的，通常假设为正态分布。这些项的均值为0，方差为${\sigma }^2$。
4. ${\theta }_1,{\theta }_2,\dots ,{\theta }_q$这些事模型的参数，每个参数的$\theta$对应一个白噪声项，它们衡量的是对应的白噪声对当前时间点的影响程度。
5. $q$是阶数，表示有多少个过去的白噪声项被纳入模型。

这个公式的含义是：在时间点$t$，观察到的值$Y_t$是由过去q个时间点的白噪声的线性组合加上一个常数$\mu$和当前时间点的白噪声决定的。各个白噪声项的权重由参数$\theta$决定。

Autoregressive Process

AR自回归模型，是一种被广泛用于分析时间序列数据的统计模型。从严格的定义上来说：在统计学中，一个自回归模型的形式是描述某个变量与其过去值的关系的。对于一个时间序列数据，自回归模型的$q$阶形式为：
$$Y_t={ϵ}_t+a_1{X}_{t-1}+\cdots +a_p{X}_{t-p}$$
其中$p$表示使用的时间滞后（lags）的数量。公式中：

1. ${ϵ}_t$表示在时间点$t$的误差项，也被称为白噪声项。
2. $a_1{X}_{t-1}+\cdots +a_p{X}_{t-p}$：这些项表示过去的$p$个时间点的时间序列值对当前时间点的影响，$a_i$是弟$i$哥滞后项的系数，代表了弟$i$个滞后项对当前时间点的相对影响力。

Evaluation Metrics for Forecast Accuracy

预测是时间序列分析中最常见的推理任务之一。为了正确衡量时序模型的性能，通常的做法是将数据集分为两部分：训练数据和测试数据。使用训练数据估计模型参数，然后使用模型生成根据测试数据进行评估的预测。
错误统计有不同的风格，每种风格都有自己的优点和缺点。

Mean Absolute Error

最常用的指标是mean absolute error或者MAE。它的主要优点是易于理解和计算。我们注意到该指标与规模相关，这意味着它无法与具有不同单位的其他时间序列进行比较。此外，最小化 MAE 的方法会生成导致中位数的预测。
$$MAE=\frac{1}{n}\sum _{t=1}^n|y_t-{\hat{y}}_t|$$

Root Mean Squared Error

另一个经常用的是指标是root mean squared error或者RMSE，虽然更难解释。最小化 RMSE 的方法生成导致平均值的预测，并且与 MAE 类似，RMSE 与规模相关。
$$RMSE=\sqrt{\frac{1}{n}\sum _{t=1}^n(y_t-{\hat{y}}_t{)}^2}$$

Mean Absolute Percentage Error

为了在具有不同单位的时间序列之间进行比较，可以使用百分比误差。一个常用的基于百分比的指标是Mean Absolute Percentage Error或 MAPE。不幸的是，百分比误差的主要缺点是，当$y_t$为零或接近零时，它们可能是无限的或未定义的。
$$MAPE=\frac{1}{n}\sum _{t=1}^n|\frac{y_t-{\hat{y}}_t}{y_t}|$$

Symmetric Mean Absolute Percentage Error

MAPE 的一个特别缺点是，与正错误相比，它对负误差的惩罚更大。这导致了另一种基于百分比的指标，称为对称平均绝对百分比误差或 SMAPE。然而，计算该指标可能依然是不稳定的如果$y_t$和${\hat{y}}_t$两者都接近于零的话。
$$SMAPE=\frac{1}{n}\sum _{t=1}^n\frac{|{\hat{y}}_t-y_t|}{|y_t|+|{\hat{y}}_t|}$$