##############################################################
#### Examples 6.9, 6.10, 6.11, 6.12
##############################################################

# This file contains complementary code for the credit default swap
# examples 6.9-6.12

# Suppose the investor invests the amount 100/k in each of k
# independent defaultable bonds. We will study VaR and ES as
# for the resulting portfolios, as functions of k. 
#
# Recall that each bond costs $97 and pays $100 in 6m if not 
# defaulted and the default probability is 0.02. 

###
# The VaR_0.05 as a function of k
###

Value.at.Risk <- c()
for(k in 1:1000){
	Value.at.Risk <- c(Value.at.Risk,100*(1-(100/97)*(1-qbinom(0.95,k,0.02)/k)))
	}

plot(Value.at.Risk,pch=20)



# Illustrate the cdf of the portfolio loss (solid), the level 0.95 (dashed), and the VaR (red dot), for different values of k
	
k <- 1
x <- c(0:k)
cdf <- pbinom(x,k,0.02)
L <- 100*(1-(100/97)*(1-x/k))
plot(c(min(L)-1, L),c(0,cdf),type="s")
points(Value.at.Risk[k],0, pch=19, col="red")
lines(c(-5,100),c(0.95,0.95),lty="dashed")


k <- 2
x <- c(0:k)
cdf <- pbinom(x,k,0.02)
L <- 100*(1-(100/97)*(1-x/k))
plot(c(min(L)-1, L),c(0,cdf),type="s")
points(Value.at.Risk[k],0, pch=19, col="red")
lines(c(-5,100),c(0.95,0.95),lty="dashed")

k <- 5
x <- c(0:k)
cdf <- pbinom(x,k,0.02)
L <- 100*(1-(100/97)*(1-x/k))
plot(c(min(L)-1, L),c(0,cdf),type="s")
points(Value.at.Risk[k],0, pch=19, col="red")
lines(c(-5,100),c(0.95,0.95),lty="dashed")

k <- 20
x <- c(0:k)
cdf <- pbinom(x,k,0.02)
L <- 100*(1-(100/97)*(1-x/k))
plot(c(min(L)-1, L),c(0,cdf),type="s")
points(Value.at.Risk[k],0, pch=19, col="red")
lines(c(-5,100),c(0.95,0.95),lty="dashed")

k <- 100
x <- c(0:k)
cdf <- pbinom(x,k,0.02)
L <- 100*(1-(100/97)*(1-x/k))
plot(c(min(L)-1, L),c(0,cdf),type="s")
points(Value.at.Risk[k],0, pch=19, col="red")
lines(c(-5,100),c(0.95,0.95),lty="dashed")

# Animation of the VaR as a function of k together with the loss cdf
# Note that the VaR sometimes jumps up when k increases, 
# illustrating the lack of subadditivity of VaR

for(k in 1:100){
	x <- c(0:k)
	cdf <- pbinom(x,k,0.02)
	L <- 100*(1-(100/97)*(1-x/k))
	plot(c(min(L)-1, L),c(0,cdf),type="s")
	points(Value.at.Risk[k],0, pch=19, col="red")
	lines(c(-5,100),c(0.95,0.95),lty="dashed")
	if(k < 10)	Sys.sleep(1)
	else if(k < 20) Sys.sleep(0.5)
	else Sys.sleep(0.3)
	}


# Illustrate the pdf of the portfolio loss (solid) and the VaR (red dot), for different values of k

k <- 2
x <- c(0:k)
pdf <- dbinom(x,k,0.02)
L <- 100*(1-(100/97)*(1-x/k))
plot(L,pdf,type="h")
points(L,pdf, pch=19)
points(Value.at.Risk[k],0, pch=19, col="red")

k <- 5
x <- c(0:k)
pdf <- dbinom(x,k,0.02)
L <- 100*(1-(100/97)*(1-x/k))
plot(L,pdf,type="h")
points(L,pdf, pch=19)
points(Value.at.Risk[k],0, pch=19, col="red")

k <- 20
x <- c(0:k)
pdf <- dbinom(x,k,0.02)
L <- 100*(1-(100/97)*(1-x/k))
plot(L,pdf,type="h")
points(L,pdf, pch=19)
points(Value.at.Risk[k],0, pch=19, col="red")

k <- 100
x <- c(0:k)
pdf <- dbinom(x,k,0.02)
L <- 100*(1-(100/97)*(1-x/k))
plot(L,pdf,type="h")
points(L,pdf, pch=19)
points(Value.at.Risk[k],0, pch=19, col="red")

# Animation of the VaR as a function of k together with the loss pdf
# Note that the VaR sometimes jumps up when k increases, 
# illustrating the lack of subadditivity of VaR


for(k in 1:100){
	x <- c(0:k)
	pdf <- dbinom(x,k,0.02)
	L <- 100*(1-(100/97)*(1-x/k))
	plot(L,pdf,type="h")
	points(L,pdf, pch=19)
	points(Value.at.Risk[k],0, pch=19, col="red")
	if(k < 10)	Sys.sleep(1)
	else if(k < 20) Sys.sleep(0.5)
	else Sys.sleep(0.3)
	}
	
###
# Normal approximation
###

# When using a Normal approximation to the Binomial the VaR
# is much smoother and the lack of subadditivity does not show up
# for the approximation

kvec <- c(1:1000)
VaR.normal.approx <- -100/97+100^2/97*sqrt(0.02*0.98/kvec)*qnorm(0.95)
plot(Value.at.Risk,xlim=c(0,1000), type="l")
lines(kvec,VaR.normal.approx, lty="dotted")	

###
# Expected Shortfall
###

# Compare VaR and ES as a function of k
#
# Note that ES behaves smoother and is subadditive


Expected.Shortfall <- c()
for(k in 1:1000){
	m <- qbinom(0.95,k,0.02)
	integral <- m*(pbinom(m,k,0.02)-0.95)+sum(c(m+1:k)*dbinom(c(m+1:k),k,0.02))
	Expected.Shortfall <- c(Expected.Shortfall, 100*(1-(100/97)*(1-(1/(0.05*k)*integral))))
	}
	
plot(Value.at.Risk,type="l")
lines(Expected.Shortfall, lty = "dotted")