使用R语言进行协整关系检验

2023-02-23 16:37:01Python082

使用R语言进行协整关系检验,第1张

使用R语言进行协整关系检验

协整检验是为了检验非平稳序列的因果关系，协整检验是解决伪回归为问题的重要方法。首先回归伪回归例子：

伪回归Spurious regression伪回归方程的拟合优度、显著性水平等指标都很好，但是其残差序列是一个非平稳序列，拟合一个伪回归:

#调用相关R包

library(lmtest)

library(tseries)

#模拟序列

set.seed(123456)

e1=rnorm(500)

e2=rnorm(500)

trd=1:500

y1=0.8*trd+cumsum(e1)

y2=0.6*trd+cumsum(e2)

sr.reg=lm(y1~y2)

#提取回归残差

error=residuals(sr.reg)

#作残差散点图

plot(error, main="Plot of error")

#对残差进行单位根检验

adf.test(error)

## Dickey-Fuller = -2.548, Lag order = 7, p-value = 0.3463

## alternative hypothesis: stationary

#伪回归结果，相关参数都显著

summary(sr.reg)

## Residuals:

## Min 1Q Median 3Q Max

## -30.654 -11.526 0.359 11.142 31.006

## Coefficients:

## Estimate Std. Error t value Pr(>|t|)

## (Intercept) -29.32697 1.36716 -21.4 <2e-16 ***

## y2 1.44079 0.00752 191.6 <2e-16 ***

## Residual standard error: 13.7 on 498 degrees of freedom

## Multiple R-squared: 0.987, Adjusted R-squared: 0.987

## F-statistic: 3.67e+04 on 1 and 498 DF, p-value: <2e-16

dwtest(sr.reg)

## DW = 0.0172, p-value <2.2e-16

恩格尔-格兰杰检验Engle-Granger第一步：建立两变量（y1,y2）的回归方程，第二部：对该回归方程的残差(resid)进行单位根检验其中，原假设两变量不存在协整关系，备择假设是两变量存在协整关系。利用最小二乘法对回归方程进行估计，从回归方程中提取残差进行检验。

set.seed(123456)

e1=rnorm(100)

e2=rnorm(100)

y1=cumsum(e1)

y2=0.6*y1+e2

# (伪)回归模型

lr.reg=lm(y2~y1)

error=residuals(lr.reg)

adf.test(error)

## Dickey-Fuller = -3.988, Lag order = 4, p-value = 0.01262

## alternative hypothesis: stationary

error.lagged=error[-c(99,100)]

#建立误差修正模型ECM.REG

dy1=diff(y1)

dy2=diff(y2)

diff.dat=data.frame(embed(cbind(dy1, dy2),2))#emed表示嵌入时间序列dy1,dy2到diff.dat

colnames(diff.dat)=c("dy1","dy2","dy1.1","dy2.1")

ecm.reg=lm(dy2~error.lagged+dy1.1+dy2.1, data=diff.dat)

summary(ecm.reg)

## Residuals:

## Min 1Q Median 3Q Max

## -2.959 -0.544 0.137 0.711 2.307

## Coefficients:

## Estimate Std. Error t value Pr(>|t|)

## (Intercept) 0.0034 0.1036 0.03 0.97

## error.lagged -0.9688 0.1585 -6.11 2.2e-08 ***

## dy1.1 0.8086 0.1120 7.22 1.4e-10 ***

## dy2.1 -1.0589 0.1084 -9.77 5.6e-16 ***

## Residual standard error: 1.03 on 94 degrees of freedom

## Multiple R-squared: 0.546, Adjusted R-squared: 0.532

## F-statistic: 37.7 on 3 and 94 DF, p-value: 4.24e-16

par(mfrow=c(2,2))

plot(ecm.reg)

Johansen-Juselius(JJ)协整检验法，该方法是一种用向量自回归(VAR)模型进行检验的方法，适用于对多重一阶单整I(1)序列进行协整检验。JJ检验有两种：特征值轨迹检验和最大特征值检验。我们可以调用urca包中的ca.jo命令完成这两种检验。其语法：

ca.jo(x, type = c("eigen", "trace"), ecdet = c("none", "const", "trend"), K = 2,spec=c("longrun", "transitory"), season = NULL, dumvar = NULL)

