rm(list=ls())
devtools::install_github('joshuaulrich/xts', force = T)
## Downloading GitHub repo joshuaulrich/xts@HEAD
##
## * checking for file ‘/tmp/Rtmp7bQz2i/remotes87e599057fd/joshuaulrich-xts-62aa765/DESCRIPTION’ ... OK
## * preparing ‘xts’:
## * checking DESCRIPTION meta-information ... OK
## * cleaning src
## * checking for LF line-endings in source and make files and shell scripts
## * checking for empty or unneeded directories
## * looking to see if a ‘data/datalist’ file should be added
## * building ‘xts_0.12.1.2.tar.gz’
## Installing package into '/home/rstudio/R/x86_64-pc-linux-gnu-library/4.1'
## (as 'lib' is unspecified)
devtools::install_github('joshuaulrich/quantmod', force = T)
## Downloading GitHub repo joshuaulrich/quantmod@HEAD
##
## * checking for file ‘/tmp/Rtmp7bQz2i/remotes87e4859404c/joshuaulrich-quantmod-f62ee26/DESCRIPTION’ ... OK
## * preparing ‘quantmod’:
## * checking DESCRIPTION meta-information ... OK
## * installing the package to process help pages
## * saving partial Rd database
## * checking for LF line-endings in source and make files and shell scripts
## * checking for empty or unneeded directories
## * building ‘quantmod_0.4.20.1.tar.gz’
## Installing package into '/home/rstudio/R/x86_64-pc-linux-gnu-library/4.1'
## (as 'lib' is unspecified)
#
library(SIT)
## Loading required package: SIT.date
## Loading required package: quantmod
## Loading required package: xts
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
##
## Attaching package: 'SIT'
## The following object is masked from 'package:TTR':
##
## DVI
## The following object is masked from 'package:base':
##
## close
library(quantmod)
library(quadprog)
#-----------------------------------------------------------------------------------------
# Here we have to load up 3 functions from original sit.gz file. However, it's in conflict
# with the SIT package. So we have to load up these 3 functions first and then load up
# SIT package then.
#-----------------------------------------------------------------------------------------
# 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 1564.1
#2012 1386.9 2060.0 1808.2 1486.3 3478.9 3291.6 1866.8
#2013 1366.0 2113.9 2301.1 1936.1 4543.2 3946.9 2403.5
#2014 1288.7 2095.4 2085.9 1738.5 4378.0 3570.9 2633.1
#2015 1187.1 1593.4 2088.3 1715.8 4894.9 0000.0 2506.8
#2016 1268.5 1993.5 2159.0 1716.0 4955.3 0000.0 2735.2
#2017 1508.5 2367.1 2749.3 2262.2 6222.8 0000.0 3212.1
#2018 1262.8 1937.9 2366.0 1755.2 5296.8 0000.0 3043.6
#2019 1487.6 2409.1 0000.0 2098.2 6191.1 0000.0 3389.1
#2020 1720.6 2541.5 0000.0 2284.1 6718.2 0000.0 4072.0
#2021 1887.4 3307.1 0000.0 0000.0 0000.0 0000.0 0000.0'
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
## Warning in if (class(val) == "try-error") return(FALSE) else return(TRUE): the
## condition has length > 1 and only the first element will be used
ef.risk.bl = portopt(ia.bl, constraints, 50, 'Black-Litterman', equally.spaced.risk = T)
## Warning in if (class(val) == "try-error") return(FALSE) else return(TRUE): the
## condition has length > 1 and only the first element will be used
# 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)
## Warning in if (class(val) == "try-error") return(FALSE) else return(TRUE): the
## condition has length > 1 and only the first element will be used
# 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)
