Portfolio Optimization, Tangency Portfolio, CML with shorts not allowed.

The functions in the previous post have been updated to put in the additional constraints that no portfolio weights can be negative, that is short selling is not allowed.

setwd("C:/Users/admin/Desktop/Data Analytics/QuantitativeFinance/QFLabs")
getwd()

library(ggplot2)
library(quadprog)

GetData20 <- function(CSVfile){
  baseData <- read.csv(CSVfile)
  N0 <- ncol(baseData)
  Close <- baseData[, 3:N0]
  logRet1 <- log(head(Close, -1) / tail(Close, -1))
  return(logRet1)
}

PlotData20 <- function(logRet1){
  Retn <- colMeans(logRet1)
  Risk <- diag(var(logRet1))
  RiskReturn <- as.data.frame(t(rbind(Retn,Risk)))
  plot1 <- ggplot(data = RiskReturn,aes(x = Risk, y = Retn) )
  plot1 <- plot1 + geom_point()
  plot1 <- plot1 + xlab("Risk / Variance") + ylab("Daily Returns") + ggtitle("Risk/Returns")
  plot1
}

rOptPort20 <- function(hRets,pRet, Shorts = TRUE){
  Dmat <- 2*cov(hRets)
  dvec <- rep(0,ncol(hRets))
  if (Shorts){
    ##message("Shorts : TRUE")    ## Shorts are allowed
    Amat <- cbind(rep(1,ncol(hRets)),colMeans(hRets))
    bvec <- c(1,pRet)
  } else {
    ##message("Shorts = FALSE")    ## Shorts are not allowed
    Amat <- cbind(rep(1,ncol(hRets)),colMeans(hRets),diag(1,nrow = ncol(hRets)))
    bvec <- c(1,pRet,rep(0,ncol(hRets)))
  }
  result = tryCatch({
    solve.QP(Dmat = Dmat, dvec = dvec, Amat = Amat, bvec = bvec, meq =2)
  }, warning = function(w){
    message("warning : ",pRet)
    return(NULL)
  }, error = function(w){
    message("infeasible : ", pRet)
    return(NULL)
  }, finally = {
  }
  )
  if (!is.null(result)){
    wP <- result$solution
    varP <- result$value
    
  }
  else{
    wP <- "error"
    varP <- "error"
  }
  retList <- list(wP,varP)
  names(retList) <- c("wP","varP")
  return(retList)
}

EFMinVar20 <- function(hRets, minRet, maxRet,Shorts = FALSE){
  smuP <- seq(minRet,maxRet,length=50)
  svarP <- sapply(smuP,function(x) rOptPort20(hRets,x,Shorts)$varP)
  EffF <- as.data.frame(cbind(smuP,svarP))
  EffF0 <- as.data.frame(EffF[EffF$svarP != "error",])
  EffF0 <- as.data.frame(apply(EffF0, 2, FUN = function(x) as.numeric(as.character(x))))
  minVar <- min(EffF0$svarP)
  L <- EffF0$svarP == minVar
  minRet <- EffF0[L,]$smuP
  minPoint <- as.data.frame(cbind(minRet,minVar))
  minVarwP <- rOptPort20(hRets,minRet,Shorts)$wP
  rList <-list(EffF0,minPoint,minVarwP)
  names(rList) <- c("EFF","minPoint","wP")
  return(rList)
}

EFMinVar20Plot <- function(list1){
  
  ealred  <- "#7D110C"
  
  plot2 <- ggplot(data = list1$EFF,aes(x = svarP, y = smuP) )
  plot2 <- plot2 + geom_point()
  plot2 <- plot2 + geom_point(data = list1$minPoint, aes(x = minVar,y = minRet),color = "red", size=3)
  plot2 <- plot2 + xlab("Variance") + ylab("Returns") + ggtitle("Efficient Frontier - MinVar")
  plot2
}

EFSharpe20 <- function(hRets, minRet, maxRet,RF, Shorts = FALSE){
  smuP <- seq(minRet,maxRet,length=50)
  svarP <- sapply(smuP,function(x) rOptPort20(hRets,x,Shorts)$varP)
  EffF <- as.data.frame(cbind(smuP,svarP))
  EffF0 <- as.data.frame(EffF[EffF$svarP != "error",])
  EffF0 <- as.data.frame(apply(EffF0, 2, FUN = function(x) as.numeric(as.character(x))))
  sharpe <- (EffF0$smuP-RF)/EffF0$svarP
  EFF <- as.data.frame(cbind(EffF0,sharpe,RF))
  L <- EFF$sharpe == max(EFF$sharpe)
  maxSharpe <- EFF[L,]
  wTP <- rOptPort20(hRets,maxSharpe$smuP,Shorts)$wP
  rList <-list(EFF,maxSharpe,wTP)
  names(rList) <- c("EFF","maxSharpe","wTP")
  return(rList)
}

