# A small set of ETF shares were selected for this experiment. To allow a rolling window
# I wanted as long a history as possible. ETFs are relatively new instruments, so there
# where a limited number that had an approximately ten year history.
# 
# LQD     iShares investment grade corporate bond
# SHY     iShares Barclays 1-3 Year Treasury Bond
# 
# IWF     Russell 1000 Large Cap Growth Index Fund
# IWB     Russell 1000 Large Cap Market Index
# IWD     Russell 1000 Large Cap Value Index Fund
# 
# IWO     Russell 2000 Small Cap Growth Index
# IWM     Russell 2000 Small Cap Index
# IWN     Russell 2000 Small Cap Value Index
# 
# IWP     Russelll Midcap Growth Index Fund
# IWR     Russelll Midcap  Index Fund
# IWS     Russelll Midcap Value Index Fund
# 
# IYM     Dow Jones US Basic Materials Sector
# IYK     Dow Jones Consumer Goods
# IYE     Dow Nones US energy Sector
#
# Proxy for the market:
#
# SPY     State Street SPDR S&P 500
#
# Author: Ian L. Kaplan
#         September 2012
#
# This code is in the public domain. This code is research code and is not intended to be used
# for actual investment. This code is published without warranty and the user accepts all risk
# in using this code.
#

library("PerformanceAnalytics")
library("tseries")
library("colorspace")

#
# Build a matrix of market quotes. The fuction returns the quotes as a zoo 
# matrix.
# Arguments:
#    symbols- a vector of one or more symbols (character strings)
#    period - d (daily), w (weekly) or m (monthly)
#    startDate - start date in ISO date format (e.g., 2011-12-31)
#    endDate - end date in ISO date format
#
buildQuoteMatrix = function( symbols, period, startDate, endDate)
{
  for (i in 1:length(symbols))
  {
    sym = symbols[i]
    print(sprintf("Reading symbol %s from yahoo", sym))
    ts = get.hist.quote( sym, startDate, endDate, provider="yahoo", compression=period, quote="AdjClose", quiet=T)
    if (i == 1) {
      quote.mat = ts
    } else {
      quote.mat = cbind(quote.mat, ts)
    }
  }
  colnames(quote.mat) = symbols
  return( quote.mat)
} # buildQuoteMatrix

#
# Calculate the weekly rate for an annual interest time series
#   symbol - the symbol for the interest rate instrument (e.g., ^TNX for the ten year t-note)
#   startDate, endDate - the start and end dates in ISO format
# 
#  Note that interest rates are reported as annual rates. Since we're dealing with weekly data,
#  the annual rate needs to be converted to the weekly rate.
# 
weeklyRate = function( symbol, startDate, endDate)
{
  # get the annual rate, as a percentage
  rf = get.hist.quote( symbol, startDate, endDate, provider="yahoo", compression="w", quote="AdjClose", quiet=T)
  rf.frac = rf / 100
  rf.week = ((1 + rf.frac)^(1/52)) - 1
  return(rf.week)
}

#
# Calculate the long-only tangency portfolio for a returns matrix.  The tangency portfolio is
# the portfolio with the highest sharp ratio.
#
longOnlyTangency = function(ret.m, rf)
{
  mu = t(apply(ret.m, 2, mean))
  rng = range(mu)
  minRet = max(rng[1], rf)
  maxRet = rng[2]
  nret = 100
  retSeq = seq(from=rng[1], to=rng[2], length=nret)
  sd.vec = c()
  targetRet = c()
  weights.l = list()
  retCov = var(ret.m)
  ix = 1
  minWeight = 0.0
  for (i in 1:length(retSeq)) {
    op = try(
         portfolio.optim(x=mu, rf=rf, pm=retSeq[i], shorts=F, covmat=retCov), silent=T
    )
    if (class(op) != "try-error") {
      targetRet = c(targetRet, retSeq[i])
      w = op$pw
      weights.l[[ix]] = round(w,4)
      sd = sqrt( w %*% retCov %*% w)
      sd.vec = c(sd.vec, sd)
      ix = ix + 1
    }
  }
  weights = matrix(unlist(weights.l), nrow=length(weights.l), ncol=ncol(ret.m), byrow=T)
  colnames(weights) = colnames(ret.m)
  sharp = (targetRet - rf) / sd.vec
  ix = which.max(sharp)
  tanWeights = weights[ix,]
  portMu = as.numeric(mu) %*% tanWeights
  portSd = sd.vec[ix]
  return(list(w = tanWeights, sd = portSd, mu = portMu))
} # longOnlyTangency

