时间序列分析的学习与应用（一）

前言
一：时间序列分析的一些基础性知识
二：相关系数和自相关函数
三：白噪声序列
四：平稳性检验原理
五：实现白噪声检验算法
六：编写单位根检验算法
七：总结

前言

接下来几篇的内容，我们系统的学习下时间序列模型相关内容，从自相关系数到检验统计量的解释，从ARIMA模型阶数截尾（拖尾）公式推导到GARCH模型的搭建等等的学习，并结合编写相关代码（或底层编写），最终应用于解决实际问题。

这里大概会进行如下几个内容的学习

基础知识的学习

1.时序概念与平稳性介绍
2.相关系数和自相关函数
3.白噪声序列与回归序列
4.单位根检验平稳的原理与实现

ARIMA（AR、MA）模型的学习

1.模型的基本原理介绍
2.模型的定阶方法与公式推导
3.模型的预测结果评价

（G）ARCH方差模型的学习

1.波动率模型的特征介绍
2.模型的基本原理介绍
3.模型的建立与预测

一：时间序列分析的一些基础性知识

从名字上就不难发现时间序列其实跟时间有关并以此产生离散数字组成的序列集合为时间序列。股价的走势，销量的走势，温度的走势等等都是可以构成时序。

时序特点

趋势性：在长期内呈现出持续向上或向下的波动
季节性：在一年内重复出现的周期性波动，常见的如天气、节假日活动等
周期性：这里的周期性不一定是固定长度的周期，主要特点是它不会朝着单一方向的持续变动，而是涨落相同的交替波动
无规则波动：除去以上三种波动，剩下的无规则随机波动可以称为无规则时序波动，一个完整的时序中尝常常伴有这种无规则小时序波动

平稳性

平稳性是时序建模能否成立的关键前提，我们有必要仔细理解一番

定义

如果给定所有的时刻 t t t 、任意正整数 k k k 和任意 k k k 个正整数， ( t 1 , t 2 , . . . t n ) (t_1,t_2,...t_n) (t1,t2,...tn) 和 ( r t 1 , r t 2 , . . . , r t k ) (r_{t_1},r_{t_2},...,r_{t_k}) (rt1,rt2,...,rtk) 的联合分布与 ( r t 1 + t , r t 2 + t , . . . , r t k + t ) (r_{t_1+t},r_{t_2+t},...,r_{t_k+t}) (rt1+t,rt2+t,...,rtk+t) 的联合分布相同，则时序 { r t } \{r_t\} {rt} 是强平稳序列。这里为什么说是强平稳序列呢，因为这里要求 ( r t 1 , r t 2 , . . . , r t k ) (r_{t_1},r_{t_2},...,r_{t_k}) (rt1,rt2,...,rtk) 的联合分布在时间的平移变换下保持不变，这一条件是非常强的，实际很难达到。所以在实际工业中，我们经常假定的平稳序列是一个较弱的方式，如下：
若时序 { r t } \{r_t\} {rt} 满足下面条件：
E ( r t ) = μ E(r_t)=\mu E(rt)=μ ( μ \mu μ 是常数)， C o v ( r t , r t − l ) = r l , r t Cov(r_t,r_{t-l})=r_l , r_t Cov(rt,rt−l)=rl,rt 只依赖于 l l l 且 l l l 是任意正数
则称时序 { r t } \{r_t\} {rt} 是弱平稳的，即该序列的均值和 r t r_{t} rt 与 r t − l r_{t-l} rt−l 的协方差不随时间而改变。

二：相关系数和自相关函数

相关系数

相关系数的公式定义应该是非常熟悉的了，假设两样本向量 X,Y 那么对应的相关系数为

ρ X Y = C o v ( X , Y ) V a r ( X ) V a r ( Y ) \rho_{XY}=\frac{Cov(X,Y)}{\sqrt{Var(X)Var(Y)}} ρXY=Var(X)Var(Y) Cov(X,Y) ，