EFSharpe20Plot <- function(list1){
  
  
  plot2 <- ggplot(data = list1$EFF,aes(x = svarP, y = smuP) )
  plot2 <- plot2 + geom_point()
  plot2 <- plot2 + geom_point(data = list1$maxSharpe, aes(x = svarP,y = smuP),colour = "red", pch =24, size=3)
  plot2 <- plot2 + geom_point(data = list1$maxSharpe, aes(x = 0,y = RF),color = "red", pch =24, size=3)
  plot2 <- plot2 + xlab("Variance") + ylab("Returns") + ggtitle("Efficient Frontier - Sharpe")
  plot2 <- plot2 + geom_abline(intercept = (list1$maxSharpe)$RF, slope = (list1$maxSharpe)$sharpe, colour = "red")
  plot2
}

#Testing

logRet <- GetData20("US5.csv")
##logRet <- GetData20("India5.csv")
PlotData20(logRet)

## Shorts = FALSE

z10 <- EFMinVar20(logRet,-0.0005,0.0002, Shorts = FALSE)
z10$wP
cminRet <- (z10$minPoint)$minRet
cminRet
z11 <- rOptPort20(logRet,cminRet,Shorts = FALSE)
z11$wP
EFMinVar20Plot(z10)
cminRet <- (z10$minPoint)$minRet
maxR <- cminRet + max(0.001,cminRet+ 2*abs(cminRet))
z10a <- EFMinVar20(logRet,cminRet, maxR, Shorts = FALSE)
EFMinVar20Plot(z10a)

z12 <- EFSharpe20(logRet,cminRet,maxR,0.0001,Shorts = FALSE)
z12$wTP
sum(z12$wTP)
(z12$maxSharpe)$smuP
(z12$maxSharpe)$svarP
maxSharpeRet <- (z12$maxSharpe)$smuP
EFSharpe20Plot(z12)

After the data is downloaded from the web and stored in  CSV file as explained in an earlier post. define the functions and test the program as follows :

> logRet <- GetData20("US5.csv")
> ##logRet <- GetData20("India5.csv")
> PlotData20(logRet)

Unlike the case were shorts are allowed, in this case, not every desired portfolio return is feasibile when shorts are not allowed. Here we search through a range of returns from -0.0005 to + 0.0002 and detect which return gives us the minimum variance. The means that are infeasible are flagged as sucn and ignored.
> z10 <- EFMinVar20(logRet,-0.0005,0.0002, Shorts = FALSE)
infeasible : -5e-04
.....
infeasible : -0.000242857142857143

Dataframe z10 contains the range of feasible returns and the corresponding variance. The portofolio weights are printed out as follows. The first one though technically negative is very, very small.
> z10$wP
[1] -1.386988e-17  4.295438e-02  1.356919e-01  4.977252e-01  3.236285e-01

we note the return achieved at the point of minimum variance
> cminRet <- (z10$minPoint)$minRet
> cminRet
[1] 0.0001285714

> EFMinVar20Plot(z10)

The efficient frontier is plotted, but this is not very useful because the range of returns that we have taken are "below" the minvariance point.

So in the next few lines, we have defined a range of returns that goes above, or higher than, the return for minimum variance

> cminRet <- (z10$minPoint)$minRet
> maxR <- cminRet + max(0.001,cminRet+ 2*abs(cminRet))
> z10a <- EFMinVar20(logRet,cminRet, maxR, Shorts = FALSE)
infeasible : 0.000761224489795919
.......
infeasible : 0.00112857142857143
> EFMinVar20Plot(z10a)

Next we identify the maximum Sharpe ratio, identify the tangency portfolio and draw the capital market line. Note that the risk free rate is 0.0001