其中：x为矩阵形式数据框；type用来设置检验方法；ecdet用于设置模型形式：none表示不带截距项，const表示带常数截距项，trend表示带趋势项。K表示自回归序列的滞后阶数；spec表示向量误差修正模型反映的序列间的长期或短期关系；season表示季节效应；dumvar表示哑变量设置。

set.seed(12345)e1=rnorm(250,0,0.5)e2=rnorm(250,0,0.5)e3=rnorm(250,0,0.5)#模拟没有移动平均的向量自回归序列；u1.ar1=arima.sim(model=list(ar=0.75), innov=e1, n=250)u2.ar1=arima.sim(model=list(ar=0.3), innov=e2, n=250)y3=cumsum(e3)y1=0.8*y3+u1.ar1y2=-0.3*y3+u2.ar1#合并y1,y2,y3构成进行JJ检验的数据库；y.mat=data.frame(y1, y2, y3)#调用urca包中cajo命令对向量自回归序列进行JJ协整检验vecm=ca.jo(y.mat)jo.results=summary(vecm)#cajorls命令可以得到限制协整阶数的向量误差修正模型的最小二乘法回归结果vecm.r2=cajorls(vecm, r=2)vecm.r2## Call:lm(formula = substitute(form1), data = data.mat)## Coefficients:## y1.d y2.d y3.d## ect1 -0.33129 0.06461 0.01268## ect2 0.09447 -0.70938 -0.00916## constant 0.16837 -0.02702 0.02526## y1.dl1-0.22768 0.02701 0.06816## y2.dl1 0.14445 -0.71561 0.04049## y3.dl1 0.12347 -0.29083 -0.07525## $beta## ect1 ect2## y1.l2 1.000e+00 0.0000## y2.l2 -3.402e-18 1.0000## y3.l2 -7.329e-01 0.2952

这里分享一下R语言实现VAR和SVAR的整个流程。

主要步骤包括：

1.单位根检验

2.确定滞后阶数

3.格兰杰因果检验

4.模型稳定性检验

5.脉冲响应

6.方差分解

(Johansen协整检验，如果需要的话)

整个过程用到的R语言的扩展包有：

library(zoo)

library(vars)

library(tseries)

首先，数据是下面的样子：

ps：数据是时间序列类型,可以通过下面方法将dataframe转成时间序列类型

data = ts(data)

1.单位根检验

#对data的第一列进行单位根检验

adf.test(data[,1])

2.滞后阶数确定

VARselect函数结果包括AIC、HQ、SC和FPE准则

#参数y为时间序列数据，lag.max为最大滞后阶数

#参数type值包括const截距，trend趋势，both同时包含截距和趋势，none不包含截距和趋势

VARselect(y=data, lag.max = 10, type = c("const"))

3.格兰杰因果检验

格兰杰因果检验有两个方法，第一个是在构造模型之前，第二个是在构造模型之后在模型的基础上进行格兰杰因果检验。

(1)构造模型之前格兰杰因果检验

#函数格式：grangertest(yt~xt)

eg：

grangertest(Value~BCI)

(2)构造模型之后格兰杰因果检验

#函数格式：causality(VARModel,cause)

var = VAR(data ,p = 2, type = "const")

causality(var,cause=c('Count','Value'))

ps:在这里如果想要构建SVAR模型的话，需要根据实际情况构建两个矩阵amat和bmat，然后使用这两个矩阵来构建SVAR模型：

svar = SVAR(var,Amat = amat,Bmat = bmat)

4.模型稳定性检验

#这里使用“OLS-CUSUM”，它给出的是残差累积和，在该检验生成的曲线图中，残差累积和曲线以时间为横坐标，

#图中绘出两条临界线，如果累积和超出了这两条临界线，则说明参数不具有稳定性。

sta = stability(var, type = c("OLS-CUSUM"), h = 0.15, dynamic = FALSE, rescale = TRUE)

plot(sta)##结果稳健

5.脉冲响应

#标题栏说明，这是BCI(或者其他变量)对各个变量（包括BCI自身）的脉冲响应

(1)VAR脉冲响应

