Independence test rugarch package
Philippe, The code is correct and returns the right results. The problem is that you are passing xts objects whereas the code cannot deal with this (ifelse and the comparison operators). The help page on VaRTest does state: "actual: A numeric vector of the actual (realized) values.? "VaR: The numeric vector of VaR." Simply transform to numeric (as.numeric(actual) and as.numeric(VaR)) and you should get the correct answer. This was already covered in a previous post to this list. Regards, Alexios
On 26 Apr 2014, at 07:43, philippe <philippe.kappeler at hotmail.com> wrote:
Hi Alexios,
Thanks for the quick response. I downloaded the package from "bitbucket" and
retested. The test now reports values for the UC test, where before I
received NaN for a confidence level 0f 0.95 (p=0.05) and rolling windows
500, 1000, 1500. The rolling window of 2000 was working fine before and
still does. So this is fine and thanks.
But now the CC test returns NaN's at all windows, whereas before window 2000
was working. I copy pasted again the function (.LR.cc) out of getAnywhere
and rerun the function in a sort of "traceback mode". I note that T11 is
still equal to N, which in my eyes is still incorrect as the
which(VaR.ind==1) results show. Please advise me if I misunderstood or
overlooked something.
Thank you very much.
The code slightly updated is again.
###############
require(devtools)
install_bitbucket("rugarch","alexiosg")
ret <- CHF
seq <- 2000 #alternativly 500, 1000, 1500
#Normal Distribution VaR
VaRdistn <- function(ret, prob=.05) {
ans <- -qnorm(prob, mean=0, sd=sd(ret))
signif(ans, digits=7)
}
### VaR
p <- c(0.05, 0.01, 0.005, 0.001) # confidence values for VaR
act <- ret[(seq+1):length(ret)] # actual return
VaRmatrix <- matrix(nrow=length(ret)-seq, ncol=length(p))
for (i in 1:nrow(VaRmatrix) ) {
VaRmatrix[i,]<- sapply(X=p, FUN=VaRdistn,ret=ret[i:(seq+i-1)] )
colnames(VaRmatrix) <- c("0.05", "0.01", "0.005", "0.001")
}
#### RUGARCH package (getAnywhere(.LR.cc))
VaR <- -VaRmatrix[,1] #change column and p value respectively to 2 for 0.01,
3 for 0.005 and 4 for 0.001
p <- 0.05
actual <- act
#function (p, actual, VaR)
#{
VaR.ind = ifelse(actual < VaR, 1, 0)
N = sum(VaR.ind)
TN = length(VaR.ind)
T00 = sum(c(0, ifelse(VaR.ind[2:TN] == 0 & VaR.ind[1:(TN -
1)] == 0, 1, 0)))
T11 = sum(c(0, ifelse(VaR.ind[2:TN] == 1 & VaR.ind[1:(TN -
1)] == 1, 1, 0)))
T01 = sum(c(0, ifelse(VaR.ind[2:TN] == 1 & VaR.ind[1:(TN -
1)] == 0, 1, 0)))
T10 = sum(c(0, ifelse(VaR.ind[2:TN] == 0 & VaR.ind[1:(TN -
1)] == 1, 1, 0)))
T0 = T00 + T01
T1 = T10 + T11
pi0 = T01/T0
pi1 = T11/T1
pe = (T01 + T11)/(T0 + T1)
stat.ind = -2 * ((T00 + T10) * log(1 - pe) + (T01 + T11) *
log(pe)) + 2 * (T00 * log(1 - pi0) + T01 * log(pi0) +
T10 * log(1 - pi1) + T11 * log(pi1))
#stat.uc = .LR.uc(p = p, TN = TN, N = N)
#stat.cc = stat.uc + stat.ind
#return(list(stat.cc = stat.cc, stat.uc = stat.uc, N = N,
# TN = TN))
#}
#### your results out of R-studio console #######
N
[1] 163
T00
[1] 3048
T00
[1] 3048
T11
[1] 163
T01
[1] 0
T10
[1] 0
which(VaR.ind == 1)
[1] 10 55 106 117 151 174 178 186 191 213 227 238 246 302 334 368 373 388 395 438 486 496 497 503 [25] 544 581 606 614 623 637 655 664 670 672 697 722 727 754 767 771 778 815 818 835 864 874 880 901 [49] 999 1000 1025 1031 1086 1100 1151 1159 1168 1184 1189 1212 1231 1272 1316 1468 1502 1526 1596 1598 1638 1639 1651 1656 [73] 1657 1675 1696 1709 1719 1723 1736 1763 1783 1785 1812 1830 1831 1840 1844 1845 1847 1848 1851 1854 1860 1873 1879 1883 [97] 1893 1902 1910 1911 1919 1927 1936 1939 1947 1954 1956 1961 1978 1981 2006 2007 2025 2032 2119 2175 2227 2230 2237 2240 [121] 2250 2260 2266 2272 2275 2292 2293 2299 2319 2328 2330 2351 2357 2370 2387 2405 2425 2443 2475 2478 2508 2516 2520 2522 [145] 2527 2528 2543 2544 2545 2573 2580 2584 2647 2757 2776 2782 2812 2954 2991 3000 3022 3072 3091
stat.ind
[1] NaN
##### my code
conf <- 0.95
p <- c(0.05, 0.01, 0.005, 0.001)
#fill vectors
uc.LR <- cc.LR <- numeric(length(p))
uc.pval <- cc.pval <- numeric(length(p))
uc.statement <- cc.statement <- numeric(length(p))
ind.LR <- ind.pval <- numeric(length(p))
VaR.ind <- matrix(nrow=length(ret)-seq, ncol=length(p))
N <- TN <- T00 <- T11 <- T01 <- T10 <- T0 <- T1 <- pi0 <- pi1 <- pe <-
numeric(length(p))
for(i in 1:length(p)){
uc.LR[i] <- VaRTest(alpha = p[i], act, -VaRmatrix[,i], conf.level =
conf)$uc.LRstat
uc.pval[i] <- VaRTest(alpha = p[i], act, -VaRmatrix[,i], conf.level =
conf)$uc.LRp
VaR.ind[,i] = ifelse(act < -VaRmatrix[,i], 1, 0)
N[i] = sum(VaR.ind[,i])
TN[i] = length(VaR.ind[,i])
T00[i] = sum(c(0, ifelse(VaR.ind[,i][2:TN[i]] == 0 & VaR.ind[,i][1:(TN[i]
- 1)] == 0, 1, 0)))
T11[i] = sum(c(0, ifelse(VaR.ind[,i][2:TN[i]] == 1 & VaR.ind[,i][1:(TN[i]
- 1)] == 1, 1, 0)))
T01[i] = sum(c(0, ifelse(VaR.ind[,i][2:TN[i]] == 1 & VaR.ind[,i][1:(TN[i]
- 1)] == 0, 1, 0)))
T10[i] = sum(c(0, ifelse(VaR.ind[,i][2:TN[i]] == 0 & VaR.ind[,i][1:(TN[i]
- 1)] == 1, 1, 0)))
T0[i] = T00[i] + T01[i]
T1[i] = T10[i] + T11[i]
pi0[i] = T01[i]/T0[i]
pi1[i] = T11[i]/T1[i]
pe[i] = (T01[i] + T11[i])/(T0[i] + T1[i])
ind.LR[i] = -2 * log((1 - pe[i])^(T00[i] + T10[i]) * pe[i]^(T01[i] +
T11[i])) +
2 * log((1 - pi0[i])^T00[i] * pi0[i]^T01[i] * (1 - pi1[i])^T10[i] *
pi1[i]^T11[i])
ind.pval[i] <- 1 - pchisq(ind.LR[i], df = 1)
cc.LR[i] <- VaRTest(alpha = p[i], act, -VaRmatrix[,i], conf.level =
conf)$cc.LRstat
cc.pval[i] <- VaRTest(alpha = p[i], act, -VaRmatrix[,i], conf.level =
conf)$cc.LRp
}
###my results out of R-studio console #####
uc.LR #from rugarch::VaRTest
[1] 0.0360193 11.4504151 10.9552213 6.9767907
ind.LR #should be the difference between uc.LR and cc.LR
[1] 2.597815 1.107278 1.044127 5.630994
cc.LR #from rugarch::VaRTest
[1] NaN NaN NaN NaN -- View this message in context: http://r.789695.n4.nabble.com/Independence-test-rugarch-package-tp4689473p4689499.html Sent from the Rmetrics mailing list archive at Nabble.com.
_______________________________________________ R-SIG-Finance at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-finance -- Subscriber-posting only. If you want to post, subscribe first. -- Also note that this is not the r-help list where general R questions should go.