# valueAtRiskDemo.ssc				Demo script to illustrate value-at-risk calculations
#
# author: Eric Zivot
# created: May 4, 2005
# updated: June 1, 2005
#
# requires Splus 7 and S+FinMetrics 2.0
#
# FM Data used:
# 1. DowJones30 - daily closing prices on Dow Jones 30 stocks from 1/2/91 - 1/2/01
#    
# references
# 1. Elements of Financial Risk Management by Peter F. Christoffersen,
#    Academic Press, 2003.
# 2. Measuring Market Risk by Kevin Dowd, Wiley Finance, 2002
# 3. Value-at-Risk by Philip Jorian, McGraw Hill, 2002
# 4. Market Models by Carol Alexander, Wiley Finance, 2002
# 5.  

# plot prices
seriesPlot(DowJones30, one.plot=F, strip=colIds(DowJones30))

# compute and plot log-returns
DowJones30.ret = getReturns(DowJones30)
seriesPlot(DowJones30.ret, one.plot=F, strip=colIds(DowJones30),
           type="l")

# compute summary statistics and normal qq-plots
summaryStats(DowJones30.ret)
qqPlot(DowJones30.ret, strip.text=colIds(DowJones30.ret), id.n = 0)

# test null hypothesis that returns are normally distributed
normalTest(DowJones30.ret, method="jb")

# compute and plot equally weighted portfolio return
w = rep(1,30)/30
port.ret.ts = rowSums(DowJones30.ret, weights=w)
colIds(port.ret.ts) = "Portfolio"
plot(port.ret.ts, reference.grid=F)
abline(h=0)

# compute summary statistics and normal qq-plot
summaryStats(port.ret.ts)
qqPlot(port.ret.ts, strip.text = colIds(port.ret.ts))

# test for normality
normalTest(port.ret.ts, method="jb")

###########################################################################
# unconditional market risk calculations
###########################################################################

#
# 1. Historical simulation (HS)
#

# compute empirical 1% quantiles (HS VaR)
# 1% VaR for equally weighted portfolio
VaR.hs.01 = quantile(port.ret.ts, probs=0.01)
VaR.hs.01
# compute $VaR for $1M investment
1000000*( exp(VaR.hs.01) - 1 )


# 1% VaR for all assets
apply(DowJones30.ret, 2 , quantile, probs=0.01)
unlist(colQuantiles(DowJones30.ret, probs=0.01))

# compute SE for 5% VaR estimate using the bootstrap
VaR.05.boot = bootstrap(port.ret.ts, statistic=quantile, 
                        args.stat=list(probs=0.05))
summary(VaR.05.boot)
plot(VaR.05.boot)

# compute 1% ETL for equally weighted portfolio
mean(port.ret.ts[port.ret.ts < VaR.hs.01])

# compute 1% ETL for all Dow Jones 30 stocks at once
ETL.hs.01 = function(x) { mean(x[x < quantile(x, probs=0.01)]) }
apply(DowJones30.ret, 2, ETL.hs.01)


#
# 2. Normal approximation
#

# function to compute estimated quantile based on normal distribution
norm.quantile = function(x, p=0.01, mu=NULL) {
	## required arguments:
	## x			timeSeries of returns (simple or continuous)
	## optional arguments:
	## p			scalar probability level
	##	mu			vector of mean values. 
	## value:
	## 			numeric value giving quantile estimate from estimated normal
	##				distribution
	if(is.null(mu))
		q = colMeans(x) + colStdevs(x)*qnorm(p)
	else
		q = mu + colStdevs(x)*qnorm(p)
	q
}

norm.ETL = function(x, VaR) {
	## required arguments:
	## x			timeSeries of returns (simple or continuous)
	## VaR			VaR estimate based on normal approximation. 
	## optional arguments:
	## value:
	## 			numeric value giving ETL estimate from estimated normal
	##			distribution
	mu = colMeans(x)
	sd = colStdevs(x)
	z = (mu - VaR)/sd)
	ETL = mu + sd*(dnorm(z)/pnorm(z))
	ETL
}

# compute 1% quantile for equally weighted portfolio based on 
# normal distribution
norm.quantile(port.ret.ts)
norm.quantile(port.ret.ts, mu=0)

# compute bootstrap confidence interval
VaR.01.boot = bootstrap(port.ret.ts, statistic=norm.quantile)
limits.emp(VaR.01.boot)

# compute 1% quantile using multivariate statistics
mu.hat = colMeans(DowJones30.ret)
Sigma.hat = var(DowJones30.ret)
w = rep(1,30)/30
t(w)%*%mu.hat + sqrt( t(w)%*%Sigma.hat%*%w )*qnorm(.01)

# compute 1% quantile for all Dow Jones 30 stocks based on 
# normal distribution

norm.quantile(DowJones30.ret, p=0.01)
norm.quantile(DowJones30.ret, p=0.01, mu=0)

# compute 1% ETL based on normal distribution

#
# 3. Extreme Value Theory
#

# compute mean excess plots for negative returns
# threshold appears to be near 0.015
par(mfrow=c(2,1))
tmp = meplot(-port.ret.ts)
shape.plot(-port.ret.ts, from=0.9, to=0.98)
par(mfrow=c(1,1))

