VaR <- function(y, tau=.01, win=250)
{
################################################################################
#                                                                              #
#             VaR Function by Breno de Andrade Pinheiro Néri                   #
#                                                                              #
#            brenoneri@gmail.com       www.fgv.br/aluno/bneri                  #
#                                                                              #
#                                                 September, 2005 - Version 13 #
#                                                                              #
################################################################################
#                                                                              #
# This code is a Monte Carlo simulation to compare seven different             #
# methodologies for estimating 1-day-ahead-forecast VaR time series, namely:   #
# 1- RiskMetrics model developed by J.P. Morgan, 1996;                         #
# 2- GARCH(1, 1) with normal errors;                                           #
# 3- APARCH(1, 1) with skewed-t errors;                                        #
# 4- ARCH(1) Quantile VaR, developed by Wu and Xiao, 2002;                     #
# 5- ARCH(0) Quantile VaR;                                                     #
# 6- QAR(1)-VaR (Quantile AutoRegressive VaR), introduced by me and my thesis  #
#    advisor, Luiz Renato Lima (see Lima and Neri, 2005);                      #
# 7- QAR(0)-VaR.                                                               #
#                                                                              #
# You enter a time series (y), a quantile (tau) and a temporal window size     #
# (win) and the code returns a matrix VaR series, one column for each          #
# methodology.                                                                 #
#                                                                              #
# Note that each forecasts are estimated based on the last 'win' observations. #
# Hence, if length(y) < win, you will have an error. If you do not have enough #
# observations, set up a lower temporal window size (win).                     #
# Note also that length(VaR) = length(y) - win + 1.                            #
#                                                                              #
# Due to the two constrained nonlinear optimization for each observation, this #
# code is extremely computational intensive!                                   #
#                                                                              #
# Defaults: tau=.01 and win=250, 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
t <- length(y)
if(t<win) error("length(y)<win")
library(quantreg)
library(skewt)
VaR <- matrix(0, nrow=t-win+1, ncol=7)
dimnames(VaR)[[2]] <- c("VaR1", "VaR2", "VaR3", "VaR4", "VaR5", "VaR6", "VaR7")
quant.norm <- qnorm(tau)

# GARCH11
garch <- function(p)
{
    sigma2 <- e[1]^2
    for(i in 2:(length(e)+1)) sigma2[i] <- p[1] + p[2]*e[i-1]^2 + p[3]*sigma2[i-1]
    sigma <- sqrt(sigma2)
    return(sigma)
}

# APARCH11
aparch <- function(p)
{
    sigmad <- abs(e[1])^p[5]
    for(i in 2:(length(e)+1)) sigmad[i] <- p[1] + p[2]*(abs(e[i-1])-p[3]*e[i-1])^p[5] + p[4]*sigmad[i-1]
    sigma <- sigmad^(1/p[5])
    return(sigma)
}

# Gaussian Joint Likelihood
L.G <- function(p)
{
    sigma <- garch(p)[1:length(e)]
    L <- -sum(log(dnorm(e/sigma)/sigma))
    return(L)
}

# Skewed Student-t Joint Likelihood
L.ST <- function(p)
{
    sigma <- aparch(p[1:5])[1:length(e)]
    L <- -sum(log(dskt(e/sigma, p[6], p[7])/sigma))
    return(L)
}

# Moving the Temporal Window
for(i in 1:(t-win+1))
{
    # Initializing
    y.win <- y[(i+1):(i+win-1)]
    lag.y.win <- y[i:(i+win-2)]
    V <- cbind(1, y.win[win-1])
    ols <- lm(y.win~lag.y.win)
    mi <- V%*%coef(ols)
    e <- residuals(ols)
    g <- constrOptim(c(.1,1e-2,.9),L.G,NULL,rbind(c(1,0,0),c(0,1,0),c(0,0,1)),c(0,0,0),control=list(reltol=1e-3))[[1]]
    skst <- constrOptim(c(.1,1e-2,0,.9,1,5,1),L.ST,NULL,rbind(c(1,0,0,0,0,0,0),c(0,1,0,0,0,0,0),c(0,0,1,0,0,0,0),c(0,0,-1,0,0,0,0),c(0,0,0,1,0,0,0),c(0,0,0,0,1,0,0),c(0,0,0,0,0,1,0)),c(0,0,-1,-1,0,0,2),control=list(reltol=1e-2))[[1]]
    sigma2.RM <- e[1]^2
    for(j in 2:win) sigma2.RM[j] <- .06*e[j-1]^2 + .94*sigma2.RM[j-1]
    sigma.g <- garch(g)
    sigma.st <- aparch(skst[1:5])

    # Estimating 1-day-ahead VaR Forecast for this temporal window
    VaR[i, 1] <- -(mi + sqrt(sigma2.RM[win])*quant.norm)
    VaR[i, 2] <- -(mi + sigma.g[length(sigma.g)]*quant.norm)
    VaR[i, 3] <- -(mi + sigma.st[length(sigma.st)]*qskt(tau, skst[6], skst[7]))
    VaR[i, 4] <- -(mi + cbind(1, abs(e[win-1]))%*%coef(rq(e[2:(win-1)]~abs(e[1:(win-2)]), tau)))
    VaR[i, 5] <- -(mi + coef(rq(e~1, tau)))
    VaR[i, 6] <- -V%*%coef(rq(y.win~lag.y.win, tau))
    VaR[i, 7] <- -coef(rq(y[i:(i+win-1)]~1, tau))
}

# Reporting
return(VaR)
}