用python处理时间序列数据，检验平稳性跟纯随机性
from statsmodels.tsa.stattools import adfuller as adf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import pandas as pd
import numpy as np
!pip install statsmodels
Requirement already satisfied: statsmodels in c:\programdata\anaconda3\lib\site-packages (0.11.0)
Requirement already satisfied: numpy>=1.14 in c:\programdata\anaconda3\lib\site-packages (from statsmodels) (1.18.1)
Requirement already satisfied: scipy>=1.0 in c:\programdata\anaconda3\lib\site-packages (from statsmodels) (1.4.1)
Requirement already satisfied: pandas>=0.21 in c:\programdata\anaconda3\lib\site-packages (from statsmodels) (1.0.1)
Requirement already satisfied: patsy>=0.5 in c:\programdata\anaconda3\lib\site-packages (from statsmodels) (0.5.1)
Requirement already satisfied: pytz>=2017.2 in c:\programdata\anaconda3\lib\site-packages (from pandas>=0.21->statsmodels) (2019.3)
Requirement already satisfied: python-dateutil>=2.6.1 in c:\programdata\anaconda3\lib\site-packages (from pandas>=0.21->statsmodels) (2.8.1)
Requirement already satisfied: six in c:\programdata\anaconda3\lib\site-packages (from patsy>=0.5->statsmodels) (1.14.0)
data=pd.read_excel('./data.xls',encoding='utf-8')
data
time
wc(误差随机项)
xt1
xt2
xt3
0
1
1.74
1.000000
1.000
1.000
1
2
-0.70
0.800000
0.800
0.800
2
3
-1.28
2.650000
1.440
1.430
3
4
0.43
7.440000
5.094
4.060
4
5
0.24
13.804000
10.204
5.645
...
...
...
...
...
...
95
96
1.55
739.086685
6908.698
146.490
96
97
0.07
748.322011
7056.498
147.345
97
98
-0.73
756.383207
7205.072
148.140
98
99
0.66
765.129924
7356.450
150.725
99
100
-0.44
773.187954
7508.420
151.515
100 rows × 5 columns
对X1做平稳性检验
xt1=data.xt1
dftest=adf(xt1)
pd.Series(dftest[0:4],index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
#p值高达0.9几
Test Statistic                  0.678947
p-value                         0.989408
#Lags Used                      1.000000
Number of Observations Used    98.000000
dtype: float64
xt1.plot()
[output_6_1.png)]
#一阶差分
xt1_1 = xt1.diff(1)
xt1_1.plot()
[(output_7_1.png)]
#一阶差分的单位根检验
dftest_1 = adf(xt1_1.dropna())
pd.Series(dftest_1[0:4],index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
#一阶差分后p值小于0.05,拒绝原假设(即不存在单位根，认为其已经平稳)
Test Statistic                -6.056515e+00
p-value                        1.241430e-07
#Lags Used                     0.000000e+00
Number of Observations Used    9.800000e+01
dtype: float64
dftest_1[1]-0.05
-0.049999875856979056
#画一阶差分之后的自相关图跟偏自相关图
plot_acf(xt1_1.dropna())
[(output_10_0.png)]
plot_pacf(xt1_1.dropna())
[(output_11_0.png)]
X2做平稳性检验
data.xt2.plot()
(output_13_1.png)]
#做二阶差分
x2_2=data.Xt2一阶差分.diff(1)
x2_2.plot()
[output_14_1.png)]
#对二阶差分后的Xt2做单位根检验
adf(x2_2.dropna())
#故拒绝原假设，该序列平稳
(-7.1910449486700525,
2.5032477359463947e-10,
4,
93,
{'1%': -3.502704609582561,
'5%': -2.8931578098779522,
'10%': -2.583636712914788},
260.6359245108364)
#画二阶差分后的Xt2的自相关图跟偏自相关图
plot_acf(x2_2.dropna())
(output_16_0.png)]
plot_pacf(x2_2.dropna())
(output_17_1.png)]
对X3做平稳性检验
data.xt3.plot()
(output_19_1.png)]
#做一阶差分
data.xt3.diff(1).plot()
(output_20_1.png)]
#对一阶差分后的序列做单位根检验
adf(data.xt3.diff(1).dropna())
(-10.661639595719135,
4.391819453885797e-19,
1,
97,
{'1%': -3.4996365338407074,
'5%': -2.8918307730370025,
'10%': -2.5829283377617176},
257.9188344687909)
adf(data.xt3.diff(1).dropna())[1]-0.05
#故拒绝原假设，该序列平稳
-0.05
#对一阶差分后的Xt3画出自相关图跟偏自相关图
plot_acf(data.xt3.diff(1).dropna())
(output_23_1.png)]
plot_pacf(data.xt3.diff(1).dropna())
(output_24_1.png)]

