- portfolio.r
| # portfolio.r # # Functions for portfolio analysis # to be used in Introduction to Computational Finance & Financial Econometrics # last updated: November 7, 2000 by Eric Zivot # Oct 15, 2003 by Tim Hesterberg # November 18, 2003 by Eric Zivot # November 9, 2004 by Eric Zivot # November 9, 2008 by Eric Zivot # August 11, 2011 by Eric Zivot # # Functions: # 1. efficient.portfolio compute minimum variance portfolio # subject to target return # 2. globalMin.portfolio compute global minimum variance portfolio # 3. tangency.portfolio compute tangency portfolio # 4. efficient.frontier compute Markowitz bullet # 5. getPortfolio create portfolio object # getPortfolio <- function(er, cov.mat, weights) { # contruct portfolio object # # inputs: # er N x 1 vector of expected returns # cov.mat N x N covariance matrix of returns # weights N x 1 vector of portfolio weights # # output is portfolio object with the following elements # call original function call # er portfolio expected return # sd portfolio standard deviation # weights N x 1 vector of portfolio weights # call <- match.call() # # check for valid inputs # asset.names <- names(er) weights <- as.vector(weights) names(weights) = names(er) er <- as.vector(er) # assign names if none exist if(length(er) != length(weights)) stop("dimensions of er and weights do not match") cov.mat <- as.matrix(cov.mat) if(length(er) != nrow(cov.mat)) stop("dimensions of er and cov.mat do not match") if(any(diag(chol(cov.mat)) <= 0)) stop("Covariance matrix not positive definite") # # create portfolio # er.port <- crossprod(er,weights) sd.port <- sqrt(weights %*% cov.mat %*% weights) ans <- list("call" = call, "er" = as.vector(er.port), "sd" = as.vector(sd.port), "weights" = weights) class(ans) <- "portfolio" ans } efficient.portfolio <- function(er, cov.mat, target.return) { # compute minimum variance portfolio subject to target return # # inputs: # er N x 1 vector of expected returns # cov.mat N x N covariance matrix of returns # target.return scalar, target expected return # # output is portfolio object with the following elements # call original function call # er portfolio expected return # sd portfolio standard deviation # weights N x 1 vector of portfolio weights # call <- match.call() # # check for valid inputs # asset.names <- names(er) er <- as.vector(er) # assign names if none exist cov.mat <- as.matrix(cov.mat) if(length(er) != nrow(cov.mat)) stop("invalid inputs") if(any(diag(chol(cov.mat)) <= 0)) stop("Covariance matrix not positive definite") # remark: could use generalized inverse if cov.mat is positive semidefinite # # compute efficient portfolio # ones <- rep(1, length(er)) top <- cbind(2*cov.mat, er, ones) bot <- cbind(rbind(er, ones), matrix(0,2,2)) A <- rbind(top, bot) b.target <- as.matrix(c(rep(0, length(er)), target.return, 1)) x <- solve(A, b.target) w <- x[1:length(er)] names(w) <- asset.names # # compute portfolio expected returns and variance # er.port <- crossprod(er,w) sd.port <- sqrt(w %*% cov.mat %*% w) ans <- list("call" = call, "er" = as.vector(er.port), "sd" = as.vector(sd.port), "weights" = w) class(ans) <- "portfolio" ans } globalMin.portfolio <- function(er, cov.mat) { # Compute global minimum variance portfolio # # inputs: # er N x 1 vector of expected returns # cov.mat N x N return covariance matrix # # output is portfolio object with the following elements # call original function call # er portfolio expected return # sd portfolio standard deviation # weights N x 1 vector of portfolio weights call <- match.call() # # check for valid inputs # asset.names <- names(er) er <- as.vector(er) # assign names if none exist cov.mat <- as.matrix(cov.mat) if(length(er) != nrow(cov.mat)) stop("invalid inputs") if(any(diag(chol(cov.mat)) <= 0)) stop("Covariance matrix not positive definite") # remark: could use generalized inverse if cov.mat is positive semi-definite # # compute global minimum portfolio # cov.mat.inv <- solve(cov.mat) one.vec <- rep(1,length(er)) # w.gmin <- cov.mat.inv %*% one.vec/as.vector(one.vec %*% cov.mat.inv %*% one.vec) w.gmin <- rowSums(cov.mat.inv) / sum(cov.mat.inv) w.gmin <- as.vector(w.gmin) names(w.gmin) <- asset.names er.gmin <- crossprod(w.gmin,er) sd.gmin <- sqrt(t(w.gmin) %*% cov.mat %*% w.gmin) gmin.port <- list("call" = call, "er" = as.vector(er.gmin), "sd" = as.vector(sd.gmin), "weights" = w.gmin) class(gmin.port) <- "portfolio" gmin.port } tangency.portfolio <- function(er,cov.mat,risk.free) { # compute tangency portfolio # # inputs: # er N x 1 vector of expected returns # cov.mat N x N return covariance matrix # risk.free scalar, risk-free rate # # output is portfolio object with the following elements # call captures function call # er portfolio expected return # sd portfolio standard deviation # weights N x 1 vector of portfolio weights call <- match.call() # # check for valid inputs # asset.names <- names(er) if(risk.free < 0) stop("Risk-free rate must be positive") er <- as.vector(er) cov.mat <- as.matrix(cov.mat) if(length(er) != nrow(cov.mat)) stop("invalid inputs") if(any(diag(chol(cov.mat)) <= 0)) stop("Covariance matrix not positive definite") # remark: could use generalized inverse if cov.mat is positive semi-definite # # compute global minimum variance portfolio # gmin.port <- globalMin.portfolio(er,cov.mat) if(gmin.port$er < risk.free) stop("Risk-free rate greater than avg return on global minimum variance portfolio") # # compute tangency portfolio # cov.mat.inv <- solve(cov.mat) w.t <- cov.mat.inv %*% (er - risk.free) # tangency portfolio w.t <- as.vector(w.t/sum(w.t)) # normalize weights names(w.t) <- asset.names er.t <- crossprod(w.t,er) sd.t <- sqrt(t(w.t) %*% cov.mat %*% w.t) tan.port <- list("call" = call, "er" = as.vector(er.t), "sd" = as.vector(sd.t), "weights" = w.t) class(tan.port) <- "portfolio" tan.port } efficient.frontier <- function(er, cov.mat, nport=20, alpha.min=-0.5, alpha.max=1.5) { # Compute efficient frontier with no short-sales constraints # # inputs: # er N x 1 vector of expected returns # cov.mat N x N return covariance matrix # nport scalar, number of efficient portfolios to compute # # output is a Markowitz object with the following elements # call captures function call # er nport x 1 vector of expected returns on efficient porfolios # sd nport x 1 vector of std deviations on efficient portfolios # weights nport x N matrix of weights on efficient portfolios call <- match.call() # # check for valid inputs # asset.names <- names(er) er <- as.vector(er) cov.mat <- as.matrix(cov.mat) if(length(er) != nrow(cov.mat)) stop("invalid inputs") if(any(diag(chol(cov.mat)) <= 0)) stop("Covariance matrix not positive definite") # # create portfolio names # port.names <- rep("port",nport) ns <- seq(1,nport) port.names <- paste(port.names,ns) # # compute global minimum variance portfolio # cov.mat.inv <- solve(cov.mat) one.vec <- rep(1,length(er)) port.gmin <- globalMin.portfolio(er,cov.mat) w.gmin <- port.gmin$weights # # compute efficient frontier as convex combinations of two efficient portfolios # 1st efficient port: global min var portfolio # 2nd efficient port: min var port with ER = max of ER for all assets # er.max <- max(er) port.max <- efficient.portfolio(er,cov.mat,er.max) w.max <- port.max$weights a <- seq(from=alpha.min,to=alpha.max,length=nport) # convex combinations we.mat <- a %o% w.gmin + (1-a) %o% w.max # rows are efficient portfolios er.e <- we.mat %*% er # expected returns of efficient portfolios er.e <- as.vector(er.e) names(er.e) <- port.names cov.e <- we.mat %*% cov.mat %*% t(we.mat) # cov mat of efficient portfolios sd.e <- sqrt(diag(cov.e)) # std devs of efficient portfolios sd.e <- as.vector(sd.e) names(sd.e) <- port.names dimnames(we.mat) <- list(port.names,asset.names) # # summarize results # ans <- list("call" = call, "er" = er.e, "sd" = sd.e, "weights" = we.mat) class(ans) <- "Markowitz" ans } # # print method for portfolio object print.portfolio <- function(x, ...) { cat("Call:\n") print(x$call, ...) cat("\nPortfolio expected return: ", format(x$er, ...), "\n") cat("Portfolio standard deviation: ", format(x$sd, ...), "\n") cat("Portfolio weights:\n") print(round(x$weights,4), ...) invisible(x) } # # summary method for portfolio object summary.portfolio <- function(object, risk.free=NULL, ...) # risk.free risk-free rate. If not null then # compute and print Sharpe ratio # { cat("Call:\n") print(object$call) cat("\nPortfolio expected return: ", format(object$er, ...), "\n") cat( "Portfolio standard deviation: ", format(object$sd, ...), "\n") if(!is.null(risk.free)) { SharpeRatio <- (object$er - risk.free)/object$sd cat("Portfolio Sharpe Ratio: ", format(SharpeRatio), "\n") } cat("Portfolio weights:\n") print(round(object$weights,4), ...) invisible(object) } # hard-coded 4 digits; prefer to let user control, via ... or other argument # # plot method for portfolio object plot.portfolio <- function(object, ...) { asset.names <- names(object$weights) barplot(object$weights, names=asset.names, xlab="Assets", ylab="Weight", main="Portfolio Weights", ...) invisible() } # # print method for Markowitz object print.Markowitz <- function(x, ...) { cat("Call:\n") print(x$call) xx <- rbind(x$er,x$sd) dimnames(xx)[[1]] <- c("ER","SD") cat("\nFrontier portfolios' expected returns and standard deviations\n") print(round(xx,4), ...) invisible(x) } # hard-coded 4, should let user control # # summary method for Markowitz object summary.Markowitz <- function(object, risk.free=NULL) { call <- object$call asset.names <- colnames(object$weights) port.names <- rownames(object$weights) if(!is.null(risk.free)) { # compute efficient portfolios with a risk-free asset nport <- length(object$er) sd.max <- object$sd[1] sd.e <- seq(from=0,to=sd.max,length=nport) names(sd.e) <- port.names # # get original er and cov.mat data from call er <- eval(object$call$er) cov.mat <- eval(object$call$cov.mat) # # compute tangency portfolio tan.port <- tangency.portfolio(er,cov.mat,risk.free) x.t <- sd.e/tan.port$sd # weights in tangency port rf <- 1 - x.t # weights in t-bills er.e <- risk.free + x.t*(tan.port$er - risk.free) names(er.e) <- port.names we.mat <- x.t %o% tan.port$weights # rows are efficient portfolios dimnames(we.mat) <- list(port.names, asset.names) we.mat <- cbind(rf,we.mat) } else { er.e <- object$er sd.e <- object$sd we.mat <- object$weights } ans <- list("call" = call, "er"=er.e, "sd"=sd.e, "weights"=we.mat) class(ans) <- "summary.Markowitz" ans } print.summary.Markowitz <- function(x, ...) { xx <- rbind(x$er,x$sd) port.names <- names(x$er) asset.names <- colnames(x$weights) dimnames(xx)[[1]] <- c("ER","SD") cat("Frontier portfolios' expected returns and standard deviations\n") print(round(xx,4), ...) cat("\nPortfolio weights:\n") print(round(x$weights,4), ...) invisible(x) } # hard-coded 4, should let user control # # plot efficient frontier # # things to add: plot original assets with names # tangency portfolio # global min portfolio # risk free asset and line connecting rf to tangency portfolio # plot.Markowitz <- function(object, plot.assets=FALSE, ...) # plot.assets logical. If true then plot asset sd and er { if (!plot.assets) { y.lim=c(0,max(object$er)) x.lim=c(0,max(object$sd)) plot(object$sd,object$er,type="b",xlim=x.lim, ylim=y.lim, xlab="Portfolio SD", ylab="Portfolio ER", main="Efficient Frontier", ...) } else { call = object$call mu.vals = eval(call$er) sd.vals = sqrt( diag( eval(call$cov.mat) ) ) y.lim = range(c(0,mu.vals,object$er)) x.lim = range(c(0,sd.vals,object$sd)) plot(object$sd,object$er,type="b", xlim=x.lim, ylim=y.lim, xlab="Portfolio SD", ylab="Portfolio ER", main="Efficient Frontier", ...) text(sd.vals, mu.vals, labels=names(mu.vals)) } invisible() } |
Source the r source codes above.
| source("H:/Quantitative Research/R/Optimization/portfolio.r") |
Sample ticker code:
SPY SPDR S&P 500 ETF TRUST
EFA ISHARES MSCI EAFE ETF
LQD ISHARES IBOXX INVESTMENT GRA
HYG ISHARES IBOXX HIGH YIELD COR
Weekly basis stats (in decimal numbers):
(1) Risk-free return
| Risk-free | 0.00018631 |
(2) Average Return and Standard Deviation of Returns
| Weekly | SPY | EFA | LQD | HYG |
| Average | 0.00260290 | 0.00120427 | 0.00110844 | 0.00132166 |
| SD | 0.02029120 | 0.02491309 | 0.00788386 | 0.01151264 |
(3) Covariance matrix
| Cov | SPY | EFA | LQD | HYG |
| SPY | 0.00041173 | 0.00045214 | -0.00000957 | 0.00017600 |
| EFA | 0.00045214 | 0.00062066 | -0.00000215 | 0.00021265 |
| LQD | -0.00000957 | -0.00000215 | 0.00006216 | 0.00002826 |
| HYG | -0.00000617 | -0.00000570 | 0.00000461 | -0.00000108 |
Annual basis stats (in decimal numbers):
(1') Geometrical-Average Risk-free return
Annual
| Risk-free | 0.00973449 |
(2') Geometrical-Average Return and Standard Deviation of Returns (SD figures of (2) are multiplied by 52^0.5.)
|
| Cov | SPY | EFA | LQD | HYG |
| SPY | 0.02141010 | 0.02351144 | -0.00049749 | 0.00915216 |
| EFA | 0.02351144 | 0.03227444 | -0.00011169 | 0.01105802 |
| LQD | -0.00049749 | -0.00011169 | 0.00323208 | 0.00146939 |
| HYG | 0.00915216 | 0.01105802 | 0.00146939 | 0.00689213 |
Stats on a annual basis:
| # annual r.free <- 0.00973449 r.free [1] 0.00973449 er <- c( "SPY" = 0.14473728, "EFA" = 0.06458429, "LQD" = 0.05929884, "EEM" = 0.07109433 ) er SPY EFA LQD EEM 0.14473728 0.06458429 0.05929884 0.07109433 covmat <- matrix(c(0.02141010, 0.02351144, -0.00049749, 0.00915216, 0.02351144, 0.03227444, -0.00011169, 0.01105802, -0.00049749, -0.00011169, 0.00323208, 0.00146939, 0.00915216, 0.01105802, 0.00146939, 0.00689213), nrow = 4, ncol = 4, byrow = TRUE) rownames(covmat) <- c("SPY", "EFA", "LQD", "HYG") colnames(covmat) <- c("SPY", "EFA", "LQD", "HYG") covmat SPY EFA LQD HYG SPY 0.02141010 0.02351144 -0.00049749 0.00915216 EFA 0.02351144 0.03227444 -0.00011169 0.01105802 LQD -0.00049749 -0.00011169 0.00323208 0.00146939 HYG 0.00915216 0.01105802 0.00146939 0.00689213 |
Equal-weighted Portfolio
| # [1 1 1 1] rep(1,4) # equally-weighted ew = rep(1,4)/4 ew [1] 0.25 0.25 0.25 0.25
equalWeight.portfolio = getPortfolio(er=er,cov.mat=covmat,weights=ew)
equalWeight.portfolio Call: getPortfolio(er = er, cov.mat = covmat, weights = ew) Portfolio expected return: 0.08492869 Portfolio standard deviation: 0.09777922 Portfolio weights: SPY EFA LQD EEM 0.25 0.25 0.25 0.25 plot(equalWeight.portfolio) # If you want to change the color, then try: # plot(equalWeight.portfolio, col="blue") |
 |
Tangency Portfolio (long-short)
tan.port <- tangency.portfolio(er, covmat, r.free)
# You can also run without a "shorts" option.
# tan.port <- tangency.portfolio(er, covmat, r.free) tan.port Call: tangency.portfolio(er = er, cov.mat = covmat, risk.free = r.free, shorts = TRUE) Portfolio expected return: 0.1371406 Portfolio standard deviation: 0.06950249 Portfolio weights: SPY EFA LQD EEM 0.9766 -0.5690 0.8118 -0.2194 print(tan.port) Call: tangency.portfolio(er = er, cov.mat = covmat, risk.free = r.free, shorts = TRUE) Portfolio expected return: 0.1371406 Portfolio standard deviation: 0.06950249 Portfolio weights: SPY EFA LQD EEM 0.9766 -0.5690 0.8118 -0.2194 plot(tan.port) |
 |
Minimum Volatility Portfolio (long-short)
gmin.port <- globalMin.portfolio(er, covmat)
gmin.port Call: globalMin.portfolio(er = er, cov.mat = covmat, shorts = TRUE) Portfolio expected return: 0.07863945 Portfolio standard deviation: 0.05111289 Portfolio weights: SPY EFA LQD EEM 0.2210 -0.1071 0.7990 0.0872 attributes(gmin.port) $names [1] "call" "er" "sd" "weights" $class [1] "portfolio" print(gmin.port) Call: globalMin.portfolio(er = er, cov.mat = covmat, shorts = TRUE) Portfolio expected return: 0.07863945 Portfolio standard deviation: 0.05111289 Portfolio weights: SPY EFA LQD EEM 0.2210 -0.1071 0.7990 0.0872 plot(gmin.port) |
 |
Target Return Portfolio (long-short)
target.return = er["SPY"]
target.return SPY 0.1447373 e.port.SPY = efficient.portfolio(er, covmat, target.return) e.port.SPY Call: efficient.portfolio(er = er, cov.mat = covmat, target.return = target.return, shorts = TRUE) Portfolio expected return: 0.1447373 Portfolio standard deviation: 0.0737838 Portfolio weights: SPY EFA LQD EEM 1.0747 -0.6290 0.8135 -0.2592 plot(e.port.SPY) |
 |
Efficient Frontier (long-short)
ef = efficient.frontier(er, covmat, alpha.min = -2, alpha.max = 1.5, nport = 20)
attributes(ef) $names [1] "call" "er" "sd" "weights" $class [1] "Markowitz" summary(ef) Frontier portfolios' expected returns and standard deviations port 1 port 2 port 3 port 4 port 5 port 6 port 7 port 8 port 9 port 10 ER 0.2769 0.2648 0.2526 0.2404 0.2282 0.2161 0.2039 0.1917 0.1795 0.1673 SD 0.1676 0.1583 0.1491 0.1399 0.1308 0.1219 0.1130 0.1044 0.0960 0.0878 port 11 port 12 port 13 port 14 port 15 port 16 port 17 port 18 port 19 ER 0.1552 0.1430 0.1308 0.1186 0.1065 0.0943 0.0821 0.0699 0.0578 SD 0.0801 0.0728 0.0662 0.0604 0.0558 0.0526 0.0512 0.0516 0.0538 port 20 ER 0.0456 SD 0.0576 Portfolio weights: SPY EFA LQD EEM port 1 2.7822 -1.6727 0.8425 -0.9520 port 2 2.6249 -1.5766 0.8399 -0.8882 port 3 2.4676 -1.4804 0.8372 -0.8244 port 4 2.3104 -1.3843 0.8345 -0.7606 port 5 2.1531 -1.2882 0.8318 -0.6968 port 6 1.9958 -1.1920 0.8292 -0.6330 port 7 1.8386 -1.0959 0.8265 -0.5692 port 8 1.6813 -0.9998 0.8238 -0.5053 port 9 1.5240 -0.9036 0.8211 -0.4415 port 10 1.3668 -0.8075 0.8185 -0.3777 port 11 1.2095 -0.7114 0.8158 -0.3139 port 12 1.0522 -0.6152 0.8131 -0.2501 port 13 0.8950 -0.5191 0.8104 -0.1863 port 14 0.7377 -0.4230 0.8078 -0.1225 port 15 0.5804 -0.3268 0.8051 -0.0587 port 16 0.4232 -0.2307 0.8024 0.0051 port 17 0.2659 -0.1346 0.7998 0.0689 port 18 0.1086 -0.0385 0.7971 0.1327 port 19 -0.0486 0.0577 0.7944 0.1966 port 20 -0.2059 0.1538 0.7917 0.2604 plot(ef, plot.assets=T, col="blue", pch=16) points(gmin.port$sd, gmin.port$er, col="green", pch=16, cex=2) points(tan.port$sd, tan.port$er, col="red", pch=16, cex=2) text(gmin.port$sd, gmin.port$er, labels="GLOBAL MIN", pos=2) text(tan.port$sd, tan.port$er, labels="TANGENCY", pos=2) sr.tan = (tan.port$er - r.free)/tan.port$sd abline(a=r.free, b=sr.tan, col="green", lwd=2) |
 |
No comments:
Post a Comment