#*****************************************************************
# First, reproduce S&P 500 chart
#*****************************************************************
library(SIT)
load.packages('quantmod')
tickers = 'SP=^GSPC'

data <- new.env()
getSymbols.extra(tickers, src = 'yahoo', from = '1970-01-01', env = data, auto.assign = T)
for(i in ls(data)) data[[i]] = adjustOHLC(data[[i]], use.Adjusted=T)

# compute Donchian channels
d.high.channel = runMax(Hi(data$SP), 55)
d.low.channel = runMin(Lo(data$SP), 55)

# compute Percentile channels
p.high.channel = runQuantile(Cl(data$SP), 55, probs=0.75)
p.low.channel = runQuantile(Cl(data$SP), 55, probs=0.25)

# compute Average
sma = SMA(Cl(data$SP), 55)

# make plot
plota(data$SP['2013:10::2014',],type='l', col='lightblue', lwd=2)
  plota.legend('SP')
plota.lines(sma, col='green', lwd=2)
plota.lines(d.high.channel, col='red', lwd=2)
plota.lines(d.low.channel, col='gray', lwd=2)
plota.lines(p.high.channel, col='yellow', lwd=2)
plota.lines(p.low.channel, col='blue', lwd=2)

#*****************************************************************
# Load historical data
#*****************************************************************
library(SIT)
load.packages('quantmod')
tickers = 'DBC+CRB'

# load saved Proxies Raw Data, data.proxy.raw
# please see http://systematicinvestor.github.io/Data-Proxy/ for more details
load('data/data.proxy.raw.Rdata')

data <- new.env()
getSymbols.extra(tickers, src = 'yahoo', from = '1970-01-01', env = data, raw.data = data.proxy.raw, auto.assign = T)
for(i in ls(data)) data[[i]] = adjustOHLC(data[[i]], use.Adjusted=T)
bt.prep(data, align='remove.na')

plota(data$DBC,type='l')
  plota.legend('DBC', 'black', data$DBC)

#*****************************************************************
# Helper functions
#*****************************************************************
donchian.channel.breakout.strategy = function(data, lockback.len, long.only=F, use.close=F) {
  if(use.close) {
    high.channel = bt.apply.matrix(data$prices, runMax, lockback.len)
    low.channel = bt.apply.matrix(data$prices, runMin, lockback.len)
   } else {
    phigh = bt.apply(data, Hi)
    plow = bt.apply(data, Lo)

    high.channel = bt.apply.matrix(phigh, runMax, lockback.len)
    low.channel = bt.apply.matrix(plow, runMin, lockback.len)
  }

  data$weight[] = NA
    if(use.close)
      data$weight[] = iif(data$prices == high.channel, 1, iif(data$prices == low.channel, -1, NA))
    else
      data$weight[] = iif(phigh == high.channel, 1, iif(plow == low.channel, iif(long.only,0,-1), NA)) 
  bt.run.share(data, clean.signal=T, silent=T)
}

percentile.channel.breakout.strategy = function(data, lockback.len, long.only=F) {
  high.channel = bt.apply.matrix(data$prices, runQuantile, lockback.len, probs=0.75)
  low.channel = bt.apply.matrix(data$prices, runQuantile, lockback.len, probs=0.25)

  data$weight[] = NA
    data$weight[] = iif(cross.up(prices, high.channel), 1, iif(cross.dn(prices, low.channel), iif(long.only,0,-1), NA))
  bt.run.share(data, clean.signal=T, silent=T)
}

#*****************************************************************
# Setup
#*****************************************************************
prices = data$prices

models = list()

#*****************************************************************
# Donchian Channel Breakout strategy
#*****************************************************************
models$Donchian55 = donchian.channel.breakout.strategy(data, 55)
models$Donchian55.Long = donchian.channel.breakout.strategy(data, 55, long.only=T)

print('#Transition Map for Donchian Channel Breakout strategy using 55 day lookback:')

plotbt.transition.map(models$Donchian55$weight)

#*****************************************************************
# Percentile Channel Breakout strategy
#*****************************************************************
models$Percentile55 = percentile.channel.breakout.strategy(data, 55)
models$Percentile55.Long = percentile.channel.breakout.strategy(data, 55, long.only=T)

print('#Transition Map for Percentile Channel Breakout strategy using 55 day lookback:')

plotbt.transition.map(models$Percentile55$weight)

models$Donchian20 = donchian.channel.breakout.strategy(data, 20)
models$Percentile20 = percentile.channel.breakout.strategy(data, 20)