var.irf<-irf(var,n.head=10)

plot(var.irf)

(2)SVAR脉冲响应

svar.irf<-irf(svar,n.ahead = 100)

plot(svar.irf)

6.方差分解

#反映了各变量的贡献率

(1)VAR方差分解

fevd1<-fevd(var, n.ahead = 10)

fevd1$Count

(2)SVAR方差分解

fevd2<-fevd(svar, n.ahead = 10)

fevd2$Value

ps：有时候需要进行Johansen协整检验

#Johansen协整检验，

#对r=0（不存在协整关系）的检验统计量大于临界值，表明拒绝原假设

yJoTest = ca.jo(data, type = c("trace"), ecdet = c("none"), K = 2)

summary(yJoTest)

网页链接

数据准备

许多实际情况中统计假设（假定观测数据抽样自正态分布或者其他性质较好的理论分布）并不一定满足，比如数据抽样于未知或混合分布、样本量过小、存在离群点、基于理论分布设计合适的统计检验过于复杂且数学上难以处理等情况，这时基于随机化和重抽样的统计方法就可派上用场。

置换检验的定义

置换检验(Permutation test)，也称随机化检验或重随机化检验，是Fisher于20世纪30年代提出的一种基于大量计算（computationally intensive），利用样本数据的全（或随机）排列，进行统计推断的方法，因其对总体分布自由，应用较为广泛，特别适用于总体分布未知的小样本资料，以及某些难以用常规方法分析资料的假设检验问题。

置换检验的原理

1、提出原假设，比如XX处理后结果没有变化

2、计算统计量，如两组的均值之差，记作t0

3、将所有样本放在一起，然后随机排序进行分组，再计算其统计量t1

4、重复第3步骤，直至所有排序可能性都齐全（比如有A组有n样本，B组有m样本，则总重复次数相当于从n+m中随机抽取n个的次数），得到一系列的统计量（t1-tn）

5、最后将这些统计量按照从小到大排序，构成抽样分布，再看t0是否落在分布的置信区间内（如95%置信区间），这时候可计算一个P值（如果抽样总体1000次统计量中大于t0的有10个，则估计的P值为10/1000=0.01），落在置信区间外则拒绝原假设

6、如果第3步骤是将所有可能性都计算了的话，则是精确检验；如果只取了计算了部分组合，则是近似结果，这时一般用蒙特卡罗模拟（Monte Carlo simulation）的方法进行置换检验

7、置换检验和参数检验都计算了统计量，但是前者是跟置换观测数据后获得的经验分布进行比较，后者则是跟理论分布进行比较。

请牢记：置换检验都是使用伪随机数来从所有可能的排列组合中进行抽样（当做近似检验时）。因此，每次检验的结果都有所不同。

coin包提供了一个进行置换检验的一般性框架。通过该包，你可以回答如下问题。

响应值与组的分配独立吗？

两个数值变量独立吗？

两个类别型变量独立吗？

表12-2列出来的每个函数都是如下形式：

function_name(formula, data, distribution=)

其中：

 formula描述的是要检验变量间的关系。示例可参见表12-2；

 data是一个数据框；

 distribution指定经验分布在零假设条件下的形式，可能值有exact，asymptotic和

approximate。

若distribution = "exact"，那么在零假设条件下，分布的计算是精确的（即依据所有可能的排列组合）。当然，也可以根据它的渐进分布（distribution = "asymptotic"）或蒙特卡洛重抽样（distribution = "approxiamate(B = #)"）来做近似计算，其中#指所需重复的次数。

distribution = "exact"当前仅可用于两样本问题。

传统t检验表明存在显著性差异（p <0.05），而精确检验却表明差异并不显著（p >0.072）。

第7章我用自己的数据进行了t检验，对比一下传统t检验和置换检验，结果如下：

两种检验方式下结果都是显著的

Wilcoxon-Mann-Whitney U检验

coin包规定所有的类别型变量都必须以因子形式编码。

wilcox.test()默认计算的也是精确分布。

K样本检验的置换检验

通过chisq_test()或cmh_test()函数，我们可用置换检验判断两类别型变量的独立性。当数据可根据第三个类别型变量进行分层时，需要使用后一个函数。若变量都是有序型，可使用lbl_test()函数来检验是否存在线性趋势。

