北京大学金融时间序列分析讲义第25章： 协整分析与向量误差修正模型

set.seed(101)
n <- 500
spureg01 <- function(n=500){
  x <- seq(n)
  eps <- rnorm(n)*50
  y <- 2 * x + eps
  lmr <- lm(y ~ x)
  summary(lmr) |> print()
  plot(x, y)
  abline(lmr, col="red")
  
  invisible(list(x=x, y=y))
}
spureg01(n)

例25.2 考虑回归误差项为AR(1)序列的情形。

set.seed(101)
spureg02 <- function(n=500){
  x <- seq(n)
  eps <- arima.sim(
    model=list(ar=c(0.7)), n=n)*36
  y <- 2 * x + eps
  lmr <- lm(y ~ x)
  summary(lmr) |> print()
  plot(x, y)
  abline(lmr, col="red")
  
  invisible(list(x=x, y=y))
}
sim2 <- spureg02(n)

arima(sim2$y, order=c(1,0,0), xreg = sim2$x)

例25.3 考虑回归误差项为I(1)序列的情形。

set.seed(101)
spureg03 <- function(n=500){
  x <- seq(n)
  eps <- cumsum(rnorm(n))*8
  y <- 2 * x + eps
  lmr <- lm(y ~ x)
  summary(lmr) |> print()
  plot(x, y)
  abline(lmr, col="red")
  
  invisible(list(x=x, y=y))
}
spureg03(n)

例25.4 如果xt和yt都是I(1)序列但是独立，
回归结果也可能很显著。

spureg04 <- function(){
  set.seed(110)
  n <- 500
  x <- cumsum(rnorm(n))
  y <- cumsum(rnorm(n))*2
  lmr <- lm(y ~ x)
  summary(lmr) |> print()
  plot(x, y)
  abline(lmr, col="red")
  
  invisible(list(x=x, y=y, lmr=lmr))
}
sim4 <- spureg04()

forecast::Acf(residuals(sim4$lmr))

fUnitRoots::adfTest(residuals(sim4$lmr), lags=1, type="c")

例25.5 如果xt和yt都是某个I(1)序列的线性变换，
则回归结果是有一定意义的，
称这两个序列有“协整关系”。

set.seed(101)
spureg05 <- function(n=500){
  z <- cumsum(rnorm(n))
  x <- z + rnorm(n)*2
  y <- 2*z + rnorm(n)*2
  lmr <- lm(y ~ x)
  summary(lmr) |> print()
  plot(x, y)
  abline(lmr, col="red")
  
  invisible(list(x=x, y=y, lmr=lmr))
}
sim5 <- spureg05(n)

forecast::Acf(residuals(sim5$lmr))

fUnitRoots::adfTest(residuals(sim5$lmr), lags=1, type="c")

library(urca, quietly = TRUE)
data(Raotbl3)
ts.Raotbl3 <- ts(Raotbl3[,1:3], start=c(1966,4), frequency=4)

library(quantmod, quietly = TRUE)
plot(as.xts(ts.Raotbl3), type="l", 
     multi.panel=FALSE, theme="white",
     main="英国消费、收入、财富季度数据",
     major.ticks="years",
     grid.ticks.on = "years")

图25.1: 英国消费、收入、财富季度数据

fUnitRoots::adfTest(Raotbl3$lc, lags=1, type="ct")

fUnitRoots::adfTest(Raotbl3$li, lags=1, type="ct")

fUnitRoots::adfTest(Raotbl3$lw, lags=1, type="ct")

fUnitRoots::adfTest(diff(Raotbl3$lc), lags=1, type="c")

fUnitRoots::adfTest(diff(Raotbl3$li), lags=1, type="c")

fUnitRoots::adfTest(diff(Raotbl3$lw), lags=1, type="c")

eglm1 <- lm(lc ~ li, data=Raotbl3); summary(eglm1)

fUnitRoots::adfTest(residuals(eglm1), lags=1, type="c")

library(urca, quietly = TRUE)
summary(ca.po(Raotbl3[-(1:2),c("lc", "li")], type="Pz", demean="constant"))

library(dse, quietly = TRUE)
#自回归部分系数矩阵
Apoly <- array(c(
  1,-0.5,0,
  0,-0.4,-0.25,
  0,-0.1,0,
  1,-0.5,0),
  c(3,2,2))
