开始于：
data("Canada")
dim(Canada) #84observations x 4 variables
VARselect(Canada) # since in small samples, AIC>BIC; VAR(3) is chosen.
现在，加拿大数据集的范围：1980.1 - 2000.4（20 年）对于建模来说已经足够长了。这长达 20 年的时期肯定包括许多危机和干预措施。因此，必须搜索数据中的结构中断。这是必要的，因为在结构破碎的序列中，SBs 的存在会改变非平稳性检验的 t 值（从而影响序列是否平稳的决定）。
由于 Narayan-Popp 2010 在多个结构断裂下的非平稳性检验在统计上对以前的检验非常有效（Lee-Strazichic2003、Zivot-Andres1992），并且自从 Joyeux 2007（在 Rao2007 中）证明了这些先前检验的不合逻辑性和 NP2013已经证明了NP2010统计能力的优越性，必须使用NP2010。由于 NP2010 的 Gauss 代码对我来说似乎很难看，我将其转换为 R 代码，并在 ggplot2 的帮助下，结果呈现得更好。
[处理结构中断也是协整检查的必要条件，因为 Osterwald-Lenum1992 CV 忽略 SB，而 Johansen-Mosconi-Nielsen2000 CV 关心 SB。]
Canada <- as.data.frame(Canada)
head(Canada)
e prod rw U
1 929.6105 405.3665 386.1361 7.53
2 929.8040 404.6398 388.1358 7.70
...................................
# Assign lexiographic row names for dates of observations
row.names(Canada) <- paste(sort(rep(seq(1980, 2000, 1), 4) ), rep(seq(1, 4, 1), 20), sep = ".")
# Insert lexiographic "date" column to the dataframe. This is necessary for creating intervention dummies.
DCanada <- data.frame(date=row.names(Canada),Canada) # dataset with obs dates in a column
head(DCanada)
date e prod rw U
1980.1 1980.1 929.6105 405.3665 386.1361 7.53
1980.2 1980.2 929.8040 404.6398 388.1358 7.70
对序列执行 Narayan-Popp 2010 非平稳性检验：
[H0：“（有 2 个结构中断）系列是非平稳的”；
H1：“（有 2 个结构中断）系列是静止的”；
“测试状态 > 临界值” => “保持 H0”； "test stat "hold H1"]
library(causfinder)
narayanpopp(DCanada[,2]) # for e
narayanpopp(DCanada[,3]) # for prod
narayanpopp(DCanada[,4]) # for rw
narayanpopp(DCanada[,5]) # for U
Narayan-Popp 2010 非平稳性测试结果（带 obs #s）：
variable t stat lag SB1 SB2 Integration Order
e -4.164 2 37:946.86 43:948.03 I(1)
prod -3.325 1 24:406.77 44:405.43 I(1)
rw -5.087 0 36:436.15 44:446.96 I(0) <trend-stationary>
U -5.737 1 43:8.169 53:11.070 I(0) <stationary pattern> (M2 computationally singular; used M1 model)
(critical values (M2): (1%,5%,10%): -5.576 -4.937 -4.596)
(critical values (M1): (1%,5%,10%): -4.958 -4.316 -3.980
由于在 VAR 结构中，所有变量都被平等对待，在系统确定结构中断时继续平等处理：
mean(c(37,24,36,43)) # 35; SB1 of system=1988.3
mean(c(43,44,44,53)) # 46; SB2 of system=1990.2
以下是克服"In Ops.factor(left, right) : >= not meaningful for factors"错误。在某些数据集中，我们需要执行以下操作：
library(readxl)
write.xlsx(Canada, file="data.xlsx", row.names=FALSE) # Take this to the below folder, add "date" column with values 1980.1,....,2000.4
mydata <- read_excel("D://eKitap//RAO 2007 Cointegration for the applied economist 2E//JoyeuxCalisma//Canada//data.xlsx")
# arrange your path accordingly in the above line.
mydata <- as.data.frame(mydata)
library(lubridate); library(zoo)
row.names(mydata) <- as.yearqtr(seq(ymd('1980-01-01'), by = '1 quarter', length.out=(84)))
Dmydata <- mydata # Hold it in a variable
使用 2 个 SB（35:1988.3 和 46:1990.2）定义干预虚拟矩阵，如下所示：
library(data.table)
DataTable <- data.table(Dmydata, keep.rownames=FALSE)
Dt <- cbind("bir"=1, # intervention dummies matrix
"D2t" = as.numeric(ifelse( DataTable[,c("date"), with=FALSE] >= "1988.3" & DataTable[,c("date"), with=FALSE] <= "1990.1", 1 , 0)),
"D3t" = as.numeric(ifelse( DataTable[,c("date"), with=FALSE] >= "1990.2" & DataTable[,c("date"), with=FALSE] <= "2000.4", 1 , 0)))
伴随干预假人的动态指标变量：
OnTheFlyIndicator <- cbind(
"I2t" = as.numeric(DataTable[, c("date"), with=FALSE] == "1988.3"),
"I3t" = as.numeric(DataTable[, c("date"), with=FALSE] == "1990.2"))
myTimeTrend <- as.matrix(cbind("TimeTrend" = as.numeric(1:nrow(Dt))))
zyDt <- Dt * as.vector(myTimeTrend) # TimeTrendDavranisDegisimleri
colnames(zyDt) <- paste(colnames(myTimeTrend), colnames(Dt), sep="*")
mydata <- mydata[,-1]
VAR订单的选择：
library(vars)
# Lag order selection with the effects of intervention dummies
VARselect(mydata, lag.max=5, "both", exogen=cbind(zyDt[drop=FALSE], Dt[drop=FALSE], OnTheFlyIndicator)) # Take VAR(3)
Joyeux2007 索引技术的滞后矩阵：
lagmatrix <- function(x, maxlag){
x <- as.matrix(x)
if(is.null(colnames(x))== TRUE){ colnames(x) <- "VarCol0" }
DondurulenDizey <- embed(c(rep(NA,maxlag),x),maxlag+1)
dimnames(DondurulenDizey)[[2]] <- c(colnames(x)[1, drop = FALSE], paste(colnames(x)[1,drop=FALSE],".",1:maxlag,"l", sep = ""))
return(DondurulenDizey)
}
分配 VAR 滞后和否。子样本数：
VARlag <- 3
Subsamples <- 3 # subsamples = no. of str breaks +1
2 个结构中断的虚拟矩阵：
dummymatrix2SB <- matrix(NA,DataTable[,.N], 10)
dummymatrix2SB <- cbind(myTimeTrend,
lagmatrix(zyDt[,c("TimeTrend*D2t"), drop=FALSE], maxlag=VARlag)[,1+VARlag, drop=FALSE],
lagmatrix(zyDt[,c("TimeTrend*D3t"), drop=FALSE], maxlag=VARlag)[,1+VARlag, drop=FALSE],
lagmatrix(Dt[,c("D2t"), drop=FALSE], maxlag=VARlag)[,1+VARlag, drop=FALSE],
lagmatrix(Dt[,c("D3t"), drop=FALSE], maxlag=VARlag)[,1+VARlag, drop=FALSE],
lagmatrix(OnTheFlyIndicator[,c("I2t"), drop=FALSE], maxlag=VARlag-1),
lagmatrix(OnTheFlyIndicator[,c("I3t"), drop=FALSE], maxlag=VARlag-1))
dummymatrix2SB[is.na(dummymatrix2SB)] <- 0 # replace NAs with 0
dummymatrix2SB # Print dummy matrix for 2 str breaks to make sure all are OK
TimeTrend TimeTrend.D2t.3l TimeTrend.D3t.3l D2t.3l D3t.3l I2t I2t.1l I2t.2l I3t I3t.1l I3t.2l
1 0 0 0 0 0 0 0 0 0 0
2 0 0 0 0 0 0 0 0 0 0
...........................................
34 0 0 0 0 0 0 0 0 0 0
35 0 0 0 0 1 0 0 0 0 0
36 0 0 0 0 0 1 0 0 0 0
37 0 0 0 0 0 0 1 0 0 0
38 35 0 1 0 0 0 0 0 0 0
39 36 0 1 0 0 0 0 0 0 0
40 37 0 1 0 0 0 0 0 0 0
41 38 0 1 0 0 0 0 0 0 0
42 39 0 1 0 0 0 0 1 0 0
43 40 0 1 0 0 0 0 0 1 0
44 41 0 1 0 0 0 0 0 0 1
45 0 42 0 1 0 0 0 0 0 0
46 0 43 0 1 0 0 0 0 0 0
............................................
83 0 80 0 1 0 0 0 0 0 0
84 0 81 0 1 0 0 0 0 0 0
VAR 的稳定性：
维克多，理论上你错了。即使在受限（协整）VAR 模型的情况下，也会从 VAR 侧检查稳定性。有关详细信息，请参阅 Joyeux2007。此外，双方的估计是相同的：
“不受限制的 VAR = 不受限制的 VECM”和
“受限 VAR = 受限 VECM”。
因此，检查不受限制的 VAR 的稳定性等于检查不受限制的 VECM 的稳定性，反之亦然。它们在数学上是相等的，只是表示不同。
此外，检查受限 VAR 的稳定性等于检查受限 VECM 的稳定性，反之亦然。它们在数学上是相等的，它们只是不同的表示。但是，您不需要对受限 VECM 情况进行此检查，因为我们正在浏览可行 VAR 的子空间。也就是说，如果restd VeCM对应的原始unr VAR稳定，则一切OK。
如果您的系列是协整的，即使在这种情况下，您也可以从 VAR 端检查稳定性！如果您想知道“是否应该检查受限 VECM 的稳定性”，答案是否定的。你不应该检查。因为，在协整情况下，您处于可行解的子空间中。也就是说，如果您坚持检查受限（协整）VECM 的稳定性，您仍然可以通过 urca::ca.jo 扩展和 vars::vec2var 扩展来做到这一点：
print(roots(VAR(mydata, p=3, "both", exogen=cbind(zyDt[drop=FALSE], Dt[drop=FALSE], OnTheFlyIndicator)), modulus=TRUE))
# [1] 0.96132524 0.77923543 0.68689517 0.68689517 0.67578368 0.67578368
[7] 0.59065419 0.59065419 0.55983617 0.55983617 0.33700725 0.09363846
print(max(roots(VAR(mydata, p=3, "both", exogen=cbind(zyDt[drop=FALSE], Dt[drop=FALSE], OnTheFlyIndicator)), modulus=TRUE)))
#0.9613252
（可选）通过 OLS-CUSUM 检查稳定性：
plot(stability(VAR(mydata, p=3, "both", exogen=cbind(zyDt[drop=FALSE], Dt[drop=FALSE], OnTheFlyIndicator)), type="OLS-CUSUM"))
VAR残差检验的非自相关：
for (j in as.integer(1:5)){
print(paste("VAR's lag no:", j))
print(serial.test(VAR(mydata, p=j, "both", exogen=cbind(zyDt[drop=FALSE], Dt[drop=FALSE], OnTheFlyIndicator)), lags.bg=4, type= c("ES")))
# lags.bg: AR order of VAR residuals
}
VAR残差检验的NORMALITY：
print(normality.test(VAR(mydata, p=3, "both", exogen=cbind(zyDt[drop=FALSE], Dt[drop=FALSE], OnTheFlyIndicator)), multivariate=TRUE))
library(normtest)
for (i in as.integer(1:4)){ # there are 4 variables
print(skewness.norm.test(resid(VAR(mydata, p=3, "both", exogen=cbind(zyDt[drop=FALSE], Dt[drop=FALSE], OnTheFlyIndicator)))[,i]))
print(kurtosis.norm.test(resid(VAR(mydata, p=3, "both", exogen=cbind(zyDt[drop=FALSE], Dt[drop=FALSE], OnTheFlyIndicator)))[,i]))
print(jb.norm.test(resid(VAR(mydata, p=3, "both", exogen=cbind(zyDt[drop=FALSE], Dt[drop=FALSE], OnTheFlyIndicator)))[,i]))
}
VAR残差检验的HOMOSCEDASTICITY：
print(arch.test(VAR(mydata, p=3, "both", exogen=cbind(zyDt[drop=FALSE], Dt[drop=FALSE], OnTheFlyIndicator))), lags.multi=6, multivariate.only=TRUE)
由于系列的积分阶不同，它们不可能是协整的。也就是说，
暂时假设所有都是 I(1) 并使用 Johansen-Mosconi-Nielsen 2000 CV 执行具有多个结构中断的协整检验：
（将 urca::cajo 扩展为 causfinder::ykJohEsbInc（即添加处理 1 SB 和 2 SB 的功能））
summary(ykJohEsbInc(mydata, type="trace", ecdet="zamanda2yk", K=3, spec="longrun", dumvar=dummymatrix2SB[,c(-1,-2,-3)]))
# summary(ykJohEsbInc(mydata, type="trace", ecdet="zamanda2yk", K=3, spec="transitory", dumvar=dummymatrix2SB[,c(-1,-2,-3)])) gives the exactly same result.
由于系统中有 2 个 SB（1988.3、1990.2），因此有 q=2+1=3 个子样本。
第一个 SB 比率：v1= (35-1)/84= 0.4047619
第二SB比率：v2= (46-1)/84= 0.5357143
因此，JMN2000 CV 用于 2 个 SB 的协整检验：
（以下是TR本地化的。可以在Giles网站上找到原始的EN-local代码）
library(gplots)
# Johansen vd. (2000) nin buldugu, yapisal kirilmalarin varliginda esbutunlesim incelemesinin degistirilmis iz sinamalarinin yanasik p degerleri ve karar degerlerini hesaplama kodu
# Ryan Godwin & David Giles (Dept. of Economics, Univesity of Victoria, Canada), 29.06.2011
# Kullanici asagidaki 4 degeri atamalidir
#======================================
degiskensayisi <- 4 # p
q<- 3 # q: verideki farkli donemlerin sayisi; q=1: 1 donem, hicbir yapisal kirilma yok demek oldugundan v1 ve v2 nin degerleri ihmal edilir
v1<- 0.4047619 # (35-1)/84 # 1.yk anı=34+1=35. Johansen et. al 2000 v1 def'n , v1: SB1 - 1
v2<- 0.5357143 # (46-1)/84 # 2nd SB moment 45+1=46.
#======================================
# iz istatistiginin biri veya her ikisi icin p degerlerinin olmasi istendiginde, sonraki 2 satirin biri veya her ikisini degistir
izZ <- 15.09 # Vz(r) istatistiginin degeri
izK <- 114.7 # Vk(r) istatistiginin degeri
#=========================================
enbuyuk_p_r<- degiskensayisi # "p-r > 10" olmasın; bkz: Johansen vd. (2000)
# "a" ve "b" nin değerleri yapısal kırılmaların sayısına (q-1) bağlıdır
# q=1 iken, hiçbir yapısal kırılma olmadığı bu durumda a=b=0 ata
# q=2 iken, 1 yapısal kırılma olduğu bu durumda a=0 (Johansen vd. 2000 4.Tabloda) ve b=min[V1 , (1-V1)] ata
# q=3 iken, 2 yapısal kırılma olduğu bu durumda a=min[V1, (V2-V1), (1-V2)] ve b=min[geriye kalan iki V ifadesi] ata
a = c(0, 0, min(v1, v2-v1, 1-v2))[q]
b = c(0, min(v1, 1-v1), median(c(v1,v2-v1,1-v2)))[q]
# YanDagOrtLog: yanaşık dağılımın ortalamasının logaritması
# YanDagDegLog: yanaşık dağılımın değişmesinin logaritması
# V(Zamanyönsemsi) veya V(Kesme) sınamalarını yansıtmak üzere adlara z veya k ekle.
# Bkz. Johansen vd. (2000) 4. Tablo.
# Önce Vz(r) sınamasının sonra Vk(r) sınamasının karar değerlerini oluştur
pr<- c(1:enbuyuk_p_r)
YanDagOrtLogZ <- 3.06+0.456*pr+1.47*a+0.993*b-0.0269*pr^2-0.0363*a*pr-0.0195*b*pr-4.21*a^2-2.35*b^2+0.000840*pr^3+6.01*a^3-1.33*a^2*b+2.04*b^3-2.05/pr-0.304*a/pr+1.06*b/pr
+9.35*a^2/pr+3.82*a*b/pr+2.12*b^2/pr-22.8*a^3/pr-7.15*a*b^2/pr-4.95*b^3/pr+0.681/pr^2-0.828*b/pr^2-5.43*a^2/pr^2+13.1*a^3/pr^2+1.5*b^3/pr^2
YanDagDegLogZ <- 3.97+0.314*pr+1.79*a+0.256*b-0.00898*pr^2-0.0688*a*pr-4.08*a^2+4.75*a^3-0.587*b^3-2.47/pr+1.62*a/pr+3.13*b/pr-4.52*a^2/pr-1.21*a*b/pr-5.87*b^2/pr+4.89*b^3/pr
+0.874/pr^2-0.865*b/pr^2
OrtalamaZ<- exp(YanDagOrtLogZ)-(3-q)*pr
DegismeZ<- exp(YanDagDegLogZ)-2*(3-q)*pr
# Sinama istatistiginin yanasik dagilimina yaklasmakta kullanilacak Gama dagiliminin sekil ve olcek degiskelerini elde etmek icin yanasik ortalama ve degismeyi kullanarak
# V0 varsayimi altinda istenen quantilelari elde et:
# quantilelar: olasilik dagiliminin araligini veya bir ornekteki gozlemleri, esit olasiliklara sahip birbirlerine bitisik araliklarla bolen kesim noktalari.
tetaZ <- DegismeZ/OrtalamaZ
kZ <- OrtalamaZ^2/DegismeZ
YanDagOrtLogK<- 2.80+0.501*pr+1.43*a+0.399*b-0.0309*pr^2-0.0600*a*pr-5.72*a^2-1.12*a*b-1.70*b^2+0.000974*pr^3+0.168*a^2*pr+6.34*a^3+1.89*a*b^2+1.85*b^3-2.19/pr-0.438*a/pr
+1.79*b/pr+6.03*a^2/pr+3.08*a*b/pr-1.97*b^2/pr-8.08*a^3/pr-5.79*a*b^2/pr+0.717/pr^2-1.29*b/pr^2-1.52*a^2/pr^2+2.87*b^2/pr^2-2.03*b^3/pr^2
YanDagDegLogK<- 3.78+0.346*pr+0.859*a-0.0106*pr^2-0.0339*a*pr-2.35*a^2+3.95*a^3-0.282*b^3-2.73/pr+0.874*a/pr+2.36*b/pr-2.88*a^2/pr-4.44*b^2/pr+4.31*b^3/pr+1.02/pr^2-0.807*b/pr^2
OrtalamaK <- exp(YanDagOrtLogK)-(3-q)*pr
DegismeK <- exp(YanDagDegLogK)-2*(3-q)*pr
# Sinama istatistiginin yanasik dagilimina yaklasmakta kullanilacak Gama dagiliminin sekil ve olcek degiskelerini elde etmek icin yanasik ortalama ve degismeyi kullanarak
# V0 varsayimi altinda istenen quantilelari elde et:
# quantilelar: olasilik dagiliminin araligini veya bir ornekteki gozlemleri, esit olasiliklara sahip birbirlerine bitisik araliklarla bolen kesim noktalari.
tetaK <- DegismeK/OrtalamaK
kK <- OrtalamaK^2/DegismeK
# (izZ veya izK den biri 0 dan farklı ise) karar değerlerini ve p değerlerini tablolaştır:
windows(6,3.8)
KararDegerleri <- cbind(sapply(c(.90,.95,.99) , function(x) sprintf("%.2f",round(c(qgamma(x, shape=kZ,scale=tetaZ)),2))),
sapply(c(.9,.95,.99) , function(x) sprintf("%.2f",round(c(qgamma(x, shape=kK,scale=tetaK)),2))))
colnames(KararDegerleri) <- rep(c(0.90,0.95,0.99),2)
# rownames(KararDegerleri) <- pr
rownames(KararDegerleri) <- c(sapply((degiskensayisi -1):1, function(i) paste(degiskensayisi - i, " ","(r<=", i, ")",sep="")), paste(degiskensayisi, " ( r=0)", sep=""))
textplot(KararDegerleri, cex=1)
text(.064,.91,"p-r",font=2)
text(.345,1,expression(paste(plain(V)[z],"(r) test")),col=2)
text(.821,1,expression(paste(plain(V)[k],"(r) test")),col=4)
title("Yanasik Karar Degerleri \n (p:duzendeki degisken sayisi; r:esbutunlesim ranki)")
if(izZ!=0){
windows(4,3.8)
pDegerleri <- matrix(sprintf("%.3f",round(1 - pgamma(izZ, shape=kZ, scale = tetaZ),3)))
# rownames(pDegerleri) <- pr
rownames(pDegerleri) <- c(sapply((degiskensayisi -1):1, function(i) paste(degiskensayisi - i, " ","(r<=", i, ")",sep="")), paste(degiskensayisi, " ( r=0)", sep=""))
textplot(pDegerleri,cex=1,show.colnames=F)
text(.69,.96,substitute(paste("Pr(",plain(V)[z],">",nn,")"),list(nn=izZ)),col=2)
text(.45,.96,"p-r",font=2)
title("Yanasik p Degerleri \n (p:duzendeki degisken sayisi; \n r:esbutunlesim ranki)")
}
if(izK!=0){
windows(3,3.8)
pDegerleri <- matrix(sprintf("%.3f",round(1 - pgamma(izK, shape=kK, scale = tetaK),3)))
#rownames(pDegerleri) <- pr
rownames(pDegerleri) <- c(sapply((degiskensayisi -1):1, function(i) paste(degiskensayisi - i, " ","(r<=", i, ")",sep="")), paste(degiskensayisi, " ( r=0)", sep=""))
textplot(pDegerleri,cex=1,show.colnames=F)
text(.78,.96,substitute(paste("Pr(",plain(V)[k],">",nn,")"),list(nn=izK)),col=4)
text(.43,.96,"p-r",font=2)
title("Yanasik p Degerleri \n (p:duzendeki degisken sayisi; \n r:esbutunlesim ranki)")
}
因此，根据 JMN2000 CV，也没有协整。因此，您对 vec2var 的使用毫无意义。因为，在协整情况下需要 vec2var。再次，假设所有系列都是协整的，让你开心（创造使用 vec2var 的需要）并继续最困难的情况（具有多个结构中断的系列的协整）；即，我们将继续“小便的人雄心勃勃地钻墙”的逻辑。
将 vars::vec2var 扩展到 causfinder::vec2var_ykJohEsbInc 以处理具有相关干预假人的“多个结构中断”情况下的转换。上面的 JMN2000 应用显示协整等级 r 不在 [1,4-1]=[1,3] 范围内。尽管为了争论，假设 JMN2000 CV 在上面导致 r=1。
因此，要将受限 VECM 转换为受限 VAR（在 multiple=2 结构中断下），请应用：
vec2var_ykJohEsbInc(ykJohEsbInc(mydata, type="trace", ecdet="zamanda2yk", K=3, spec="longrun", dumvar=dummymatrix2SB[,c(-1,-2,-3)]),r=1)
这些结果：
Deterministic coefficients (detcoeffs):
e prod rw U
kesme 22.6612871 -0.215892151 32.0610121 -9.26649249 #(const)
zyonsemesi 0.2505164 -0.009900004 0.3503561 -0.10494714 #(trend)
zy*D2t_3 0.2238060 -0.008844454 0.3130007 -0.09375756
zy*D3t_3 -0.1234803 0.004879743 -0.1726916 0.05172878
$deterministic
kesme zyonsemesi zy*D2t_3 zy*D3t_3 D2t.3l D3t.3l
e 22.6612871 0.250516390 0.223806048 -0.123480327 -8.8012612 5.3052074
prod -0.2158922 -0.009900004 -0.008844454 0.004879743 -0.1157137 -0.3396206
rw 32.0610121 0.350356063 0.313000702 -0.172691620 -12.5838458 7.2201840
U -9.2664925 -0.104947142 -0.093757559 0.051728781 3.5836119 -2.2921099
I2t I2t.1l I2t.2l I3t I3t.1l I3t.2l
e -0.2584379 0.08470453 0.2102661 -0.51366831 -1.0110891 -2.08728944
prod 0.3013044 0.25103445 -0.8640467 0.08804425 -0.2362783 -0.05606892
rw -0.5838161 0.28400182 1.2073483 -0.67760848 -2.2650094 -0.70586316
U 0.1305258 0.03559119 0.1476985 0.14614290 0.6847273 1.27469940
$A
$A$A1
e.1g prod.1g rw.1g U.1g
e 1.4817704 0.1771082 -0.2274936 0.2332402
prod -0.1605790 1.1846699 0.0406294 -0.9398689
rw -0.8366449 -0.1910611 0.9774874 0.4667430
U -0.4245817 -0.1498295 0.1226085 0.7557885
$A$A2
e.2g prod.2g rw.2g U.2g
e -0.8441175 -0.04277845 0.01128282 -0.01896916
prod -0.3909984 -0.25960184 -0.20426749 0.79420691
rw 1.4181448 -0.03659278 -0.12240211 -0.06579174
U 0.4299422 0.09070905 0.04935195 -0.12691817
$A$A3
e.3g prod.3g rw.3g U.3g
e 0.40149641+0i -0.07067529+0i -0.008175418-0i 0.2286283+0i
prod 0.55003024+0i 0.07241639+0i 0.172505474-0i 0.1281593+0i
rw -0.52674826+0i 0.31667695+0i -0.168897398-0i 0.2184591+0i
U -0.02176108-0i 0.03245409-0i -0.077959841+0i 0.1855889-0i
所以，现在，检查根：
print(roots(vec2var_ykJohEsbInc(ykJohEsbInc(mydata, type="trace", ecdet="zamanda2yk", K=3, spec="longrun", dumvar=dummymatrix2SB[,c(-1,-2,-3)]),r=1), modulus=TRUE))
这个结果是“Please provide an object of class 'varest', generated by 'VAR()'.”，因为vars::roots没有被扩展，因为：我们不需要这个扩展！正如我之前所说，即使在 VECM 受限的情况下，也会从 VAR 端检查稳定性。您必须逐行阅读 Joyeux2007 才能看到这一点。
我将彻底提供上述功能的输出（打印屏幕）以进一步澄清。
出于教学原因，我还将向vars::root 写扩展名。