rm(list=ls())
# Load some reuired library
library(SIT)
library(quantmod)
library(quadprog)
# 1. Load up aa.test.hist.capitalization()
aa.test.hist.capitalization <- function()
{
symbols = spl('Australia Canada France Germany Japan United Kingdom United States', '\t')
data =
'1988 138.0 242.0 245.0 252.0 3910.0 771.0 2790.0
1989 141.0 291.0 365.0 365.0 4390.0 827.0 3510.0
1990 109.0 242.0 314.0 355.0 2920.0 849.0 3060.0
1991 149.0 267.0 348.0 393.0 3130.0 988.0 4090.0
1992 145.0 243.0 351.0 348.0 2400.0 927.0 4490.0
1993 204.9 326.5 456.1 463.5 2999.8 1151.6 5136.2
1994 218.9 315.0 451.3 470.5 3719.9 1210.2 5067.0
1995 245.2 366.3 522.1 577.4 3667.3 1407.7 6857.6
1996 312.0 486.3 591.1 671.0 3088.9 1740.2 8484.4
1997 295.8 567.6 674.4 825.2 2216.7 1996.2 11308.8
1998 328.9 543.4 991.5 1094.0 2495.8 2374.3 13451.4
1999 427.7 800.9 1475.5 1432.2 4546.9 2933.3 16635.1
2000 372.8 841.4 1446.6 1270.2 3157.2 2577.0 15104.0
2001 375.1 700.8 1174.4 1071.7 2251.8 2164.7 13854.6
2002 378.8 575.3 967.0 691.1 2126.1 1864.3 11098.1
2003 585.5 894.0 1355.9 1079.0 3040.7 2460.1 14266.3
2004 776.4 1177.5 1559.1 1194.5 3678.3 2815.9 16323.7
2005 804.1 1480.9 1758.7 1221.3 4736.5 3058.2 16970.9
2006 1095.9 1700.7 2428.6 1637.8 4726.3 3794.3 19425.9
2007 1298.4 2186.6 2771.2 2105.5 4453.5 3858.5 19947.3
2008 675.6 1002.2 1492.3 1108.0 3220.5 1852.0 11737.6
2009 1258.5 1681.0 1972.0 1297.6 3377.9 2796.4 15077.3
2010 1454.5 2160.2 1926.5 1429.7 4099.6 3107.0 17139.0'
#2011 1198.2 1912.9 1554.0 1184.5 3325.4 2932.2 15640.7
#2012 1386.9 2060.0 1808.2 1486.3 3478.8 3291.6 18668.3
#2013 1366.0 2113.8 2301.1 1936.1 4543.2 3946.9 24034.9
#2014 1288.7 2095.4 2085.9 1738.5 4378.0 3570.9 26330.6
#2015 1187.1 1593.4 2088.3 1715.8 4894.9 0 25067.5
#2016 1268.5 1993.5 2159.0 1716.0 4955.3 0 27352.2
#2017 1508.5 2367.1 2749.3 2262.2 6222.8 0 32120.7
#2018 1262.8 1937.9 2366.0 1755.2 5296.8 0 30436.3
#2019 1487.6 2409.1 0 2098.2 6191.1 0 33890.8
#2020 1720.6 2641.5 0 2284.1 6718.2 0 40719.7'
hist.caps = matrix( as.double(spl( gsub('\n', '\t', data), '\t')),
nrow = len(spl(data, '\n')), byrow=TRUE)
load.packages('quantmod')
symbol.names = symbols
hist.caps = as.xts( hist.caps[,-1] ,
as.Date(paste('1/1/', hist.caps[,1], sep=''), '%d/%m/%Y')
)
colnames(hist.caps) = symbols
return(hist.caps)
}
# 2. Load up aa.test.create.ia.country()
aa.test.create.ia.country <- function(dates = '1990::2010')
{
# load.packages('quantmod,quadprog')
symbols = spl('EWA,EWC,EWQ,EWG,EWJ,EWU,SPY')
symbol.names = spl('Australia, Canada, France, Germany, Japan, UK, USA')
getSymbols(symbols, from = '1980-01-01', auto.assign = TRUE)
hist.prices = merge(EWA,EWC,EWQ,EWG,EWJ,EWU,SPY)
period.ends = endpoints(hist.prices, 'months')
hist.prices = Ad(hist.prices)[period.ends, ]
colnames(hist.prices) = symbol.names
annual.factor = 12
hist.prices = na.omit(hist.prices[dates])
hist.returns = na.omit( ROC(hist.prices, type = 'discrete') )
ia = create.historical.ia(hist.returns, annual.factor)
return(ia)
}
# 3. Load up efficient frontier plotting function:
plot.ef <- function(
ia,
efs,
portfolio.risk.fn = portfolio.risk,
transition.map = TRUE,
layout = NULL
)
{
risk.label = as.character(substitute(portfolio.risk.fn))
n = ia$n
x = match.fun(portfolio.risk.fn)(diag(n), ia)
y = ia$expected.return
xlim = range(c(0, x,
max( sapply(efs, function(x) max(match.fun(portfolio.risk.fn)(x$weight,ia))) )
), na.rm = T)
ylim = range(c(0, y,
min( sapply(efs, function(x) min(portfolio.return(x$weight,ia))) ),
max( sapply(efs, function(x) max(portfolio.return(x$weight,ia))) )
), na.rm = T)
x = 100 * x
y = 100 * y
xlim = 100 * xlim
ylim = 100 * ylim
if( !transition.map ) layout = T
if( is.null(layout) ) layout(1:2)
par(mar = c(4,3,2,1), cex = 0.8)
plot(x, y, xlim = xlim, ylim = ylim,
xlab='', ylab='', main=paste(risk.label, 'vs Return'), col='black')
mtext('Return', side = 2,line = 2, cex = par('cex'))
mtext(risk.label, side = 1,line = 2, cex = par('cex'))
grid();
text(x, y, ia$symbols, col = 'blue', adj = c(1,1), cex = 0.8)
for(i in len(efs):1) {
ef = efs[[ i ]]
x = 100 * match.fun(portfolio.risk.fn)(ef$weight, ia)
y = 100 * ef$return
lines(x, y, col=i)
}
plota.legend(sapply(efs, function(x) x$name), 1:len(efs))
if(transition.map) {
plot.transition.map(efs[[i]]$weight, x, risk.label, efs[[i]]$name)
}
}
#--------------------------------------------------------------------------
# Visualize Market Capitalization History
#--------------------------------------------------------------------------
hist.caps = aa.test.hist.capitalization()
hist.caps.weight = hist.caps/rowSums(hist.caps)
# Plot Transition of Market Cap Weights in time
plot.transition.map(hist.caps.weight, index(hist.caps.weight), xlab='', name='Market Capitalization Weight History')
# Plot History for each Country's Market Cap
layout( matrix(1:9, nrow = 3, byrow=T) )
col = plota.colors(ncol(hist.caps))
for(i in 1:ncol(hist.caps)) {
plota(hist.caps[,i], type='l', lwd=5, col=col[i], main=colnames(hist.caps)[i])
}
# Use reverse optimization to compute the vector of equilibrium returns
bl.compute.eqret <- function(
risk.aversion, # Risk Aversion
cov, # Covariance matrix
cap.weight, # Market Capitalization Weights
risk.free = 0 # Rsik Free Interest Rate
)
{
return( risk.aversion * cov %*% cap.weight + risk.free)
}
#--------------------------------------------------------------------------
# Compute Risk Aversion, prepare Black-Litterman input assumptions
#--------------------------------------------------------------------------
ia = aa.test.create.ia.country()
# compute Risk Aversion
risk.aversion = bl.compute.risk.aversion( ia$hist.returns$` USA` )
# the latest market capitalization weights
cap.weight = last(hist.caps.weight)
# create Black-Litterman input assumptions
ia.bl = ia
ia.bl$expected.return = bl.compute.eqret( risk.aversion, ia$cov, as.vector(cap.weight) )
# Plot market capitalization weights and implied equilibrium returns
layout( matrix(c(1,1,2,3), nrow=2, byrow=T) )
pie(coredata(cap.weight), paste(colnames(cap.weight), round(100*cap.weight), '%'),
main = paste('Country Market Capitalization Weights for', format(index(cap.weight),'%b %Y'))
, col=plota.colors(ia$n))
plot.ia(ia.bl, T)
#
#--------------------------------------------------------------------------
# Create Efficient Frontier(s)
#--------------------------------------------------------------------------
n = ia$n
# -1 <= x.i <= 1
constraints = new.constraints(n, lb = 0, ub = 1)
# SUM x.i = 1
constraints = add.constraints(rep(1, n), 1, type = '=', constraints)
# create efficient frontier(s)
ef.risk = portopt(ia, constraints, 50, 'Historical', equally.spaced.risk = T)
##
## Attaching package: 'corpcor'
## The following object is masked from 'package:SIT':
##
## cov.shrink
##
## Attaching package: 'kernlab'
## The following object is masked from 'package:SIT':
##
## cross
ef.risk.bl = portopt(ia.bl, constraints, 50, 'Black-Litterman', equally.spaced.risk = T)
# Plot multiple Efficient Frontiers and Transition Maps
layout( matrix(1:4, nrow = 2) )
plot.ef(ia, list(ef.risk), portfolio.risk, T, T)
plot.ef(ia.bl, list(ef.risk.bl), portfolio.risk, T, T)
##
bl.compute.posterior <- function(
mu, # Equilibrium returns
cov, # Covariance matrix
pmat=NULL, # Views pick matrix
qmat=NULL, # Views mean vector
tau=0.025 # Measure of uncertainty of the prior estimate of the mean returns
)
{
out = list()
omega = diag(c(1,diag(tau * pmat %*% cov %*% t(pmat))))[-1,-1]
temp = solve(solve(tau * cov) + t(pmat) %*% solve(omega) %*% pmat)
out$cov = cov + temp
out$expected.return = temp %*% (solve(tau * cov) %*% mu + t(pmat) %*% solve(omega) %*% qmat)
return(out)
}
#--------------------------------------------------------------------------
# Create Views
#--------------------------------------------------------------------------
temp = matrix(rep(0, n), nrow = 1)
colnames(temp) = ia$symbols
# Relative View
# Japan will outperform UK by 2%
temp[,' Japan'] = 1
temp[,' UK'] = -1
pmat = temp
qmat = c(0.02)
# Absolute View
# Australia's expected return is 12%
temp[] = 0
# temp[,'Australia'] = 1
temp[,1] = 1
pmat = rbind(pmat, temp)
qmat = c(qmat, 0.12)
# compute posterior distribution parameters
post = bl.compute.posterior(ia.bl$expected.return, ia$cov, pmat, qmat, tau = 0.025 )
# create Black-Litterman input assumptions with Views
ia.bl.view = ia.bl
ia.bl.view$expected.return = post$expected.return
ia.bl.view$cov = post$cov
ia.bl.view$risk = sqrt(diag(ia.bl.view$cov))
# create efficient frontier(s)
ef.risk.bl.view = portopt(ia.bl.view, constraints, 50, 'Black-Litterman + View(s)', equally.spaced.risk = T)
# Plot multiple Efficient Frontiers and Transition Maps
layout( matrix(1:4, nrow = 2) )
plot.ef(ia.bl, list(ef.risk.bl), portfolio.risk, T, T)
plot.ef(ia.bl.view, list(ef.risk.bl.view), portfolio.risk, T, T)
