BackTest <- function(y, VaR, tau=.01)
{
################################################################################
#                                                                              #
#           BackTest Function by Breno de Andrade Pinheiro Néri                #
#                                                                              #
#            brenoneri@gmail.com        www.fgv.br/aluno/bneri                 #
#                                                                              #
#                                                     August, 2005 - Version 2 #
#                                                                              #
################################################################################
#                                                                              #
# This code computes five backtests statistics and p-values, namely:           #
# 1- Kupiec (1995), Unconditional Coverage;                                    #
# 2- Christoffersen (1997), Independence;                                      #
# 3- Christoffersen (1997), Conditional Coverage;                              #
# 4- Engle and Manganelli (2002), Dynamic Quantile;                            #
# 5- Lopez (1998), Magnitude Loss Function.                                    #
#                                                                              #
# You enter a vector of returns (y), a vector of Value-at-Risk (VaR) and a     #
# quantile (tau). VaR variable may also receive a matrix with several VaR      #
# estimates (different methodologies, for example), one VaR for each column    #
# (i.e. observations in the rows).                                             #
#                                                                              #
# Defaults: tau=.01, the most common value, due to regulatory obligation (1996 #
# Amendment to Basle Capital Accord).                                          #
#                                                                              #
################################################################################
#                                                                              #
# You may use this code only if you accept the conditions:                     #
# (1) I am not liable for any problem caused by this code (i.e. you use it at  #
# your own risk);                                                              #
# (2) You must give me credit in your papers where this code has been used.    #
#                                                                              #
################################################################################

# Initializing
VaR <- as.matrix(VaR)
T <- nrow(VaR)-1
n <- ncol(VaR)
if(n-1) VaR <- as.matrix(VaR[1:T,]) else VaR <- as.matrix(VaR[1:T])
y <- y[(length(y)-T+1):length(y)]

# Kupiec (1995), Unconditional Coverage
I <- y< -VaR
T1 <- apply(I, 2, sum)
T0 <- T-T1
P1 <- T1/T
P0 <- 1-P1
LR.uc <- 2*log((P1/tau)^T1*(P0/(1-tau))^T0)
p.value.uc <- 1-pchisq(LR.uc, df=1)

# Christoffersen (1997), Independence
lag.I <- as.matrix(I[1:(T-1),])
I1 <- as.matrix(I[2:T,])
T11 <- apply(lag.I&I1, 2, sum)
T00 <- apply(!lag.I&!I1, 2, sum)
T01 <- T1-T11
T10 <- T0-T00
P01 <- T01/T0
P11 <- T11/T1
P00 <- 1-P01
P10 <- 1-P11
LR.i <- 2*log(P00^T00*P01^T01*P10^T10+P11^T11/(P1^T1*P0^T0))
p.value.i <- 1-pchisq(LR.i, df=1)

# Christoffersen (1997), Conditional Coverage
LR.cc <- LR.uc + LR.i
p.value.cc <- 1-pchisq(LR.cc, df=2)

# Engle and Manganelli (2002), Dynamic Quantile
q <- 3 # modify this parameter to alter the maximum order of captured dependence
t <- T-q
laghit <- array(0, c(t, q))
DQ <- 0
for(j in 1:n)
{
    hit <- I[,j]-tau
    for(i in q:1) laghit[, i] <- hit[(q+1-i):(T-i)]
    hit <- hit[(q+1):T]
    var <- VaR[(q+1):T, j]
    x <- cbind(rep(1, t), laghit, var)
    l <- as.matrix(coef(lm(hit~x-1)))
    DQ[j] <- t(l)%*%t(x)%*%x%*%l/(tau*(1-tau))
}
p.value.DQ <- 1-pchisq(DQ, df=ncol(x))

# Lopez (1998), Magnitude Loss Function
C <- apply((1+(y+VaR)^2)*I, 2, sum)

# Reporting
Result <- list(P1, T1, T, LR.uc, LR.i, LR.cc, DQ, p.value.uc, p.value.i, p.value.cc, p.value.DQ, C)
names(Result) <- c("P1", "T1", "T", "LR.uc", "LR.i", "LR.cc", "DQ", "p.value.uc", "p.value.i", "p.value.cc", "p.value.DQ", "C")
return(Result)
}