#  Calculate a long-only, efficient portfolio (using mean variance) and calculate the return
#  and standard deviation of the forward, out of sample portfolio.
# 
#  ret.mat : a matrix of returns
#  rf.vec : a vector of risk free rates (where length = nrow(ref.mat))
#  width : the width of the portfolio window
#  step : the number of steps to take into the future
# 
#   Return
#     zoo vectors with a date index
#       predMu  : the predicted portfolio return
#       actMu   : the actual portfolio return
#       predSd  : the predicted standard deviation
#       actSd   : the actual standard deviation
#       eqWt    : the mean return of an equally weighted portfolio
#     
#       weights : a matrix (with date row names) of portfolio weights
# 
portfolioInOutSample = function( ret.mat, rf.vec, width, step)
{
  numRows = nrow(ret.mat) 
  start = 1
  end = width
  ix = 1
  
  predMu = list()
  predSd = list()
  actMu = list()
  actSd = list()
  eqWeight = list()
  weights = list()
  
  n = ncol(ret.mat)
  eqW = rep(1/n, n)

  while((end + step) < numRows) {
    window = ret.mat[start:end,]
    rf = mean(rf.vec[start:end])
    t = longOnlyTangency(window, rf)
    len = nrow(window)
    y = as.Date(names(window[,1]))[len]
    predMu[[ix]] = t$mu
    predSd[[ix]] = t$sd
    names( predMu[[ix]] ) = y
    names( predSd[[ix]]) = y
    w = t$w
    weights[[ix]] = w
    outStart = end + 1
    outEnd = outStart + step
    outOfSample = ret.mat[outStart:outEnd,]
    actRet = as.numeric(outOfSample %*% w)
    actMu[[ix]] = mean( actRet )
    actSd[[ix]] = sd( actRet )
    eqWeight[[ix]] = mean(as.numeric(outOfSample %*% eqW))
    len = nrow(outOfSample)
    y = as.Date(names(outOfSample[,1]))[len]
    names( actMu[[ix]] ) = y
    names( actSd[[ix]]) = y
    names( eqWeight[[ix]]) = y
    start = start + step
    end = end + step
    ix = ix + 1
  }
  
  predMu = unlist(predMu)
  predSd = unlist(predSd)
  actMu = unlist(actMu)
  actSd = unlist(actSd)
  eqWeight = unlist(eqWeight)
  
  predMu.z = as.zoo(predMu)
  index(predMu.z) = as.Date(names(predMu))
  predSd.z = as.zoo(predSd)
  index(predSd.z) = as.Date(names(predSd))
  actMu.z = as.zoo(actMu)
  index(actMu.z) = as.Date(names(actMu))
  actSd.z = as.zoo(actSd)
  index(actSd.z) = as.Date(names(actSd))
  eqWeight.z = as.zoo(eqWeight)
  index(eqWeight.z) = as.Date(names(eqWeight))
  
  numRows = length(weights)
  numCols = ncol(ret.mat)
  weights =matrix(data = unlist(weights), ncol=numCols, nrow=numRows, byrow=T)
  colnames(weights) = colnames(ret.mat)
  rownames(weights) = names(predMu)
  
  return(list(predMu = predMu.z, actMu = actMu.z, 
              predSd = predSd.z, actSd = actSd.z, 
              eqWt = eqWeight.z, 
              weights = weights))
}  # portfolioInOutSample

