常用的一元时间序列方法 单位根及协整检验 VARX模型与状态空间模型 多元时间序列数据 常用的一元时间序列方法 单位根及协整检验 VARX模型与状态空间模型

时间序列是指数据按照一定时间间隔收集的一系列数据 时间序列是指数据按照一定时间间隔收集的一系列数据. 时间序列可以是一维或一元的, 也就是说只有一个按照时间记录的变量, 比如一个气象站点获取的降水量数据. 也有些是多维的, 即变量是多维向量, 比如一个气象站按照同样时间间隔收集的气压、气温、降水量、风力等多个指标. 研究时间序列的一个主要目的是做对同样变量的未来值的预测. 这就意味着下面的假定必须成立: 这个未来值能够完全由同样变量的现在和过去值预测, 而不受任何其他变量的影响. 这个假定很强, 往往不能满足, 但人们往往又有意无意地无视这个假定. 当然也有加入其他变量的模型. 中国通常的时间序列教科书内容的主要部分都是讲的一元时间序列, 这是因为其数学推导和结论被研究得比较透彻, 比较容易讲. 到了高维时间序列, 一切都不那么清晰和漂亮, 一元时间序列的一些漂亮数学结论和公式很难推广到高维情况.

多元时间序列在预测上是否就比一元强些呢. 这不见得 多元时间序列在预测上是否就比一元强些呢? 这不见得. 时间序列, 特别是经济领域的时间序列受到大量其他因素的影响, 比如疾病(如SARS), 自然灾害(如地震海啸), 法律政策的改变, 全球的金融危机等等. 这些经济特征的时间序列严重受制于那些大的环境因素的变化, 而后者是几乎无法用数学方法预料的.

常用的一元时间序列方法 时间序列的组成和分解, 差分及平滑 例7.1 货币基数(AMBNS.csv). 这是在圣路易的美国联邦储备银行发布的货币基数. 该数据是月度数据, 从1918年1月1日到2012年1月1日, 单位是10亿美元. w=read.csv("AMBNS.csv") w1=ts(w[,2],start = c(1918,1), frequency = 12,) par(mfrow=c(1,2)) ts.plot(w1);lines(diff(w1),lty=3,col=4) w2=ts(w1[949:1129],start = c(1997,1), frequency = 12,) ts.plot(w2,ylim=c(0,max(w2)));lines(diff(w2),lty=3,col=4)

Holt-Winters滤波函数(可做指数平滑) b2 <- HoltWinters(w2, gamma = FALSE, beta = FALSE) par(mfrow=c(1,2)) plot(b2) #原序列及拟合曲线图 ee2=b2$x-b2$fit[,1] plot(ee2) #拟合误差图

Loess方法对橙汁的冰冻橙汁厂家价格做季节趋势分解 data(Orange,package="Ecdat") ts.plot(Orange[,1]) #原始变量点图 a=stl(Orange[,1], "period") #对第一个变量(Orange[,1])分解 ts.plot(a$time.series[,1:3])#画出分解出来的三部分

ARIMA模型 例7.1的货币基数数据, 截取货币基数的一部分数据(从1918年1月到2007年12月), 试试ARIMA(0,1,8)(1,0,0)模型. w0=ts(w1[1:1080],start = c(1918,1), frequency = 12,)#到2007-12-01 854.997 b0=arima(w0,order=c(0,1,6),seasonal = list(order = c(1, 0, 0))) bp72=predict(b0,72);plot(w2);lines(bp72$pre,lty=2);tsdiag(b0)#图7.4

Ljung-Box检验的零假设为序列独立(对于某个滞后) B=NULL;for( i in 1:30) B=c(B,Box.test(b0$resi, lag = i, type = "Ljung-Box")$p.value) plot(B,main="Ljung-Box tests", ylab="p-value", xlab="lag",pch=16) abline(h=.05,lty=2)

7.3 尼罗河(Nile.txt). 这是是在阿斯旺(Ashwan)所测量的1871–1970年尼罗河的年度流量 利用下面语句来点出尼罗河数据的acf和pacf函数图(图7.6) data(Nile,package="datasets") layout(matrix(c(1,1,2,3),2,2,byrow=TRUE)) plot(Nile); acf(Nile) ;pacf(Nile) ; ar(Nile) # 选择了二阶AR模型 rn=arima(Nile, c(2, 0, 0)) 函数自动根据AIC选择AR模型的阶数

