季节性单位根
正如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我们甚至有这样的感觉，在循环系列中也有一个单位根。。。
所以在这里，在短时间内，我们否决了单位根假设，即使是对我们温度序列中的季节性因素。或许有一天我们依然会有许多高频问题需要处理（虽然我不认为我在这个课程中有时间来介绍远距离依赖。。。）