# 
#  - Plot the mean for the predicted, actual and equally weighted portfolio returns, over time
#  - Plot the predicted and actual standard deviations
#  - Plot a bar chart showning the portfolio weights over time
# 
#  This function takes a list containing the following zoo objects:
#    predMu : the predicted mean return for the portfolio
#    actMu  : the actual mean return for the portfolio
#    predSd : the predicted standard deviation for the portfolio
#    actSd  : the actual standard deviation
#    eqWt   : the mean return for the equally weighted portfolio
#    
#    weights : this is not a zoo object, but a matrix of portfolio weights over time
# 
plotPortfolios = function( portData, title )
{
  predMu.z = portData$predMu
  actMu.z = portData$actMu
  predSd.z = portData$predSd
  actSd.z = portData$actSd
  eqWeight.z = portData$eqWt
  weights = portData$weights
  zero = as.zoo(rep(0, length(actMu.z)))
  index(zero) = index(actMu.z)
  
  # Annualize the mean and standard deviation
  predMu.z = ((1 + predMu.z) ^ 52) - 1
  actMu.z = ((1 + actMu.z) ^ 52) -1
  eqWeight.z = ((1 + eqWeight.z) ^ 52) - 1
  predSd.z = predSd.z * sqrt(52)
  actSd.z = actSd.z * sqrt(52)
  
  minDate = index(predMu.z[1])
  maxDate = index(actMu.z[length(actMu.z)])
  
  par(mfrow=c(2,1), mar=c(2,4.5,1.5,1))
  ymin = min(predMu.z, actMu.z, eqWeight.z)
  ymin = ymin - abs(ymin * 0.2)
  ymax = max(predMu.z, actMu.z, eqWeight.z)
  ymax = ymax + (ymax * 0.5)
  plot(predMu.z, type="l", col="blue", main=title, ylab="Annualized Mean Return", 
       ylim=c(ymin, ymax), xlim=c(minDate, maxDate))
  lines(zero, type="l", col="black")
  points(predMu.z, pch=15, col="blue")
  lines(actMu.z, type="l", col="green")
  points(actMu.z, pch=17, col="green")
  lines(eqWeight.z, type="l", col="black")
  legend(x="bottomleft", legend=c("predicted", "actual", "eq. wt"), col=c("blue", "green", "black"), 
         lwd=4, cex=0.8, horiz=T)
  
  minDate = index(predSd.z[1])
  maxDate = index(actSd.z[length(actMu.z)])
  ymax = max(predSd.z, actSd.z)
  ymax = ymax + (ymax * 0.2)
  plot(predSd.z, type="l", col="blue", cex.lab=1.3, ylab=expression(paste("Annualized ", sigma)), 
       xlim=c(minDate, maxDate), ylim=c(0, ymax))
  points(predSd.z, pch=15, col="blue")
  lines(actSd.z, type="l", col="green")
  points(actSd.z, pch=17, col="green")
  legend(x="topleft", legend=c("predicted", "actual"), col=c("blue", "green"), lwd=4, cex=0.8, horiz=F)
  par(mfrow=c(1,1))
  
  par(mfrow=c(1,1))
  chart.StackedBar(w = weights, 
                   #                  colorset = rainbow(ncol(weights))[c(seq(from=1, to=ncol(weights), by=2), 
                   #                                                      seq(from=2, to=ncol(weights), by=2))], 
                   colorset = rainbow(ncol(weights)),
                   ylab="Weight",
                   space=0.1, main=title)
  par(mfrow=c(1,1))
}  # plotPortfolios

# 
#   Generate cumulative return plots for the portfolio predicted return, the actual return,
#   the return on a bond ETF and the market return. The portfolio returns are returns over
#   the number of weeks defined in the weeks argument. The bond and market returns are 
#   weekly returns, which must be aggregated to match the portfolio return
# 
plotPortReturn = function(predMu, actMu, bondMu, mktRet, offset, weeks )
{
  pMu = coredata(predMu)
  aMu = coredata(actMu)
  if (class(mktRet) == "zoo") {
    mktRet = coredata(mktRet)
  }
  # For the bond fund and the market, calculate the mean return over teh weeks window
  bMu = round(rollapply(bondMu[(offset+1):length(bondMu)], width=weeks, by=weeks, FUN=mean, align="right"), 6)
  mMu = round(rollapply(mktRet[(offset+1):length(mktRet)], width=weeks, by=weeks, FUN=mean, align="right"), 6)
  
  # The returns are the mean return over the weeks period. However, for cumulative returns to work we
  # need to calculate the compound return over that period.
  adjPredMu = ((1 + pMu)^weeks) - 1
  adjActMu = ((1 + aMu)^weeks) - 1
  adjBndMu = ((1 + bMu)^weeks) - 1
  adjMktMu = ((1 + mMu)^weeks) - 1

  r = cbind(adjPredMu, adjActMu, adjBndMu, adjMktMu)
  colnames(r) = c("Predicted", "Actual", "Bonds", "S&P 500")
  par(mfrow=c(1,1))
  chart.CumReturns(r, wealth.index=T, main=paste("Cumulative Mean Returns (", weeks, " week period)", sep=""),
                   col=c("blue", "green", "red", "black"), 
                   legend.loc="topleft", grid.color="slategrey")
  par(mfrow=c(1,1))
} # plotPortReturn


