季节性单位根

正如MAT8181课程中所讨论的那样,至少有两种非平稳的时间序列:存在趋势的和存在单位根(这种类型被称为 单整的)。单位根测试不能用来评估一个时间序列是否平稳,这种方法只能检测整合的时间序列。 季节性单位根也是如此。

之前的文章中,我们已经看到,对小时温度建模很困难,因为大多数测试不会否决单位根假设。我们依然使用蒙特利尔的数据, 来考虑平均每月温度。

> montreal=read.table("http://freakonometrics.free.fr/temp-montreal-monthly.txt")
> M=as.matrix(montreal[,2:13])
> X=as.numeric(t(M))
> tsm=ts(X,start=1948,freq=12)
> plot(tsm)

对于那些不了解蒙特利尔的人来说,冬季夏季非常不同。我们可以使用可视化每个月份差异:

> month=rep(1:12,length(tsm)/12)
> plot(month,as.numeric(tsm))
> lines(1:12,apply(M,2,mean),col="red",type="b",pch=19)

或者,如果安装从CRAN仓库中被移除的uroot软件包,我们可以使用:

> library(uroot)
> bbplot(tsm)

或者:

> bb3D(tsm)
Loading required package: tclt

因为每年的季节模式,这个时间序列看起来是是循环的。下面是自相关函数:

> acf(tsm,lag=120)

同样的,我们可以可视化这个循环:

> persp(1948:2013,1:12,M,theta=-50,col="yellow",shade=TRUE,
+ xlab="Year",ylab="Month",zlab="Temperature",ticktype="detailed")

现在的问题是这里有季节性单位根吗?这意味着我们的模型应该是如下的形式:

如果我们忘记了自回归和移动平均分量,我们可以估计:

如果这里有季节单位根\alpha无论如何都应该接近1。

> arima(tsm,order=c(0,0,0),seasonal=list(order=c(1,0,0),period=12))

Call:
arima(x = tsm, order = c(0, 0, 0), seasonal = list(order = c(1, 0, 0), period = 12))

Coefficients:
        sar1  intercept
      0.9702     6.4566
s.e.  0.0071     2.1515

从上面的结果我们可以看出它的值离1并不远。但是实际上,它并不能过于接近1.如果它过于接近了,我们就会收到一条错误消息...

为了展示一些有趣的模型,让我们也来考虑一下季度温度,

> N=cbind(apply(montreal[,2:4],1,sum),apply(montreal[,5:7],1,sum),apply(montreal[,8:10],1,sum),apply(montreal[,11:13],1,sum))
> X=as.numeric(t(N))
> tsq=ts(X,start=1948,freq=4)
> persp(1948:2013,1:4,N,theta=-50,col="yellow",shade=TRUE,
+ xlab="Year",ylab="Quarter",zlab="Temperature",ticktype="detailed")

(如果需要的话,再做一次,当然目的只是为了能够写出一些方程式)

我们为什么不考虑对季度温度使用V AR(1) 模型呢?就像是下面这样:

即,

A 是一个4*4矩阵时,这个模型很容易被预测。

> library(vars)
> df=data.frame(N)
> names(df)=paste("y",1:4,sep="")
> model=VAR(df)
> model

VAR Estimation Results:
======================= 

Estimated coefficients for equation y1: 
======================================= 
Call:
y1 = y1.l1 + y2.l1 + y3.l1 + y4.l1 + const 

       y1.l1        y2.l1        y3.l1        y4.l1        const 
 -0.13943065   0.21451118   0.08921237   0.30362065 -34.74793931 

Estimated coefficients for equation y2: 
======================================= 
Call:
y2 = y1.l1 + y2.l1 + y3.l1 + y4.l1 + const 

      y1.l1       y2.l1       y3.l1       y4.l1       const 
 0.02520938  0.05288958 -0.13277377  0.05134148 40.68955266 

Estimated coefficients for equation y3: 
======================================= 
Call:
y3 = y1.l1 + y2.l1 + y3.l1 + y4.l1 + const 

      y1.l1       y2.l1       y3.l1       y4.l1       const 
 0.07740824 -0.21142726  0.11180518  0.12963931 56.81087283 

Estimated coefficients for equation y4: 
======================================= 
Call:
y4 = y1.l1 + y2.l1 + y3.l1 + y4.l1 + const 

      y1.l1       y2.l1       y3.l1       y4.l1       const 
 0.18842863 -0.31964579  0.25099508 -0.04452577  5.73228873

然后矩阵A 在这里:

> A=rbind(
+ coefficients(model$varresult$y1)[1:4],
+ coefficients(model$varresult$y2)[1:4],
+ coefficients(model$varresult$y3)[1:4],
+ coefficients(model$varresult$y4)[1:4])
> A
           y1.l1       y2.l1       y3.l1       y4.l1
