R/rgarch-riskreporting.R

In rgarch: Flexible GARCH modelling in R

#################################################################################
##
##   R package rgarch by Alexios Ghalanos Copyright (C) 2008, 2009, 2010, 2011
##   This file is part of the R package rgarch.
##
##   The R package rgarch is free software: you can redistribute it and/or modify
##   it under the terms of the GNU General Public License as published by
##   the Free Software Foundation, either version 3 of the License, or
##   (at your option) any later version.
##
##   The R package rgarch is distributed in the hope that it will be useful,
##   but WITHOUT ANY WARRANTY; without even the implied warranty of
##   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##   GNU General Public License for more details.
##
#################################################################################


.VaRplot = function(varname , p, actual, dates, VaR)
{
	
	y.actual = actual
	y.VaR = VaR
	x.dates = dates
	zd = .makedate(as.character(dates))
	if(zd$status == 0) dates = 1:length(dates) else dates = zd$dates
	title = paste("Daily Returns & Value-at-Risk Exceedances\n","(Series: ", varname,", alpha=", p,")",sep="")
	if(zd$status){
		plot(x = as.Date(dates, format = zd$dformat), y = y.actual, type = "n",
				main = title, ylab = "Daily Log Return", xlab = "time", 
				ylim = c(min(y.actual, y.VaR), max(y.actual, y.VaR)))
	} else{
		plot(as.numeric(dates), y.actual, type = "n",
				main = title, ylab = "Daily Log Return", xlab = "time", 
				ylim = c(min(y.actual, y.VaR), max(y.actual, y.VaR)))
	}
	abline(h = 0, col = 2, lty = 2)
	sel  =  which(y.actual>0)
	points(dates[sel], y.actual[sel], pch = 18, col = "green")
	sel  =  which(y.actual<0)
	points(dates[sel], y.actual[sel], pch = 18, col = "orange")
	sel  =  which(y.actual<y.VaR)
	points(dates[sel], y.actual[sel], pch = 18, col = "red")
	if(zd$status){
		lines(x = as.Date(dates, format = "%Y-%m-%d"), y = y.VaR, lwd = 2, col = "black")
	} else{
		lines(as.numeric(dates), y.VaR, lwd = 2, col = "black")
	}
	legend("topleft", max(actual),c("return >= 0","VaR <= return < 0","return < VaR","VaR"),
			col=c("green","orange","red","black"), cex=0.75,
			pch = c(18,18,18,-1), lty=c(0,0,0,1), lwd=c(0,0,0,2), bty = "n")
}



.VaRreport = function(varname, garchmodel, distribution, p, actual, VaR, conf.level = 0.95)
{
	actual<-as.matrix(actual)
	VaR<-as.matrix(VaR)
	result <- LR.cc.test(p=p,actual=actual,VaR=VaR,conf.level=conf.level)
	cat("VaR Backtest Report\n")
	cat("===========================================\n")
	cat(paste("Model:\t\t\t","",garchmodel,"-",distribution,"\n",sep=""))
	cat(paste("Backtest Length:\t",result$TN,"\n",sep=""))
	cat(paste("Data:\t\t\t","",varname,"\n",sep=""))
	cat("\n==========================================\n")		
	cat(paste("alpha:\t\t\t\t",round(100*p,1),"%\n",sep=""))
	cat(paste("Expected Exceed:\t",round(p*result$TN,1),"\n",sep=""))
	cat(paste("Actual VaR Exceed:\t",result$N,"\n",sep=""))
	cat(paste("Actual %:\t\t\t",round(100*result$N/result$TN,1),"%\n",sep=""))
	cat("\nUnconditional Coverage (Kupiec)")
	cat("\nNull-Hypothesis:\tCorrect Exceedances\n")
	cat(paste("LR.uc Statistic:\t",round(result$stat.uc,3),"\n",sep=""))
	cat(paste("LR.uc Critical:\t\t",round(result$crit.val.uc,3),"\n",sep=""))
	cat(paste("LR.uc p-value:\t\t",round(result$p.value.uc,3),"\n",sep=""))
	cat(paste("Reject Null:\t\t",ifelse(round(result$p.value.uc,3)< (1-conf.level), "YES","NO"),"\n",sep=""))
	cat("\nConditional Coverage (Christoffersen)")
	cat("\nNull-Hypothesis:\tCorrect Exceedances & Independence of Failures\n")
	cat(paste("LR.cc Statistic:\t",round(result$stat.cc,3),"\n",sep=""))
	cat(paste("LR.cc Critical:\t\t",round(result$crit.val.cc,3),"\n",sep=""))
	cat(paste("LR.cc p-value:\t\t",round(result$p.value.cc,3),"\n",sep=""))
	cat(paste("Reject Null:\t\t",ifelse(result$reject,"YES","NO"),"\n",sep=""))
}


########################################################################
# Available in a number of locations/textbooks.
# This code originally presented in a webcast by Guy Yollin Feb 2006 (from Insightful)
# Functions to perform Hypothesis test
# on VaR models based on # of exceedances
# calc LR.uc statistic

LR.uc = function(p, TN, N)
{
	stat.uc = -2*log((1-p)^(TN-N)*p^N )+2*log((1-N/TN)^(TN-N)*(N/TN)^N)
	
	return(stat.uc)
}

# calc LR.cc statistic

LR.cc = 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*log((1-pe)^(T00+T10)*pe^(T01+T11))+2*log((1-pi0)^T00*pi0^T01*(1-pi1)^T10*pi1^T11)
	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))
}

# perform LR.cc confidence test
LR.cc.test = function(p, actual, VaR, conf.level = 0.95)
{
	result = LR.cc(p = p, actual = actual,VaR = VaR)
	crit.val.uc = qchisq(conf.level, df = 1)
	crit.val.cc = qchisq(conf.level, df = 2)
	p.value.cc = 1-pchisq(result$stat.cc, df = 2)
	p.value.uc = 1-pchisq(result$stat.uc, df = 1)
	reject = ifelse(p.value.cc<1-conf.level, TRUE, FALSE)
	return(list(stat.cc = result$stat.cc, stat.uc = result$stat.uc,
					p.value.cc = p.value.cc, p.value.uc = p.value.uc,
					conf.level=conf.level, reject = reject,
					N = result$N, TN = result$TN, crit.val.uc = crit.val.uc,
					crit.val.cc = crit.val.cc))
}