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