Python学习——时间序列分析

发表于 | 分类于

面向乱搞的程序设计？

序列的分析

序列的平稳性

如果时间序列 \(x_t\) 在某一常数附近波动并且范围不大，方差和均值为常数，且延迟 \(k\) 期的子自协方差和自相关系数是相等的，则称 \(x_t\) 是平稳序列. 1

自协方差和自相关系数的定义如下：

  • 自协方差：\(\gamma(t,s) = \mathbb{E}[(x_t-\mu_t)(x_s-\mu_s)]\)
  • 自相关系数：\(\rho(t,s) = \frac{cov(x_t, x_s)}{\sigma_t\sigma_s}\)

平稳序列的判定：平稳性检验

一般有以下两种方法：

  1. 时序图检验：根据平稳序列的均值和方差都为常数的性质，平稳时间序列的时序图应当显示该序列始终在一个常数附近随机波动，且波动范围有限；
  2. 自相关检验：平稳时间序列具有短期相关性，故平稳时间序列通常只有近期的项才对当前值的影响较为明显。随着间隔 \(k\) 的增大，平稳序列的自相关系数会不断衰减至 \(0\)，而非平稳序列的衰减较慢。

白噪声序列

  • 白噪声序列指的是纯随机序列，序列的各项之间毫无关系，从中无法获取任何有用的信息。

白噪声序列的判定：纯随机性检验

纯随机序列满足 \(\gamma(k) = 0\)，其样本自相关系数很接近零，并且在零附近随机扰动。一般构造 \(Q\) 统计量和 \(LB\) 统计量来进行检验.

  • 在 Python 中可以使用 acorr_ljungbox 方法来检验.
  • 直接 print(acorr_ljungbox(ts, lags = 1)).

常用的模型原理

时间序列趋势主要受到总体趋势 \(T\)、季节性因素 \(S\)、周期性波动 \(C\) 和随机扰动 \(\varepsilon\) 的影响，一般来说有加法和乘法两种模型。

  1. 加法模型(additive model): \(x_t = T_t+S_t+C_t+\varepsilon_t\);
  2. 乘法模型(multiplicative model): \(x_t = T_t\times S_t\times C_t\times\varepsilon_t\).

在作分析时，常常采用的是AR-MA模型：

  1. AR模型 \(AR(p)\)：每一项都与前 \(p\) 项线性相关, 与随机扰动项无关. \[x_t=\varphi_0 + \varphi_1 x_{t-1}+\varphi_2 x_{t-2}+\cdots+\varphi_p x_{t-p}+\varepsilon_t = \varphi_0+ \sum\limits_{k=1}^p \varphi_k x_{t-k}+\varepsilon_t\] 其中随机扰动项 \(\varepsilon_t\) 为均值为 \(0\) 的白噪声序列. 均值方差均为常数，ACF拖尾，PACF \(q\) 阶截尾.
  2. MA模型 \(MA(q)\)：每一项都与前 \(q\) 项的随机扰动线性相关, 与真实值无关. \[x_t = \mu + \varepsilon_t-\theta_1 \varepsilon_{t-1}-\theta_2\varepsilon_{t-2}-\cdots-\theta_q \varepsilon_{t-q} = \mu + \varepsilon_t - \sum\limits_{k=1}^q \theta_k\varepsilon_{t-k}\] 其中随机扰动项 \(\varepsilon_t\) 为均值为 \(0\) 的白噪声序列，\(\mu\) 为序列均值. 均值方差为常数，ACF为 \(q\) 阶截尾，PACF拖尾.
  3. ARMA模型 \(ARMA(p,q)\)：每一项都与前 \(p\) 项线性相关, 同时与前 \(q\) 项的随机扰动有关. \[ \begin{align} x_t &= \varphi_0 + \varphi_1 x_{t-1}+\varphi_2 x_{t-2}+\cdots+\varphi_p x_{t-p} + \varepsilon_t-\theta_1 \varepsilon_{t-1}-\theta_2\varepsilon_{t-2}-\cdots-\theta_q \varepsilon_{t-q} \\ &= \varphi_0 + \sum\limits_{i=1}^p \varphi_i x_{t-i} + \varepsilon_t - \sum\limits_{j=1}^q \theta_j\varepsilon_{t-j} \end{align} \] 其中随机扰动项 \(\varepsilon_t\) 为均值为 \(0\) 的白噪声序列，\(\mu\) 为序列均值. 均值方差为常数，ACF为 \(q\) 阶截尾，PACF \(p\) 阶截尾.

容易看出AR模型和MA模型分别是ARMA在 \(p=0\) 和 \(q=0\) 时的特例.

一般在分析时分析的多是平稳序列，可以考虑用ARMA模型。有一些序列本身不是平稳序列，但做一阶差分之后成为平稳序列，这种序列被称为差分平稳序列，可以使用ARIMA模型进行拟合.

Python中的AMRA模型

导入库

