#======================================================
# Use systematic investors toolbox (SIT)
# Load Systematic Investor Toolbox (SIT)
# http://www.r-bloggers.com/backtesting-minimum-variance-portfolios/
# https://systematicinvestor.wordpress.com/2011/12/13/backtesting-minimum-variance-portfolios/
# https://systematicinvestor.wordpress.com/2013/03/22/maximum-sharpe-portfolio/
#=========================================================
#setInternet2(TRUE)
# Load Systematic Investor Toolbox (SIT)
rm(list=ls())
#setInternet2(TRUE)
con = gzcon(url('https://github.com/systematicinvestor/SIT/raw/master/sit.gz', 'rb'))
source(con)
close(con)

#*****************************************************************
# Load historical data
#****************************************************************** 
load.packages('quantmod,quadprog,lpSolve')
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Attaching package: 'TTR'
## The following object is masked _by_ '.GlobalEnv':
## 
##     DVI
## Version 0.4-0 included new data defaults. See ?getSymbols.
tickers = spl('EA,LMT,TSLA')


data <- new.env()
getSymbols(tickers, src = 'yahoo', from = '1980-01-01', env = data, auto.assign = T)
## '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.
## [1] "EA"   "LMT"  "TSLA"
for(i in ls(data)) data[[i]] = adjustOHLC(data[[i]], use.Adjusted=T)

names(data)
## [1] "TSLA"        "LMT"         "EA"          ".getSymbols"
data.weekly <- new.env()
for(i in tickers) data.weekly[[i]] = to.weekly(data[[i]], indexAt='endof')
data.monthly <- new.env()
for(i in tickers) data.monthly[[i]] = to.monthly(data[[i]], indexAt='endof')

bt.prep(data, align='remove.na', dates='2010::2018')
bt.prep(data.monthly, align='remove.na', dates='2010::2018')
names(data)
## [1] "prices"          "TSLA"            "LMT"             "dates"          
## [5] "EA"              "weight"          ".getSymbols"     "symbolnames"    
## [9] "execution.price"
#*****************************************************************
# Code Strategies
#****************************************************************** 
prices = data$prices   
n = ncol(prices)
n
## [1] 3
# find week ends
week.ends = endpoints(prices, 'weeks')
week.ends = week.ends[week.ends > 0]  
month.ends = endpoints(prices, 'months')
month.ends = month.ends[month.ends > 0]  


# Equal Weight 1/N Benchmark
data$weight[] = NA
data$weight[week.ends,] = ntop(prices[week.ends,], n)       
data$weight[month.ends,] = ntop(prices[month.ends,], n) 
#capital = 100000
#data$weight[] = (capital / prices) * data$weight

equal.weight = bt.run(data)
## Latest weights :
##                        EA   LMT  TSLA
## 2018-12-31 08:00:00 33.33 33.33 33.33
## 
## Performance summary :
##  CAGR    Best    Worst   
##  30.7    8.1 -9  
names(equal.weight)
## [1] "weight"      "type"        "ret"         "best"        "worst"      
## [6] "equity"      "cagr"        "dates.index"
head(equal.weight$equity)
##                            EA
## 2010-06-29 08:00:00 1.0000000
## 2010-06-30 08:00:00 1.0000000
## 2010-07-01 08:00:00 0.9720431
## 2010-07-02 08:00:00 0.9373379
## 2010-07-06 08:00:00 0.8829540
## 2010-07-07 08:00:00 0.8921729
#*****************************************************************
# Create Constraints
#*****************************************************************
constraints = new.constraints(n, lb = -Inf, ub = +Inf)

# SUM x.i = 1
constraints = add.constraints(rep(1, n), 1, type = '=', constraints)        


ret = prices / mlag(prices) - 1
weight = coredata(prices)
weight[] = NA

for( i in month.ends[month.ends >= (756 + 1)] ) {
  # one quarter is 63 days
  hist = ret[ (i- 756 +1):i, ]
  
  # create historical input assumptions
  ia = create.historical.ia(hist, 756)
  #s0 = apply(coredata(hist),2,sd)     
  #ia$cov = cor(coredata(hist), use='complete.obs',method='pearson') * (s0 %*% t(s0))
  
  weight[i,] = min.risk.portfolio(ia, constraints)
}

tail(weight,10)
##                EA       LMT       TSLA
## [2133,]        NA        NA         NA
## [2134,]        NA        NA         NA
## [2135,]        NA        NA         NA
## [2136,]        NA        NA         NA
## [2137,]        NA        NA         NA
## [2138,]        NA        NA         NA
## [2139,]        NA        NA         NA
## [2140,]        NA        NA         NA
## [2141,]        NA        NA         NA
## [2142,] 0.1918519 0.7383225 0.06982564
# Minimum Variance
data$weight[] = weight      
#capital = 100000
#data$weight[] = (capital / prices) * data$weight
min.var.daily = bt.run(data)
## Latest weights :
##                        EA   LMT TSLA
## 2018-12-31 08:00:00 17.15 75.51 7.34
## 
## Performance summary :
##  CAGR    Best    Worst   
##  12.4    6.1 -5.4    
names(min.var.daily)
## [1] "weight"      "type"        "ret"         "best"        "worst"      
## [6] "equity"      "cagr"        "dates.index"
#min.var.daily$equity
#min.var.daily$ret
#min.var.daily$weight
#
#*****************************************************************
# Code Strategies: Weekly
#******************************************************************     
retw = data.weekly$prices / mlag(data.weekly$prices) - 1
weightw = coredata(prices)