那么把两样本向量按照展开为 ρ X Y = ∑ t = 1 T ( x t − x ˉ ) ( y t − y ˉ ) ∑ t = 1 T ( x t − x ˉ ) 2 ∑ t = 1 T ( y t − y ˉ ) 2 \rho_{XY}=\frac{\sum_{t=1}^{T}{(x_t-\bar{x})(y_t-\bar{y})}}{\sqrt{\sum_{t=1}^{T}{(x_t-\bar{x})^2\sum_{t=1}^{T}{(y_t-\bar{y})^2}}}} ρXY=∑t=1T(xt−xˉ)2∑t=1T(yt−yˉ)2 ∑t=1T(xt−xˉ)(yt−yˉ) ，

那么自相关函数又如何呢？

自相关函数

自相关函数系数度量了两个向量的线性相关性，而在平稳时序 { r t } \{r_t\} {rt} 中，我们想知道 r t r_t rt 与 r t − l r_{t-l} rt−l 的相关系数称为 r t r_t rt 的间隔为 l l l 的自相关系数 ρ l \rho_l ρl ，具体定义如下：

ρ l = C o v ( r t , r t − l ) V a r ( r t ) V a r ( r t − l ) = C o v ( r t , r t − l ) V a r ( r t ) \rho_l=\frac{Cov(r_t,r_{t-l})}{\sqrt{Var(r_t)Var(r_{t-l})}}=\frac{Cov(r_t,r_{t-l})}{Var(r_t)} ρl=Var(rt)Var(rt−l) Cov(rt,rt−l)=Var(rt)Cov(rt,rt−l) ，

这里等式右边成立主要是因为

C o v ( r t , r t − l ) = r l Cov(r_t,r_{t-l})=r_l Cov(rt,rt−l)=rl ，当 l = 0 l=0 l=0 时， C o v ( r t , r t ) = r 0 Cov(r_t,r_t)=r_0 Cov(rt,rt)=r0 为一常数，那么同理有：

ρ l = ∑ t = l + 1 T ( r t − r ˉ ) ( r t − l − r ˉ ) ∑ t = 1 T ( r t − r ˉ ) 2 , 0 ≤ l ≤ T − 1 \rho_{l}=\frac{\sum_{t=l+1}^{T}{(r_t-\bar{r})(r_{t-l}-\bar{r})}}{\sum_{t=1}^{T}{(r_t-\bar{r})^2}},0\leq l\leq T-1 ρl=∑t=1T(rt−rˉ)2∑t=l+1T(rt−rˉ)(rt−l−rˉ),0≤l≤T−1 。

三：白噪声序列

若随机过程 X t ( t = 1 , 2 , 3 , . . . ) X_t(t=1,2,3,...) Xt(t=1,2,3,...) 是由不相关的随机变量的序列构成的，即对于所有 s ≠ t s\ne t s=t ,随机变量 X t X_t Xt 与 X s X_s Xs 的协方差为0，则称 X t X_t Xt 为纯随机过程（白噪声）。
在真正建立时序模型时候，我们开始是有必要判断此时序是否是纯随机过程，如果是纯随机序列，那么就终止建模，因为纯随机序列没有任何可以提取的相关性信息。

检验方法

纯随机序列的检验方法，一般是自相关函数结合相关统计量来检验。我们认为自相关函数中，所有的值都为0时（实际上不可能完全为0，只是非常接近于0），则认为序列是完全不相关的，因此需要检验多个自相关系数，我们建立熟悉的Ljung-Box(LB)统计量来检验上述过程：
原假设： H 0 ： ρ 1 = ρ 2 = . . . = ρ m = 0 H_0：\rho_1=\rho_2=...=\rho_m=0 H0：ρ1=ρ2=...=ρm=0 ,

备择假设： H 1 : ∃ i ∈ 1 , 2 , . . . , m , ρ i ≠ 0 H_1:\exists i\in1,2,...,m,\rho_i\ne0 H1:∃i∈1,2,...,m,ρi=0 ,

LB检验统计量公式如下：
Q ( m ) = T ( T + 2 ) ∑ l = 1 m ρ l 2 T − l Q(m)=T(T+2)\sum_{l=1}^{m}{\frac{\rho_l^2}{T-l}} Q(m)=T(T+2)∑l=1mT−lρl2 ，这里注意 Q ( m ) Q(m) Q(m) 渐进服从自由度为 m m m 的 χ 2 \chi^2 χ2 分布。