1
2
3
4
5
6
7
8
9
10
11
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime
from statsmodels.graphics.tsaplots import plot_acf # 自相关图
from statsmodels.graphics.tsaplots import plot_pacf # 偏自相关图
from statsmodels.stats.diagnostic import acorr_ljungbox # 白噪声检验
from statsmodels.tsa.stattools import adfuller # 平稳性检测
from statsmodels.tsa.seasonal import seasonal_decompose # 季节性分解
from statsmodels.tsa.arima_model import ARMA
from statsmodels.tsa.arima_model import ARIMA

乱码问题处理

1
2
plt.rcParams['font.sans-serif'] = ['SimHei'] # 显示中文标签
plt.rcParams['axes.unicode_minus'] = False # 显示负号

数据读取与处理

按照我的习惯会将日期一列命名为 date.

数据读取

一般数据的格式为 XLS 或者 CSV.

  1. 数据为CSV格式 

    1
    df = pd.read_csv('文件地址', encoding = 'utf-8', index_col = 'date')

  2. 数据为XLS格式 

    1
    df = pd.read_excel('文件地址', encoding = 'utf-8', index_col = 'data')

若数据为数据库形式(如 MySQL), 则可以直接用 SQL 语句从数据库中提取.

数据格式的处理

以上读取数据的函数只会将读入的数据保存为 DataFrame 对象，但我们做预测所需要的数据类型是 Series，故需要将其转换为该格式，并将 date 指定为索引.2

1
2
df.index = pd.to_datetime(df.index)  # 指定 date 为索引
ts = df['x']  # 生成Series对象
注意此处数据文件中的日期必须为日期格式. 如 (yyyy/mm/dd) 是合法的，但 “20080808” 就是不合法的，必须进行处理.3

如此得到的 ts 即为一个可用的时间序列，不仅可以用类似于 2014-1-1 的字符串进行访问，也可以用时间对象 datetime(2014, 1, 1) 进行访问. 即以下两种访问方式都是可行的：

1
2
3
ts['2014-1-1']
ts[datetime(2014, 1, 1)]
ts['2014-1-1' : '2014-5-1'] # 切片操作. 注意此处的切片为闭区间, 端点都包括.

关于 Series 的更多内容参见 https://blog.csdn.net/zutsoft/article/details/51482573

数据的趋势分析

平稳性检验

可以进行 ADF 检验. 弱鸡还不理解其原理…暂且放在这里. 直接 print(adfuller(ts)) ，输出有一大串. 第一个值比后面带百分号的都小说明比较平稳.

平稳性处理：季节性分解

为了使序列更平稳，在做计算时经常对其取对数: 

1
ts_log = np.log(ts)

若原始数据有很强的周期性，则可以考虑对其进行季节性的分解，即将元素数据分离为总体趋势，季节性趋势以及残差三部分.

1
2
3
4
decomposition = seasonal_decompose(ts_log, model="additive")
trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid

将三种因素分离后，即可分别做预测，再相加. statsmodels 中还提供了乘法模型，只需将上述代码中的 additive 改为 multiplicative 即可.

模型识别

只需要对所要拟合的数据集绘制其 ACF(自相关) 与 PACF(偏自相关) 图，观察其截尾与拖尾性质即可确定模型阶数.
相关图的绘制代码如下：

1
2
plot_acf(ts_log).show()
plot.pacf(ts_log).show()

例如最终识别为 \(p=q=1\)，则预测模型为 \(ARMA(1,1)\)，实现代码如下：

1
2
model = ARMA(ts_log, order = (1,1))
result = model.fit(disp = -1, method = 'css')

随后用以下代码即可预测：

1
2
predict = result.predict(begin, end) # begin 和 end 表示所需预测的起点与终点
predict = np.exp(predict) # 若之前取了对数 or 做了差分, 预测之后应还原

Python中的ARIMA模型

没有什么卵用的模型… 多了个参数表示差分的阶数. 手动差分效果应当完全没差.

  1. 参见 https://zhuanlan.zhihu.com/p/35128342↩

  2. python时间序列分析: https://www.cnblogs.com/foley/p/5582358.html↩

  3. 将 “yyyymmdd” 转换为 “yyyy/mm/dd” 的方法很简单，以下给出 C++ 代码：

    date_format.cpp 处理日期格式
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    #include <bits/stdc++.h>
    using namespace std;

    int f(char c)
    {
        return c - '0';
    }

    int main()
    {
        freopen("out.txt", "w", stdout); // 输入到文件里, 方便使用
        std::string s; // 输入格式为 yyyymmdd 字符串
        while(std::cin >> s)
        {
            int y = 1000 * f(s[0]) + 100 * f(s[1]) + 10 * f(s[2]) + f(s[3]);
            int m = 10 * f(s[4]) + f(s[5]);
            int d = 10 * f(s[6]) + f(s[7]);
            printf("%d/%d/%d\n", y, m, d);
        }

        return 0;
    }
    ↩