> z12 <- EFSharpe20(logRet,cminRet,maxR,0.0001,Shorts = FALSE)
infeasible : 0.000761224489795919
....
infeasible : 0.00112857142857143
> z12$wTP
[1]  2.515328e-01  0.000000e+00 -1.456853e-17  7.484672e-01 -6.938894e-18
> sum(z12$wTP)
[1] 1
> (z12$maxSharpe)$smuP
[1] 0.0005163265
> (z12$maxSharpe)$svarP
[1] 0.0002621069
> maxSharpeRet <- (z12$maxSharpe)$smuP
> EFSharpe20Plot(z12)

Please note that by setting SHORTS = TRUE, this set of functions should be able to recreate the results given in the previous post.

-------------------------------------------------------------------------------------------
These programs can be run on StatAce, the free hosted R environment and get the same result.

Portfolio Optimisation, Tangency Portfolio and CML -- shorts allowed

I wish I could document this better, but let me get the code out first. The data used by the program shown in this post is available in a CSV file that can be created in the manner shown in an earlier post.

setwd("C:/Users/admin/Desktop/Data Analytics/QuantitativeFinance/QFLabs")
getwd()
##baseData <- read.csv("US5.csv")
baseData <- read.csv("India5.csv")
N0 <- ncol(baseData)
Close <- baseData[, 3:N0]
logRet <- log(head(Close, -1) / tail(Close, -1))
Retn <- colMeans(logRet)
Risk <- diag(var(logRet))
RiskReturn <- as.data.frame(t(rbind(Retn,Risk)))

library(ggplot2)

plot1 <- ggplot(data = RiskReturn,aes(x = Risk, y = Retn) )
plot1 <- plot1 + geom_point()
plot1 <- plot1 + xlab("Risk / Variance") + ylab("Daily Returns") + ggtitle("Risk/Returns")
plot1

# using portfolio.optim 
library(tseries)
w0 <- portfolio.optim(as.matrix(logRet),pm = 0.005,shorts = TRUE,riskless= FALSE)
w0$pw
sum(w0$pw)

# understanding the code

library(quadprog)

# the basic function that calculates the min variance portfolio with NO SHORTS

rOptPort10 <- function(hRets,pRet){
  Dmat <- 2*cov(hRets)
  dvec <- rep(0,ncol(hRets))
  Amat <- cbind(rep(1,ncol(hRets)),colMeans(hRets))
  bvec <- c(1,pRet)
  result <- solve.QP(Dmat = Dmat, dvec = dvec, Amat = Amat, bvec = bvec, meq =2)
  wP <- result$solution
  varP <- result$value
  retList <- list(wP,varP)
  names(retList) <- c("wP","varP")
  return(retList)
}

# testing out the function with expected return

z <- rOptPort10(logRet,0.005)

z$wP
z$varP
sum(z$wP)

# here we create the Efficient Frontier for a given range of returns

EFMinVar10 <- function(hRets, minRet, maxRet){
  smuP <- seq(minRet,maxRet,length=50)
  svarP <- sapply(smuP,function(x) rOptPort10(hRets,x)$varP)
  EffF <- as.data.frame(cbind(smuP,svarP))  
  minVar <- min(EffF$svarP)
  L <- EffF$svarP == minVar
  minRet <- EffF[L,]$smuP
  minPoint <- as.data.frame(cbind(minRet,minVar))
  minVarwP <- rOptPort10(hRets,minRet)$wP
  rList <-list(EffF,minPoint,minVarwP)
  names(rList) <- c("EFF","minPoint","wP")
  return(rList)
}

# We use the above function to get the points of the Efficient Frontier
# and numerically detect the point of minimum variance

z10 <- EFMinVar10(logRet,-0.005,.005)
z10$wp
cminRet <- (z10$minPoint)$minRet
z11 <- rOptPort10(logRet,cminRet)
z11$wP

# This function plots the Efficient Frontier and showing the point of minimum variance

EFMinVar10Plot <- function(list1){
  
  plot2 <- ggplot(data = list1$EFF,aes(x = svarP, y = smuP) )
  plot2 <- plot2 + geom_point()
  plot2 <- plot2 + geom_point(data = list1$minPoint, aes(x = minVar,y = minRet),color = "red", size=3)
  plot2 <- plot2 + xlab("Variance") + ylab("Returns") + ggtitle("Efficient Frontier - MinVar")
  plot2
}

EFMinVar10Plot(z10)
cminRet <- (z10$minPoint)$minRet 

z10a <- EFMinVar10(logRet,min(0,cminRet),0.001)
EFMinVar10Plot(z10a)


# This function calculates the Max Sharpe Ratio
# and the Tangency Portfolio weights