我们在拿到一个观察值序列之后，首先是判断它的平稳性。通过平稳性的检验，序列可以分为平稳序列和非平稳序列两大类。（通过平稳性检验判断序列的平稳性）

对于非平稳序列，由于它不具有二阶平稳的性质，所以对它的统计分析要费一些周折，通常要通过进一步的检验、变换或处理，才能确定适当的拟合模型。

对于平稳序列，情况就简单多了。我们有一套非常成熟的平稳序列建模方法。但是，并不是所有的平稳序列都值得建模。只有那些序列值之间具有密切的相关关系，历史数据对未来的发展有一定影响的序列，才值得我们花时间去挖掘历史数据中的有效信息，用来预测序列未来的发展。

如果序列值彼此之间没有任何相关性，那就意味着该序列是一个没有记忆的序列，过去的行为对将来的发展没有丝毫影响，这种序列称为纯随机序列。从统计角度而言，纯随机序列是没有任何分析价值的序列。

确定了序列是平稳序列后，还要进行随机性检验来确定平稳序列是否值得继续分析下去。
纯随机序列，也称为白噪声序列。之所以称为白噪声序列，是因为人们最初发现白光具有这种特性。容易证明白噪声序列一定是平稳序列，而且是最简单的平稳序列。
1.随机产生1000个服从标准正态分布的白噪声序列观察值，并绘制时序图。

正态分布随机数生成函数是rnorm。rnorm函数的命令格式为：

rnorm(n=, mean=, sd=)

式中：

n：随机数个数

mean：均值，缺省值默认为0

sd：标准差，缺省值默认为1

注：如果要产生n个服从标准正太分布的随机数，可以简写为rnorm(n)
所以本次命令与输出结果如下。
> white_noise<-rnorm(1000)
> white_noise<-ts(white_noise)
> plot(white_noise)



白噪声序列的性质

（1）随机性

这说明白噪声序列的各项之间没有任何相关关系，这种“没有记忆”的序列就是纯随机序列。一旦某个随机事件呈现出纯随机波动的特征，就认为该随机事件没有包含任何值得提取的有用信息，我们就应该终止分析。

（2）方差齐性

在时间序列分析中，方差齐性是一个非常重要的限制条件。因为根据马尔可夫定理，只有方差齐性假定成立时，我们用最小二乘法得到的未知参数估计值才是准确、有效的。如果假定不成立，最小二乘法就不是方差最小线性无偏估计，拟合模型的精度会受到很大影响。

所以我们在进行模型拟合时，检验内容之一就是要检验拟合模型的残差是否满足方差齐性假定。如果不满足，那就说明残差序列还不是白噪声序列，即拟合模型没有提取随机序列中的相关信息，这时拟合模型的精度是值得怀疑的。在这种场合下，我们通常需要使用适当的条件异方差模型来处理异方差信息。
纯随机性检验

纯随机性检验也称白噪声检验，是专门用来检验序列是否为纯随机序列的一种方法。我们知道如果一个序列是纯随机序列，那么它的序列值之间应该没有任何相关关系。但这是一种理论上上才会出现的理想状况。实际上，由于观察值序列的有限性，导致纯随机序列的样本自相关系数不会绝对为零。
绘制上述随机产生1000个服从标准正态分布的白噪声序列的样本自相关图。
> acf(white_noise)