[1,] -0.13943065  0.21451118  0.08921237  0.30362065
[2,]  0.02520938  0.05288958 -0.13277377  0.05134148
[3,]  0.07740824 -0.21142726  0.11180518  0.12963931
[4,]  0.18842863 -0.31964579  0.25099508 -0.04452577

由于这个多时间序列的平稳性与这个矩阵的特征值密切相关,让我们看看它们,

> eigen(A)$values
[1]  0.35834830 -0.32824657 -0.14042175  0.09105836
> Mod(eigen(A)$values)
[1] 0.35834830 0.32824657 0.14042175 0.09105836

所以在这里看起来并没有平稳性问题。Paap和Franses在这里讨论了一个被称之为PAR(1)周期性自回归模型,

当:

并且:

我们可以发现这是一个VAR(1) 模型,因为,

这个模型可以使用专门的包来进行估计(也可以看一下相关介绍,以更好地理解这里的意思)

> library(partsm)
> detcomp <- list(regular=c(0,0,0), seasonal=c(1,0), regvar=0)
> model=fit.ar.par(wts=tsq, detcomp=detcomp, type="PAR", p=1)
> PAR.MVrepr(model)
----
    Multivariate representation of a PAR model.

  Phi0:

  1.000  0.000  0.000 0
 -0.242  1.000  0.000 0
  0.000 -0.261  1.000 0
  0.000  0.000 -0.492 1

  Phi1:

 0 0 0 0.314
 0 0 0 0.000
 0 0 0 0.000
 0 0 0 0.000

  Eigen values of Gamma = Phi0^{-1} %*% Phi1:
0.01 0 0 0 

  Time varing accumulation of shocks:

 0.010 0.040 0.155 0.314
 0.002 0.010 0.037 0.076
 0.001 0.003 0.010 0.020
 0.000 0.001 0.005 0.010

在这里,特征方程为:

当满足下面的情况时,具有(季节性)单位根,

对于那些明显不符合这种情况的,可以执行卡诺瓦 - 汉森测试

> CH.test(tsq)

  ------ - ------ ----
  Canova & Hansen test
  ------ - ------ ----

  Null hypothesis: Stationarity.
  Alternative hypothesis: Unit root.
  Frequency of the tested cycles: pi/2 , pi , 

  L-statistic: 1.122  
  Lag truncation parameter: 5 

  Critical values:

  0.10 0.05 0.025 0.01
 0.846 1.01  1.16 1.35

这个想法是多项式(1-L^4)C有四个根:

因为,

如果我们回到按月份统计的数据(1-L^{12}) , 则会有十二个根,

每个根都有不同的几何解析,

在这里,我们每年可以按照1个周期(每12个月一组),每年2个周期(每6个月一组),每年3个周期(每4个月一组),每年4个周期(每3个月一组),甚至每年6个周期(每2个月一组)。这将分别取决于不同的根的参数,

如下是测试的输出,

> CH.test(tsm)

  ------ - ------ ----
  Canova & Hansen test
  ------ - ------ ----

  Null hypothesis: Stationarity.
  Alternative hypothesis: Unit root.
  Frequency of the tested cycles: pi/6 , pi/3 , pi/2 , 2pi/3 , 5pi/6 , pi , 

  L-statistic: 1.964  
  Lag truncation parameter: 20 

  Critical values:

 0.10 0.05 0.025 0.01
 2.49 2.75  2.99 3.27

看起来我们否决了季节单位根的假设。我们可以使用下面的测试程序,

> library(forecast)
> nsdiffs(tsm, test="ch")
[1] 0

当输出为“1”时,表示存在季节性单位根,反之输出为“0”则表示没有季节性单位根。很简单易懂对不对?如果我们考虑基于按月更新的周期性自回归模型,则输出为,

> model=fit.ar.par(wts=tsm, detcomp=detcomp, type="PAR", p=1)
> model
----
  PAR model of order 1 .

  y_t = alpha_{1,s}*y_{t-1} + alpha_{2,s}*y_{t-2} + ... + alpha_{p,s}*y_{t-p} + coeffs*detcomp + epsilon_t,  for s=1,2,...,12
----
  Autoregressive coefficients. 

          s=1  s=2  s=3  s=4  s=5 s=6 s=7  s=8  s=9 s=10 s=11 s=12
alpha_1s 0.15 0.05 0.07 0.33 0.11   0 0.3 0.38 0.31 0.19 0.15 0.37

因此,无论如何测试,我们总是可以否决假设存在季节性单位根。这并不意味着我们没有一个明显的循环!实际上,这个系列几乎是周期性的,但是没有单位根!所以这一切都讲得通了(我很难想象季节性序列有单位根,而温度序列没有)。

为了确保我们做到准确无误,考虑下面两个灵感来自前一个的时间序列。第一个是周期性序列(噪声非常小,为了避免非正定矩阵的问题),第二个明显是单整的。

