时间序列模型简介

平稳序列
判断样本的平稳性
时间序列模型
参数选择
实验代码

1. 平稳序列

时间序列是一列观测值 $X_t$ 的集合, 其中每个观测值是在时段 $t$ 观测所得( $t$ 是自然数 ). 给定时间序列 $\{X_t\}_{t=1}^n$ , 如果对任意的 $t=1,\ldots,n$ , 它满足下列条件:
i. $\text{var}(X_t)<\infty$
ii. $\mathbb{E}(X_t) = \mu$
iii. $\text{cov}(X_r, X_s) = \text{cov}(X_{r+t}, X_{s+t}), \forall r,s = 1, ..., n$
我们把它叫做(弱)平稳(weakly stationary)序列.(下文我们简称平稳序列.)

通俗地讲, 平稳序列的期望, 方差, 协方差不随时间变化. 例如, $X_t$ 服从同一个分布时, 它是平稳的.

例1 下图中的时间序列由 $X_t\sim N(2, 1)$ 生成. 从直观上看, 这个序列是"平稳的".

Fig1. 平稳的时间序列

例2 下图的中的时间序列由 $X_t= 2 + 0.9X_{t-1}+W_t$ 生成, 其中 $X_0=0$ , $W_t= N(0, 1)$ . 它起初有明显地增长, 然后趋于平稳. 利用ADF检验(详情见下文), 我们发现该序列是平稳的(p-value < 0.01).

Fig2. 平稳的时间序列

Remark 弱平稳性的"弱"主要体现在时间序列在全局上是平稳的, 即，时间序列局部是波动的，但整体上看是平稳的, 或者随着时间的变化其样本的均值收敛.

2. 判断样本的平稳性

我们用统计学中假设检验的方法来判断样本的平稳性. 常用的是Augmented Dickey-Fuller(ADF)检验^[1].

Null Hypothesis $(H_0)$ : 样本中存在unit root^[2]. 如果接受 $H_0$ , 则意味着样本是非平稳的.
Alternate Hypothesis $(H_1)$ : 样本中不存在unit root. 如果拒绝 $H_0$ , 则意味着样本是平稳的.

在显著水平 $\alpha=0.05$ 的条件下, 我们可以通过计算p-value来接受或者拒绝 $H_0$ :

$\text{p-value} > 0.05$ : 接受 $H_0$ .
$\text{p-value} \leq 0.05$ : 拒绝 $H_0$ .

Python3中statsmodels.tas.stattools中的adfuller函数^[3]实现了ADF检验. 使用方法如下所示.

from statsmodels.tsa.stattools import adfuller


def test_stationarity(data, alpha=0.05, print_detail=True):
    """ Test stationarity of time series data.
    :param data: time series data, formatted as list
    :param alpha: significance level.
    :param print_detail: if True print additional information. 
    """
    result = adfuller(data)
    is_stationary = True if result[1] <= alpha else False
    if print_detail:
        print('ADF statistic: %f' % result[0])
        print('p-value: %f' % result[1])
        print('critical values:')
        for key, value in result[4].items():
            print('\t%s: %.3f' % (key, value))
    return is_stationary

3. 时间序列模型

前面之所以介绍平稳序列的概念及检验方法, 是因为它是很多基础的时间序列模型的前提假设. 在本节我们介绍一些常见的时间序列模型(更多内容可以参考^[4], ^[5]).

AR(p)

AR代表自回归(Autoregression). 假设时间序列 $\{X_t\}$ 是平稳的, 它可以被表示成如下形式:
$X_t = \delta + \sum_{i=1}^p \phi_iX_{t-i} + W_t.$

$\delta$ 是常数.
$W_t\overset{\text{iid}}{\sim} N(0, \sigma^2)$ , 它表示时段 $t$ 的误差(随机变量).
$p$ 代表自回归阶数.

MA(q)

MA代表移动平均(Moving Average). 假设时间序列 $\{X_t\}$ 是平稳的, 它可以被表示成如下形式:
$X_t = \delta + \sum_{i=1}^q \theta_iW_{t-i} + W_t.$

$\delta$ 是常数.
$W_t\overset{\text{iid}}{\sim} N(0, \sigma^2)$ , 它表示时段 $t$ 的误差(随机变量).
$q$ 代表移动平均阶数.

ARMA(p,q)

ARMA模型是AR和MA的组合. 假设同上. 它可以被表示为如下形式:
$X_t = \delta + \sum_{i=1}^p \phi_iX_{t-i} + \sum_{i=1}^q \theta_iW_{t-i} + W_t.$

$p$ 是自回归阶数
$q$ 是移动平均阶数

ARIMA(p,d,q)

ARIMA模型是ARMA模型的推广, 全称是Autoregressive Integrated Moving Average. 当时间序列 $\{X_t\}$ 不满足平稳性时, 我们通常使用差分的技巧把序列变得平稳, 然后再应用ARMA模型.

参数 $d$ 代表差分的阶数. 下面是差分的计算公式( $\Delta$ 为差分算子):