Apoly[2,,]
#滑动平均部分系数矩阵，此处有变换
B <- matrix(c(3,0,0,2),ncol=2)*0.01
#生成模型
var2 <- ARMA(A=Apoly, B=B)
print(var2)
#模拟数据
varsim <- simulate(
  var2, sampleT=250, 
  noise=list(w=matrix(rnorm(500),nrow=250,ncol=2)),
  rng=list(seed=c(123456)))
vardat <- matrix(varsim$output, nrow=250, ncol=2)

25.4.1 VARMA模型协整关系

这里将VAR和VARMA模型略推广，
不要求模型一定满足平稳性条件。
如果不然，
因为多元弱平稳其分量一定也是一元弱平稳的，
就不能容许分量存在单位根。

考虑如下的二元VARMA(1,1)模型的例子。

=(x1tx2t)−(0.5−0.25−1.00.5)(x1,t−1×2,t−1)(a1ta2t)+(0.2−0.1−0.40.2)(a1,t−1a2,t−1)(25.2)

其中的扰动at的协方差阵Σ为正定阵。

这个模型不是弱平稳的，
因为|P(z)|=1−z有一个单位根。

将模型改写为

(1−0.5B0.25BB1−0.5B)(x1tx2t)=(1−0.2B0.1B0.4B1−0.2B)(a1ta2t)

利用伴随矩阵，
因为A⋅adj(A)=|A|I，
模型两边都左乘P(B)的伴随矩阵

(1−0.5B−0.25B−B1−0.5B)

而|P(B)|=1−B，
模型转化成

(1−B001−B)(x1tx2t)=(1−0.7B−0.15B−0.6B1−0.7B)(a1ta2t)

即

x1t−x1,t−1=x2t−x2,t−1=a1t−0.7a1,t−1−0.6a2,t−1a2t−0.7a2,t−1−0.15a1,t−1

这两个模型左边是一阶差分，
右边是MA(1)序列，
所以x1t和x2t都是一元单位根过程。
这两个分量并不独立，
因为模型右边包含共同的扰动成分。

设法对两个分量进行线性变换，
使得单位根只出现在一个分量中。
令

L=yt=bt=(1.00.5−2.01.0)(y1ty2t)=L(x1tx2t)Lat

原来的xt模型可以写成

xt=Φxt−1+at+Θat−1

两边左乘L得到关于yt的模型：

yt===LxtLΦxt−1+Lat+LΘat−1LΦL−1yt−1+bt+LΘL−1bt−1

计算其中的系数矩阵，得到xt线性变换得到的yt序列的模型为

=(y1ty2t)−(1000)(y1,t−1y2,t−1)(b1tb2t)+(−0.4000)(b1,t−1b2,t−1)

即

y1t−y1,t−1=y2t=b1t−0.4b1,t−1b2t

其中

y1t=y2t=x1t−2x2t0.5x1t+x2t

y1t仍是一元单位根过程，
但y2t已经是弱平稳列，
所以x1t和x2t之间存在协整关系。

一般地，
对于k元时间序列xt建立了VAR或者VARMA模型后，
如果分量都是一元单位根过程，
但AR部分的特征多项式|P(z)|的根1的重数小于k，
分量间就存在协整关系。
如果特征多项式的根1的重数为h<k，
称k−h为协整因子个数，
这是能找到的线性组合后平稳的线性独立的线性组合系数的个数。
这种线性组合称为协整向量。
在上面的例子中因为y2t平稳所以(0.5,1)T是一个协整向量。

25.5.1 误差修正模型

因为在协整系统中，
单位根非平稳分量的个数多于单位根的个数
（通过线性组合可以使得单位根非平稳的分量减少），
所以如果对每个单位根非平稳分量计算差分，
虽然使得分量都平稳了，
但是会造成过度差分，
使得部分分量的ARMA模型的MA部分有单位根，
这样的模型平稳但不可逆，
不可逆的模型在估计和预测上比较困难。

(Engle and Granger 1987)讨论了协整系统的误差修正表示，
称为向量误差修正模型(VECM, vector error correction model)。

以(25.2)的协整系统为例。
令Δxt=xt−xt−1，
模型可以化为

Δxt===xt−xt−1=(Δx1tΔx2t)(−0.5−0.25−1−0.5)xt−1+at+(−0.20.10.4−0.2)at−1(−1−0.5)(0.51.0)xt−1+at+(−0.20.10.4−0.2)at−1