#
# This function attempts to simulate the investment of the MV portfolio by buying at the
# start of a period and selling at the end.
#
# The portfolio is a long-only tangency portfolio. The investment of this portfolio is
# simulated by calculating the out-of-sample return over a perios of weeks. The sample used 
# to calculate is a sliding window of width weeks.
#
# returns : a matrix of log returns
# prices : a matrix of weekly prices
# rf : a vector of the weekly risk free rate
# width : the width of the window used to estimate the portfolio
# weeks : the investment period
# W : the initial "wealth"
#
# The returns used here are defined as
#
#     R = (P[t] - P[t-1])/P[t-1]
#
# The portfolio is calculated using weekly log returns, so the mean return from the
# portfolio calculation is the weekly mean return.
#
# The out-of-sample return is calculated as if the portfolio were invested at the start
# of the period and sold at the end of the period.
#
# To allow comparison of the out-of-sample return to the portfolio return, the portfolio return
# is converted to the return over the period of "weeks".
#
investmentSimulation = function( returns, prices, rf.vec, width, weeks )
{
  numRows = nrow(returns) 
  start = 1
  end = width
  ix = 1
  
  predMu = list()
  portRet = list()
  eqWtRet = list()
  
  n = ncol(returns)
  eqW = rep(1/n, n)
  
  while((end + weeks) < numRows) {
    window = returns[start:end,]
    rf = mean(rf.vec[start:end])
    t = longOnlyTangency(window, rf)
    len = nrow(window)
    y = as.character(rownames(window)[len])
    weeklyMu = t$mu
    periodMu = ((1 + weeklyMu)^weeks) - 1
    predMu[[ix]] = periodMu
    names( predMu[[ix]] ) = y
    w = t$w
    outStart = end + 1
    outEnd = outStart + weeks
    # we assume that we buy at the start of the out of sample period (which is just
    # after the window used to calculate the portfolio) and sell at the end of the
    # out of sample period
    R = (prices[outEnd,] - prices[outStart,])/prices[outStart,]
    portRet[[ix]] = sum( R * w)
    eqWtRet[[ix]] = sum(R * eqW)
    y = rownames(prices)[outEnd]
    names( portRet[[ix]] ) = y
    names( eqWtRet[[ix]]) = y
    start = start + weeks
    end = end + weeks
    ix = ix + 1
  }
  
  predMu = unlist(predMu)
  portRet = unlist(portRet)
  eqWtRet = unlist(eqWtRet)
  
  predMu.z = as.zoo(predMu)
  index(predMu.z) = as.Date(names(predMu))
  portRet.z = as.zoo(portRet)
  index(portRet.z) = as.Date(names(portRet))
  eqWtRet.z = as.zoo(eqWtRet)
  index(eqWtRet.z) = as.Date(names(eqWtRet))
  
  return(list(predMu = predMu.z, 
              portRet = portRet.z, 
              eqWtRet = eqWtRet.z))  
}  # investmentSimulation

# The returns used here are defined as
#
#     R = (P[t] - P[t-1])/P[t-1]
#
regionReturn = function( region )
{
  len = length(region)
  R = (region[len] - region[1])/region[1]
  return( R )
} # regionReturn


#
# Plot the simulated investment portfolios, along with the bond and market returns.
# Here the return is the return that would be obtained by buying at the start of
# a period and selling at the end of the period.
#
plotInvest = function( predMu, actRet, eqWtRet, bond, market, window, weeks)
{
  bondRet = rollapply(bond[(window+1):length(bond)], width=weeks, by=weeks, FUN=regionReturn, align="right")
  mktRet = rollapply(market[(window+1):length(market)], width=weeks, by=weeks, FUN=regionReturn, align="right")
  retData = cbind(coredata(predMu), coredata(actRet), coredata(eqWtRet), coredata(bondRet), coredata(mktRet))
  rownames(retData) = as.character(index(actRet))
  colnames(retData) = c("MV Predicted", "Actual", "eqWtRet", "bond", "market")
  par(mfrow=c(1,1))
  chart.CumReturns(retData, wealth.index=T, main=paste("Cumulative Period Returns (", weeks, " week period)", sep=""),
                   col=c("blue", "green", "magenta", "red", "black"), 
                   legend.loc="topleft", grid.color="slategrey")
  par(mfrow=c(1,1))
}