决策规则：当 Q ( m ) > χ α 2 Q(m)>\chi^2_\alpha Q(m)>χα2,拒绝原假设 H 0 H_0 H0 ，即 Q ( m ) Q(m) Q(m) 的值大于自由度为 m m m 的 χ 2 \chi^2 χ2 分布 100 ( 1 − α ) 100(1-\alpha) 100(1−α) 分位点时，拒绝 H 0 H_0 H0 。

以上就是时序建模的相关预知识点，当然检验的方法不止LB一种，原理都是不是大相径同。

四：平稳性检验原理

一般时序模型的平稳性检验用的是单位根检验方法，这里，我们详细讲解一下单位根原理以及特征方程的运用。

非平稳介绍

我们令时序 { r t } \{r_t\} {rt} 满足表达式 r t = r t − 1 + ϵ t r_t=r_{t-1}+\epsilon_t rt=rt−1+ϵt ，其中 x 0 = 0 , ϵ t ∼ N ( 0 , σ ϵ 2 ) x_0=0,\epsilon_t\sim N(0,\sigma_\epsilon^2) x0=0,ϵt∼N(0,σϵ2) ，

那么通过简单递推公式我们有 r t = r t − 1 + ϵ t = r t − 2 + ϵ t − 1 + ϵ t + . . . = ∑ i = 1 t ϵ i r_t=r_{t-1}+\epsilon_t =r_{t-2}+\epsilon_{t-1}+\epsilon_t +...=\sum_{i=1}^{t}{\epsilon_i} rt=rt−1+ϵt=rt−2+ϵt−1+ϵt+...=∑i=1tϵi ，

那么很容易得到 E ( r t ) = 0 , V a r ( r t ) = t σ ϵ 2 E(r_{t})=0,Var(r_t)=t\sigma_\epsilon^2 E(rt)=0,Var(rt)=tσϵ2 ，

其中方差跟时间 t t t 有关，很明显这不满足平稳时序条件。那么从简单的一个例子我们发现，当在时序中，一直存在随机的扰动，且扰动没有减弱（如果在无穷处减弱，那么此序列还是弱平稳序列），则此序列为非平稳状态。

滞后算子

我们定义滞后算子 L L L 可得 L r t = r t − 1 , L 2 r t = r t − 2 , . . . ， Lr_t=r_{t-1},L^2r_t=r_{t-2},...， Lrt=rt−1,L2rt=rt−2,...，

那么对一般时序方程 r t = ϕ 0 + ϕ 1 r t − 1 + . . . + ϕ p r t − p + a 0 + θ 1 a t − 1 + . . . + θ q a t − q r_t=\phi_0+\phi_1r_{t-1}+...+\phi_pr_{t-p}+a_0+\theta_1a_{t-1}+...+\theta_qa_{t-q} rt=ϕ0+ϕ1rt−1+...+ϕprt−p+a0+θ1at−1+...+θqat−q ，

这里的 a t a_t at 为白噪声序列，总共为 q q q 阶， θ q \theta_q θq 为常数；这里的 r t r_t rt 即为 t t t 时刻预测的数据， a 0 、 ϕ 0 a_0、\phi_0 a0、ϕ0 都为一常数，然后根据上述算子改写为

r t = ∑ i = 1 q θ i a t − i + a 0 + ϕ 0 + ∑ j = 1 p ϕ j L j r t r_t=\sum_{i=1}^{q}\theta_i{a_{t-i}}+a_0+\phi_0+\sum_{j=1}^{p}{\phi_jL^jr_t} rt=∑i=1qθiat−i+a0+ϕ0+∑j=1pϕjLjrt ，

那么 r t = ∑ i = 1 q θ i a t − i + a 0 + ϕ 0 1 − ∑ j = 1 p ϕ j L j r_t=\frac{\sum_{i=1}^{q}\theta_i{a_{t-i}}+a_0+\phi_0}{1-\sum_{j=1}^{p}{\phi_jL^j}} rt=1−∑j=1pϕjLj∑i=1qθiat−i+a0+ϕ0 ，即改写为了相应的关于算子 L 的特征方程。