B=c(B,Box.test(rn$resi, lag = i, type = "Ljung-Box")$p.value) 打印出来的结果为: Call: arima(x = Nile, order = c(2, 0, 0)) Coefficients: ar1 ar2 intercept 0.4096 0.1987 919.8397 s.e. 0.0974 0.0990 35.6410 sigma^2 estimated as 20291: log likelihood=-637.98,aic=1283.96 用下面语句点出Ljung-Box检验的p值(虚线为0.05水平线), 残差的acf和pacf函数图(图7.7), 看来拟合虽然不是那么完美, 但也还过得去. B=NULL;for( i in 1:30) B=c(B,Box.test(rn$resi, lag = i, type = "Ljung-Box")$p.value) layout(matrix(c(1,1,2,3),2,2,byrow=TRUE)) plot(B,main="Ljung-Box tests", ylab="p-value", xlab="lag",pch=16,ylim=c(0,1)) abline(h=.05,lty=2) acf(rn$res) ;pacf(rn$res)

结构时间序列模型/一元状态空间模型 测量方程 动态方程(状态方程) data(Nile,package="datasets") (fit <- StructTS(Nile, type = "level")) plot(Nile) lines(fitted(fit), lty = 2) # 同期平滑得到的状态 lines(tsSmooth(fit), lty = 3, col = 4) # 固定区间平滑得到的状态 legend("bottomleft", lty=1:3, c("Nile", "contemporaneous smoothing","fixed-interval smoothing")) par(mfrow=c(1,2));acf(residuals(fit));pacf(residuals(fit))

单位根及协整检验 本节的方法在计量经济学中是比较精彩的方法. 很多人都在使用. 但是, 这些数学公式都是人们用有限的数学语言对现实问题的一种近似. 首先, 这里所有的方法都要求时间序列具有线性表示, 甚至如Johansen方法中的更强的VAR表示, 没有这种假定, 所有的方法, 所有的渐近性质和结论都毫无意义. 想想看, 宇宙间有多少规律是严格线性的? 人们之所以处处使用线性模型, 或用线性模型近似, 仅仅是因为线性关系是我们或多或少可以控制的少数几个关系之一. 其次, 所有的方法的结论都严重依赖于大样本性质, 这是因为独立同正态分布的条件几乎不可能满足, 也不可能验证. 而没有任何统计学家(无论他们证明了多少大样本定理)可以断定在实际中样本量多少才算是大样本. 而在小样本情况下,可能所有这些结论根本不成立, 显著性检验也没意义.

单位根及协整检验 一些时间序列可能是相关的, 他们在一定时期可能显示出较强的相关性. 但是, 用它们自己来预测自己的未来或者互相预测, 则风险很大. 就拿宏观经济时间序列来说, 它们都受到经济危机、金融危机、政策变化、政权更迭、天灾人祸疾病等偶然事件、法律和规则的改变以及人们心态改变的影响, 而这些无法预料的因素根本无法加入到这几个只包含了若干可测量的时间序列的模型之中, 某些短期预测往往不用时间序列也可以猜出来, 而长期预测则往往是极端不可靠的. 单纯利用时间序列来做出推断, 而且还给出``政策建议'', 不但是不科学的, 也是很不负责任的. 我们还会遇到不同方法得到的结论不相同的问题. 这么多检验方法(这里只介绍几个)产生出不同的结论, 原因在于这些方法的数学模型不同, 这些模型及假定与事实不符. 不仅仅不同方法会造成不同结论, 同一个方法的不同角度也会造成不同的结论. 比如Engle-Granger检验, 很可能轮换回归得到的结论就不同, 不仅如此, 该方法有两个步骤, 第二个步骤肯定继承了第一个步骤的误差, 最终误差可能达到不可承受的地步.