注意到y2t=(0.5,1)xt−1是平稳的，
所以上式两边都是平稳的，
MA部分没有出现单位根，不引入不可逆性。
上式描述了平稳序列。
这样的表示称为VECM。

此表示可以推广到一般的VARMA模型。
对(25.1)的VARMA模型，
如果其中有m个协整因子，m<k，
单位根的个数大于m，
必存在如下的误差修正形式(VECM):

Δxt=αβTxt−1+∑j=1p−1Φ∗jΔxt−j+at+∑j=1qΘjat−j

其中α和β都是k×m列满秩矩阵，
MA部分没有单位根，
m维时间序列yt=βTxt是平稳列（没有单位根），
β的每一列都是xt的一个协整系数。
Φ∗j和α, β都依赖于原来的AR部分的系数矩阵
Φj，关系为：

Φ∗j=αβT=−∑i=j+1pΦi, j=1,2,…,p−1Φp+⋯+Φ1−I=−P(1)

可以对−P(1)作奇异值分解得到α和β，
α和β并不唯一。
当模型参数从样本估计时，
不会有严格的单位根，
也不能通过矩阵秩得到−P(1)=αβT这样的分解，
需要进行检验。

25.5.2 VAR模型的误差修正形式

为简单起见仅考虑VAR模型（不要求平稳性条件）。
VAR模型的估计与解释都更容易。

考虑xt的如下的带有非随机线性趋势项的k维VAR(p)模型

xt=μt+Φ1xt−1+⋯+Φpxt−p+at(25.3)

其中at∼Nk(0,Σ)，
μt=μ0+μ1t。

记

P(z)=I−Φ1z−⋯−Φpzp

当用行列式|P(z)|表示z的多项式的所有根都在复平面的单位圆外时{xt}
是弱平稳的（简称平稳的），
在讨论单位根和协整问题时记弱平稳列为I(0)序列。

为简单起见，
设xt的所有分量至多需要一次差分就可以变平稳，
即至多为I(1)序列。
这就是说，
如果分量xit不平稳，
Δxit=xit−xi,t−1就一定是平稳的。
这也是实际的大部分序列都能满足的。

上面的VAR(p)模型有如下的误差修正形式（VECM）

Δxt=μt+Πxt−1+Φ∗1Δxt−1+⋯+Φ∗p−1Δxt−p+1+at(25.4)

其中

Φ∗j=Π=−∑i=j+1pΦi, j=1,2,…,p−1αβT=−P(1)

α和β是满秩的k×m矩阵，
m是协整向量个数。
容易看出Φj可以从VECM的参数中递推地反解为

Φ1=Φj=I+Π+Φ∗1Φ∗j−Φ∗j−1, j=2,…,p, Φ∗p=0

因为假定每个分量至多I(1)
所以VECM式(25.4)中的{Δxt}序列是I(0)序列。

如果xt有单位根，
则|P(1)|=0，
从而Π=−P(1)不满秩。
根据Π的秩，
将(25.4)的VECM分为如下三种情况：

情况1：
rank(Π)=0。
这时Π=0，
xt有k个单位根，
不存在协整关系，
xt不是协整的，
(25.4)的VECM变成了如下的关于Δxt的VAR(p−1)模型

Δxt=μt+Φ∗1Δxt−1+⋯+Φ∗p−1Δxt−p+1+at

只要关于差分序列建模即可。

情况2：
rank(Π)=k满秩。
这时|P(1)|≠0，
没有单位根，
xt本身是平稳列，
只要直接关于xt建立平稳的VAR模型即可。

情况3：
0<m=rank(Π)<k，
这时xt存在协整关系。
下面主要考虑情况3。

Π可以分解为

Πk×k=αβT

其中α和β都是k×m的列满秩矩阵，
β的每一列都是xt的一个协整向量，
m元时间序列ωt=βTxt是I(0)序列。
α称为载荷矩阵(loadings matrix)。
xt有k−m个单位根，
这些单位根给出了xt的k−m个共同趋势，
相当于自由变化的趋势；
而m个协整分量ωt，
则可以看成是xt的k个分量的m个线性约束，
或长期线性依赖关系。
本身就是I(0)的分量可以计入这m个当中。
(25.4)的VECM可以写成