一阶差分 $\Delta X_t = X_t - X_{t-1}$
二阶差分 $\Delta^{(2)} X_t = \Delta(X_t-X_{t-1}) = X_t - 2X_{t-1}+X_{t-2} = (1-\Delta)^2X_t$
$d$ 阶差分 $\Delta^{(d)} X_t = (1-\Delta)^dX_t$

例3 下图是原始的时间序列. 通过观察, 它的均值有明显的上升趋势且不收敛, 因此不是平稳序列(ADF检验的p-value为0.94).

对该序列进行一阶差分后, 我们得到如下平稳的时间序列(p-value为0.00).

ARIMA(p,d,q) $\times$ (P,D,Q,s)

该记号代表季节性(或周期性)ARIMA模型, 详细的表达式可以参考^[4](4.1 Seasonal ARIMA models), 其中

$p,d,q$ 的意义同上.
$P$ 代表周期性自回归阶数(前 $P$ 个周期对应观测值的自回归).
$D$ 代表周期性差分阶数.
$Q$ 代表周期性移动平均阶数(前 $Q$ 个周期对应的移动平均).
$s$ 代表一个周期的长度.

我们可以把它看成两阶段模型: 第一阶段在全局使用ARIMA(p,d,q); 第二阶段通过指定周期长度 $s$ , 再利用ARIMA(P,Q,D)模型考虑周期之间的关系.

例4 考虑如下周期性的平稳时间序列( $s=18$ ).

season.png

对序列进行周期性差分: $Y_t = X_t - X_{t-18}$ 得到新的时间序列 $\{Y_t\}$ 如下图所示(红色部分)

deseason.png

通过使用周期性差分, 我们可以把原有时间序列的周期性移除. 同理, 通过采用周期性的自回归和移动平均系数, 我们可以把周期之间的依赖关系考虑进模型.

例5 考虑周期s=18的数据(蓝色曲线). 用 $\text{ARIMA}(1,0, 0)$ 和 $\text{ARIMA}(1, 0, 0)\times(0, 1, 0, 18)$ 分别进行预测的结果如下.

不考虑周期性的ARIMA模型的预测结果(灰色曲线)逐渐收敛到时间序列的均值. 由于序列是平稳的, 这样的预测结果符合我们的期望. 考虑到该时间序列有比较强的周期性, 且通过观察发现周期 $s=18$ . 在本例中, 我们仅使用周期差分, 最终得到了如图所示(红色曲线)的周期性预测结果.

ARCH(p)

ARCH的全称是Autoregressive Conditionally Heteroscedasticity, 它可以用来考虑样本的方差随着时间变化(或震荡)的时间序列. 设时间序列 $\{X_t\}$ 是平稳的, $\text{ARCH}(p)$ 模型可以被表示成如下形式:

$X_t = \sigma_tW_t,\quad W_t \overset{\text{iid}}{\sim} N(0,1)$
其中
$\sigma_t^2 = \alpha_0 + \sum_{i=1}^p \alpha_iX_{t-i}^2.$

$p$ 代表 $X_t^2$ 的自回归阶数.

GARCH(p,q)

GARCH即Generalized ARCH, 是ARCH模型的推广^[6]. 设时间序列 $\{X_t\}$ 是平稳的, $\text{GARCH}(p,q)$ 模型可以被表示成如下形式:
$X_t = \sigma_tW_t,\quad W_t \overset{\text{iid}}{\sim} N(0,1)$
其中
$\sigma_t^2 = \alpha_0 + \sum_{i=1}^p \alpha_iX_{t-i}^2 + \sum_{i=1}^q \beta_i\sigma^2_{t-i}.$

$p$ 代表 $X_t^2$ 的自回归阶数.
$q$ 代表 $\sigma_t^2$ 的移动平均阶数.

Remark ARCH/GARCH随机过程产生的数据是什么样的? 前面提到它们允许样本的方差随时间变化, 但是由于 $\{X_t\}$ 必须满足平稳性(前提假设), 因此样本的方差从局部看是变化(震荡)的, 但从整体看应该是"平稳的"序列. 例如下图是一个 $\text{GARCH}(1,1)$ 过程生成的时间序列( $\alpha_0=5, \alpha_1=\beta_1=0.5$ ).

VAR(p)

VAR即Vector Autoregression, 它是多变量的自回归模型. 类似地, 我们有 $\text{VARMA}(p, q)$ , 它是 $\text{ARMA}(p,q)$ 的向量版本. 需要注意的是, VARMA模型处理的时间序列可以有趋势. 我们不做详细的展开, 感兴趣的读者可以参考^[4]章节11.2: Vector Autoregressive models VAR(p) models.

4. 参数选择

给定时间序列的观测样本, 选定预测模型之后如何确定模型的参数? 本节我们介绍两种常用的方法: 1. 画出ACF/PACF图, 然后观察出 $p,q$ 的值; 2. 通过计算相关的统计指标, 自动化地选择参数.

4.1 观察ACF/PACF

ACF