# estimate parameters of GPD with u = 0.015
gpd.fit = gpd(-port.ret.ts, threshold=0.015)
gpd.fit
tailplot(gpd.fit)

# compute 1% VaR and ETL from fitted GPD
riskmeasures(gpd.fit, p=0.99)
gpd.q(pp=0.99, ci.p = 0.95, plot=F)
gpd.sfall(pp=0.99, ci.p = 0.95, plot=F)

###########################################################################
# conditional market risk calculations
###########################################################################

#
# Risk-Metrics approach (EWMA)
#

# compute EWMA estimates of daily variances and VaR

var.ewma = EWMA(port.ret.ts^2, lambda = 0.94)
plot(var.ewma)

par(mfrow=c(2,1))
	plot(port.ret.ts, reference.grid=F)
	abline(h=0)
	plot(sqrt(var.ewma), reference.grid=F)
par(mfrow=c(1,1))

# conditional Normal VaR based on EWMA estimate of variances
# and zero means
VaR.01.ewma = var.ewma*qnorm(0.01)

#
# GARCH Models
#

#
# GARCH + EVT Models
#

###########################################################################
# Backtesting
###########################################################################

# 
# 6. backtest VaR estimates using simulated normal data
#


sim.norm.ret = rnorm(numRows(port.ret.ts), 
                     mean = colMeans(port.ret.ts),
                     sd = colStdevs(port.ret.ts))
tsplot(sim.norm.ret)
sim.norm.ret.ts = timeSeries(pos=positions(port.ret.ts), 
                             data = sim.norm.ret)
plot(sim.norm.ret.ts, reference.grid=F)
abline(h=0)

# 
# 5. backtest VaR estimates using equally weighted portfolio data
#

# use aggregateSeries() to computing rolling VaR estimates using a 
# a 250 day rolling window

# compute rolling VaR.05 based on normal distribution
VaR.norm.ts = aggregateSeries(sim.norm.ret.ts, FUN=norm.quantile, moving=250,
                              adj=1, colnames="VaR.normal")

# compute rolling VaR.05 based on historical quantiles
hs.quantile = function(x) { quantile(x, probs=0.05) }
VaR.hs.ts = aggregateSeries(sim.norm.ret.ts, FUN=hs.quantile, moving=250,
                            adj=1, colnames="VaR.hs")

# determine VaR violations
VaR.ts = seriesMerge(-sim.norm.ret.ts,-VaR.norm.ts)
colIds(VaR.ts) = c("Portfolio","VaR.normal")

violations.norm.idx = VaR.ts[,1] > VaR.ts[,2] 
violations.norm.idx = ifelse(violations.norm.idx == F, NA, T)
violations.norm.ret = violations.norm.idx*VaR.ts[,1]
n.trials = numRows(VaR.norm.ts)
n.violations.norm = sum(violations.norm.idx, na.rm=T)

# plot VaR violations along with returns
plot(VaR.ts, violations.norm.ret, plot.args=list(type=c("l","l","b")))

plot(-sim.norm.ret.ts, -VaR.norm.ts, main="5% VaR Estimates and Actual Returns",
      reference.grid=F, plot.args=list(col=c(2,3), lwd=c(1,2)))
legend(0,0.07,legend=c("Negative returns","5% Normal VaR"),
       lty=c(1,1),lwd=c(1,2),col=c(2,3))

# Test Null hypothesis that true proportion of VaR violations is 5%
prop.test(n.violations.norm, n.trials, p=0.05)




# 
# 5. backtest VaR estimates using equally weighted portfolio data
#

# use aggregateSeries() to computing rolling VaR estimates using a 
# a 250 day rolling window

# compute rolling VaR.05 based on normal distribution
VaR.norm.ts = aggregateSeries(port.ret.ts, FUN=norm.quantile, moving=250,
                              adj=1, colnames="VaR.normal")

# compute rolling VaR.05 based on historical quantiles
hs.quantile = function(x) { quantile(x, probs=0.05) }
VaR.hs.ts = aggregateSeries(port.ret.ts, FUN=hs.quantile, moving=250,
                            adj=1, colnames="VaR.hs")

# determine VaR violations
VaR.ts = seriesMerge(-port.ret.ts,-VaR.norm.ts)
colIds(VaR.ts) = c("Portfolio","VaR.normal")

violations.norm.idx = VaR.ts[,1] > VaR.ts[,2] 
violations.norm.idx = ifelse(violations.norm.idx == F, NA, T)
violations.norm.ret = violations.norm.idx*VaR.ts[,1]
n.trials = numRows(VaR.norm.ts)
n.violations.norm = sum(violations.norm.idx, na.rm=T)

# plot VaR violations along with returns
plot(VaR.ts, violations.norm.ret, plot.args=list(type=c("l","l","b")))

plot(-port.ret.ts, -VaR.norm.ts, main="5% VaR Estimates and Actual Returns",
      reference.grid=F, plot.args=list(col=c(2,3), lwd=c(1,2)))
legend(0,0.07,legend=c("Negative returns","5% Normal VaR"),
       lty=c(1,1),lwd=c(1,2),col=c(2,3))

# Test Null hypothesis that true proportion of VaR violations is 5%
prop.test(n.violations.norm, n.trials, p=0.05)