weightw[] = NA

week.ends<-week.ends[-length(week.ends)]
for( i in month.ends[month.ends >= (756 + 1)] ) {   
  # map
  j = which(index(ret[i,]) == index(retw))
  
  # one quarter = 13 weeks
  #hist = retw[ (j- 156 +1):j, ]
  
  # create historical input assumptions
  ia = create.historical.ia(hist, n)
  s0 = apply(coredata(hist),2,sd)     
  ia$cov = cor(coredata(hist), use='complete.obs',method='pearson') * (s0 %*% t(s0))
  
  weightw[i,] = min.risk.portfolio(ia, constraints)
}   

data$weight[] = weightw     
capital = 100000
data$weight[] = (capital / prices) * data$weight
min.var.weekly = bt.run(data, type='share', capital=capital)
## Latest weights :
##                        EA   LMT TSLA
## 2018-12-31 08:00:00 20.34 72.18 7.48
## 
## Performance summary :
##  CAGR    Best    Worst   
##  13.1    5.1 -5.2    
#*****************************************************************
# Code Strategies: Monthly
#******************************************************************    
retm = data.monthly$prices / mlag(data.monthly$prices) - 1
weightm = coredata(prices)
weightm[] = NA

month.ends<-month.ends[-length(month.ends)]
for( i in month.ends[month.ends >= (756 + 1)] ) {   
  # map
  j = which(index(ret[i,]) == index(retm))
  
  # one quarter = 13 weeks
  hist = retm[ (j- 36 +1):j, ]
  
  # create historical input assumptions
  ia = create.historical.ia(hist, 36)
  s0 = apply(coredata(hist),2,sd)     
  ia$cov = cor(coredata(hist), use='complete.obs',method='pearson') * (s0 %*% t(s0))
  
  weightm[i,] = min.risk.portfolio(ia, constraints)
}   

data$weight[] = weightm     
capital = 100000
data$weight[] = (capital / prices) * data$weight
min.var.monthly = bt.run(data, type='share', capital=capital)
## Latest weights :
##                       EA   LMT  TSLA
## 2018-12-31 08:00:00 3.47 72.78 23.76
## 
## Performance summary :
##  CAGR    Best    Worst   
##  13.2    5.7 -5.1    
#*****************************************************************
# Create Report
#****************************************************************** 
plotbt.custom.report.part3(min.var.daily, min.var.weekly, min.var.monthly, equal.weight)
#
models<-list("Min.var.daily" = min.var.daily, 
             "Min.var.weekly" = min.var.weekly,
             "Min.var.monthly" = min.var.monthly,
             "Equal.weight" = equal.weight)
#
strategy.performance.snapshoot(models, T)

## NULL
strategy.performance.snapshoot(models, control=list(comparison=T), 
                               sort.performance=T)

plotbt.strategy.sidebyside(models, return.table=T)
##            Min.var.daily       Min.var.weekly      Min.var.monthly    
## Period     "Jun2010 - Dec2018" "Jun2010 - Dec2018" "Jun2010 - Dec2018"
## Cagr       "12.43"             "13.14"             "13.23"            
## Sharpe     "0.95"              "1"                 "0.97"             
## DVR        "0.84"              "0.89"              "0.86"             
## Volatility "13.27"             "13.25"             "13.75"            
## MaxDD      "-30.59"            "-30.05"            "-28.22"           
## AvgDD      "-2.46"             "-2.34"             "-2.57"            
## VaR        "-1.33"             "-1.33"             "-1.41"            
## CVaR       "-2.13"             "-2.11"             "-2.18"            
## Exposure   "63.68"             "63.68"             "63.68"            
##            Equal.weight       
## Period     "Jun2010 - Dec2018"
## Cagr       "30.71"            
## Sharpe     "1.2"              
## DVR        "1.13"             
## Volatility "24.93"            
## MaxDD      "-28.83"           
## AvgDD      "-3.6"             
## VaR        "-2.42"            
## CVaR       "-3.39"            
## Exposure   "99.91"
plotbt.strategy.sidebyside(min.var.daily, return.table=T)
##            min.var.daily      
## Period     "Jun2010 - Dec2018"
## Cagr       "12.43"            
## Sharpe     "0.95"             
## DVR        "0.84"             
## Volatility "13.27"            
## MaxDD      "-30.59"           
## AvgDD      "-2.46"            
## VaR        "-1.33"            
## CVaR       "-2.13"            
## Exposure   "63.68"
# plot Daily and Weekly transition maps
layout(1:3)

plotbt.transition.map(min.var.daily$weight)
legend('topright', legend = 'min.var.daily', bty = 'n')
plotbt.transition.map(min.var.weekly$weight)
legend('topright', legend = 'min.var.weekly', bty = 'n')
plotbt.transition.map(min.var.monthly$weight)
legend('topright', legend = 'min.var.monthy', bty = 'n')