#*****************************************************************
# Report
#strategy.performance.snapshoot(models, T)
plotbt(models, plotX = T, log = 'y', LeftMargin = 3, main = NULL)	    	
	mtext('Cumulative Performance', side = 2, line = 1)

print(plotbt.strategy.sidebyside(models, make.plot=F, return.table=T))

#*****************************************************************
# Load historical data
#*****************************************************************
library(SIT)
load.packages('quantmod')

# load saved Proxies Raw Data, data.proxy.raw
# please see http://systematicinvestor.github.io/Data-Proxy/ for more details
load('data/data.proxy.raw.Rdata')

tickers = '
EQ = VTI +VTSMX # (or SPY)
RE = IYR + VGSIX # (or ICF)
CORP.FI = LQD + VWESX
COM = DBC + CRB
CASH = SHY + TB3Y
'

data <- new.env()
getSymbols.extra(tickers, src = 'yahoo', from = '1970-01-01', env = data, raw.data = data.proxy.raw, set.symbolnames = T, auto.assign = T)
for(i in data$symbolnames) data[[i]] = adjustOHLC(data[[i]], use.Adjusted=T)

#print(bt.start.dates(data))
bt.prep(data, align='remove.na', fill.gaps = T)

# Check data
plota.matplot(scale.one(data$prices),main='Asset Perfromance')

#*****************************************************************
# Setup
#*****************************************************************
data$universe = data$prices > 0
  # do not allocate to CASH, or BENCH
  data$universe$CASH = NA

prices = data$prices * data$universe
  n = ncol(prices)
  nperiods = nrow(prices)


frequency = 'months'
# find period ends, can be 'weeks', 'months', 'quarters', 'years'
period.ends = endpoints(prices, frequency)
  period.ends = period.ends[period.ends > 0]

models = list()

commission = list(cps = 0.01, fixed = 10.0, percentage = 0.0)

# lag prices by 1 day
#prices = mlag(prices)

#*****************************************************************
# Equal Weight each re-balancing period
#******************************************************************
data$weight[] = NA
  data$weight[period.ends,] = ntop(prices[period.ends,], n)
models$ew = bt.run.share(data, clean.signal=F, commission = commission, trade.summary=T, silent=T)



#*****************************************************************
# Risk Parity each re-balancing period
#******************************************************************
ret = diff(log(prices))
hist.vol = sqrt(252) * bt.apply.matrix(ret, runSD, n = 20)

# risk-parity
weight = 1 / hist.vol
rp.weight = weight / rowSums(weight, na.rm=T)

data$weight[] = NA
  data$weight[period.ends,] = rp.weight[period.ends,]
models$rp = bt.run.share(data, clean.signal=F, commission = commission, trade.summary=T, silent=T)




#*****************************************************************
# Helper functions
#*****************************************************************
reallocate = function(allocation, data, period.ends,	lookback.len,
	prefix = '',
	min.risk.fns = 'min.var.portfolio',
	silent = F
) {
	allocation = ifna(allocation, 0)[period.ends,]

	# only allocate to selected assets
	obj = portfolio.allocation.helper(
   data$prices, 
   period.ends=period.ends,
   lookback.len = lookback.len,
   universe = allocation > 0,
   prefix = prefix,
   min.risk.fns = min.risk.fns,
   silent = silent
  )

  # rescale weights to be proportionate allocation
  for(i in names(obj$weights)) {
   weight = allocation * obj$weights[[i]]
   weight = ifna(rowSums(allocation) * weight / rowSums(weight), 0)
   weight$CASH = 1 - rowSums(weight)
   obj$weights[[i]] = weight
  }

	create.strategies(obj, data, silent = silent)$models
}



#*****************************************************************
# Strategy:
#
# 1) Use 60,120,180, 252-day percentile channels
# - corresponding to 3,6,9 and 12 months in the momentum literature- 
# (4 separate systems) with a .75 long entry and .25 exit threshold with 
# long triggered above .75 and holding through until exiting below .25 
# (just like in the previous post) - no shorts!!!
#
# 2) If the indicator shows that you should be in cash, hold SHY
#
# 3) Use 20-day historical volatility for risk parity position-sizing 
# among active assets (no leverage is used). This is 1/volatility (asset A) 
# divided by the sum of 1/volatility for all assets to determine the position size.
#******************************************************************