Δxt=μt+αβTxt−1+Φ∗1Δxt−1+⋯+Φ∗p−1Δxt−p+1+at

记m=rank(Π)。
在0<m<k时，
得到这k−m个公共趋势的方法是先计算α
的正交补矩阵α⊥，
α⊥是k×(k−m)矩阵，
(α,α⊥)为满秩k阶方阵，
且αT⊥α=0。
令yt=αT⊥xt,
则yt是m−k元的I(1)序列，
yt不再有协整关系，
在(25.4)的VECM两边同时左乘αT⊥可以去掉
Πxt−1项，
变成了

Δyt=αT⊥μt+αT⊥Φ∗1Δxt−1+⋯+αT⊥Φ∗p−1Δxt−p+1+αT⊥at

将Δyt表示成了右侧的k−m元平稳列的形式。

容易看出Π=αβT的分解是不唯一的，
对任意m阶正交阵Ω，
都有

αβT=αΩΩTβT=(αΩ)(βΩ)T=α∗βT∗

其中α∗和β∗也都是k×m列满秩矩阵。
为此增加约束条件，
要求

β=(Imβ1)

β1为(k−m)×m矩阵。
将xt分解为

xt=(utvt)

ut为m元，
vt为k−m元，
则
βTxt=ut+β1vt。
实际中这可能需要将xt的各分量重新排序使得有单位根的分量（至少有m个）都排在前面，
将没有单位根的分量都排在后面。

为了使得ωt=βTxt平稳，
还需要其它的参数条件，
就像没有协整的I(0)的VAR(1)需要平稳性条件那样。
比如，
考虑有一个协整向量的二元VAR(1)模型，
这里k=2，
m=1，
两个分量都为单位根过程，
VECM为

Δxt=μt+(α1α2)(1β1)xt−1+at

两边同时左乘βT，
注意ωt=βTxt
则Δωt=βTΔxt，所以

Δωt=βTμt+(α1+α2β1)ωt−1+βTat

记bt=βTat，
则{bt}为白噪声列，
将方程左边的ωt−1移项到右边，得

ωt=βTμt+(1+α1+α2β1)ωt−1+bt

为使得模型平稳需要α和β满足

|1+α1+α2β1|<1.

因为协整的条件是0<rank(Π)<k，
所以从数据中估计得到VECM模型后对协整的检验就是对rank(Π)的检验。
(Johansen 1988),
(Johansen 1995),
Reinsel和Ahn(1992)正是这样进行协整检验的。

25.5.3 关于VECM中的非随机趋势

为了基于VECM模型进行协整检验，
需要规定μt的形式，
μt的形式会影响到检验统计量的极限分布。

设xt的所有分量都为I(1)序列。
文献中已有的μt的形式包括：

（1）
μt=0。
这时xt的分量为不带漂移的I(1)过程，
平稳序列ωt=βTxt
均值为零。

（2）
μt=μ0=αc0，
即μt是有特殊形式的非零常数。
VECM变成

Δxt=α(βTxt−1+c0)+Φ∗1Δxt−1+⋯+Φ∗p−1Δxt−p+1+at

这时xt的分量是不带漂移的I(1)过程，
但是平稳列ωt的均值为非零的−c0。

（3）
μt=μ0≠0，
即μt是不带特殊结构的非零常数。
这时xt的分量是带有漂移项μ0的I(1)序列，
且ωt非零均值。

（4）
μt=μ0+αc1t，
即μt是线性项t的系数有特殊结构的线性函数，
VECM变成

Δxt=μ0+α(βTxt−1+c1t)+Φ∗1Δxt−1+⋯+Φ∗p−1Δxt−p+1+at

这时xt的分量是带有漂移项μ0的I(1)过程，
且ωt有非随机的c1t线性趋势项。

（5）
μt=μ0+μ1t，
μ0, μ1非零。
这是没有约束的线性趋势。
这时xt的分量都是带有二次时间趋势的I(1)过程，
ωt是线性趋势项加平稳列。

最后一种情况不常见；
第一种情况在经济序列中也不常见。
第三种情况在许多对数价格序列中常见。

25.5.4 VECM的最大似然估计

设xt有观测数据
{xt,t=1,2,…,T}。
记μt=μdt，
其中dt=(1,t)T。
设Π的秩为m，
0<m<k。
VECM为