EFSharpe10 <- function(hRets, minRet, maxRet,RF){
  smuP <- seq(minRet,maxRet,length=50)
  svarP <- sapply(smuP,function(x) rOptPort10(hRets,x)$varP)
  sharpe <- (smuP-RF)/svarP
  EFF <- as.data.frame(cbind(smuP,svarP,sharpe,RF))
  L <- EFF$sharpe == max(EFF$sharpe)
  maxSharpe <- EFF[L,]
  wTP <- rOptPort10(hRets,maxSharpe$smuP)$wP
  rList <-list(EFF,maxSharpe,wTP)
  names(rList) <- c("EFF","maxSharpe","wTP")
  return(rList)
}

z11 <- EFSharpe10(logRet,min(0,cminRet),0.001,0.0001)
z11$wTP
sum(z11$wTP)
(z11$maxSharpe)$smuP
(z11$maxSharpe)$svarP
maxSharpeRet <- (z11$maxSharpe)$smuP

# This function plots the Efficient Fronter and shows
# The Tangency Point and the Capital Market Line

EFSharpe10Plot <- function(list1){
 
  
  plot2 <- ggplot(data = list1$EFF,aes(x = svarP, y = smuP) )
  plot2 <- plot2 + geom_point()
  plot2 <- plot2 + geom_point(data = list1$maxSharpe, aes(x = svarP,y = smuP),colour = "red", pch =24, size=3)
  plot2 <- plot2 + geom_point(data = list1$maxSharpe, aes(x = 0,y = RF),color = "red", pch =24, size=3)
  plot2 <- plot2 + xlab("Variance") + ylab("Returns") + ggtitle("Efficient Frontier - Sharpe")
  plot2 <- plot2 + geom_abline(intercept = (list1$maxSharpe)$RF, slope = (list1$maxSharpe)$sharpe, colour = "red")
  plot2
}

EFSharpe10Plot(z11)

---

DONT PANIC ! This 20 minute video will explain what this code is doing


What does negative weights mean ? It means short selling. So what is short selling ? Let us explain.

If you have a Rs 100 you can invest Rs 80 in Stock A and Rs 20 in Stock B and so the weights are (0.8, 0.2). Note that the weights add up to 1.

 But through a process known as short selling, you could sell Rs 80 worth of Stock B [ that you do not have, at the moment ] get Rs 80. To this you add the Rs 100 that you have and invest Rs 180 [ Rs 100 + Rs 80 ] in Stock A. In this case your weights are (180,-80) or (1.8 and -0.8) and you would note that weights add up to 1 again.

Why would you do short selling ? There are many reasons. For example, you may believe that Stock A will give you very high returns. Also if you sell "short", that is sell shares that you do not have, you will have to, at some point in future, or in the next "period" buy the shares from the market and deliver it to the person to whom you had sold the shares.

Beta, SCL with bulk NSE data

As a continuation of the previous post on calculating beta, we test the model with the closing share prices of 44 companies that are a part of the CNX Nifty index.

Instead of following the Yahoo Finance based data collection strategy shown in the sample program of the text book, we have written our own data extraction routine based on Quandl. Note that
niftStks04_* is a set of four lists that have names of NIFTY companies for which we have data from 2004
niftStks09 is a list of names of NIFTY companies for which we have data from 2009
niftStks13 is a list of names of NIFTY companies for whch we have data from 2013

For the purpose of our exercise we have used data from 2009 onwards so that we have excluded the last list.

As a sanity check, we have extracted the CNX Nifty data, twice, once along with the company share price data and once as a standalone data with the corresponding dates. By laying out this data side-by-side in the Close01 data.frame through a cbind, we can see the dates of the closing prices of the company data.

========================================
   setwd("C:/Users/admin/Desktop/Data Analytics/QuantitativeFinance/QFLabs")
getwd()
library(Quandl)

Quandl.auth("xxxxxxxxxxxxx")

Cstartdate = '2009-06-01'
Cenddate = '2014-06-01'
IIR00 <- Quandl('BCB/17901', trim_start = Cstartdate, trim_end = Cenddate)
IIR00$Value <- IIR00$Value/1200
NIFTY <- Quandl("NSE/CNX_NIFTY", trim_start = Cstartdate, trim_end = Cenddate, collapse="monthly")
NIFTY01 <- NIFTY[,c('Date','Close')]