样本自相关图显示这个纯随机序列没有一个样本自相关系数严格等于零。但这些自相关系数确实都非常小，都在零值附近以一个很小的幅度随机波动。这提醒我们应该考虑样本自相关系数分分布性质，从统计意义上来判断序列的纯随机性质。
R语言中使用Box.test函数进行纯随机性检验（白噪声检验）。该函数的命令格式为：

Box.test(x, type=, lag= )

式中：

x：变量名

type：检验统计量类型

（1）type=“Box-Pierce”，输出白噪声检验的Q统计量，该统计量为系统默认输出结果。

（2）TYPE=“Ljung-Box”，输出白噪声检验的LQ统计量。

lag：延迟阶数，lag=表示输出滞后n阶的白噪声检验统计量，忽略该选项时，默认输出滞后1阶的检验统计量结果。
计算上述随机产生1000个服从标准正态分布的白噪声序列延迟6期、延迟12期的QLB统计量的值，并判断该序列的随机性（α=0.05）



R是最小结果输出软件，所以每次调出Box.test函数，只能给出一个检验结果。上面为了得到延迟6阶和延迟12阶两个LB统计量的结果，我们调用两次Box.test函数。如果希望能够像其他统计软件一样，一次得到多个白噪声检验结果，可以编写一个循环命令。最简单的循环函数是for函数，它的命令格式如下：

for(x, in, n1:n2) state

式中：

x：循环变量名

n1：n2：给出的循环取值区间

state：需要循环执行的命令

于是等价命令和输出结果为：



由于P值显著性水平α，所以该序列不能拒绝纯随机的原假设，换言之，我们可以认为该序列的波动没有任何统计规律可循，因而可以停止对该序列的统计分析。

还需要解释的一点是，为什么上述检验了前6期和前12期延迟的Q统计量和LB统计量就直接判断该序列是白噪声序列？为什么不进行全部999期延迟检验？

这是因为平稳序列通常具有短期相关性，如果序列值之间存在显著的相关关系，通常只存在于延迟时期比较短的序列值之间。所以，如果一个平稳序列短期延迟的序列值之间都不存在显著的相关关系，通常长期延迟之间就更不会存在显著的相关关系了。

另一方面，假如一个平稳序列显示出显著的短期相关性，那么该序列就一定不是白噪声序列，我们就可以对序列值之间存在的相关性进行分析。因为平稳序列只要延迟时期足够长，自相关系数都会收敛于零。
以1949-1998年北京市每年最高气温序列为例进行平稳性与纯随机性检验。

1.绘制1949-1998年北京市每年最高气温序列时序图。
> b<-read.table("file6.csv",sep=",",header=T)
> temp<-ts(b$temp,start=1949)
> plot(temp)



上图显示北京市每年的最高气温始终围绕在37℃附近随机波动，没有明显趋势或周期，基本可以视为平稳序列。但为了稳妥起见，我们还需要利用自相关图进一步辅助识别。
2.绘制1949-1998年北京市每年最高气温序列自相关图
> acf(temp)



自相关图显示该序列的自相关系数一直都比较小，始终控制在2倍标准差范围以内，可以认为该序列自始至终都在零轴附近波动，这是随机性非常强的平稳时间序列通常具有的自相关特征。
3.对1949-1998年北京市最高气温序列做白噪声检验（α=0.05）
> for(i in 1:2)print(Box.test(temp,lag=6*i))



根据这个检验结果，不能拒绝序列纯随机的原假设。因而可以认为北京市最高气温的变动属于纯随机波动。这说明我们很难根据历史信息预测未来年份的最高气温。至此，对该序列的分析也就结束了。

