R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)
##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

# https://systematicinvestor.wordpress.com/2011/11/16/black-litterman-model/
# load Systematic Investor Toolbox
# setInternet2(TRUE)
#
# source(gzcon(url('https://github.com/systematicinvestor/SIT/raw/master/sit.gz', 'rb')))

rm(list=ls())
devtools::install_github('joshuaulrich/xts', force = T)
## Downloading GitHub repo joshuaulrich/xts@HEAD
## 
## * checking for file ‘/tmp/RtmpajXL42/remotes2ae58daec00/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 '/cloud/lib/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/RtmpajXL42/remotes2ae79ce6c0e/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 '/cloud/lib/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'
  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)       
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.1'
## (as 'lib' is unspecified)
## 
## Attaching package: 'corpcor'
## The following object is masked from 'package:SIT':
## 
##     cov.shrink
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.1'
## (as 'lib' is unspecified)
## 
## 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) 

plot(pressure)

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.