nifStks04_1 <- c("ACC","ASIANPAINT","BANKBARODA","BHEL","CIPLA","DRREDDY","GAIL","GRASIM","HCLTECH","HDFCBANK","HDFC","ITC","ICICIBANK","IDFC","INDUSINDBK","JPASSOCIAT")
nifStks04_2 <- c("JINDALSTEL","KOTAKBANK","LT","LUPIN","MM","MARUTI","NTPC","ONGC","PNB","RANBAXY","RELIANCE")
nifStks04_3 <- c("SBIN","SUNPHARMA","TCS","TATAMOTORS","TATAPOWER","TATASTEEL","ULTRACEMCO","WIPRO")
nifStks09_1 <- c("AMBUJACEM","AXISBANK","BHARTIARTL","CAIRN","HINDALCO","HINDUNILVR","BAJAJ_AUTO","DLF","NMDC")
nifStks13_1 <- c("COALINDIA","HEROMOTOCO","SSLT","INFY")

#comps <- c("CNX_NIFTY","TATAELXSI", "TCS", "TATASTEEL","TATACOMM","TATAMOTORS","TATAMETALI","TATASPONGE","TATACHEM","TATAPOWER","TATACOFFEE")

comps <- c("CNX_NIFTY",nifStks04_1,nifStks04_2,nifStks04_3,nifStks09_1)

Close <- as.data.frame(sapply(comps, function(x) {
  NSEx <- paste("NSE", x, sep = "/")
  Quandl(NSEx, trim_start = Cstartdate, trim_end = Cenddate, collapse="monthly")$Close
}))

Close01 <- cbind(NIFTY01,Close)
write.csv(Close01,"NIFTY44Close.csv")
head(Close01)


IIR01 <- IIR00$Value

logreturn1 <- function(x) log(head(x,-1)/tail(x,-1))
##riskpremium <- function(x) logreturn1(x) - head(IIR01,-1)
riskpremium <- function(x) logreturn1(x) - IIR00$Value

Returns <- as.data.frame(apply(Close,2,riskpremium))

head(Returns)

r <-t(sapply(comps, function(comp)
  c(beta = lm(Returns[,comp] ~ riskpremium(NIFTY01$Close))$coefficients[[2]],
    mean = mean(Returns[,comp]))
  ))
r1 <- as.data.frame(r)
plot(r1$beta, r1$mean)
abline(lm(r1$mean ~ r1$beta), col = 'red')
========================================

Once the program is run, the following graph is generated

Calculating Beta for a specific stock

The Capital Asset Pricing Model, explained in Chapter 3 of our reference Resource 1 [ Introduction to R for Quantitative Finance ], leads us to the concept of Beta that is a measure of the risk of a particular security when compared to the risk of the market in general.

Beta is calculated on the basis of three input parameters, namely

• The returns on a given security, based on the logreturn formula, for example Tata Steel
• The returns of the market, which in our case will be returns on the CNX Nifty
• The risk free returns, represented in our case by bank interest in India.

Now we can run the the following program

========================================
setwd("C:/Users/admin/Desktop/Data Analytics/QuantitativeFinance/QFLabs")
getwd()
library(Quandl)
Quandl.auth("xxxxxxxxxxxxxxxxxxxx")

Cstartdate = '2007-06-01'
Cenddate = '2014-06-01'
##LIBOR <- Quandl('FED/RILSPDEPM01_N_B', start_date = Cstartdate, end_date = Cenddate)
IIR00 <- Quandl('BCB/17901', trim_start = Cstartdate, trim_end = Cenddate)
IIR00$Value <- IIR00$Value/1200
NIFTY <- Quandl("NSE/CNX_NIFTY", trim_start = Cstartdate, trim_end = Cenddate, collapse="monthly")
TISCO <- Quandl("NSE/TATASTEEL", trim_start = Cstartdate, trim_end = Cenddate, collapse="monthly")

cdates <- Reduce(intersect,list(IIR00$Date,NIFTY$Date,TISCO$Date))

IIR01 <- IIR00[IIR00$Date %in% cdates,'Value']
NIFTY01 <- NIFTY[NIFTY$Date %in% cdates,'Close']
TISCO01 <- TISCO[TISCO$Date %in% cdates,'Close']


## this is perhaps clear accurate
logreturn1 <- function(x) log(head(x,-1)/tail(x,-1))
riskpremium <- function(x) logreturn1(x) - head(IIR01,-1)

lrTIS1 <- logreturn1(TISCO01)