ACF的全称是Autocorrelation Function. 对变量 $h=1,2, ...$ , ACF的值代表 $X_t$ 与 $X_{t-h}$ 之间的相关性.
$\text{ACF}(h) = \frac{\text{cov}(X_t, X_{t-h})}{\text{var}(X_t)}$

PACF

PACF的全称是Partial Autocorrelation Function. 对变量 $h=1,2,...$ , PACF的值代表已知 $X_{t-1}, ... X_{t-h+1}$ 的条件下, $X_t$ 与 $X_{t-h}$ 之间的相关性.
$\text{PACF}(h) = \frac{\text{cov}(X_t,X_{t-h}|x_{t-1}, ... x_{t-h+1})}{\sqrt{\text{var}(X_t|x_{t-1}, ... x_{t-h+1})\text{var}(X_{t-h}|x_{t-1}, ... x_{t-h+1})}}$

例6 设 $W_t\sim N(0, 1)$ . 考虑下面三个模型生成的时间序列, 并计算相应的ACF/PACF.

$\text{AR}(1): X_t=2 + 0.5X_{t-1}+W_t$
$\text{MA}(1): X_t=2 + 0.5W_{t-1}+W_t$
$\text{ARMA}(1,1): X_t=2 + 0.5X_{t-1}+ 0.5W_{t-1}+W_t$

指导原则(参考^[4])

模型	ACF	PACF	说明
AR(p)	逐渐趋近于0或像正弦曲线一样收敛到0	前 $p$ 个值非常显著, 其余的值不显著	$p$ 值主要参考PACF
MA(q)	前 $q$ 个值非常显著, 其余的值不显著	逐渐趋近于0或像正弦曲线一样收敛到0	$q$ 值主要参考ACF
ARMA(p,q)	逐渐趋近于0或像正弦曲线一样收敛到0	逐渐趋近于0或像正弦曲线一样收敛到0	$p,q$ 值靠猜

4.2 自动化地决定参数

基本思想是通过计算一些指标, 并选择参数使得相关的指标值尽可能小. 下面我们介绍一些常用的指标.

为方便描述, 我们先定义一些记号.

$n$ = 样本的大小
$k$ = 模型中需要拟合的参数数量(例如正态分布有的数量是2: $\mu$ 和 $\sigma$ )
$L_{\max}$ = 通过最大似然估计得到的最大Likelihood

AIC(Akaike Information Criterion)^[7]

$\text{AIC} = 2k - 2\ln(L_{\max})$

AICc^[8]

(AIC的改良版, 解决小样本过拟合的问题)
$\text{AICc} = \text{AIC} + \frac{2k^2 + 2k}{n-k-1}$

BIC(Bayesian Information Criterion)^[9]

(也称为Schwartz Criterion, SBC, SBIC)

$\text{BIC} = \ln(n)k - 2\ln(L_{\max})$

HQIC(Hannan–Quinn Information Criterion)^[10]

$\text{HQIC} = 2k\ln(\ln(n)) -2\ln(L_{\max})$

Remark 建议在实际中综合考虑这些指标.

5. 实验代码

Python3 code on Github

参考文献

Wikipedia. Augmented Dickey-Fuller test. ↩
Dickey, D. A.; Fuller, W. A. Distribution of the Estimators for Autoregressive Time Series with a Unit Root. Journal of the American Statistical Association. 74(366): 427–431, 1979. ↩
https://www.statsmodels.org/dev/generated/statsmodels.tsa.stattools.adfuller.html ↩
Lecture notes of Applied Time Series Analysis (SATA 510). The Pennsylvania State University. ↩ ↩ ↩ ↩
Jan Grandell. Time series analysis (lecture notes). ↩
Bollerslev, T. Generalized autoregressive conditional heteroscedasticity. Journal of Econometrics, 31, 307–327, 1986. ↩
Hirotugu Akaike. A new look at the statistical model identification, IEEE Transactions on Automatic Control, 19 (6): 716–723, 1974. ↩
Wikipedia. https://en.wikipedia.org/wiki/Akaike_information_criterion#AICc. ↩
Schwartz E.S. The Stochastic Behavior of Commodity Prices: Implications for Valuation and Hedging'. J Finance 52(3) Papers and Proceedings Fifty-Seventh Annual Meeting, American Finance Association, New Orleans, Louisiana, 923-973, 1997. ↩
Hannan, E. J. and B. G. Quinn. The Determination of the order of an autoregression, Journal of the Royal Statistical Society, Series B, 41: 190–195, 1979. ↩

最后编辑于：2018.10.20 10:46:09

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 159,015评论 4赞 362
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 67,262评论 1赞 292
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 108,727评论 0赞 243
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 43,986评论 0赞 205
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 52,363评论 3赞 287
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 40,610评论 1赞 219
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 31,871评论 2赞 312
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 30,582评论 0赞 198
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 34,297评论 1赞 242
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 30,551评论 2赞 246
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 32,053评论 1赞 260
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 28,385评论 2赞 253
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 33,035评论 3赞 236
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 26,079评论 0赞 8
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 26,841评论 0赞 195
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 35,648评论 2赞 274
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 35,550评论 2赞 270