data609project
Portfolio Optimization
jim lung
November 26, 2018
Data 609 Project proposal - Mathematical Modeling Techniques for Data Analytics
Authors: Jim Lung
Aim:
** Use mathematical models to make a decision for portfolio optimization Investigating portfolio optimization with expected return on investiment in risk control **
Data source:
** Historical stock price data are readily accessible using functions in “Quantmod” package The filtered data for this application selects total 9 stock cases from 2015 to August 2018 The datasets include the date, daily market close price, market volumes, the closing price will be used to make for portfolio optimization **
In the following sections, we use a variety of mathematical tools to perform the following tasks:
* Data loading * Graphical Exploration * Compute daily, monthly and yearly return * Calculate the Mean Variance model by specific stocks * Use minimax model to optimize portfolio * Use linear programming techniques to compare the log return * Use a quadratic programming approach to determine appropriate portfolio weights
Load Libraries
suppressMessages(suppressWarnings(library(knitr)))
suppressMessages(suppressWarnings(library(ggplot2)))
suppressMessages(suppressWarnings(library(ggthemes)))
suppressMessages(suppressWarnings(library(dplyr)))
suppressMessages(suppressWarnings(library(tidyr)))
suppressMessages(suppressWarnings(library(stringr)))
suppressMessages(suppressWarnings(library(lpSolveAPI)))
suppressMessages(suppressWarnings(library(xts)))
suppressMessages(suppressWarnings(library(matrixcalc)))
suppressMessages(suppressWarnings(library(dygraphs)))
suppressMessages(suppressWarnings(library(ggthemes)))
suppressMessages(suppressWarnings(library(highcharter)))
suppressMessages(suppressWarnings(library(viridis)))
suppressMessages(suppressWarnings(library(tibbletime)))
suppressMessages(suppressWarnings(library(timetk)))
suppressMessages(suppressWarnings(library(tidyquant)))
suppressMessages(suppressWarnings(library(tidyverse)))1. Data
# input and print available stocks
symbols <- c("AAPL", "AMD", "ADI", "ABBV", "AET", "A", "APD", "AA","CF")
stock <- c("Apple Inc.", "Advanced Micro Devices, Inc.", "Analog Devices, Inc.", "AbbVie Inc","AETNA INC", "Agilent Technologies Inc","Air Products & Chemicals, Inc.","Alcoa Corp","CF Industries Holdings, Inc.")
df <- data.frame(ticker=symbols, stocks=stock)
knitr::kable(df[,1:2], col.names = c("Ticker","Stocks")) | Ticker | Stocks |
|---|---|
| AAPL | Apple Inc. |
| AMD | Advanced Micro Devices, Inc. |
| ADI | Analog Devices, Inc. |
| ABBV | AbbVie Inc |
| AET | AETNA INC |
| A | Agilent Technologies Inc |
| APD | Air Products & Chemicals, Inc. |
| AA | Alcoa Corp |
| CF | CF Industries Holdings, Inc. |
Construct a vector of tickers and gather prices for them using the getSymbols function within quantmod. We will next calculate returns and convert the data to a time series object.
# The symbols vector holds our tickers.
tickers <- c("AAPL", "AMD", "ADI", "ABBV", "AET", "A", "APD", "AA","CF")
# The prices object will hold our raw price data
close_price<-
getSymbols(tickers,src="yahoo", from = "2017-01-01",
auto.assign = TRUE, warnings = FALSE) %>%
map(~Ad(get(.))) %>% #Extract (transformed) data from a suitable OHLC object.
reduce(merge) %>% #reduce() combines from the left, reduce_right() combines from the right
`colnames<-`(tickers )## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
##
## This message is shown once per session and may be disabled by setting
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.##
## WARNING: There have been significant changes to Yahoo Finance data.
## Please see the Warning section of '?getSymbols.yahoo' for details.
##
## This message is shown once per session and may be disabled by setting
## options("getSymbols.yahoo.warning"=FALSE).## pausing 1 second between requests for more than 5 symbols
## pausing 1 second between requests for more than 5 symbols
## pausing 1 second between requests for more than 5 symbols
## pausing 1 second between requests for more than 5 symbols
## pausing 1 second between requests for more than 5 symbolsclose_price<-tk_zoo(close_price)%>%na.omit()
# Calculate returns
prices <- close_price
returns.data <- sapply(close_price, CalculateReturns)
returns.data <- na.omit(returns.data)
returns.data %>%tk_tbl()%>%tail()## # A tibble: 6 x 10
## index AAPL AMD ADI ABBV AET A APD
## <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2018-11-20 -0.0478 0.00523 0.0409 -0.0155 -0.0168 0.0784 -0.0146
## 2 2018-11-21 -0.00113 -0.0250 0.00741 -0.0220 0.00147 0.00296 0.0103
## 3 2018-11-23 -0.0254 0.0347 -0.0116 -0.00604 0.00519 0.00724 -0.0106
## 4 2018-11-26 0.0135 0.0361 0.0161 0.00712 0.0231 0.0152 0.00692
## 5 2018-11-27 -0.00218 0.0483 -0.00544 0.0167 0.00904 0.00159 -0.0146
## 6 2018-11-28 0.0385 0.0138 0.0182 0.0144 0.00330 0.0306 0.0153
## # ... with 2 more variables: AA <dbl>, CF <dbl>x = close_price
date = index(close_price)
close_price = data.frame(date,x)Graphical Exploration
Use adjusted closed price to plot graph from 2017 until now:
start="2017-01-01"
end=today()
dateWindow <-cat(paste(sQuote("2017-01-01"),sQuote(end),sep=","), "\n")## '2017-01-01','2018-12-04'dygraph(close_price, main = "Value", group = "stock") %>%
dyRebase(value = 100) %>%
dyRangeSelector(dateWindow = dateWindow)According to the historial price movements, it indicate the future trends and it concern the maximum of expected return, it is difficult to find out any good investment stratege in plot graph. So we should calculate the return by daily, weekly or yearly to decide the short , long term investment.
2. Compute daily, monthly and yearly return
#daily return
daily_return <- as.data.frame(apply(close_price[,-1] , 2 , diff )/close_price[-479,-1])
daily_return <- cbind(Date = close_price[-1,1],daily_return)
#monthly returns
temp <- lapply(daily_return[,-1], function(x) x+1)
temp$year <- format(daily_return[,1], "%Y")
temp$month <- format(daily_return[,1], "%m")
temp<- as.data.frame(temp)
monthly_return <- aggregate(temp[-c(10,11)] ,by=list(temp$month,temp$year), FUN = prod, na.rm=TRUE)
monthly_return <-cbind(monthly_return[,c(1,2)],monthly_return[,-c(1,2)]-1)
colnames(monthly_return)[1] <-'Month'
colnames(monthly_return)[2] <-'Year'
monthly_return$Date <- as.Date(paste0(monthly_return$Year,'-',monthly_return$Month,"-01"))
monthly_return <- monthly_return[,c(12,1:11)]
#yearly returns
annual_return <- aggregate(temp[-c(10,11)] ,list(temp$year), FUN = prod, na.rm=TRUE)
annual_return <-cbind(annual_return[,1],annual_return[,-1]-1)
colnames(annual_return)[1] <-'Year'
#$Year <- as.Date(ISOdate(plot_dt_Y$Year, 1, 1))
# time period: year
n <- tail(as.numeric(annual_return$Year), n=1) - as.numeric(as.numeric(annual_return$Year))[1]
# 5 years average annual return
ave_annual_return <- apply(annual_return[,-1]+1,2,prod)
ave_annual_return <- ave_annual_return^(1/n)-1Compare the annual return:
annualreturn <- as.data.frame(ave_annual_return)
colnames(annualreturn)<-'ave_yearly_return'
kable(annualreturn)| ave_yearly_return | |
|---|---|
| AAPL | 0.6044175 |
| AMD | 0.8627081 |
| ADI | 0.3042133 |
| ABBV | 0.5339709 |
| AET | 0.7718944 |
| A | 0.5605331 |
| APD | 0.1653859 |
| AA | 0.1571870 |
| CF | 0.4239632 |
| To comp | are the average yearly return, AMD and AET are the most increasing return yearly. |
3.Linear Programming - Mean Variance model
Investors are risk averse in that they prefer higher return for a given level of risk (variance, standard deviation), or they want to minimize risk for a given level of returns, so we go to minimize the variance and maximize the return.
Suppose data are observed for N securities, over T time periods. Let
\(y_{jt}\) = Return on one dollar invested in security j in time period t.
\(\bar{y_j}\)=Average Return on security j
\(w_j=\) Portfolio allocation to security j.
\(y_{pt}\) = Return on portfolio in time period t
\(E_p\) =Average Return on portfolio
\(M_p\) = Minimum return on portfolio
The objective function: \[min \displaystyle\sum_{j=k}^{N}\displaystyle\sum_{j=1}^{N}w_jw_ks_{jk}\]
subject to: \[\displaystyle\sum_{j=1}^{N}w_j\bar{y_j} \geq G \]
with: \[s_{jk} = \frac{1}{T-N}\displaystyle\sum_{t=1}^{T}(y_{jt} - \bar{y_j})(y_{kt} -\bar{y_k})\]
#The subset of the data as candidate
sub_return <- monthly_return[,c('Date','Month','Year',tickers)]
# convert to matrix and vector
cov_matrix <- cov(sub_return[sub_return$Date > as.Date('2018-01-01'),][,-c(1:3)])
pstart <- rep(1/12,8)
# arithmetic mean of monthly return
mu <- apply(sub_return[sub_return$Date > as.Date('2018-01-01'),][,-c(1:3)],2,mean)
# standard deviation
std <- apply(sub_return[sub_return$Date > as.Date('2018-01-01'),][,-c(1:3)],2,sd)
#sharp ratio
sharpe <- mu /std
############################################
#Minimize the variance
############################################
SSRE <- function(parx) {
par<- c(parx,1-sum(parx))
if(all(par > 0 & par < 1)) { # parameters meet requirements
sqrt(t(par)%*% cov_matrix%*% par) # this is a linear algebra version of your objective without the division by xi
} else 1e7 # penalty for parameters not satisfying constraints
}
SSRE(pstart)## [,1]
## [1,] 0.06851034#par<- c(pstart,1-sum(pstart))
#vr <- sqrt(t(par)%*% cov_matrix%*% par)
opt_1 <- optim(pstart,SSRE)
opt_1## $par
## [1] 2.831649e-02 4.470701e-06 2.141530e-01 1.270327e-01 9.970038e-05
## [6] 2.654679e-01 4.358655e-02 1.120192e-01
##
## $value
## [1] 0.04817881
##
## $counts
## function gradient
## 501 NA
##
## $convergence
## [1] 1
##
## $message
## NULLpara <- c(opt_1$par, 1-sum(opt_1$par)) # final contributions
weight<- as.data.frame(t(round(para,4)))
colnames(weight)<- tickers
display <- cbind(t(weight),'Ave.Return'=mu,'Stdev'=std,'Sharp'=sharpe)
colnames(display)[1]<-'Weight'
data.frame(display)## Weight Ave.Return Stdev Sharp
## AAPL 0.0283 0.0138927012 0.10329584 0.134494296
## AMD 0.0000 0.0729485380 0.23835453 0.306050566
## ADI 0.2142 0.0025407452 0.06428978 0.039520204
## ABBV 0.1270 -0.0153324794 0.09888394 -0.155055300
## AET 0.0001 0.0148118327 0.04606748 0.321524677
## A 0.2655 -0.0004756023 0.06037735 -0.007877164
## APD 0.0436 -0.0032581455 0.03793639 -0.085884435
## AA 0.1120 -0.0403532553 0.08293752 -0.486550045
## CF 0.2093 0.0068655083 0.09323691 0.073635089cat('\n','The monthly return of the portfolio at the minimum global variance is:',max_value <-sprintf("%1.3f%%", para%*%mu*100),'\n')##
## The monthly return of the portfolio at the minimum global variance is: -0.436%cat('\n','The minimum global variance is:',min_vr <-opt_1$value,'\n')##
## The minimum global variance is: 0.04817881#or
#opt_1 <-nlminb(pstart, SSRE)
############################################
#Maximize the return
############################################
MaxReturn <- function(parx) {
par<- c(parx,1-sum(parx))
if(all(par > 0 & par < 1)) { # parameters meet requirements
-t(par)%*%mu # this is a linear algebra version of your objective without the division by xi
} else 1e7 # penalty for parameters not satisfying constraints
}
-MaxReturn (pstart)## [,1]
## [1,] 0.006019697opt_2 <-nlminb(pstart, MaxReturn)
opt_2## $par
## [1] 9.573517e-02 1.999590e-01 7.570084e-02 4.415754e-02 9.735729e-02
## [6] 7.037748e-02 6.546675e-02 1.574947e-15
##
## $objective
## [1] -0.01903878
##
## $convergence
## [1] 1
##
## $iterations
## [1] 43
##
## $evaluations
## function gradient
## 83 352
##
## $message
## [1] "false convergence (8)"para_2 <- c(opt_2$par, 1-sum(opt_2$par)) # final contributions
weight<- as.data.frame(t(round(para_2,4)))
colnames(weight)<- tickers
data.frame('Weight'=t(weight))## Weight
## AAPL 0.0957
## AMD 0.2000
## ADI 0.0757
## ABBV 0.0442
## AET 0.0974
## A 0.0704
## APD 0.0655
## AA 0.0000
## CF 0.3512cat('\n','The maximized monthly return of the portfolio is:',max_value <-sprintf("%1.3f%%", -opt_2$objective*100),'\n')##
## The maximized monthly return of the portfolio is: 1.904%cat('\n','The global variance is:',vr_mxreturn <- sqrt(t(para_2)%*%cov_matrix%*%para_2),'\n')##
## The global variance is: 0.09151696The average monthly return of the portofolio at the evenly distributed allocation is 6.8 %. After optimization, the average monthly return of portfolio is -0.436 % when the global variance is at minimum 0.048. The maximized monthly return of portfolio is 1.904 % when the global variance is 0.0915.
4.Linear Programming - Minimax Model
The minimax model will maximize return with respect to one of these prior distributions providing valuable insight regarding an investor’s risk attitude and decision behavior.
Suppose data are observed for N securities, over T time periods. Let
\(y_{jt}\) = Return on one dollar invested in security j in time period t.
\(\bar{y_j}\)=Average Return on security j
\(w_j=\) Portfolio allocation to security j.
\(y_{pt}\) = Return on portfolio in time period t
\(E_p\) =Average Return on portfolio
\(M_p\) = Minimum return on portfolio
The objective function:
\[max\ M_p\]
subject to:
\[\displaystyle\sum_{j=1}^{N}w_jy_{jt} \geq M_p \] t = 1,…,T
\[\displaystyle\sum_{j=1}^{N}w_j\bar{y_j} \geq G \]
#the worst monthly return was calculated by subtract average return by the standard deviation
worst_return<-mu-std
#The subset of the data as candidate
sub_return <- monthly_return[,c('Date','Month','Year',tickers)]
# convert to matrix and vector for later use
cov_matrix <- cov(sub_return[sub_return$Date > as.Date('2018-01-01'),][,-c(1:3)])
pstart <- rep(1/12,8)
# arithmetic mean of monthly return
mu <- apply(sub_return[sub_return$Date > as.Date('2018-01-01'),][,-c(1:3)],2,mean)
# standard deviation
std <- apply(sub_return[sub_return$Date > as.Date('2018-01-01'),][,-c(1:3)],2,sd)
#sharp ratio
sharpe <- mu /std
worst_return <- mu-std
############################################
#Minimize the variance
############################################
portfolio_return <- function(parx) {
par<- c(parx,1-sum(parx))
if(all(par > 0 & par < 1)) { # parameters meet requirements
-par%*%worst_return # this is a linear algebra version of your objective without the division by xi
} else 1e7 # penalty for parameters not satisfying constraints
}
portfolio_return(pstart)## [,1]
## [1,] 0.08607118opt_3 <- optim(pstart, portfolio_return)
opt_3## $par
## [1] 4.354848e-02 1.180378e-05 2.235999e-01 1.010767e-04 1.569046e-01
## [6] 2.356440e-01 3.175448e-02 4.985266e-04
##
## $value
## [1] 0.06492427
##
## $counts
## function gradient
## 501 NA
##
## $convergence
## [1] 1
##
## $message
## NULLpara_3 <- c(opt_3$par, 1-sum(opt_3$par)) # final contributions
weight<- as.data.frame(t(round(para_3,4)))
colnames(weight)<- tickers
display <- cbind(t(weight),'Worst.Return'=worst_return,'Ave.Return'=mu,'Stdev'=std,'Sharp'=sharpe)
colnames(display)[1]<-'Weight'
data.frame(display)## Weight Worst.Return Ave.Return Stdev Sharp
## AAPL 0.0435 -0.08940314 0.0138927012 0.10329584 0.134494296
## AMD 0.0000 -0.16540599 0.0729485380 0.23835453 0.306050566
## ADI 0.2236 -0.06174904 0.0025407452 0.06428978 0.039520204
## ABBV 0.0001 -0.11421642 -0.0153324794 0.09888394 -0.155055300
## AET 0.1569 -0.03125565 0.0148118327 0.04606748 0.321524677
## A 0.2356 -0.06085295 -0.0004756023 0.06037735 -0.007877164
## APD 0.0318 -0.04119454 -0.0032581455 0.03793639 -0.085884435
## AA 0.0005 -0.12329078 -0.0403532553 0.08293752 -0.486550045
## CF 0.3079 -0.08637140 0.0068655083 0.09323691 0.073635089cat('\n','The average minimum losses is:',min_losses <-sprintf("%1.3f%%", opt_3$value*100,'\n'))##
## The average minimum losses is: 6.492%cat('\n','The variance is:',SSRE(para_3[-12]))##
## The variance is: 1e+07alocation <- rbind('worst_return'=worst_return,'weight'=weight)
colnames(alocation) <- tickers
kable(alocation)| AAPL | AMD | ADI | ABBV | AET | A | APD | AA | CF | |
|---|---|---|---|---|---|---|---|---|---|
| worst_return | -0.0894031 | -0.165406 | -0.061749 | -0.1142164 | -0.0312557 | -0.060853 | -0.0411945 | -0.1232908 | -0.0863714 |
| weight | 0.0435000 | 0.000000 | 0.223600 | 0.0001000 | 0.1569000 | 0.235600 | 0.0318000 | 0.0005000 | 0.3079000 |
| Average monthly | return is 8. | 6%, After op | timization, | mininum avera | ge loss is 6. | 49 % when va | riance is 1e+ | 07. |
5. linear programming vs log returns:
Modeling linear vs log returns: Now we are ready to obtain the sample estimates from the returns\(\mathbf{x}_t\) \[ \begin{align} \hat{\boldsymbol{\mu}} & = \frac{1}{T}\sum_{t=1}^T \mathbf{x}_t\\ \hat{\boldsymbol{\Sigma}} & = \frac{1}{T-1}\sum_{t=1}^T (\mathbf{x}_t - \hat{\boldsymbol{\mu}})(\mathbf{x}_t - \hat{\boldsymbol{\mu}})^T \end{align} \]
Daily rebalancing We will start with a daily rebalancing since we already have the daily returns readily available.
# compute log-returns and linear returns
X_log <- CalculateReturns(prices, "log")[-1]
X_lin <- CalculateReturns(prices)[-1]
N <- ncol(X_log) # number of stocks
T <- nrow(X_log) # number of days
# split data into training and test data
T_trn <- round(0.7*T) # 70% of data
X_log_trn <- X_log[1:T_trn, ]
X_log_tst <- X_log[(T_trn+1):T, ]
X_lin_trn <- X_lin[1:T_trn, ]
X_lin_tst <- X_lin[(T_trn+1):T, ]
# Method 1: directly from linear returns
mu_lin <- colMeans(X_lin_trn)
X_ <- X_lin_trn - matrix(mu_lin, T_trn, N, byrow = TRUE) #remove mean
Sigma_lin <- 1/(T_trn-1) * t(X_) %*% X_
# Method 2: directly from log returns
mu_log <- colMeans(X_log_trn)
Sigma_log <- cov(X_log_trn)
# Method 3: from log returns plus transformation
momentsReturnLog2Lin <- function(mu, Sigma) {
K <- ncol(Sigma)
mu_ <- exp(mu + 0.5*diag(Sigma)) - 1
Sigma_ <- matrix(NA, nrow=K, ncol=K)
for(ii in 1:K)
for(jj in 1:K)
Sigma_[ii,jj] <- exp( mu[ii] + mu[jj] + 0.5*(Sigma[ii,ii]+Sigma[jj,jj]) ) * (exp(Sigma[ii,jj])-1)
return( list(mu=mu_, Sigma=Sigma_) )
}
tmp <- momentsReturnLog2Lin(mu_log, Sigma_log)
mu_log_trans <- tmp$mu
Sigma_log_trans <- tmp$SigmaNow let compute the three corresponding GMV portfolios:
# create function for GMVP
portolioGMVP <- function(Sigma) {
ones <- rep(1, nrow(Sigma))
Sigma_inv_1 <- solve(Sigma, ones) #same as: inv(Sigma) %*% ones
w <- (1/as.numeric(ones %*% Sigma_inv_1)) * Sigma_inv_1
return(w)
}
# compute the three versions of GMVP
w_lin <- portolioGMVP(Sigma_lin)
w_log <- portolioGMVP(Sigma_log)
w_log_trans <- portolioGMVP(Sigma_log_trans)
w_all <- cbind(w_lin, w_log, w_log_trans)
w_all## w_lin w_log w_log_trans
## AAPL 0.133266875 0.133093599 0.132663160
## AMD -0.011305201 -0.010945478 -0.010885562
## ADI 0.075478176 0.074269407 0.074486872
## ABBV 0.070315726 0.065447418 0.065214751
## AET 0.302006519 0.313954359 0.313825748
## A 0.098051485 0.094692532 0.094554219
## APD 0.308204376 0.306635367 0.307459326
## AA 0.004160171 0.003564821 0.003389748
## CF 0.019821873 0.019287975 0.019291738# plot to compare the allocations
barplot(t(w_all), col = c("darkblue","darkcyan", "darkgoldenrod"), legend = colnames(w_all),
main = "Portfolio allocation", xlab = "stocks", ylab = "dollars", beside = TRUE)
By portfolio allocation, AAPL, AET and APD are shown the most positive in investing value, but it is not significate in difference between log and tranformation.
6. Qudratic programming
Modelling portfolio optimization as quadratic program
mu <- colMeans(X_log_trn)
Sigma <- cov(X_log_trn)Uniform portfolio
The uniform portfolio allocates equal weight to each stock: \(\mathbf{w} = \frac{1}{N}\mathbf{1}\)
w_unif <- rep(1/N, N)GMVP
The Global Minimum Variance Portfolio (GMVP) is formulated as \[ \begin{array}{ll} \underset{\mathbf{w}}{\textsf{minimize}} & \mathbf{w}^T\mathbf{\Sigma}\mathbf{w}\\ {\textsf{subject to}} & \mathbf{1}^T\mathbf{w} = 1\\ & \mathbf{w}\ge\mathbf{0} \end{array}\] We can compute it with a solver (the closed-form solution does not exist for w???0):
library(CVXR)##
## Attaching package: 'CVXR'## The following object is masked from 'package:purrr':
##
## is_vector## The following object is masked from 'package:matrixcalc':
##
## vec## The following object is masked from 'package:dplyr':
##
## id## The following object is masked from 'package:stats':
##
## powerportolioGMVP <- function(Sigma) {
w <- Variable(nrow(Sigma))
prob <- Problem(Minimize(quad_form(w, Sigma)),
constraints = list(w >= 0, sum(w) == 1))
result <- solve(prob)
return(as.vector(result$getValue(w)))
}
w_GMVP <- portolioGMVP(Sigma)
w_GMVP## [1] 1.296308e-01 1.290853e-07 6.611604e-02 6.360544e-02 3.146232e-01
## [6] 9.372526e-02 3.118962e-01 2.233130e-03 1.816985e-02Markowitz portfolio
The mean-variance Markowitz portfolio with no shorting is formulated as:
\[\begin{array}{ll} \underset{\mathbf{w}}{\textsf{maximize}} & \boldsymbol{\mu}^T\mathbf{w} -\lambda\mathbf{w}^T\mathbf{\Sigma}\mathbf{w}\\ {\textsf{subject to}} & \mathbf{1}^T\mathbf{w} = 1\\ & \mathbf{w}\ge\mathbf{0} \end{array}\]
portolioMarkowitz <- function(mu, Sigma, lmd = 0.5) {
w <- Variable(nrow(Sigma))
prob <- Problem(Maximize(t(mu) %*% w - lmd*quad_form(w, Sigma)),
constraints = list(w >= 0, sum(w) == 1))
result <- solve(prob)
return(as.vector(result$getValue(w)))
}
w_Markowitz <- portolioMarkowitz(mu, Sigma, lmd = 2)Return-risk tradeoff for all portfolios
We can now compare the allocations of the portfolios:
# put together all portfolios
w_all <- cbind(w_unif, w_GMVP, w_Markowitz)
rownames(w_all) <- colnames(X_lin)
colnames(w_all) <- c("uniform", "GMVP", "Markowitz")
round(w_all, digits = 2)## uniform GMVP Markowitz
## AAPL 0.11 0.13 0.38
## AMD 0.11 0.00 0.00
## ADI 0.11 0.07 0.00
## ABBV 0.11 0.06 0.35
## AET 0.11 0.31 0.00
## A 0.11 0.09 0.00
## APD 0.11 0.31 0.00
## AA 0.11 0.00 0.27
## CF 0.11 0.02 0.00barplot(t(w_all),
main = "Portfolio allocation", xlab = "stocks", ylab = "dollars", beside = TRUE,
legend = c(colnames(w_all)),
col = c("darkblue", "brown", "darkcyan"))
Then we can compare the performance:
# compute returns of all portfolios
ret_all <- xts(X_lin %*% w_all, index(X_lin))
ret_all_trn <- ret_all[1:T_trn, ]
ret_all_tst <- ret_all[-c(1:T_trn), ]
# performance
table.AnnualizedReturns(ret_all_trn)## uniform GMVP Markowitz
## Annualized Return 0.3301 0.2820 0.5328
## Annualized Std Dev 0.1653 0.1279 0.1921
## Annualized Sharpe (Rf=0%) 1.9970 2.2045 2.7741table.AnnualizedReturns(ret_all_tst)## uniform GMVP Markowitz
## Annualized Return 0.1430 0.1373 -0.2311
## Annualized Std Dev 0.1950 0.1329 0.2239
## Annualized Sharpe (Rf=0%) 0.7333 1.0329 -1.0321{ chart.CumReturns(ret_all, main = "Performance of different portfolios",
wealth.index = TRUE, legend.loc = "topleft", colorset = rich8equal)
addEventLines(xts("training", index(X_lin[T_trn])), srt=90, pos=2, lwd = 2, col = "darkblue") }
Finally, we can plot the expected return vs the standard deviation along with the efficient frontier:
# first, compute the efficient frontier
w_frontier_trn <- NULL
lmd_sweep <- exp(seq(-5, 5, by = 0.5))
for (lmd in lmd_sweep)
w_frontier_trn <- cbind(w_frontier_trn, portolioMarkowitz(mu, Sigma, lmd))
ret_frontier_trn <- xts(X_lin_trn %*% w_frontier_trn, index(X_lin_trn))
mu_frontier_trn <- table.AnnualizedReturns(ret_frontier_trn)[1, ]
sd_frontier_trn <- table.AnnualizedReturns(ret_frontier_trn)[2, ]
# plot in-sample sd-mu scatter plot
#maxSR <- table.AnnualizedReturns(ret_all_trn[, "maxSR"])[3, ]
chart.RiskReturnScatter(ret_all_trn,
main = "Annualized Return and Risk",
symbolset = c(rep(21, 9), rep(22, 4+3)),
colorset = c(rep("red", 9), rep("blue", 4), rep("green", 3)),
bg = "black")
#add.sharpe = maxSR)
lines(sd_frontier_trn, mu_frontier_trn)
From the graph of annualized Return and Risk, Markowitz optimization has larger annualized return and risk.
Conclusion
We can conclude with the following points:
To compare the average yearly return, AMD and AET are the most increasing return yearly. It seems that using linear returns or log-returns does not make any significant difference for daily, weekly, and monthly returns. Mean-variance Markowitz portfolio has larger annualized return and risk, and Global Minimum Variance Portfolio (GMVP) has lower annualized return and risk.