################################################################# # Create some data ################################################################# require(quantmod) # not used in abbreviated example, but useful for reporting startYear = 1928 endYear = 2010 YEARS =startYear:endYear # nominal returns SP500 = c(0.4381,-0.083,-0.2512,-0.4384,-0.0864,0.4998,-0.0119,0.4674,0.3194,-0.3534,0.2928, -0.011,-0.1067,-0.1277,0.1917,0.2506,0.1903,0.3582,-0.0843,0.052,0.057,0.183,0.3081,0.2368, 0.1815,-0.0121,0.5256,0.326,0.0744,-0.1046,0.4372,0.1206,0.0034,0.2664,-0.0881,0.2261,0.1642, 0.124,-0.0997,0.238,0.1081,-0.0824,0.0356,0.1422,0.1876,-0.1431,-0.259,0.370,0.2383,-0.0698, 0.0651,0.1852,0.3174,-0.047,0.2042,0.2234,0.0615,0.3124,0.1849,0.0581,0.1654,0.3148,-0.0306, 0.3023,0.0749,0.0997,0.0133,0.372,0.2382,0.3186,0.2834,0.2089,-0.0903,-0.1185,-0.2197,0.2836, 0.1074,0.0483,0.1561,0.0548,-0.3658,0.2592,0.148600) BILLS = c(0.0308,0.0316,0.0455,0.0231,0.0107,0.0096,0.0032,0.0018,0.0017,0.003,0.0008,0.0004, 0.0003,0.0008,0.0034,0.0038,0.0038,0.0038,0.0038,0.0057,0.0102,0.011,0.0117,0.0148,0.0167, 0.0189,0.0096,0.0166,0.0256,0.0323,0.0178,0.0326,0.0305,0.0227,0.0278,0.0311,0.0351,0.039, 0.0484,0.0433,0.0526,0.0656,0.0669,0.0454,0.0395,0.0673,0.0778,0.0599,0.0497,0.0513,0.0693, 0.0994,0.1122,0.143,0.1101,0.0845,0.0961,0.0749,0.0604,0.0572,0.0645,0.0811,0.0755,0.0561, 0.0341,0.0298,0.0399,0.0552,0.0502,0.0505,0.0473,0.0451,0.0576,0.0367,0.0166,0.0103,0.0123, 0.0301,0.0468,0.0464,0.0159,0.0014,0.001300) BONDS=c(0.0084,0.042,0.0454,-0.0256,0.0879,0.0186,0.0796,0.0447,0.0502,0.0138,0.0421,0.0441, 0.054,-0.0202,0.0229,0.0249,0.0258,0.038,0.0313,0.0092,0.0195,0.0466,0.0043,-0.003,0.0227, 0.0414,0.0329,-0.0134,-0.0226,0.068,-0.021,-0.0265,0.1164,0.0206,0.0569,0.0168,0.0373,0.0072, 0.0291,-0.0158,0.0327,-0.0501,0.1675,0.0979,0.0282,0.0366,0.0199,0.0361,0.1598,0.0129,-0.0078, 0.0067,-0.0299,0.082,0.3281,0.032,0.1373,0.2571,0.2428,-0.0496,0.0822,0.1769,0.0624,0.150, 0.0936,0.1421,-0.0804,0.2348,0.0143,0.0994,0.1492,-0.0825,0.1666,0.0557,0.1512,0.0038,0.0449, 0.0287,0.0196,0.1021,0.201,-0.1112,0.084600) CPI=c(-0.0115607,0.005848,-0.0639535,-0.0931677,-0.1027397,0.0076336,0.0151515,0.0298507, 0.0144928,0.0285714,-0.0277778,0.0000,0.0071429,0.0992908,0.0903226,0.0295858,0.0229885, 0.0224719,0.1813187,0.0883721,0.0299145,-0.0207469,0.059322,0.0600,0.0075472,0.0074906, -0.0074349,0.0037453,0.0298507,0.0289855,0.0176056,0.017301,0.0136054,0.0067114,0.0133333, 0.0164474,0.0097087,0.0192308,0.0345912,0.0303951,0.0471976,0.0619718,0.0557029,0.0326633, 0.0340633,0.0870588,0.1233766,0.0693642,0.0486486,0.0670103,0.0901771,0.1329394,0.125163, 0.0892236,0.0382979,0.0379098,0.0394867,0.0379867,0.010979,0.0443439,0.0441941,0.046473, 0.0610626,0.0306428,0.0290065,0.0274841,0.026749,0.0253841,0.0332248,0.017024,0.016119, 0.0268456,0.0338681,0.0155172,0.0237691,0.0187949,0.0325556,0.0341566,0.0254065,0.0408127, 0.0009141,0.0272133,0.0149572) GOLD = c(0.000000000,0.000000000,0.000000000,0.000000000,0.000000000,0.447002860,0.079661828, 0.000000000,0.000000000,0.000000000,0.000000000,0.000000000,-0.014388737,0.028573372,0.000000000, 0.027779564,-0.006872879,0.027212564,0.026491615,0.117056555,-0.023530497,-0.036367644, -0.006191970,-0.006230550,-0.033039854,-0.086306904,-0.007067167,-0.002840911,0.001421464, 0.001419447,0.000000000,0.000000000,0.034846731,-0.027779564,-0.004234304,-0.002832863, 0.002832863,0.004234304,-0.002820876,0.002820876,0.203228242,-0.059188871,-0.052577816, 0.136739608,0.358646094,0.511577221,0.545727802,-0.132280611,-0.253090628,0.168898536, 0.265477915,0.464157559,0.689884535,-0.285505793,-0.201637346,0.120144312,-0.163629424, -0.127202258,0.149181164,0.192236014,-0.020385757,-0.134512586,0.005221944,-0.058998341, -0.051002554,0.045462374,0.064538521,0.002600782,0.010336009,-0.155436854,-0.121167134, -0.052185753,0.000000000,-0.028987537,0.133990846,0.157360955,0.119003292,0.084161792, 0.308209839,0.142551544,0.220855221,0.110501915,0.228258652) # put into a data frame (matrix) fnominal=data.frame(stocks=SP500, bills=BILLS, bonds=BONDS, gold=GOLD, CPI=CPI) freal=data.frame(stocks=SP500-CPI, bills=BILLS-CPI, bonds=BONDS-CPI, gold=GOLD-CPI) # compute real returns realreturns = apply(freal, 2, mean) realreturnspct = realreturns*100 # print them realreturnspct # compute real volatility (standard deviation of real returns) realsds = apply(freal, 2, sd) realsdspct = realsds*100 # print them realsdspct ################################################################# # Run efficient frontier and transition map with Systematic Investor code ################################################################# # put input assumption in suitable format ia = list() ia$n = length(freal) ia$annual.factor = 1 ia$symbols = names(freal) ia$symbol.names = names(freal) ia$hist.returns = freal ia$arithmetic.return = apply(ia$hist.returns, 2, mean, na.rm = T) ia$arithmetic.return = (1 + ia$arithmetic.return)^ia$annual.factor - 1 ia$geometric.return = apply(ia$hist.returns, 2, function(x) prod(1+x)^(1/length(x))-1 ) ia$geometric.return = (1 + ia$geometric.return)^ia$annual.factor - 1 ia$risk = apply(ia$hist.returns, 2, sd, na.rm = T) ia$risk = sqrt(ia$annual.factor) * ia$risk ia$correlation = cor(ia$hist.returns, use = 'complete.obs', method = 'pearson') ia$cov = cov(ia$hist.returns, use = 'complete.obs', method = 'pearson') ia$expected.return = ia$arithmetic.return ############################################################################### # Load Systematic Investor Toolbox (SIT) # http://systematicinvestor.wordpress.com/systematic-investor-toolbox/ ############################################################################### setInternet2(TRUE) con = gzcon(url('http://www.systematicportfolio.com/sit.gz', 'rb')) source(con) close(con) # if the above doesn't work, download and unpack the URL above in a local dir and try # source('C:\\Temp\\sit.r') # do optimizations 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) ef.risk = portopt(ia, constraints, 50, 'Historical', equally.spaced.risk = T) # chart layout( matrix(1:2, nrow = 2) ) # can skip this in RStudio, use back button plot.ef(ia, list(ef.risk), portfolio.risk, T, T) # you should see 2 plots ################################################################# # use linear programming to find maximum return portfolio (all stocks) ################################################################# install.packages('lpSolve') require(lpSolve) # find maximum return portfolio (rightmost point of efficient frontier # will be 100% of highest return asset # maximize # w1 * stocks +w2 *bills +w3*bonds + w4 * gold # subject to 0 <= w <= 1 for each W # will pick highest return asset with w=1 # skipping >0 constraint, no negative return assets, so not binding opt.objective <- realreturns opt.constraints <- matrix (c(1, 1, 1, 1, # constrain sum of weights to 1 1, 0, 0, 0, # constrain w1 < 1 0, 1, 0, 0, # constrain w2 < 1 0, 0, 1, 0, # constrain w3 < 1 0, 0, 0, 1) # constrain w4 < 1 , nrow=5, byrow=TRUE) opt.operator <- c("=", "<=", "<=", "<=", "<=") opt.rhs <- c(1, 1, 1, 1, 1) opt.dir="max" solution.maxret = lp (direction = opt.dir, opt.objective, opt.constraints, opt.operator, opt.rhs) # portfolio weights for max return portfolio wts.maxret=solution.maxret$solution # return for max return portfolio ret.maxret=solution.maxret$objval # compute return covariance matrix to determine volatility of this portfolio covmatrix = cov(freal, use = 'complete.obs', method = 'pearson') # multiply weights x covariances x weights, gives variance var.maxret = wts.maxret %*% covmatrix %*% wts.maxret # square root gives standard deviation (volatility) vol.maxret = sqrt(var.maxret) ################################################################# # find minimum volatility portfolio ################################################################# install.packages('quadprog') require(quadprog) # minimize variance: w %*% covmatrix %*% t(w) # subject to sum of ws = 1 # subject to each w > 0 # solution.minvol <- solve.QP(covmatrix, zeros, t(opt.constraints), opt.rhs, meq = opt.meq) # first 2 parameters covmatrix, zeros define function to be minimized # if zeros is all 0s, the function minimized ends up equal to port variance / 2 # opt.constraints is the left hand side of the constraints, ie the cs in # c1 w1 + c2 w2 ... + cn wn = K # opt.rhs is the Ks in the above equation # meq means the first meq rows are 'equals' constraints, remainder are >= constraints # if you want to do a <= constraint, multiply by -1 to make it a >= constraint # does not appear to accept 0 RHS, so we make it a tiny number> 0 # compute covariance matrix covmatrix = cov(freal, use = 'complete.obs', method = 'pearson') nObs <- length(freal$stocks) nAssets <- length(freal) # 1 x numassets array of 1s opt.constraints <- matrix (c(1, 1, 1, 1, # sum of weights =1 1, 0, 0, 0, # w1 >= 0 0, 1, 0, 0, # w1 >= 0 0, 0, 1, 0, # w1 >= 0 0, 0, 0, 1) # w2 >= 0 , nrow=5, byrow=TRUE) opt.rhs <- matrix(c(1, 0.000001, 0.000001, 0.000001, 0.000001)) opt.meq <- 1 # first constraint is '=', rest are '>=' # numassets x 1 array of 0s zeros <- array(0, dim = c(nAssets,1)) solution.minvol <- solve.QP(covmatrix, zeros, t(opt.constraints), opt.rhs, meq = opt.meq) wts.minvol = solution.minvol$solution var.minvol = solution.minvol$value *2 ret.minvol = realreturns %*% wts.minvol vol.minvol = sqrt(var.minvol) # generate a sequence of 50 evenly spaced returns between min var return and max return lowreturn=ret.minvol highreturn=ret.maxret minreturns=seq(lowreturn, highreturn, length.out=50) # add a return constraint: sum of weight * return >= x retconst= rbind(opt.constraints, realreturns) retrhs=rbind(opt.rhs, ret.minvol) # create vectors for the returns, vols, and weights along the frontier, # starting with the minvol portfolio out.ret=c(ret.minvol) out.vol=c(vol.minvol) out.stocks=c(wts.minvol[1]) out.bills=c(wts.minvol[2]) out.bonds=c(wts.minvol[3]) out.gold=c(wts.minvol[4]) ################################################################# # loop and run a minimum volatility optimization for each return level from 2-49 ################################################################# for(i in 2:(length(minreturns) - 1)) { print(i) # start with existing constraints, no return constraint tmp.constraints = retconst tmp.rhs=retrhs # set return constraint tmp.rhs[6] = minreturns[i] tmpsol <- solve.QP(covmatrix, zeros, t(tmp.constraints), tmp.rhs, meq = opt.meq) tmp.wts = tmpsol$solution tmp.var = tmpsol$value *2 out.ret[i] = realreturns %*% tmp.wts out.vol[i] = sqrt(tmp.var) out.stocks[i]=tmp.wts[1] out.bills[i]=tmp.wts[2] out.bonds[i]=tmp.wts[3] out.gold[i]=tmp.wts[4] } # put maxreturn portfolio in return series for max return, index =50 out.ret[50]=c(ret.maxret) out.vol[50]=c(vol.maxret) out.stocks[50]=c(wts.maxret[1]) out.bills[50]=c(wts.maxret[2]) out.bonds[50]=c(wts.maxret[3]) out.gold[50]=c(wts.maxret[4]) # organize in a data frame efrontier=data.frame(out.ret*100) efrontier$vol=out.vol*100 efrontier$stocks=out.stocks*100 efrontier$bills=out.bills*100 efrontier$bonds=out.bonds*100 efrontier$gold=out.gold*100 names(efrontier) = c("Return", "Risk", "%Stocks", "%Bills", "%Bonds", "%Gold") ############################################################ # charts ############################################################ install.packages('ggplot2') require(ggplot2) # define colors dvblue = "#000099" dvred = "#e41a1c" dvgreen = "#4daf4a" dvpurple = "#984ea3" dvorange = "#ff7f00" dvyellow = "#ffff33" dvgray="#666666" apoints=data.frame(realsdspct) apoints$realreturns = realreturnspct ggplot(data=efrontier, aes(x=Risk, y=Return)) + opts(title="Efficient Frontier") + theme_bw() + geom_line(size=1.4) + geom_point(aes(x=apoints$realsdspct, y=apoints$realreturns)) + scale_x_continuous(limits=c(1,24)) + annotate("text", apoints[1,1], apoints[1,2],label=" stocks", hjust=0) + annotate("text", apoints[2,1], apoints[2,2],label=" bills", hjust=0) + annotate("text", apoints[3,1], apoints[3,2],label=" bonds", hjust=0) + annotate("text", apoints[4,1], apoints[4,2],label=" gold", hjust=0) + annotate("text", 19,0.3,label="streeteye.com", hjust=0, alpha=0.5) keep=c("Risk", "%Stocks","%Bills","%Bonds","%Gold") efrontier.tmp = efrontier[keep] efrontier.m = melt(efrontier.tmp, id ='Risk') ggplot(data=efrontier.m, aes(x=Risk, y=value, colour=variable, fill=variable)) + theme_bw() + opts(title="Transition Map", legend.position="top", legend.direction="horizontal") + ylab('% Portfolio') + geom_area() + scale_colour_manual("", breaks=c("%Stocks", "%Bills", "%Bonds","%Gold"), values = c(dvblue,dvgreen,dvred,dvyellow), labels=c('%Stocks', '%Bills','%Bonds','%Gold')) + scale_fill_manual("", breaks=c("%Stocks", "%Bills", "%Bonds","%Gold"), values = c(dvblue,dvgreen,dvred,dvyellow), labels=c('%Stocks', '%Bills','%Bonds','%Gold')) + annotate("text", 16,-2.5,label="streeteye.com", hjust=0, alpha=0.5)