allocation = 0 * ifna(prices, 0)
for(lockback.len in c(60,120,180, 252)) {
  high.channel = bt.apply.matrix(data$prices, runQuantile, lockback.len, probs=0.75)
  low.channel = bt.apply.matrix(data$prices, runQuantile, lockback.len, probs=0.25)
  signal = iif(cross.up(prices, high.channel), 1, iif(cross.dn(prices, low.channel), 0, NA))
  allocation = allocation + ifna( bt.apply.matrix(signal, ifna.prev), 0) 
}

# convert to weights, i.e. 4 assets times 4 systems
allocation = allocation / (4 * 4)

# equal-weight
weight = allocation
weight$CASH = 1 - rowSums(weight)

data$weight[] = NA
  data$weight[period.ends,] = weight[period.ends,]
models$strategy.ew = bt.run.share(data, clean.signal=F, commission = commission, trade.summary=T, silent=T)


# risk-parity
weight = allocation * rp.weight
weight = rowSums(allocation, na.rm=T) * weight / rowSums(weight, na.rm=T)
weight$CASH = 1 - rowSums(weight, na.rm=T)

data$weight[] = NA
  data$weight[period.ends,] = ifna(weight[period.ends,], 0)
models$strategy.rp = bt.run.share(data, clean.signal=F, commission = commission, trade.summary=T, silent=T)


# let's apply leverage
weight1 = ifna(weight, 0)
data$weight[] = NA
	data$weight[period.ends,] = target.vol.strategy(models$strategy.rp, weight1, 
		target=6/100, lookback.len=2*20, max.portfolio.leverage=150/100)[period.ends,]
	rs = rowSums(data$weight, na.rm=T)[period.ends]
	data$weight$CASH[period.ends] = data$weight$CASH[period.ends] + iif(rs > 1, 0, 1-rs)
   	
models$strategy.rp.L = bt.run.share(data, clean.signal=F, commission = commission, trade.summary=T, silent=T)




#*****************************************************************
# Alternative Way
#*****************************************************************
models = c(models, 
  reallocate(allocation, data, period.ends,	20,
   prefix = 'REP.STG.',
   min.risk.fns = list(
    EW=equal.weight.portfolio,
    RP=risk.parity.portfolio()
   ),
   silent = T
  )
)

#*****************************************************************
#The [Quantitative Approach To Tactical Asset Allocation Strategy(QATAA) by Mebane T. Faber](http://mebfaber.com/timing-model/)
#[SSRN paper](http://papers.ssrn.com/sol3/papers.cfm?abstract_id=962461)
#******************************************************************
# compute 10 month moving average
sma = bt.apply.matrix(prices, SMA, 200)

# go to cash if prices falls below 10 month moving average
go2cash = prices < sma
  go2cash = ifna(go2cash, T)

# equal weight target allocation
target.allocation = ntop(prices,n)

# If asset is above it's 10 month moving average it gets allocation
qt.weight = iif(go2cash, 0, target.allocation)

# otherwise, it's weight is allocated to cash
weight = qt.weight
weight$CASH = 1 - rowSums(weight)

data$weight[] = NA
  data$weight[period.ends,] = weight[period.ends,]
models$QATAA.EW = bt.run.share(data, clean.signal=F, commission = commission, trade.summary=T, silent=T)


# same, but risk-parity
weight = qt.weight * rp.weight
weight = rowSums(qt.weight, na.rm=T) * weight / rowSums(weight, na.rm=T)
weight$CASH = 1 - rowSums(weight, na.rm=T)

data$weight[] = NA
  data$weight[period.ends,] = ifna(weight[period.ends,], 0)
models$QATAA.RP = bt.run.share(data, clean.signal=F, commission = commission, trade.summary=T, silent=T)


#*****************************************************************
# Alternative Way
#*****************************************************************
models = c(models, 
  reallocate(qt.weight, data, period.ends,	20,
   prefix = 'REP.QATAA.',
   min.risk.fns = list(
    EW=equal.weight.portfolio,
    RP=risk.parity.portfolio()
   ),
   silent = T
  )
)


#*****************************************************************
# Report
#*****************************************************************
#strategy.performance.snapshoot(models, T)
plotbt(models, plotX = T, log = 'y', LeftMargin = 3, main = NULL)
	mtext('Cumulative Performance', side = 2, line = 1)

print(plotbt.strategy.sidebyside(models, make.plot=F, return.table=T))

for(m in names(models)) {
  print('#', m)
  plotbt.transition.map(models[[m]]$weight, name=m)
    legend('topright', legend = m, bty = 'n')
  
  print('Last Trades:')  
  print(last.trades(models[m],make.plot=F, return.table=T))
                                              
  print('Current Allocation:')
  print(to.percent(last(models[[m]]$weight)))
}