VaR <- function(y)
{
################################################################################
#                                                                              #
#             VaR Function by Breno de Andrade Pinheiro Néri                   #
#                                                                              #
#            brenoneri@gmail.com       www.fgv.br/aluno/bneri                  #
#                                                                              #
#                                                  September, 2005 - Version 3 #
#                                                                              #
################################################################################
#                                                                              #
# 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.                                                               #
#                                                                              #
# Observe that the 4th and the 5th methodology are estimated by an Ox code,    #
# called by this R code!                                                       #
#                                                                              #
# You enter a time series (y) 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
tau=.01
win=250
t <- length(y)
if(t<win) error("Catastrophic Error: length(y)<win!")
write(y, "vY.mat", ncolumns=1)
system("oxl VaR.ox")
library(quantreg)
VaR <- cbind(as.matrix(read.table("mVaR.mat")), 0, 0, 0, 0)
dimnames(VaR)[[2]] <- c("VaR1", "VaR2", "VaR3", "VaR4", "VaR5", "VaR6", "VaR7")

# Moving the Temporal Window
for(i in 1:(t-win+1))
{
    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)
    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))
}
return(VaR)
}