引言
根据Wold分解定理，任意一个离散平稳时间序列都可以分解为一个确定性平稳序列和一个随机性平稳序列之和。且确定性序列可以表达为历史序列值的线性函数，而随机性序列可以表达为历史新息（历史纯随机波动）的线性组合，即
上式在统计上被称为自回归移动平均模型 (auto-regression moving average), 简称为 ARMA模型。
Wold分解定理保证了平稳序列一定可以用某个 ARMA 模型等价表达, 所以 ARMA模型是目前最常用的平稳序列拟合与预测模型。
ARMA 模型实际上是一个模型族, 它又可细分为 AR 模型, MA 模型和 ARMA 模型三大类。
AR模型的定义
具有如下结构的模型称为p阶自回归模型，简记为  
特别当  
 时，称为中心化
模型
令                  
，称 
是
的中心化序列
自回归系数多项式
引进延迟算子，中心化 
 模型又可以简记为 
称下式为p阶自回归系数多项式
延迟算子
延迟算子类似于一个时间指针，当前序列值乘以一个延迟算子，就相当于把当前序列值的时间向过去拨了一个时刻 
记B为延迟算子，有 
所以
模型的简写形式如下导出
延迟算子的性质
AR模型平稳性判别 
判别原因
要拟合一个平稳序列的发展, 用来拟合的模型显然也应该是平稳的。AR 模型是常用的平稳序列的拟合模型之一, 但并非所有的 AR 模型都是平稳的。 
判别方法
特征根判别法
平稳域判别法
例1
考察如下四个模型的平稳性
平稳特征
非平稳特征
特征根判别
p阶自回归序列平稳，要求p个非零特征根都在单位圆内，即
在引入延迟算子之后, 我们还可以推导出跟特征根判别等价的性质: p阶自回归序列平稳的条件是自回归系数多项式的p 个根都在单位圆外
平稳域判别
对于一个
模型而言，如果没有平稳性的要求，实际上也就意味着对参数向量没有任何限制，它们可以取遍维欧氏空间的任意一点
如果加上了平稳性限制，参数向量就只能取维欧氏空间的一个子集，使得特征根都在单位圆内的系数集合
对于低阶自回归模型用平稳域的方法判别模型的平稳性通常更为简便
 AR(1)模型平稳条件
方程结构
特征根
平稳域
AR(2)模型的平稳条件
方程结构
特征根
平稳域
AR(2)的平稳域
例1续
分别用特征根判别法和平稳域判别法检验例3-1中四个 AR 模型的平稳性
平稳性判别
平稳AR模型的统计性质——均值 
如果AR(p)模型满足平稳性条件，则有
根据平稳序列均值为常数，且
为白噪声序列，有
推导出
平稳AR模型的统计性质——方差
要得到平稳 AR(p)模型的方差, 需要借助 Green函数的帮助
Green函数的定义
假设 
 为任意p阶的平稳 AR 模型, 那么一定存在一个常数序列 
  使得
可以等价表达为纯随机序列
的线性组合, 即
这个常数序列
就称为 Green函数 
基于Green函数，可以求出AR(p)模型的方差为
Green函数的递推公式
原理
方法：待定系数法
例2
求平稳 AR(1) 模型 
 的 Green函数的递推公式, 并基于 Green函数求解 AR(1) 模型的方差。
平稳AR(1)模型的Green函数递推公式为
平稳AR(1)模型的方差为
平稳AR模型的统计性质——协方差函数
在平稳AR(p)模型两边同乘 
，再求期望
根据
得协方差函数的递推公式
例3
求平稳AR(1)模型的协方差
递推公式
因为平稳AR(1)模型的方差为
所以协方差函数的递推公式为
例4
求平稳AR(2)模型的协方差
平稳AR(2)模型的协方差函数递推公式为
平稳AR模型的统计性质——自相关系数
自相关系数的定义
平稳AR(P)模型的自相关系数递推公式
常用AR模型自相关系数递推公式
AR(1)模型
AR(2)模型
AR模型自相关系数的性质
AR模型的自相关系数的表达式实际上是一个齐次差分方程，它的通解形式为
根据自相关系数的通解形式，可以判断AR模型的自相关系数具有如下特征
呈指数衰减
拖尾性
例5
考察如下AR模型的自相关图
图示解释
从上图中可以看到, 这四个平稳 AR 模型, 不论它们是 AR(1)模型还是 AR(2)模型, 不论它们的特征根是实根还是复根, 是正根还是负根, 它们的自相关系数都呈现出拖尾性和呈指数衰减到零值附近的性质。
但由于特征根不同, 它们的自相关系数衰减的方式也不一样 
有的是按负指数单调衰减 (如模型(1))
有的是正负相间地衰减 (如模型(2)) 
有的呈现出类似于周期性的余弦衰减, 即具有 “伪周期” 特征 (如模型(3))
有的是不规则衰减（如模型(4)）
偏自相关系数(partial autocorrelation function,PACF )
时间序列过程的偏自相关函数就是时间序列在两个时间随机变量之间，排除了其间各个时间随机变量影响的相关系数。
偏自相关系数的定义
对于平稳 
序列，所谓滞后k偏自相关系数就是指在给定中间 
 个随机变量   