Δxt=μdt+αβTxt−1+Φ∗1Δxt−1+⋯+Φ∗p−1Δxt−p+1+at(25.5)

为了求各个参数的最大似然估计，
首先求解如下的多对多的线性回归问题

Δxt=xt−1=γ0dt+Ω1Δxt−1+⋯+Ωp−1Δxt−p+1+utγ1dt+Ξ1Δxt−1+⋯+Ξp−1Δxt−p+1+vt

令û t和v̂ t为最小二乘估计的残差，
定义如下的样本协方差阵

S00=S01=S11=1T−p∑t=p+1Tû tû Tt,1T−p∑t=p+1Tû tv̂ Tt,S10=ST01,1T−p∑t=p+1Tv̂ tv̂ Tt

然后，计算矩阵
S10S−100S01关于S11的广义特征值和特征向量，
即求λ和ξ使得

S10S−100S01ξ=λS11ξ

设k个特征值、特征向量为(λ̂ i,ei)，
其中λ̂ 1>λ̂ 2>⋯>λ̂ k，
记e=(e1,e2,…,ek)，
标准化使得
eTS11e=Ik。

协整向量组成的矩阵β的未标准化的最大似然估计为
β̂ =(e1,…,em)。
由此可以得到β的满足可辨识条件和标准化条件的MLE，
将得到的估计记为β̂ c，
其中下标c表示需要满足一定的限制条件。

其它参数的MLE可以从下式用最小二乘方法得到：

Δxt=μdt+αβ̂ cxt−1+Φ∗1Δxt−1+⋯+Φ∗p−1Δxt−p+1+at

基于m个协整向量的似然函数的最大值满足

L−2/Tmax∝|S00|∏i=1m(1−λ̂ i)

在计算关于rank(Π)=m的似然比统计量时需要用到这个似然函数最大值。

α和β的正交补可以用如下公式计算：

α̂ ⊥=β̂ ⊥=S−100S11(em+1,…,ek)S11(em+1,…,ek)

25.5.5 基于VECM最大似然估计的Johansen协整检验

在给定了μt形式后，
基于VECM的最大似然估计对rank(Π)进行系列的似然比检验。
记H(m)表示零假设rank(Π)=m，
在序贯地进行检验时可以看成是rank(Π)≤m，
有

H(0)⊂H(1)⊂⋯⊂H(k)

将(25.5)变成

Δxt=μdt+Πxt−1+Φ∗1Δxt−1+⋯+Φ∗p−1Δxt−p+1+at,t=p+1,…,T(25.6)

目标是检验Π的秩，
这等于Π的非零奇异值的个数，
如果能得到Π的相合估计就可以估计rank(Π)。
从(25.6)可以看出，
Π是排除了μt和Δxt−i,
i=1,…,p−1的线性影响后，
Δxt对xt−1的多对多的线性回归系数。
从上面的MLE步骤可知，
调整后的Δxt和xt−1
分别为û t和v̂ t，
所以

û t=Πv̂ t+at

在正态性假定下，
可以根据û t与v̂ t
的典型相关分析的方法估计rank(Π)，
λ̂ i是û t与v̂ t
的典型相关系数平方。

Johansen提出了两种检验方法。
首先，
考虑假设

H0:rank(Π)=m⟷Ha:rank(Π)>m

(Johansen 1988)提出如下的似然比统计量：

LRtr(m)=−(T−p)∑i=m+1kln(1−λ̂ i)

如果rank(Π)≤m，
则i>m时λ̂ i应该很小，
因此LRtr(m)也应该很小，
当统计量值超过一个右侧临界值时拒绝零假设。
这个检验称为Johansen迹协整检验。
因为单位根的存在，
似然比统计量在零假设下的渐近分布不再服从卡方分布，
其临界值需要用模拟方法。
从H(0)开始到H(1)序贯地进行检验，
在第一次被接受时，
令被接受的H(m)的m为选定的秩。

(Johansen 1988)提出的另一种似然比检验方法针对如下假设

H0:rank(Π)=m⟷Ha:rank(Π)=m+1

取最大特征值统计量

LRmax(m)=−(T−p)ln(1−λ̂ m+1)

当统计量值超过一个右侧临界值时拒绝零假设。
临界值也需要通过模拟获得。
从H(0)开始到H(1)序贯地进行检验，
在第一次被接受时，
令被接受的H(m)的m为选定的秩。