beta1_TISCO <- cov(logreturn1(TISCO01)-head(IIR01,-1),logreturn1(NIFTY01)-head(IIR01,-1))/var(logreturn1(NIFTY01)-head(IIR01,-1))
beta2_TISCO <- cov(riskpremium(TISCO01),riskpremium(NIFTY01))/var(riskpremium(NIFTY01))


(fit <- lm(riskpremium(TISCO01) ~ riskpremium(NIFTY01)))

plot(riskpremium(NIFTY01),riskpremium(TISCO01))
abline(fit,col = 'red')

(fit1 <- lm(riskpremium(TISCO01) ~ -1 +riskpremium(NIFTY01)))

summary(fit)
summary(fit1
========================================

Here what we are doing is as follows
We are getting the India Interest Rates, the NIFTY CNX closing prices and the Tata Steel Closing Prices for the same set of data. Note that for the 7 year data requested, we get 84 data points for the interest rates and 85 data points for the NSE data, which is good because in calculating returns from closing prices, we will be using data for two consecutive points and so we will end up with 84 data points

Also we note that unlike Google Data used in the book, the India data comes in descending order of dates [ later dates come earlier in the list ] and so the logreturn is calculated as
logreturn1 <- function(x) log(head(x,-1)/tail(x,-1) that is different from what is given in the book.

The program runs and we note that the beta calculated from ratio of covariances [ 1.7685] is the same as that obtained from the regression coeffient [1.768501] of the risk premium.

The Security Characteristic Line (SCL) is obtained as follows

---

---

The same analysis can be carried out using a spreadsheet and this is demonstrated in this Google Docs spreadsheet. The beta values calculated using the ratio of Covariance to Variance is very close to the value obtained using R (1.7691 vs 1.7685 ). However slope obtained in the spreadsheet (1.7905) is a little different.

A Risk Free Asset and the Tangency Portfolio

In many cases, we have a risk free asset ( like a Government Bond ) that offers zero risk and very small return and we would like to introduce the same into our portfolio.

Once again, we will see how we can create this portfolio that consists of a basket of risky assets ( the tangency portfolio) plus the risk free security. In the R program, shown below, we have taken the previous function ( created in an earlier post ) and enhanced the covariance matrix by forcibly adding a row and column of 0s and then further enhancing it with a row of returns that now includes the risk free return of the extra risk free asset.

========================================
   setwd("C:/Users/admin/Desktop/Data Analytics/QuantitativeFinance/QFLabs")
getwd()
India5 <- read.csv("India5.csv")
assets <- India5[,3:7]
histRets <- log(tail(assets, -1) / head(assets, -1))

##-----------------------------------------------------
## modifying the previous function

OptWeightsRF <- function(rets,mu,rf){
  ## we need one more column for riskless assets
  n <- ncol(rets)+1  
  ## we add a column of 0s to the right of covariance matrix
  Q <- cbind(cov(rets), rep(0, n - 1)) 
  ## we add a rows of 0s to the bottom covariance matrix
  Q <- rbind(Q, rep(0, n))
  ## we create a new row consisting of existing mean returns + risk free
  r <- c(colMeans(rets), rf)
  ## the cov matrix is enhanced by adding this extra rwo at bottom
  Q <- rbind(Q, rep(1, n), r)
  ## here the matrix is made into a square matrix
  Q <- cbind(Q, rbind(t(tail(Q, 2)), matrix(0, 2, 2)))
  b <- c(rep(0, n), 1, mu)
  w0 <- solve(Q, b)
  head(w0, -2)
}

OptWeightsRF(histRets,0.005,0.0001)
sum(OptWeightsRF(histRets,0.005,0.0001))

library("tseries")

p2 <- portfolio.optim(as.matrix(histRets),pm=0.005,riskless=TRUE,shorts=TRUE,rf=0.0001)
p2$pw
sum(p2$pw)
p2 <- portfolio.optim(as.matrix(histRets),riskless=TRUE,shorts=FALSE,rf=0.0001)
p2$pw
sum(p2$pw)
p2$pm

========================================

on running the program, the output from the OptWeightsRF function is as follows

> OptWeightsRF(histRets,0.005,0.0001)
      Rel      HDFC     Tisco       LnT       HUL        
 4.166645  4.136603  1.622548 -2.219220  3.104350 -9.810924
> sum(OptWeightsRF(histRets,0.005,0.0001))
[1] 1
where we see the weights for the 5 risky securities and that of the risk free security and we note that the sum adds up to 1. However this includes shorts.

When we run the same data through the portfolio.optim ( introduced in an earlier post)  function from the package tseries, we get the identical answer for the risky assets

> p2 <- portfolio.optim(as.matrix(histRets),pm=0.005,riskless=TRUE,shorts=TRUE,rf=0.0001)
> p2$pw
[1]  4.166645  4.136603  1.622548 -2.219220  3.104350
> sum(p2$pw)
[1] 10.81092

since the weight of the risk-free asset is not shown in the output but it can be calculated so that the sum adds up to 1. This means that the weight of the risk free asset should be -9.81 which is what was determined in the OptWeightsRF function.

finally, we run portfolio.optim so that SHORTs are not allowed and the answer is as follows

> p2 <- portfolio.optim(as.matrix(histRets),riskless=TRUE,shorts=FALSE,rf=0.0001)
> p2$pw
[1] 0.25225852 0.23839309 0.08740829 0.00000000 0.22403593
> sum(p2$pw)
[1] 0.8020958
> p2$pm
[1] 0.0003773413

here we note that all weights are positive including the riskfree asset = 1 - 0.802 = 0.198

Portfolio Optimisation - II

The analysis in the previous post created portfolios that had negative weights, or "short" positions. To overcome we search for options to block short positions and arrived at a couple of other websites that demonstrated portfolio optimization with more complex models. Both the Economist at Large site and the Alphaism site give good inputs on how to address this problem and this leads us to the function portfolio.optim available in the {tseries} that addresses this problem.

In the following program, this portfolio.optim function has been used instead of the custom built function shown in the previous post.

========================================
setwd("C:/Users/admin/Desktop/Data Analytics/QuantitativeFinance/QFLabs")
getwd()
India5 <- read.csv("India5.csv")
assets <- India5[,3:7]
returns <- log(tail(assets, -1) / head(assets, -1))

library("tseries")

## Create an optimum portfolio with expected means defined and shorts allowed

w2 <- portfolio.optim(as.matrix(returns),pm = 0.005,shorts = TRUE,riskless= FALSE)
w2$pw
sum(w2$pw)

## Create an optimum portfolio with expected means defined and shorts NOT ALLOWED

w2 <- portfolio.optim(as.matrix(returns),shorts = FALSE,riskless = FALSE)
w2$pw
sum(w2$pw)
w2$pm

# Create a function to create the frontier using portfolio.optim 

frontier2 <- function(return,minRet,maxRet){
  rbase <- seq(minRet,maxRet,length=100)
  s <- sapply(rbase,function(x){
    p2 <- portfolio.optim(as.matrix(returns),pm = x, shorts = TRUE)
    p2$ps^2
  })
  plot(s,rbase,xlab="Variance",ylab="Return", main="w/Portfolio.Optim")
}

frontier2(returns,-0.0005,0.05)
========================================

We observe that with Short positions allowed ( shorts=true), the portfolio weights are the same as what was obtained in the previous program ( previous post ), and as expected they sum up to 1

> w2$pw
[1]   4.869750   5.651428   3.375111 -10.711613  -2.184675
> sum(w2$pw)
[1] 1

however when shorts are not allowed ( shorts = false ) a different portfolio is created

> w2$pw
[1] 0.22797760 0.19765858 0.03686214 0.16357332 0.37392835
> sum(w2$pw)
[1] 1
> w2$pm
[1] 0.0003971317

and they sum up to 1 as well.

However, when short positions are not allowed, it is not possible to specify the portfolio returns before hand ( as it was, in the shorts allowed case, where pm = 0.005). With this constraint, the portfolio returns has dropped to 0.000397.

The efficient frontier remains the same as before

The program given in Alphaism uses portfolio.optim and then extends it further.

Portfolio Optimisation - I

Chapter 2 in the book Introduction to R for Quantitative Finance, gives a good example of how to extract data on data on returns on 5 US companies and use to to create an optimum portfolio and draw the efficient frontier. The actual working code for this is available here.

While the program works, one thing that is a little confusing is how the daily returns on the 5 stocks is calculated by taking the log of ratios of two subsequent daily closing prices! return <- log(tail(assets, -1) / head(assets, -1))

To understand how this is arrived at, see Zivot's notes on calculating returns or his slide deck.

We next tried to see how the same technique can be applied to 5 India companies that are listed in the NSE. The following R program shows how this data can be extracted from Quandl and converted into a CSV file for repeated usage.

========================================
setwd("C:/Users/admin/Desktop/Data Analytics/QuantitativeFinance/QFLabs")
getwd()
library(Quandl)
## Quandl Key for academics@yantrajaal.com
Quandl.auth("xxxxxxxxxxxxxxxx")
FRel1314 <- Quandl("NSE/RELIANCE", trim_start="2013-01-01", trim_end="2014-06-30")
FTisco1314 <- Quandl("NSE/TATASTEEL", trim_start="2013-01-01", trim_end="2014-06-30")
FLT1314 <- Quandl("NSE/LT", trim_start="2013-01-01", trim_end="2014-06-30")
FHDFC1314 <- Quandl("GOOG/BOM_500180", trim_start="2013-01-01", trim_end="2014-06-30")
FHUL1314 <- Quandl("GOOG/NSE_HINDUNILVR", trim_start="2013-01-01", trim_end="2014-06-30")
library(plyr)
CREL <- FRel1314[,c('Date','Close')]
CTIS <- FTisco1314[,c('Date','Close')]
CLNT <- FLT1314[,c('Date','Close')]
CHUL <- FHUL1314[,c('Date','Close')]
CHDF <- FHDFC1314[,c('Date','Close')]
Stocks <- merge(CREL,CHDF,by.x='Date',by.y='Date')
Stocks <- rename(Stocks,c("Close.x"="Rel","Close.y"="HDFC"))
Stocks <- merge(Stocks,CTIS,by.x='Date',by.y='Date')
Stocks <- rename(Stocks,c("Close"="Tisco"))
Stocks <- merge(Stocks,CLNT,by.x='Date',by.y='Date')
Stocks <- rename(Stocks,c("Close"="LnT"))
Stocks <- merge(Stocks,CHUL,by.x='Date',by.y='Date')
Stocks <- rename(Stocks,c("Close"="HUL"))
write.csv(Stocks,"India5.csv")
========================================

Please note that  you need to get your own Quandl authorisation code by registering ( free of course) at Quandl. Also note that the data comes from two different sources so in each case you need to chose the fields that are required, in this case, Date and Closing Prices. However not all the data has closing prices for all days and so we need to merge the data (something like an SQL join) to make sure that your data frame has closing prices for all 5 companies for the same dates and they are in the same row. Once the data is ready, they are stored in a CSV file that will be accessed in the next program.

========================================
setwd("C:/Users/admin/Desktop/Data Analytics/QuantitativeFinance/QFLabs")
getwd()
## read the data
India5 <- read.csv("India5.csv")

## select the data for the 5 companies in columns 3 to 7
assets <- India5[,3:7]

## calulate the returns
returns <- log(tail(assets, -1) / head(assets, -1))

## Define a function to calculate the optimum portfolio weights
OptWeights <- function(return, mu = 0.005) {
  Q <- rbind(cov(return), rep(1, ncol(assets)), colMeans(return))
  Q <- cbind(Q, rbind(t(tail(Q, 2)), matrix(0, 2, 2)))
  b <- c(rep(0, ncol(assets)), 1, mu)
  head(solve(Q, b),-2)
} 

## Call the function with data and note that the weights add up to 1
OptWeights(returns)
sum(OptWeights(returns))

## Define a function to create the graph of the Efficient frontier
frontier <- function(return,minRet,maxRet){
  Q <- cov(return)
  n <- ncol(assets)
  r <- colMeans(return)
  Q1 <- rbind(Q,rep(1,n),r)
  Q1 <- cbind(Q1,rbind(t(tail(Q1,2)),matrix(0,2,2)))
  rbase <- seq(minRet,maxRet,length=100)
  s <- sapply(rbase,function(x){
    y <- head(solve(Q1,c(rep(0,n),1,x)),n)
    y %*% Q %*% y
  })
  plot(s,rbase,xlab="Variance",ylab="Return", main = "Custom")
}

## Call the function with the data and plot the frontier between two values
## of returns

frontier(returns,-0.0005,0.05)
========================================

The optimum weights are calculated as follows and they do add up to 1 as expected

> OptWeights(returns)
       Rel       HDFC      Tisco        LnT        HUL
  4.869750   5.651428   3.375111 -10.711613  -2.184675
> sum(OptWeights(returns))
[1] 1

However we note that some of the weights are negative, which means that the portfolio allows "short" positions. It is possible to bar short positions. This will demonstrated in the next post.