的条件下，或者说，在剔除了中间k-1个随机变量的干扰之后， 
 对
 影响的相关度量。用数学语言描述就是
偏自相关系数的计算
基于Yule-Walker方程组计算偏自相关系数
在方程      
等号两边同时乘以 
，等号两边求期望再除以方差，得
取前k个方程构成的方程组即Yule-Walker方程组
解Yule-Walker方程组可以得到参数
的解，最后一个参数的解即为延迟K偏自相关系数
AR(1)模型偏自相关系数的计算
AR(1)模型
Jule-Walker方程
偏自相关系数的解
AR(2)模型偏自相关系数的计算
基于矩阵结构计算偏自相关系数
证明AR(p)模型偏自相关系数p阶截尾
所谓p 阶截尾, 是指 
  。 要证明这一点, 实际上只要能证明当 
时,       
 即可。
例5续
求如下AR模型的偏自相关系数，并考察它们的偏自相关图特征
1、已知AR(1)模型为: 
解：
2、已知某AR(2)模型为:
解：
3、已知某AR(2)模型为:
解：
4、已知AR(2)序列为:
解：
5、证明对任意常数c,如下定义的AR(3)序列-定是非平稳序列:
解：
6、对于AR(1)模型,
判断如下命题是否正确

这一部分是时间序列预处理R语言的实现。

目标是将课本和上课知识点整合。

老师是用一节课讲完的，本篇文章只做了平稳性检验~~~

下一篇再写纯随机性检验

全部代码
#input data
yield <- c(15.2,16.9,15.3,14.9,15.7,15.1,16.7)

#对比向量和时间序列画图的差别
par(mfrow=c(1,2))
plot(yield)

yield <- ts(yield,start = 1884)
plot(yield)

#介绍plot函数
help(plot)
plot(yield,type = "o",pch = 5,lty=1,lwd=3,col=1,
     main="yield",xlab="Year",ylab="yield",xlim=c(1884,1890),ylim=c(14.5,17))
abline(v=c(1889),lty=2)#添加一条线/多条线
abline(h=16)




#计算自相关系数
acf(yield,lag.max = 24)
AC <- acf(yield)
AC$acf
help(acf)
abline(h=1.96*1/sqrt(7),col=3)#α=0.05时，α/2 的分位点————假设检验。超出了线就是拒绝原假设，说明自相关系数不等于0.
abline(h=- 1.96*1/sqrt(7),col=2)

T <- length(yield)
cor(yield[1:T-1],yield[2:T])


#平稳性检验
sha <- read.csv()
output <- ts(sha$output,start=1994)
plot(output)
acf(output,lag.max=24,ylim=c(-1,1))

a <- read.csv()
milk <- ts(a$milk,start =c(1962,1) ,frequency = 12)
plot(milk)
acf(milk,lag.max=24,ylim=c(-1,1))

b <- read.csv()
temp <- ts(b$temp,start = 1949)
plot(temp)
acf(temp,ylim=c(-1,1))

平稳性检验的R语言的实现

1 对比时序图和向量区别
#input data
yield <- c(15.2,16.9,15.3,14.9,15.7,15.1,16.7)