单位根及协整检验 此外, 在实际工作者使用这些方法时, 往往把``不能拒绝零假设‘’当成``零假设正确‘’或``接受零假设‘’的同义词, 这在在逻辑上是完全错误的, 在统计上是绝对不允许的. ``不能拒绝零假设‘’意味着拒绝的证据不足, 并不意味着零假设正确. 在拒绝零假设时, 至少给出了犯错误的大致概率, 即p值, 但任何人都给不出``接受零假设‘’时犯错误的概率, 做出结论而又不给出该结论所包含的风险是极端不负责的行为. 但实际工作者往往需要在显著性检验中无法拒绝零假设时做出决策, 合乎逻辑的说法是: “在所有关于模型形式(比如某种线性表示以及对数据的各种假定都成立的前提下, 我没有足够证据否定零假设, 因此把零假设当成另一个附加的假定.”

单位根及检验 如果一个时间序列是平稳的, 则没有任何预测价值, 因为平稳序列的均值不变, 任何预测都不会有什么有价值的结果. 因此人们只对非平稳序列感兴趣, 而对平稳序列的研究, 也是因为对于非平稳序列, 总希望可以通过差分等方式转换成平稳的, 而后者是可以通过数学方式予以解释的. 如果由$d$次差分可以将一个非平稳序列转换成平稳的, 则称其为$d$阶单整的, 记为$I(d)$. 换句话说, 如果$\Delta^d X_t$为平稳的, 则序列$X_t$称为有$d$个单位根. 最简单的情况为(这里带有截距)随机游走 随机游走为$I(1)$单整的, 它有1个单位根.

例7.2 芬兰数据(finland.csv). 该数据来自Johansen and Juselius(1990), 有4个时间序列变量, 它们是从1958年第2季度到1984年第3季度的货币供应量M1的对数(lrm1), 实际收入的对数(lny), 边际利率(lnmr), 通货膨胀率(difp).

lrm1.df=ur.df(lrm1,lags=5,type='trend');summary(lrm1.df) 下面就是对于例7.2数据的ADF检验, 看其4个变量有没有单位根. library(urca) data(finland) attach(finland) lrm1.df=ur.df(lrm1,lags=5,type='trend');summary(lrm1.df) lny.df=ur.df(lny,lags=5,type='trend');summary(lny.df) lnmr.df=ur.df(lnmr,lags=5,type='trend');summary(lnmr.df) difp.df=ur.df(difp,lags=5,type='trend');summary(difp.df) ADF显著性水平为0.1的临界值为-3.13(值越小就越显著), 而这4个检验统计量的值分别为-2.2702, -1.5932, -3.1016, -2.7569, 因此没有足够证据拒绝零假设, 也就是说, 没有证据说它们是平稳的. 因此假定它们都至少是$I(1)$的.

lrm12=diff(lrm1);lrm12.df=ur.df(lrm12,lags=5,type='trend') 为了确定这些系列的单整阶数, 再对它们的差分数做ADF单位根检验, 零假设是差分序列有单位根. 下面是代码: lrm12=diff(lrm1);lrm12.df=ur.df(lrm12,lags=5,type='trend') summary(lrm12.df)#.01显著I(1) lny2=diff(lny); lny2.df=ur.df( lny2,lags=5,type='trend') summary( lny2.df)#.05显著I(1) lnmr2=diff(lnmr);lnmr2.df=ur.df(lnmr2,lags=5,type='trend') summary(lnmr2.df)#.01显著I(1) difp2=diff(difp);difp2.df=ur.df(difp2,lags=5,type='trend') summary(difp2.df)#.01显著I(1) 这次, ADF显著性水平为0.01, 0.05及0.1的临界值分别为-3.99, -3.45, -3.13(值越小就越显著), 而上述差分的4个检验统计量分别为 -4.89, -3.7416, -6.4769, -5.7403, 分别在0.01, 0.05, 0.01, 0.05水平上显著, 这就拒绝了差分有单位根的零假设, 即他们都不大会是$I(2),$ 有可能是$I(1)$单整的.

协整检验 如果若干时间序的线性组合的单整系数阶数小于其成分的单整阶数, 就称这些序列间存在协整(cointegration). 那些线性组合系数称为协整向量(cointegrating vector), 协整向量个数小于变量数目. 协整概念类似于线性代数中的线性相关的概念. 一个时间序列向量, 只有当它们是协整的, 才有同时研究的价值. 协整意味着向量分量之间存在长期关系. 在短期中, 这些向量可能关系不那么显著, 各自由不同的动态过程所支配, 然而, 长远来说, 协整把变量绑在一起.

对于多元时间序列 带有输入项(哑元$D_t$)的向量自回归模型$VAR(p)$为

Egle-Granger协整检验 这个检验的思想很简单: 先用这些时间序列变量互相做通常最小二乘(OLS)回归, 再通过单位根检验它们的残差是否为$I(0)$, 如果是, 则这些变量可能存在协整关系. 我们通过例7.2来说明. 首先对这些系列互相做OLS回归: lrm1=ts(lrm1,start=c(1958,2),end=c(1984,3),frequency=4) lny=ts( lny,start=c(1958,2),end=c(1984,3),frequency=4) lnmr=ts(lnmr,start=c(1958,2),end=c(1984,3),frequency=4) difp=ts(difp,start=c(1958,2),end=c(1984,3),frequency=4) flcons=window(cbind(lrm1,lny,lnmr,difp),start=c(1958,2), end=c(1984,3)) #各个变量轮流做因变量: lrm1.eq=summary(lm(lrm1~lny+lnmr+difp,data=flcons));lrm1.eq lny.eq=summary(lm(lny~lrm1+lnmr+difp,data=flcons));lny.eq lnmr.eq=summary(lm(lnmr~lrm1+lny+difp,data=flcons));lnmr.eq difp.eq=summary(lm(difp~lrm1+lny+lnmr,data=flcons));difp.eq #下面是残差序列: error.lrm1=ts(resid(lrm1.eq),start=c(1958,3), end=c(1984,3),frequency=4) error.lny=ts(resid( lny.eq),start=c(1958,3), error.lnmr=ts(resid(lnmr.eq),start=c(1958,3), error.difp=ts(resid(difp.eq),start=c(1958,3),

df.lrm1=ur.df(error.lrm1,lags=0,type='none');summary(df.lrm1) 通过这些回归的输出, 可以看出有些变量之间是有些关系, 但由于OLS的检验条件不一定符合, 不足为据. 必须对残差做进一步研究. 下面对这些残差做ADF单位根检验, 这里零假设是残差序列存在单位根, 即它们之间不存在协整. df.lrm1=ur.df(error.lrm1,lags=0,type='none');summary(df.lrm1) df.lny=ur.df( error.lny,lags=0,type='none'); summary( df.lny) df.lnmr=ur.df(error.lnmr,lags=0,type='none');summary(df.lnmr) df.difp=ur.df(error.difp,lags=0,type='none');summary(df.difp) 这里ADF检验的显著性为0.01, 0.05, 0.1的临界值分别为-2.58 -1.95 -1.62, 而这四个残差的检验统计量的值为-3.5256, -3.904, -5.1339, -8.7288, 都在0.01的水平上显著, 这就是说, 可以拒绝零假设, 而得出结论说残差不存在单位根, 为$I(0)$序列, 因而这些序列之间很可能存在协整关系.

有了协整关系, 我们希望得到短期的误差修正模型, 对于例7.2, 它可以是这样的(非矩阵)形式: 这些系数可以用简单线性模型估计: lrm12=diff(lrm1) lny2=diff( lny) lnmr2=diff(lnmr) difp2=diff(difp) leq2=lag(error.lrm1) ecm=summary(lm(lrm12~lny2+lnmr2+difp2+leq2))#OK ecm$coef[,1] #得到ECM系数 当然, 这些估计是比较粗糙的. 但可以看出, 最后一项, 即误差项的系数为负数, 这就是``误差修正''的含义, 描述了长期效应收敛的快慢.

Pillips-Ouliaris协整检验

flcons=window(cbind(lrm1,lny,lnmr,difp), 对于例7.2数据, 方差率检验为 flcons=window(cbind(lrm1,lny,lnmr,difp), start=c(1958,2), end=c(1984,3)) pu.test=summary(ca.po(flcons,demean='const',type='Pu')) pu.test #.1不显著 统计量的值为26.0919, 而显著性水平为0.1的临界值为39.6949, 因此不显著. 多元迹检验代码为: pz.test=summary(ca.po(flcons,demean='const',type='Pz')) pz.test #.01显著 统计量的值为184.3248, 而显著性水平为0.01的临界值为153.4504, 因此很显著. 按照多元迹检验, 可能存在协整, 而根据方差率检验, 没有证据说存在协整. 这就产生了矛盾结论.

Johansen 方法

summary(ca.jo(data.frame(lrm1,lny,lnmr,difp), 对于例7.2, 使用Johansen方法的迹检验的代码如下: summary(ca.jo(data.frame(lrm1,lny,lnmr,difp), type="trace",ecdet="const")) #迹检验 程序输出下面检验结果: Values of teststatistic and critical values of test: test 10pct 5pct 1pct r <= 3 | 7.76 7.52 9.24 12.97 r <= 2 | 18.90 17.85 19.96 24.60 r <= 1 | 51.67 32.00 34.91 41.07 r = 0 | 95.78 49.65 53.12 60.16 按照这个结果, 有$r=3$个协整向量.

summary(ca.jo(data.frame(lrm1,lny,lnmr,difp), 对于例7.2, 使用Johansen方法的特征值检验的代码如下: summary(ca.jo(data.frame(lrm1,lny,lnmr,difp), type="eigen",ecdet="const"))#特征值检验 程序输出下面检验结果: Values of teststatistic and critical values of test: test 10pct 5pct 1pct r <= 3 | 7.76 7.52 9.24 12.97 r <= 2 | 11.15 13.75 15.67 20.20 r <= 1 | 32.77 19.77 22.00 26.81 r = 0 | 44.11 25.56 28.14 33.24 按照这个结果, 有可能有$r=2$个协整向量.

VARX模型与状态空间模型 前面一节讨论了单位根检验及协整等问题. 其主要目的是要确定这些时间序列之间存在长期的关联关系, 但究竟是怎样的关系, 则需要用一些模型来拟合数据. 本章前面已经介绍了带有输入的向量自回归模型(即VARX模型), 本节还要简单介绍状态空间(state space)模型, 然后用这两个模型拟合例7.2的数据.

VARX模型拟合 data(finland,package="urca") #finlandêâ VAR(finland, p = 2, type = "none") VAR(finland, p = 2, type = "const") VAR(finland, p = 2, type = "trend") VAR(finland, p = 2, type = "both")

VARX模型拟合

data(finland,package="urca") #finland数据 attach(finland) library(dse) #把第一个变量(M1)作为输出, 另外三个作为输入: fld <- TSdata(input= finland[,2:4],output= finland[, 1]) #把数据标为时间序列: fld <-tframed(fld,list(start=c(1958,2), frequency=4)) #把变量名标上,似乎无法取代默认名字: Series 1, Series 2等: seriesNamesInput(fld) <- c("lny","lnmr","difp") seriesNamesOutput(fld) <- "lrm1“ #拟合VARX模型: fld.ls <- estVARXls(fld,max.lag=2)#默认max.lag=6 print(fld.ls) stability(fld.ls) rr=checkResiduals(fld.ls) par(mfrow=c(1,2));acf(rr$re);pacf(rr$re)

输出的拟合矩阵为 neg. log likelihood= -81.20899 A(L) = 1-0.733L1-0.173L2 B(L) = 1 C(L) = 0.487-0.410L1 0.245-0.329L1 -0.770+0.091L1 这意味着拟合的(7.6)为 或者

状态空间模型拟合

fld.ss <- estSSfromVARX(fld,max.lag=3) print(fld.ss) stability(fld.ss) rs=checkResiduals(fld.ss) par(mfrow=c(1,2));acf(rr$re);pacf(rr$re) 输出的拟合矩阵为 {\begin{verbatim} neg. log likelihood= -58.54483 F = [,1] [,2] [,3] [1,] 0 0 -0.0678 [2,] 1 0 0.2391 [3,] 0 1 0.7257 G = [,1] [,2] [,3] [1,] 0.0721 -0.162 0.3296 [2,] -0.4238 -0.191 0.0517 [3,] 0.4371 0.204 -0.9037 H = [,1] [,2] [,3] [1,] 0 0 1 K = [,1] [1,] -0.0678 [2,] 0.2391 [3,] 0.7257 \end{verbatim}}

模型的比较 fld.bb <- estBlackBox(fld, max.lag=6)#?estBlackBox 到底选择哪个模型来拟合我们的数据呢? 软件包dse有一个自动在各种VARX模型及状态空间模型中选择最好模型的函数estBlackBox() (黑匣子). 下面来运行一下. neg. log likelihood= -81.20899 F = [,1] [,2] [1,] 0.542 -0.682 [2,] -0.406 0.192 G = [,1] [,2] [,3] [1,] -0.471 -0.253 0.6713 [2,] -0.456 -0.396 -0.0366 H = [,1] [,2] [1,] -1.14 0.11 K = [,1] [1,] -0.605 [2,] 0.387 fld.bb <- estBlackBox(fld, max.lag=6)#?estBlackBox print(fld.bb) stability(fld.bb) rb=checkResiduals(fld.bb) par(mfrow=c(1,2));acf(rb$re);pacf(rr$re)

对三个模型的比较, 可以用下面语句: informationTests(fld.ls, fld.ss,fld.bb) tfplot(fld.ls, fld.ss,fld.bb,start=c(1960,1))