1## 2## bh5lrtest 3## 4bh5lrtest <- function (z, H, r) 5{ 6 if (!(class(z) == "ca.jo")) { 7 stop("\nPlease, provide object of class 'ca.jo' as 'z'.\n") 8 } 9 if (r >= z@P || r < 1) { 10 stop("\nCount of cointegrating relationships is out of allowable range.\n") 11 } 12 if (z@ecdet == "none") { 13 P <- z@P 14 } else { 15 P <- z@P + 1 16 } 17 r <- as.integer(r) 18 H <- as.matrix(H) 19 if (!(nrow(H) == P)) { 20 stop("\nRow number of 'H' is unequal to VAR order.\n") 21 } 22 if(ncol(H) > r - 1){ 23 stop("\nToo many columns in H for provided r.\n") 24 } 25 r1 <- ncol(H) 26 r2 <- r - r1 27 type <- "Estimation and testing under partly known beta" 28 N <- nrow(z@Z0) 29 M00 <- crossprod(z@Z0)/N 30 M11 <- crossprod(z@Z1)/N 31 MKK <- crossprod(z@ZK)/N 32 M01 <- crossprod(z@Z0, z@Z1)/N 33 M0K <- crossprod(z@Z0, z@ZK)/N 34 MK0 <- crossprod(z@ZK, z@Z0)/N 35 M10 <- crossprod(z@Z1, z@Z0)/N 36 M1K <- crossprod(z@Z1, z@ZK)/N 37 MK1 <- crossprod(z@ZK, z@Z1)/N 38 M11inv <- solve(M11) 39 S00 <- M00 - M01 %*% M11inv %*% M10 40 S0K <- M0K - M01 %*% M11inv %*% M1K 41 SK0 <- MK0 - MK1 %*% M11inv %*% M10 42 SKK <- MKK - MK1 %*% M11inv %*% M1K 43 S00.h <- S00 - S0K%*%H%*%solve(t(H)%*%SKK%*%H)%*%t(H)%*%SK0 44 S0K.h <- S0K - S0K%*%H%*%solve(t(H)%*%SKK%*%H)%*%t(H)%*%SKK 45 SK0.h <- SK0 - SKK%*%H%*%solve(t(H)%*%SKK%*%H)%*%t(H)%*%SK0 46 SKK.h <- SKK - SKK%*%H%*%solve(t(H)%*%SKK%*%H)%*%t(H)%*%SKK 47 valeigen <- eigen(SKK.h) 48 C <- valeigen$vectors[ ,1:(P-r1)]%*%diag(1/sqrt(valeigen$values[1:(P-r1)])) 49 valeigen <- eigen(t(C)%*%SK0.h%*%solve(S00.h)%*%S0K.h%*%C) 50 PSI <- C%*%valeigen$vectors[,1:r2] 51 Dtemp <- chol(t(H)%*%SKK%*%H, pivot = TRUE) 52 pivot <- attr(Dtemp, "pivot") 53 oo <- order(pivot) 54 D <- t(Dtemp[, oo]) 55 Dinv <- solve(D) 56 rho <- eigen(Dinv %*% t(H) %*% SK0 %*% solve(S00) %*% S0K %*% H %*% t(Dinv)) 57 Vorg <- cbind(H, PSI) 58 idx <- ncol(PSI) 59 PSI <- sapply(1:idx, function(x) PSI[, x]/PSI[1, x]) 60 V <- cbind(H, PSI) 61 W <- S0K %*% V %*% solve(t(V) %*% SKK %*% V) 62 PI <- W %*% t(V) 63 DELTA <- S00 - S0K %*% V %*% solve(t(V) %*% SKK %*% V) %*% t(V) %*% SK0 64 GAMMA <- M01 %*% M11inv - PI %*% MK1 %*% M11inv 65 lambda.res <- valeigen$values 66 lambda <- z@lambda 67 teststat <- N *(sum(log(1 - rho$values[1:r1])) + sum(log(1-lambda.res[1:r2])) - sum(log(1 - lambda[1:r]))) 68 df <- (P - r)*r1 69 pval <- c(1 - pchisq(teststat, df), df) 70 new("cajo.test", Z0 = z@Z0, Z1 = z@Z1, ZK = z@ZK, ecdet = z@ecdet, H = H, A = NULL, B = NULL, type = type, teststat = teststat, pval = pval, lambda = lambda.res, Vorg = Vorg, V = V, W = W, PI = PI, DELTA = DELTA, DELTA.bb = NULL, DELTA.ab = NULL, DELTA.aa.b = NULL, GAMMA = GAMMA, test.name = "Johansen-Procedure") 71} 72