#对比向量和时间序列画图的差别
par(mfrow=c(1,2))
plot(yield)
yield <- ts(yield,start = 1884)
plot(yield)

画两个图的时候运用par这个函数



可以看到两者的区别一个是散点图，另一个是连起来的。

x轴一个是时间，一个是第几个数。


2 基本画图设置
#介绍plot函数
help(plot)
plot(yield,type = "o",pch = 5,lty=1,lwd=3,col=1   main="yield",xlab="Year",ylab="yield",xlim=c(1884,1890),ylim=c(14.5,17))
abline(v=c(1889),lty=2)#添加一条线/多条线
abline(h=16)


plot函数在{graphics}程序包里
type:

“p” for points,

“l” for lines,

“b” for both, 线不穿过点

“c” for the lines part alone of “b”,

“o” for both ‘overplotted’, 线穿过点

“h” for ‘histogram’ like (or ‘high-density’) vertical lines,

“s” for stair steps,

“S” for other steps, see ‘Details’ below,

“n” for no plotting.
main：给主图加标题
sub：给子图加标题
xlab：给横坐标加标题
ylab：给纵坐标加标题
pch：有25种形式——形状


lty：
lwd：
color：
xlim ylim：x，y轴的取值



abline()：在已有的图上表示添加一条线/多条线

abline(v=16)：横着

abline(h=16)：竖着


3 计算自相关系数
#计算自相关系数
acf(yield,lag.max = 24)
AC <- acf(yield)
AC$acf
abline(h=1.96*1/sqrt(7),col=3)#α=0.05时，α/2 的分位点————假设检验。超出了线就是拒绝原假设，说明自相关系数不等于0.
abline(h=- 1.96*1/sqrt(7),col=2)



自相关图：


acf(yield,lag.max=···,type，na.action) ：

yield是数据
lag.max是滞后阶数
type = c(“correlation”, “covariance”, “partial”)，默认自相关函数，第二个自协方差，第三个偏自相关
na.action：对缺失值的处理




#计算自相关系数

acf(yield,lag.max = 24)

AC <- acf(yield)

AC$acf

，  ，  1

[,1]

[1,]  1.00000000

[2,] -0.35565808       #滞后一阶自相关

[3,] -0.07731467

[4,] -0.04470251

[5,] -0.28825962

[6,]  0.39394421

[7,] -0.12800933


滞后自相关：$\rho_0=1$$X_t和X_{t-1}\to\rho_1$$X_t和X_{t-2}\to \rho_2$$X_t和X_{t-3}\to \rho_3$$\dots和\dots\to\dots$$X_t和X_{t-6}\to \rho_6$

4 用相关函数的方法计算自相关函数
#滞后一阶
T <- length(yield)
cor(yield[1:T-1],yield[2:T])

T <- length(yield)

cor(yield[1:T-1],yield[2:T])

[1] -0.4182568

R语言自相关函数是：acf(yield),也可以运用教材23页公式计算自相关函数，因为是约等，所以计算结果会有一点差别，当样本量较大时差别会很小：

cov（yield[1:T-1],yield[2:T]）/var(yield[1:T])

cor(yield[1:T-1],yield[2:T])

5 平稳性检验的三个例子
#平稳性检验
sha <- read.csv()
output <- ts(sha$output,start=1994)
plot(output)
acf(output,lag.max=24,ylim=c(-1,1))
时序图：有趋势

自相关系数图：先正的递减，后负的递增。而平稳的是呈现衰减的趋势，所以不是平稳时间序列。


a <- read.csv()
milk <- ts(a$milk,start =c(1962,1) ,frequency = 12)
plot(milk)
acf(milk,lag.max=24,ylim=c(-1,1))
时序图：有周期

自相关图：未衰减


b <- read.csv()
temp <- ts(b$temp,start = 1949)
plot(temp)
acf(temp,ylim=c(-1,1))
时序图：无递增递减趋势，无周期趋势，初步判断是平稳的。

自相关图：前几阶比较大，后衰减接近0