R - Efficient Frontier and Optimal Risky Portfolio of Several Assets
1. Introduction
We are still following up on our Risk and Return demo, CAL and Correlation tutorials. Remember,
we had a three-stock portfolio composed of Mc Donalds (MCD), General Electric (GE) and Apple
(AAPL) with the following weights: 25%, 25%, 50%. We have added three more securities :
Microsoft (MSFT), Bank of America (BAC) and an ETF on gold (GLD). The weights are no longer
relevant.
This six-stock portfolio is going to be our risky asset. We have already defined our risk-free asset in
the Risk and Return demo: the average (i.e. mean) T-bill rate over the 2010-01-01 to 2020-04-28
period.
The weights that we have been using so far to build our risky portfolio were chosen randomly. We
will now see how to find the optimal weight of each security so that we can build the optimal risky
portfolio.
2. Opening your Script
If you did not save the script that you wrote for the Risk and Return demo, don’t worry, you can do
it again… or use mine (note that we actually do not need the first part on Netflix, I included it just in
case some of you would want to have it):
—————————————————
library(quantmod)
#It is good practice to write a line about the code you want to write
below. Here we are going to download stock prices of Netflix
getSymbols("NFLX")
#Let's see what the last few trading days look like
tail(NFLX)
#Generate summary statistics
summary(NFLX)
#Create a chart of NFLX price time series
chartSeries(NFLX)
#Get annual returns of NFLX
annualReturn(NFLX)
1
#Get monthly arithmetic returns of NFLX in %
annualReturn(NFLX)*100
#For ease of writing codes from now on, let's assign a name to the values
returned by annualReturn(NFLX)*100
NFLX_a_retruns <- annualReturn(NFLX)*100
#Get the variance and the standard deviation of the annual returns of
NFLX
var(NFLX_a_retruns)
sd(NFLX_a_retruns)
#Create a graph of the annual returns of NFLX
chartSeries(NFLX_a_retruns)
#Now let's move to daily returns
dailyReturn(NFLX)
#let's assign a name to the values returned by dailyReturn(NFLX)
NFLX_d_retruns <- dailyReturn(NFLX)
#Create a graph of the daily returns of NFLX using a different layout
plot(NFLX_d_retruns)
#Create a histogram of daily returns of NFLX to see their distribution
hist(NFLX_d_retruns)
#Get a better visual of the histogram by increasing the number of breaks
(bins) and reducing the limits of the x-axis
hist(NFLX_d_retruns, breaks=200, xlim=c(-0.2,0.2))
#Transform the y-axis in Density (not frequency)
hist(NFLX_d_retruns, freq=FALSE, breaks=200, xlim=c(-0.2,0.2))
#Add a normal distribution curve with same mean and SD as NFLX daily
returns
curve(dnorm(x, mean=mean(NFLX_d_retruns), sd=sd(NFLX_d_retruns)),
add=TRUE)
#Kurtosis and Skewness values
install.packages("moments")
library(moments)
kurtosis(NFLX_d_retruns)
skewness(NFLX_d_retruns)
#Create a list of ticker symbols to download
tickers=c("GE","MCD","AAPL")
2
#Download data from "2010-01-01" to "2020-04-28"
getSymbols(tickers, from = "2010-01-01", to = "2020-04-28")
#Get monthly returns
mr_MCD <- monthlyReturn(MCD)
mr_GE <- monthlyReturn(GE)
mr_AAPL <- monthlyReturn(AAPL)
#Merge the monthly returns into one object
merged_returns <- merge.xts(mr_MCD, mr_GE, mr_AAPL)
#Assess weight of each stock
w <- c(.25, .25, .50)
library(PerformanceAnalytics)
#Calculate portfolio returns and plot them
portfolio_monthly_returns <- Return.portfolio(merged_returns, weights =
w)
chartSeries(portfolio_monthly_returns)
#Calculate the growth of 1 dollar invested in this portfolio at the
beginning till now and plot it
dollar_growth <- Return.portfolio(merged_returns, weights = w,
wealth.index = TRUE)
plot(dollar_growth)
#Download 13-week T-bills rates from "2010-01-01" to "2020-04-28"
getSymbols("^IRX",from = "2010-01-01", to = "2020-04-28")
#Calculate average risk-free rate, using the "Close" column", removing
NAs and transforming it into a decimal value
riskfree <- mean(IRX$IRX.Close, na.rm = TRUE)/100
#Get Sharpe ratio with risk-free rate = riskfree
SharpeRatio(portfolio_monthly_returns, Rf = riskfree)
#Get Slope of CAL
slope <- CAPM.CML.slope(portfolio_monthly_returns, Rf = riskfree)
#Get Curve of CAL
curve(slope*x+riskfree*100, from=0, to=10, xlab="Standard Deviation",
ylab="Return")
#Calculate average portfolio return and standard deviation
exp_m_return <- mean(portfolio_monthly_returns)*100
sd_m_return <- StdDev(portfolio_monthly_returns)*100
3
#Add the risky portfolio's point on the CAL
points(sd_m_return, exp_m_return, col="red")
#Calculate expected return and SD of a 70% risky/30% risk-free portfolio
assuming SD of riskfree=0
exp_return70_30 <- 0.7*exp_m_return + 0.3*riskfree*100
sd70_30 <- 0.7*sd_m_return
#Add the 70/30 portfolio's point on the CAL
points(sd70_30, exp_return70_30, col="blue")
#THIS IS THE NON-MANDATORY PART: Calculate SD of riskfree so that we have
the coordinates of the 2nd point on the new CAL: risk-free asset
sd_riskfree <- StdDev(IRX$IRX.Close)
#Calculate slope of new CAL
new_slope <- (exp_m_return - riskfree*100)/(sd_m_return - sd_riskfree)
#New CAL equation: y = new_slope*x + b. We want to figure out b. Let's
use x and y from the risky asset
#so that exp_m_retrun = new_slope*sd_m_return + b or
# b = exp_m_return - new_slope*sd_m_return
b_intercept <- exp_m_return - new_slope*sd_m_return
#Get Curve of new CAL
curve(new_slope*x+b_intercept, from=0, to=10, xlab="Standard Deviation",
ylab="Return")
#Add the risky portfolio's point on the CAL
points(sd_m_return, exp_m_return, col="red")
#Calculate expected return and SD of a 70% risky/30% risk-free portfolio
with actual SD of riskfree
exp_return70_30 <- 0.7*exp_m_return + 0.3*riskfree*100
sd70_30 <- 0.7*sd_m_return + 0.3*sd_riskfree
#Add the 70/30 portfolio's point on the CAL
points(sd70_30, exp_return70_30, col=“blue")
#Calculate correlation coefficient between Monthly returns of AAPL and
MCD
cor(mr_AAPL, mr_MCD)
#Calculate correlation coefficient between Monthly returns of GE and MCD
cor(mr_GE, mr_MCD)
4
#Calculate correlation coefficient between Monthly returns of AAPL and GE
cor(mr_AAPL, mr_GE)
#Download Microsoft (MSFT) price data with the same time frame as the one
we used before
getSymbols("MSFT", from = "2010-01-01", to = "2020-04-28")
#Calculate monthly returns of MSFT
mr_MSFT <- monthlyReturn(MSFT)
#Calculate correlation coefficient between Monthly returns of AAPL and
MSFT
cor(mr_AAPL, mr_MSFT)
#Download Bank of America (BAC) price data with the same time frame as
the one we used before
getSymbols("BAC", from = "2010-01-01", to = "2020-04-28")
#Calculate monthly returns of MSFT
mr_BAC <- monthlyReturn(BAC)
#Download Gold ETF (GLD) price data with the same time frame as the one
we used before
getSymbols("GLD", from = "2010-01-01", to = "2020-04-28")
#Calculate monthly returns of MSFT
mr_GLD <- monthlyReturn(GLD)
#Install and load package corrplot to create correlation matrix
install.packages("corrplot")
library(corrplot)
#We want to create a unique dataset with all the monthly returns of each
of the stock in our portfolio.
#We are going to merge mr_AAPL, mr_GE, mr_MCD, mr_MSFT, and BAC. But
first we need to rename the column in each "mr_" object.
#otherwise, we will lose their name in the data set.
colnames(mr_AAPL) <- "AAPL"
colnames(mr_GE) <- "GE"
colnames(mr_MCD) <- "MCD"
colnames(mr_MSFT) <- "MSFT"
colnames(mr_BAC) <- "BAC"
colnames(mr_GLD) <- "GLD"
#Create our dataset
risky_portfolio <- Reduce(merge, list(mr_AAPL, mr_GE, mr_MCD, mr_MSFT,
mr_BAC, mr_GLD))
5
#Calculate correlation matrix
risky_portfolio.cor <- cor(risky_portfolio)
#Get correlogram of the stocks in our portfolio
corrplot(risky_portfolio.cor)
—————————————————
Highlight what you just wrote (or pasted) and click “run” or cmd + enter key on a mac or ctrl +
enter key on a PC.
3. Building the Efficient Frontier
It is rather difficult (and long) to build the efficient frontier because it requires for us to compute
several “random” portfolios (i.e. with random weights) and select only the best one for every level
of risk (i.e. select the one with the highest return for every level of standard deviation).
Fortunately, there is an R package that deals with that: “PortfolioAnalytics” (https://
www.rdocumentation.org/packages/PortfolioAnalytics).
Install it and load it. Depending on your version of R and RStudio, you may have to install the
following packages along with “PortfolioAnaytics”:
install.packages('foreach')
install.packages("Rglpk")
install.packages("ROI.plugin.quadprog")
install.packages(“ROI.plugin.glpk")
Here are the steps to build the efficient frontier:
a) Set the portfolio specification
First we have to let “PortfolioAnalytics” knows the specifications of our portfolio: which assets to
use, the constraints and the objectives of our portfolio.
We are going to use the “portfolio.spec”, “add.constraint” and “add.objective” functions from
“PortfolioAnalytics”.
#Set Portfolio Specification: Which assets?
port_spec <- portfolio.spec(asset = colnames(risky_portfolio))
6
We are going to store the portfolio specifications in an object named port_spec.
Question 1: Using port_spec <- add.constraint(…) and port_spec <- add.objective(…), add the
following specifications to our portfolio:
a) type of constraint: “long_only” (we want no short-selling),
b) type of objective #1: ‘risk’, defined by ‘var’ (we want to minimize the risk)
c) type of objective #2: ‘return’, defined by ‘mean’ (we want to maximize the return)
Use ?add.constraint and ?add.objective to get some information on those functions or see some
examples.
TO INCLUDE IN YOUR EXAM COPY: Copy and paste your code in a .doc file along withe the
output of print(port_spec).
b) Computing the Efficient Frontier
Now that “PortfolioAnalytics” knows the specifications of our portfolio, we can calculate the
efficient frontier:
#Compute efficient frontier
efficient.frontier <- create.EfficientFrontier(
R = risky_portfolio,
portfolio = port_spec,
type = "mean-var",
n.portfolios = 30)
Note that we are using 30 portfolios to plot on the efficient frontier because it is a number large
enough to get a precise frontier but yet small enough so that we can still identify the portfolios on
the frontier.
c) Plotting the Efficient Frontier
We can now create the graph of the efficient frontier and the CAL optimal using our risk-free asset
(we are going to assume 0 standard deviation for ou risk-free asset).
#Plot efficient frontier
chart.EfficientFrontier(efficient.frontier,
match.col="StdDev",
RAR.text="Sharpe Ratio",
rf = riskfree,
tangent.line = TRUE,
chart.assets=TRUE,
labels.assets=TRUE,
xlim=c(0,0.15),
ylim=c(-0.01,0.05))
7
I have used the xlim and ylim (i.e. scales of the axis) so that every single stock can be displayed).
Note that GE has a negative expected return while having one of the highest standard deviation
amongst our stock sample. Most likely, GE will not be part of the optimal portfolio (i.e. its weight
will be = 0).
8
d) Finding the Optimal Risky Portfolio
There exist different optimisers that will return slightly different weights for the optimal portfolio.
We are going to rely on the one embedded in the chart.EfficientFrontier function.
#Get weights of portfolios along the frontier
summary(efficient.frontier, digits = 4)
We can see the weight for each stock in the 30 portfolios along the efficient frontier as well as their
corresponding mean return and standard deviation. For instance, we learn that the last portfolio
(30th) is composed entirely of AAPL, which means that the single stock of AAPL lies right on the
efficient frontier.
Although chart.EfficientFrontier is capable of identifying the tangency portfolio (i.e. the optimal
risky portfolio), I haven’t found a quick way to extract its corresponding weights directly from this
function… So for the lack of a better method, we can “manually” identify the tangency portfolio by
noting its position along the frontier. Here, we see that it is in 9th position, starting from the tail of
the frontier. Remark that choosing 4 digits was a bit of an overkill. We don’t need that sort of
precision (2 would have been enough). So portfolio #22 seems to be a good approximation of the
optimal risky portfolio. It is composed of 49% AAPL, 20% MCD, and 31% MSFT and yields a
mean monthly return of 1.7% with a 4.9% standard deviation. Its Sharpe Ratio (as displayed on the
graph) is 0.224.
**************************************************
PortfolioAnalytics Efficient Frontier
**************************************************
Call:
create.EfficientFrontier(R = risky_portfolio, portfolio = port_spec,
type = "mean-var", n.portfolios = 30)
Efficient Frontier Points: 30
Weights along the efficient frontier:
AAPL GE MCD MSFT BAC GLD
1 0.0020 0.0870 0.3503 0.1209 0.0248 0.4150
2 0.0179 0.0726 0.3552 0.1280 0.0239 0.4023
3 0.0338 0.0582 0.3602 0.1352 0.0230 0.3895
4 0.0497 0.0439 0.3652 0.1424 0.0220 0.3768
5 0.0656 0.0295 0.3702 0.1495 0.0211 0.3641
6 0.0815 0.0151 0.3752 0.1567 0.0202 0.3513
7 0.0974 0.0007 0.3801 0.1638 0.0193 0.3386
9
8 0.1195 0.0000 0.3799 0.1761 0.0091 0.3153
9 0.1422 0.0000 0.3792 0.1879 0.0000 0.2908
10 0.1668 0.0000 0.3768 0.1947 0.0000 0.2617
11 0.1914 0.0000 0.3743 0.2015 0.0000 0.2327
12 0.2161 0.0000 0.3719 0.2083 0.0000 0.2037
13 0.2407 0.0000 0.3695 0.2151 0.0000 0.1747
14 0.2653 0.0000 0.3671 0.2220 0.0000 0.1457
15 0.2899 0.0000 0.3646 0.2288 0.0000 0.1166
16 0.3146 0.0000 0.3622 0.2356 0.0000 0.0876
17 0.3392 0.0000 0.3598 0.2424 0.0000 0.0586
18 0.3638 0.0000 0.3574 0.2493 0.0000 0.0296
19 0.3884 0.0000 0.3550 0.2561 0.0000 0.0005
20 0.4211 0.0000 0.3054 0.2736 0.0000 0.0000
21 0.4538 0.0000 0.2549 0.2913 0.0000 0.0000
22 0.4866 0.0000 0.2044 0.3090 0.0000 0.0000
23 0.5194 0.0000 0.1539 0.3267 0.0000 0.0000
24 0.5522 0.0000 0.1034 0.3444 0.0000 0.0000
25 0.5850 0.0000 0.0529 0.3621 0.0000 0.0000
26 0.6178 0.0000 0.0024 0.3798 0.0000 0.0000
27 0.7111 0.0000 0.0000 0.2889 0.0000 0.0000
28 0.8074 0.0000 0.0000 0.1926 0.0000 0.0000
29 0.9037 0.0000 0.0000 0.0963 0.0000 0.0000
30 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000
Risk and return metrics along the efficient frontier:
mean StdDev out
1 0.0070 0.0300 0.0009
2 0.0075 0.0300 0.0009
3 0.0080 0.0301 0.0009
4 0.0085 0.0303 0.0009
5 0.0090 0.0306 0.0009
6 0.0094 0.0309 0.0010
7 0.0099 0.0313 0.0010
8 0.0104 0.0318 0.0010
9 0.0109 0.0325 0.0011
10 0.0114 0.0333 0.0011
11 0.0118 0.0342 0.0012
12 0.0123 0.0352 0.0012
13 0.0128 0.0364 0.0013
14 0.0133 0.0377 0.0014
15 0.0138 0.0391 0.0015
16 0.0142 0.0405 0.0016
17 0.0147 0.0421 0.0018
18 0.0152 0.0437 0.0019
19 0.0157 0.0454 0.0021
20 0.0162 0.0472 0.0022
21 0.0166 0.0491 0.0024
22 0.0171 0.0512 0.0026
10
23 0.0176 0.0534 0.0029
24 0.0181 0.0557 0.0031
25 0.0185 0.0581 0.0034
26 0.0190 0.0606 0.0037
27 0.0195 0.0635 0.0040
28 0.0200 0.0669 0.0045
29 0.0205 0.0709 0.0050
30 0.0209 0.0754 0.0057
3. Comparing our Optimal Risky Portfolio to a Benchmark
Now that we have found the optima risky portfolio we want to see how it has performed compared
to a given benchmark. Since our portfolio is composed of large US stocks, it makes sense to use the
S&P 500 index as a benchmark.
Question 2:
a) Download the index time series of the S&P 500 with the same time frame as the one used to
construct our optimal risky portfolio.
b) Calculate the monthly returns of the S&P 500 over the time period and assign them the name
“mr_SP500”
c) Calculate the growth of 1 dollar invested in our optimal risky portfolio at the beginning of our
time frame till the end of it and plot it (we used a similar approach in a previous tutorial).
d) Calculate the growth of 1 dollar invested in a S&P500 portfolio at the beginning of our time
frame till the end of it and plot it on the same plot as the one used in c) (the “line” function will help
you).
TO INCLUDE IN YOUR EXAM COPY: Copy and paste your code in your .doc file along withe
the plot.
4. Using a Single-Factor Model
With the help of a powerful software/language like R (or Python, or Matlab), building efficient
frontiers and hence finding the optimal risky portfolio has become relatively easy. But the SingleFactor
Model can provide with a different way of getting to the optimal risky portfolio.
First, let’s plot the Security Characteristic Line of each individual stock. In the package
“PerformanceAnalytics” (we have used it before), there is a function for that:
11
library(PerformanceAnalytics)
#Plot the SCL of AAPL
chart.Regression(
mr_SP500,
mr_AAPL,
Rf = riskfree,
excess.returns = TRUE,fit = "linear",
reference.grid = TRUE,
main = "AAPL's SCL")
12
Question 3:
Do the same for MCD and MSFT
TO INCLUDE IN YOUR EXAM COPY: Copy and paste your code in your .doc file along withe
the 2 plots.
Then, we can retrieve the betas, alphas, and epsilons of the stocks within the risky_portfolio by:
#Get betas, alphas, and epsilons
betas <- CAPM.beta(risky_portfolio, mr_SP500, Rf = riskfree)
alphas <- CAPM.alpha(risky_portfolio, mr_SP500, Rf = riskfree)
epsilons <- CAPM.epsilon(risky_portfolio, mr_SP500, Rf = riskfree)
NOT MANDATORY: Please find below the rest of the methodology to get to the optimal risky
portfolio weights. Note that the model is far from perfect.
#Calculate variance of mr_SP500
var_SP500 <- var(mr_SP500)
#Calculate optimal weights
w_GE <- alphas[,"GE"] / (var(mr_GE) - var_SP500)
w_MCD <- alphas[,"MCD"] / (var(mr_MCD) - var_SP500)
w_MSFT <- alphas[,"MSFT"] / (var(mr_MSFT) - var_SP500)
w_BAC <- alphas[,"BAC"] / (var(mr_BAC) - var_SP500)
w_GLD <- alphas[,"GLD"] / (var(mr_GLD) - var_SP500)
#Scale initial weights to 1 (only positive ones because of longonly
constraint)
weight_AAPL <- w_AAPL / (w_AAPL + w_MCD + w_MSFT)
weight_MCD <- w_MCD / (w_AAPL + w_MCD + w_MSFT)
weight_MSFT <- w_MSFT / (w_AAPL + w_MCD + w_MSFT)
———————————————————-