卡方独立性检验

卡方独立性检验的置换检验

你可能会有疑问，为什么需要把变量Improved从一个有序因子变成一个分类因子？（好问题！）这是因为，如果你用有序因子，coin()将会生成一个线性与线性趋势检验，而不是卡方检验。

结果解读：两种检验下p值都是小于0.05，说明Treatment和Improved之间相互不独立

自己数据的演示

结果解读：p值均为1，表明nitrogen和variety相互独立。

spearman_test()函数提供了两数值变量的独立性置换检验。

当处于不同组的观测已经被分配得当，或者使用了重复测量时，样本相关检验便可派上用场。

对于两配对组的置换检验，可使用wilcoxsign_test()函数；多于两组时，使用friedman_test()函数。

自己数据演示

lmPerm包可做线性模型的置换检验。比如lmp()和aovp()函数即lm()和aov()函数的修改版，能够进行置换检验，而非正态理论检验。

lmp()和aovp()函数的参数与lm()和aov()函数类似，只额外添加了perm =参数。

perm =选项的可选值有"Exact"、"Prob"或"SPR"。Exact根据所有可能的排列组合生成精确检验。Prob从所有可能的排列中不断抽样，直至估计的标准差在估计的p值0.1之下，判停准则由可选的Ca参数控制。SPR使用贯序概率比检验来判断何时停止抽样。注意，若观测数大于10，perm = "Exact"将自动默认转为perm = "Prob"，因为精确检验只适用于小样本问题。

简单线性回归的置换检验

R语言实战的例子：

多项式回归的置换检验

R语言实战的例子：

自己数据集的例子：

R语言实战的例子：

自己数据集的例子：

当两种方法所得结果不一致时，你需要更加谨慎地审视数据，这很可能是因为违反了正态性假设或者存在离群点。

R语言实战的例子：

自己数据集的例子：

R语言实战的例子：

自己数据集的例子：

R语言实战的例子：

自己数据集的例子：

值得注意的是，当将aovp()应用到方差分析设计中时，它默认使用唯一平方和法（SAS也称为类型III平方和）。每种效应都会依据其他效应做相应调整。R中默认的参数化方差分析设计使用的是序贯平方和（SAS是类型I平方和）。每种效应依据模型中先出现的效应做相应调整。对于平衡设计，两种方法结果相同，但是对于每个单元格观测数不同的不平衡设计，两种方法结果则不同。不平衡性越大，结果分歧越大。若在aovp()函数中设定seqs = TRUE，可以生成你想要的序贯平方和。

你可能已经注意到，基于正态理论的检验与上面置换检验的结果非常接近。在这些问题中数据表现非常好，两种方法结果的一致性也验证了正态理论方法适用于上述示例。

当然，置换检验真正发挥功用的地方是处理非正态数据（如分布偏倚很大）、存在离群点、样本很小或无法做参数检验等情况。不过，如果初始样本对感兴趣的总体情况代表性很差，即使是置换检验也无法提高推断效果。

置换检验主要用于生成检验零假设的p值，它有助于回答“效应是否存在”这样的问题。不过，置换方法对于获取置信区间和估计测量精度是比较困难的。幸运的是，这正是自助法大显神通的地方。

所谓自助法，即从初始样本重复随机替换抽样，生成一个或一系列待检验统计量的经验分布。无需假设一个特定的理论分布，便可生成统计量的置信区间，并能检验统计假设。

倘若你假设均值的样本分布不是正态分布，该怎么办呢？可使用自助法。

(1)从样本中随机选择10个观测，抽样后再放回。有些观测可能会被选择多次，有些可能一直都不会被选中。

(2)计算并记录样本均值。

(3)重复1和2一千次。

(4)将1000个样本均值从小到大排序。

(5)找出样本均值2.5%和97.5%的分位点。此时即初始位置和最末位置的第25个数，它们就限定了95%的置信区间。

boot包扩展了自助法和重抽样的相关用途。你可以对一个统计量(如中位数)或一个统计量向量(如一列回归系数)使用自助法。

一般来说，自助法有三个主要步骤。