25.5.6 VECM预测

估计的VECM模型可以用于预测。
首先可以从模型得到Δxt序列的预测值，
然后可以从Δxt反解得到xt的预测值。
VECM预测与VAR预测的区别在于VECM允许有单位根和协整关系，
VAR预测不允许有单位根。

25.5.7 Johansen方法进行协整检验的例子

R扩展包urca的ca.jo()函数可以进行计算Johansen的两种检验。
参见(Pfaff 2008)。

例4.2
考虑urca包的UKpppuip数据，
此数据用来检验英国经济的PPP和UIP假设。
季度数据，从1971年第1季度到1987年第2季度。
分量p1是英国批发价格指数，
p2是外贸加权批发价格指数，
e12是英镑实际汇率，
i1是三月期国债利率，
i2是三月期欧元利率，
dpoil0是即期国际油价，
dpoil1是上期国际油价。

载入数据：

library(urca, quietly = TRUE)
data(UKpppuip)
ts.UKpppuip <- ts(UKpppuip, start=c(1971,1), frequency=4)

数据表见附录。

时间序列图：

library(quantmod, quietly = TRUE)
plot(as.xts(ts.UKpppuip), type="l", 
     multi.panel=TRUE, theme="white",
     main="英国PPP-UIP数据",
     major.ticks="years",
     grid.ticks.on = "years")

图25.2: 英国PPP-UIP数据

下面对p1, p2, e12, i1, i2执行Johansen的两种协整检验。
ADF结果表明部分序列有单位根。
以油价作为外生变量。
选项type="trace"使用LRtr统计量，
选项type="eigen"使用LRmax统计量。
选项ecdet指定要考虑的非随机趋势，
"none"表示没有非随机趋势，
"constant"表示常数趋势，
"trend"表示线性趋势。
选项K表示VAR模型的阶。
选项season=4表示有与季度数据相关联的哑变量存在，
用dumvar=输入保存了哑变量(外生变量)的数据框，
其行数必须与用来检验的输入数据行数相同。

dat1 <- UKpppuip[,c("p1", "p2", "e12", "i1", "i2")]
dat2 <- UKpppuip[,c("doilp0", "doilp1")]
summary(ca.jo(
  dat1, type="trace", K=2, season=4, dumvar=dat2))

## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: trace statistic , with linear trend 
## 
## Eigenvalues (lambda):
## [1] 0.40672818 0.28538240 0.25415335 0.10230406 0.08287097
## 
## Values of teststatistic and critical values of test:
## 
##           test 10pct  5pct  1pct
## r <= 4 |  5.19  6.50  8.18 11.65
## r <= 3 | 11.67 15.66 17.95 23.52
## r <= 2 | 29.26 28.71 31.52 37.22
## r <= 1 | 49.42 45.23 48.28 55.43
## r = 0  | 80.75 66.49 70.60 78.87
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##             p1.l2      p2.l2    e12.l2      i1.l2      i2.l2
## p1.l2   1.0000000   1.000000  1.000000  1.0000000  1.0000000
## p2.l2  -0.9086265  -1.143047 -1.272628 -2.4001444 -1.4528820
## e12.l2 -0.9321133  -3.363042  1.113631  1.1221619 -0.4805235
## i1.l2  -3.3746393  35.243576  2.746828 -0.4088865  2.2775510
## i2.l2  -1.8906210 -32.917370 -2.835714  2.9863624  0.7628011
## 
## Weights W:
## (This is the loading matrix)
## 
##             p1.l2         p2.l2       e12.l2        i1.l2       i2.l2
## p1.d  -0.06816507  0.0011795779 -0.002790218  0.001373599 -0.01333013
## p2.d  -0.01773477  0.0001220008 -0.014159241  0.013178503  0.00755575
## e12.d  0.10065321 -0.0001432122 -0.055628059 -0.035400025 -0.04707585
## i1.d   0.03434737 -0.0041631581 -0.010363374  0.012309982 -0.02394672
## i2.d   0.05766426  0.0082830953  0.004821036  0.026984801 -0.01006765

summary(ca.jo(
  dat1, type="eigen", K=2, season=4, dumvar=dat2))

## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: maximal eigenvalue statistic (lambda max) , with linear trend 
## 
## Eigenvalues (lambda):
## [1] 0.40672818 0.28538240 0.25415335 0.10230406 0.08287097
## 
## Values of teststatistic and critical values of test:
## 
##           test 10pct  5pct  1pct
## r <= 4 |  5.19  6.50  8.18 11.65
## r <= 3 |  6.48 12.91 14.90 19.19
## r <= 2 | 17.59 18.90 21.07 25.75
## r <= 1 | 20.16 24.78 27.14 32.14
## r = 0  | 31.33 30.84 33.32 38.78
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##             p1.l2      p2.l2    e12.l2      i1.l2      i2.l2
## p1.l2   1.0000000   1.000000  1.000000  1.0000000  1.0000000
## p2.l2  -0.9086265  -1.143047 -1.272628 -2.4001444 -1.4528820
## e12.l2 -0.9321133  -3.363042  1.113631  1.1221619 -0.4805235
## i1.l2  -3.3746393  35.243576  2.746828 -0.4088865  2.2775510
## i2.l2  -1.8906210 -32.917370 -2.835714  2.9863624  0.7628011
## 
## Weights W:
## (This is the loading matrix)
## 
##             p1.l2         p2.l2       e12.l2        i1.l2       i2.l2
## p1.d  -0.06816507  0.0011795779 -0.002790218  0.001373599 -0.01333013
## p2.d  -0.01773477  0.0001220008 -0.014159241  0.013178503  0.00755575
## e12.d  0.10065321 -0.0001432122 -0.055628059 -0.035400025 -0.04707585
## i1.d   0.03434737 -0.0041631581 -0.010363374  0.012309982 -0.02394672
## i2.d   0.05766426  0.0082830953  0.004821036  0.026984801 -0.01006765

Johansen检验是序贯进行的。
从迹方法的检验结果看出，
在0.05水平下，
H0:r=0（这里r就是前面的m，Π的秩）
统计量为80.75，右侧临界值是70.60，拒绝零假设，
至少有一个协整关系；
H0:r=1的检验统计量为49.42，
右侧临界值为48.28，
拒绝零假设，认为至少有两个协整关系；
H0:r=2的检验统计量为29.26,
右侧临界值为31.52，
不拒绝零假设，
认为恰好有两个协整关系。
如果使用最大特征值方法，
检验结论则有所不同。

结果中包含了所有的协整变量，
前两个协整变量为

⎛⎝⎜⎜⎜⎜⎜⎜p1p2e12i1i21−0.91−0.93−3.37−1.891−1.14−3.3635.2432.91⎞⎠⎟⎟⎟⎟⎟⎟

人为地计算两个协整变量，
检验其平稳性：

dat1 <- UKpppuip[,c("p1", "p2", "e12", "i1", "i2")]
f <- function(){
  x1 <- as.matrix(dat1) %*% c(1, -0.91, -0.93, -3.37, -1.89)
  x2 <- as.matrix(dat1) %*% c(1, -1.14, -3.36, -35.24, -32.91)
  print(fUnitRoots::adfTest(x1, lags=1, type="c"))
  print(fUnitRoots::adfTest(x2, lags=1, type="c"))
  invisible()
}
f()

## Warning in if (class(x) == "timeSeries") x = series(x): 条件的长度大于一，因此只
## 能用其第一元素

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -2.5995
##   P VALUE:
##     0.0994 
## 
## Description:
##  Tue Feb 22 12:01:46 2022 by user: Lenovo

## Warning in if (class(x) == "timeSeries") x = series(x): 条件的长度大于一，因此只
## 能用其第一元素

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -3.0532
##   P VALUE:
##     0.03853 
## 
## Description:
##  Tue Feb 22 12:01:46 2022 by user: Lenovo

在0.05水平下，
第一个协整向量有单位根，
第二个协整向量没有单位根。
但是，因为估计这些权重消耗了自由度，
所以ADF检验的p值是不可用的。

定理25.1 (最小二乘估计的弱相合性) 记Sn=XTnXn,
设n充分大时Xn列满秩，
则β̂ n是β的弱相合估计，
当且仅当limn→∞S−1n=0。

定理25.2 (最小二乘估计的强相合性) 若e1,e2,… 独立同N(0, σ2)分布，
limn→∞S−1n→0，
则β̂ n是β的强相合估计。

25.7.1 英国消费研究数据

25.7.2 英国PPP研究数据