plotDrawdown = function( predMu, actRet, eqWtRet, bond, market, window, weeks)
{
  bondRet = rollapply(bond[(window+1):length(bond)], width=weeks, by=weeks, FUN=regionReturn, align="right")
  mktRet = rollapply(market[(window+1):length(market)], width=weeks, by=weeks, FUN=regionReturn, align="right")
  retData = cbind(coredata(predMu), coredata(actRet), coredata(eqWtRet), coredata(bondRet), coredata(mktRet))
  rownames(retData) = as.character(index(actRet))
  colnames(retData) = c("MV Predicted", "Actual", "eqWtRet", "bond", "market")
  par(mfrow=c(1,1))
  chart.Drawdown(retData, main=paste("Drawdown, (rebalance period = ", weeks, " weeks)", sep=""), cex.main=0.8,
                   colorset=c("blue", "green", "magenta", "red", "black"), 
                   legend.loc="bottomleft" )
  par(mfrow=c(1,1))
}


startDate = "2002-08-01"
endDate = "2012-08-01"

# The universe of Exchange Traded Funds (ETFs)
symbols = c("LQD", "SHY","IWF","IWB", "IWD","IWO","IWM","IWN", "IWP", "IWR","IWS", "IYM","IYK", "IYE")
# The ten year T-note annual rate
rfSym = "^TNX"
# Market ETF
mktSym = "SPY"

returnFile = "etfPrices.R"

# Load the ETF returns object, if the save file exists or write out the save file
if (file.exists( returnFile ) == T)
{
  load( returnFile )
} else {
  prices.mat = buildQuoteMatrix(symbols, period="w", startDate, endDate)
  rfTenYear = weeklyRate(rfSym, startDate, endDate)
  mkt = get.hist.quote( mktSym, startDate, endDate, provider="yahoo", compression="w", quote="AdjClose", quiet=T)
  save(prices.mat, mkt, rfTenYear, file=returnFile)
}
logPrices.mat = apply(prices.mat, 2, log)
ret.mat = apply(logPrices.mat, 2, diff)
mktRet = diff(log(mkt))

yearPort = portfolioInOutSample(ret.mat, rfTenYear, 52, 52)
plotPortfolios( yearPort, "52-week period" )
quarterPort = portfolioInOutSample(ret.mat, rfTenYear, 52, 12)
plotPortfolios( quarterPort, "12-week period" )

# Get the weekly bond fund returns, for comparision
lqd = ret.mat[,"LQD"]

# plotPortReturn(yearPort$predMu, yearPort$actMu, lqd, mktRet, 52, 52)
# plotPortReturn(quarterPort$predMu, quarterPort$actMu, lqd, mktRet, 52, 12)

prices = coredata(prices.mat)
rownames(prices) = as.character(index(prices.mat))

yearInvest = investmentSimulation( ret.mat, prices, rfTenYear, 52, 52)
quarterInvest = investmentSimulation( ret.mat, prices, rfTenYear, 52, 12)
# weekly prices for the bond fund LQD
lqdPrices = prices.mat[,"LQD"]
plotInvest(yearInvest$predMu, yearInvest$portRet, yearInvest$eqWtRet, lqdPrices, mkt, 52, 52)
plotInvest( quarterInvest$predMu, quarterInvest$portRet, quarterInvest$eqWtRet, lqdPrices, mkt, 52, 12)

#
# Generate the drawdown plot for the quarterly rebalanced portfolio
#
plotDrawdown( quarterInvest$predMu, quarterInvest$portRet, quarterInvest$eqWtRet, lqdPrices, mkt, 52, 12)

#
# Generate the efficient portfolio graph (quick and dirty, since it uses the fPortfolioPackage)
#
cov = function( mat )
{
  return(var(mat))
}

Constraints = "LongOnly"
Spec = portfolioSpec()
setNFrontierPoints(Spec) = 80
setRiskFreeRate(Spec) = mean( rfTenYear)
ret.mat.ts = as.timeSeries(ret.mat)
frontier = portfolioFrontier(data=ret.mat.ts, spec=Spec, 
                             constraints=Constraints, description="Long-Only Efficient Frontier")
par(mfrow=c(1,1), mar=c(2,4.5,1.5,1))
plot(frontier)
par(mfrow=c(1,1))