> Xp1=Xp2=as.numeric(t(M))
> for(t in 13:length(M)){
+ Xp1[t]=Xp1[t-12]
+ Xp2[t]=Xp2[t-12]+rnorm(1,0,2)
+ }
> Xp1=Xp1+rnorm(length(Xp1),0,.02)
> tsp1=ts(Xp1,start=1948,freq=12)
> tsp2=ts(Xp2,start=1948,freq=12)
> par(mfrow=c(2,1))
> plot(tsp1)
> plot(tsp2)

同样,我们也可以把它们可视化,

> par(mfrow=c(1,2))
> bb3D(tsp1)
> bb3D(tsp2)

如果我们快速观察这些时间序列,我会说第一个没有单位根 - 即使它不是平稳的,这是因为该系列是周期性的 - 而第二个是具有单位根的。如果我们使用卡诺瓦 - 汉森测试,我们会得到,

> CH.test(tsp1)

  ------ - ------ ----
  Canova & Hansen test
  ------ - ------ ----

  Null hypothesis: Stationarity.
  Alternative hypothesis: Unit root.
  Frequency of the tested cycles: pi/6 , pi/3 , pi/2 , 2pi/3 , 5pi/6 , pi , 

  L-statistic: 2.234
  Lag truncation parameter: 20 

  Critical values:

 0.10 0.05 0.025 0.01
 2.49 2.75  2.99 3.27

> CH.test(tsp2)

  ------ - ------ ----
  Canova & Hansen test
  ------ - ------ ----

  Null hypothesis: Stationarity.
  Alternative hypothesis: Unit root.
  Frequency of the tested cycles: pi/6 , pi/3 , pi/2 , 2pi/3 , 5pi/6 , pi , 

  L-statistic: 5.448  
  Lag truncation parameter: 20 

  Critical values:

 0.10 0.05 0.025 0.01
 2.49 2.75  2.99 3.27

我知道这个软件包已经被删除了,或许我不应该使用它。而是应该改为使用,

> nsdiffs(tsp1, 12,test="ch")
[1] 0
> nsdiffs(tsp2, 12,test="ch")
[1] 1

我们得到了同样的结论。第一个没有单位根,但第二个有。不过我们要留意Osborn-Chui-Smith-Birchenhall测试

> nsdiffs(tsp1, 12,test="ocsb")
[1] 1
> nsdiffs(tsp2, 12,test="ocsb")
[1] 1

我们甚至有这样的感觉,在循环系列中也有一个单位根。。。

所以在这里,在短时间内,我们否决了单位根假设,即使是对我们温度序列中的季节性因素。或许有一天我们依然会有许多高频问题需要处理(虽然我不认为我在这个课程中有时间来介绍远距离依赖。。。)

本文的版权归 Erikdu 所有,如需转载请联系作者。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏专知

【论文推荐】最新六篇聊天机器人相关论文—弱监督信息、内容驱动、对话管理系统、可扩展情感序列到序列、自主性

2372
来自专栏生信技能树

叫一声lncRNA你敢答应么[男女不限]

lncRNA 的全称是long noncoding RNA。即又长又表达且还不能编码翻译成蛋白质的一类RNA。

2251
来自专栏null的专栏

《数学之美》拾遗——TF-IDF

开篇序     在学习机器学习的过程中,我写了简单易学的机器学习算法的专题,依然还有很多的算法会陆续写出来。网上已经有很多人分享过类似的材料,我只是通过自己的理...

31810
来自专栏机器学习人工学weekly

机器学习人工学weekly-2018/9/23

Rosetta: Understanding text in images and videos with machine learning

1085
来自专栏AI研习社

Github 项目推荐 | 中文突发事件语料库

https://github.com/shijiebei2009/CEC-Corpus

2424
来自专栏专知

【读书笔记】基于知识库的问答:生成查询图进行语义分析

【导读】将DBPedia和Freebase这样的大规模知识库组织并存储在一个结构化的数据库,这已成为支持开放领域问题问答的重要资源。 KB-QA的大多数方法基于...

5087
来自专栏林欣哲

自然语言处理--文本处理

自然语言处理的目的是让机器试图理解和处理人类的文字。通常来说,人的语言是冗余的,含有歧义的,而机器是准确的,无歧义的,要让机器理解,这之间存在一个转换的问题。 ...

3908
来自专栏iOSDevLog

机器学习研究和开发所需的组件列表

Here is a list of components that are needed for the successful machine learning...

1022
来自专栏机器学习人工学weekly

机器学习人工学weekly-2018/5/27

Prefrontal cortex as a meta-reinforcement learning system

1454
来自专栏iOSDevLog

Scikit-Learn教程:棒球分析 (一)

一个scikit-learn教程,通过将数据建模到KMeans聚类模型和线性回归模型来预测MLB每赛季的胜利。

2122

扫码关注云+社区

领取腾讯云代金券