这里要注意，当 i → ∞ i\rightarrow\infty i→∞ 时，我们假定 θ i → 0 \theta_i\rightarrow0 θi→0 。

利用级数性质对平稳性判定

由上式最终的特征方程表达式我们知道，当分母不为0时，时序 { r t } \{r_t\} {rt} 要平稳存在，则分母表达式要满足等比级数收敛，那么对应的特征根的绝对值 ∣ L ∣ < 1 \left| L \right|<1 ∣L∣<1 的，也就是特征值的解的模要在单位圆中（否则当 ∣ L ∣ ≥ 1 \left| L \right|\geq1 ∣L∣≥1 时，时序 { r t } \{r_t\} {rt} 为一发散序列）。

五：实现白噪声检验算法

先是自相关函数的编写

在import pandas as pd
import numpy as npdata = pd.read_excel(r'...\demo_data.xls')##自编自相关系数
lag = 11
def Autocorr_function(rt,lag):corr_ls = []rtmean = np.mean(rt)Var = sum(np.power((rt - rtmean),2))##自相关系数分母的运算值for i in range(lag):res_temp = (np.array(rt[i:len(rt)]) - rtmean) * (np.array(rt[:len(rt)-i]) - rtmean)##自相关系数分子表达式求解corr_ls.append(sum(res_temp) / Var)return corr_ls
x = data['销量']
resacf = Autocorr_function(x,lag)
print('ACF_value：',resacf[1:])这里插入代码片

再是LB统计量的编写

###ljung_box统计量值的编写
def lbvalue(rt,lag,pl):T = len(rt)lbv = []for i in range(lag):l = np.array([j for j in range(i)])res_temp = np.power(pl[1:i+1],2) / (T - l - 1)##这里的编写逻辑要注意res_fin = T * (T + 2) * sum(res_temp)lbv.append(res_fin)return lbv[1:]
lb = lbvalue(x,lag,resacf)
print('lb：',lb)

最后是借用卡方累积函数求解 p p p值

##p值的求解
from scipy.stats import chi2
print('p_value：',1.0 - chi2.cdf(lb,[l for l in range(1,lag)]))

当然Python是提供求解p值的相关函数的，如果我们直接调用第三方模块，代码如下：

import statsmodels.api as sm
import pandas as pd
import numpy as np
def stat_che(data):m = 10# 检验10个自相关系数acf,q,p = sm.tsa.acf(data,nlags=m,qstat=True)out = np.c_[range(1,m+1),acf[1:],q,p]output = pd.DataFrame(out,columns=['lag','自相关系数','统计量Q值','p_values'])output = output.set_index('lag')# 设置第一列索引名称,可省略重复索引列1return output

六：编写单位根检验算法

由于干扰项存在对平稳性影响可以忽略，我们将基于自回归（AR）模型进行平稳性检验。

import statsmodels.api as sm
model = sm.tsa.AR(x)
result_AR = model.fit()# print('逆特征方程的解：',len(result_AR.roots),result_AR.roots)
nu_root = np.array(result_AR.roots)
nu_root_spe = np.power(nu_root,-1)ls = []
for i in nu_root_spe:res1 = i.realres2 = i.imagres = np.sqrt(res1 ** 2 + res2 ** 2)ls.append(res)
print('特征根的模：',ls)
for i in ls:if i <= 1:passelse:print('\033[1;31m由于存在模大于1的特征根，此序列不平稳！！\033[0m')break

可视化单位圆看看效果

import matplotlib.pyplot as plt
pi,sin,cos = np.pi,np.sin,np.cos
r = 1
theta = np.linspace(0,2 * pi,360)
x = r * cos(theta)
y = r * sin(theta)plt.plot(x,y,'r')
X = [0 for i in range(len(ls))]
plt.scatter(X,ls,marker='*')
plt.show()

以上就是通过单位根算法检验的结果，当然也可直接调用ADF函数进行平稳性检验。

七：总结

这一章我们主要介绍了时序模型的各个基本概念，从理论解释到算法实现，后面关于时间序列分析的学习主要涉及一些结论的公式推导与模型的建立，有阶数截尾与拖尾的由来，ARIMA与GARCH模型等等，望对有帮助有兴趣的朋友，点个赞收藏支持一下，在此表示感谢！