Applied Portfolio VaR Decomposition

#*****************************************************************
# Load historical data
#*****************************************************************
library(SIT)
load.packages('quantmod')
tickers = spl('AAPL,DIS,IBM,JNJ,KO,NKE,TXN')

data <- new.env()
	data$symbolnames = tickers
getSymbols(tickers, src = 'yahoo', from = '1970-01-01', env = data, auto.assign = T)
for(i in data$symbolnames) data[[i]] = adjustOHLC(data[[i]], use.Adjusted=T)
bt.prep(data, align='remove.na')

#*****************************************************************
# Setup
#*****************************************************************
prices = last(data$prices, 3*252)
	n = ncol(prices)
 last.price = last(prices)
 last.date = index(last(prices))
rets = prices / mlag(prices) - 1
#rets = diff(log(prices))
		rets = coredata(na.omit(rets))

shares = list(AAPL = 530, DIS = 1050, IBM = 150, JNJ = 12700, KO = 3100, NKE = 3500, TXN = 19650)
	shares = unlist(shares[colnames(prices)])

capital = sum(last.price * shares)
exposure = last.price * shares
weight = exposure / capital

#*****************************************************************
# Compute VaR
#*****************************************************************
compute.VaR = function(rets = NULL, cov.mat = cov(rets), 
	weight = rep(1/nrow(cov.mat), nrow(cov.mat)), 
 	capital = 1, con.levl = 0.95, date = Sys.Date()) 
 {	
	z.stat = qnorm(con.levl)
 	weight = as.vector(weight)
 
	asset.vol = abs(weight) * sqrt(diag(cov.mat))
	asset.VaR = z.stat * asset.vol
	undiversified.VaR = sum(asset.VaR)
 
	port.vol = sqrt(t(weight) %*% cov.mat %*% weight)[1]
 	port.VaR = z.stat * port.vol
 	beta = (cov.mat %*% weight) / port.vol^2
 
 	marginal.VaR = z.stat * (cov.mat %*% weight) / port.vol
	component.VaR = marginal.VaR * weight
 
 list(con.levl = con.levl,
 	date = date,
  	n = nrow(cov.mat),
  	capital = capital,
  	weight = weight,
  	tickers = colnames(cov.mat),

		 undiversified.VaR = capital * undiversified.VaR,
  	port.vol = sqrt(252) * port.vol,
  	port.VaR = capital * port.VaR,
  	asset.VaR = capital * asset.VaR,
  	component.VaR = capital * component.VaR,
  	component.VaR.percent = component.VaR / port.VaR,
  	marginal.VaR = marginal.VaR
 )
}

#*****************************************************************
# Helper function to display VaR
#*****************************************************************
summary.VaR = function(data.VaR, capital = 1) {
	plot.data = list(
	'Portfolio Risk Report as of' = format(data.VaR$date, '%d-%b-%Y'),
 	'VaR estimation at' = paste(to.percent(data.VaR$con.levl), 'confidence level'),
 	'Number of assets:' = data.VaR$n,
 	'Current exposure:' = to.cash(capital * data.VaR$capital),
 	'Portfolio VaR (undiversified):' = to.cash(capital * data.VaR$undiversified.VaR),
 	'Portfolio VaR (diversified):' = to.cash(capital * data.VaR$port.VaR),
 	'Component VaR [individual VaR]' = ''
 	)
 
	plot.data1 = cbind( 
 	to.percent(data.VaR$component.VaR.percent),
  	to.cash(capital * data.VaR$component.VaR),
  	to.cash(capital * data.VaR$asset.VaR)
  	)
 
	rbind(cbind(sapply(plot.data,identity),'',''), plot.data1)
}


#*****************************************************************
# Compute VaR
#*****************************************************************
base.VaR = compute.VaR(rets, weight = weight, capital = capital, date = last.date)
print(summary.VaR(base.VaR))

# compute min risk portfolio
ia = create.ia(rets)
	ia$cov = cov(rets)
load.packages('quadprog')
weight = min.var.portfolio(ia, create.basic.constraints(n))

min.risk.VaR = compute.VaR(rets, weight = weight, capital = capital, date = last.date)
print(summary.VaR(min.risk.VaR))

#*****************************************************************
# Helper function to display VaR comparison for two portfolios
#*****************************************************************
compare.VaR = function(base.VaR, test.VaR) {
	plot.data = list(
	'Original position' = to.percent(base.VaR$weight),
	'Marginal VaR' = to.cash(base.VaR$marginal.VaR,5),
	'New position' = to.percent(test.VaR$weight),
	'Marginal VaR' = to.cash(test.VaR$marginal.VaR,5)
	)

	plot.data1 = list(
	'Portfolio VaR' = to.cash(c(base.VaR$port.VaR, test.VaR$port.VaR)),
	' ' = c('', to.percent(test.VaR$port.VaR / base.VaR$port.VaR - 1)),
	'Annulized Vol' = to.percent(c(base.VaR$port.vol, test.VaR$port.vol))
	)
 
	plot.data1 = sapply(plot.data1,identity)
 	plot.data1 = col.name2row(plot.data1)
	
	rbind(cbind(base.VaR$ticker, sapply(plot.data,identity)),
		cbind(plot.data1[1,], plot.data1[2,],'',plot.data1[3,],'')
		)
}


print('Risk-Minimising Position scenario:')

print(compare.VaR(base.VaR, min.risk.VaR))