Sortino vs Omega ratios

library(SIT)
load.packages('quantmod')

#*****************************************************************
# Helper functions
#*****************************************************************
sortino.ratio = function(R,MAR) {
	r = R[R < MAR]
	(mean(R) - MAR) / sqrt(sum((MAR - r)^2) / length(r) ) 
}

omega.ratio = function(R,MAR) {
	#sum(pmax(R - MAR, 0)) / sum(pmax(MAR - R, 0))
	r = R - MAR
	- sum(r[r > 0]) / sum(r[r < 0])
}

#*****************************************************************
# check: Omega has a value of 1 at the mean of the distribution
#*****************************************************************
R = rnorm(100)
MAR = mean(R)

print(list(Sortino = sortino.ratio(R, MAR), Omega = omega.ratio(R,MAR)))

#*****************************************************************
# Plot Sortino vs Omega
# based on [R graph with two y-axes](http://robjhyndman.com/hyndsight/r-graph-with-two-y-axes/)
#*****************************************************************
plot.sortino.omega = function(R, main='') {
  x = quantile(R, c(0.05, 0.95))
  x = seq(x[1],x[2], length.out=100)

  par(mar=c(5,4,4,5))
  plot(x, sapply(x, function(y) sortino.ratio(R,y)), xlab='MAR', ylab='Sortino', 
    main=paste('Sortino vs Omega', main), type='l', col='red')

  par(new=TRUE)
  plot(x, sapply(x, function(y) omega.ratio(R,y)), xaxt="n",yaxt="n",xlab="", ylab='', type='l', col='blue')

  axis(4)
  mtext('Omega',side=4,line=3)
  legend("topright",col=c('red','blue'),lty=1,legend=c('Sortino','Omega'), bty='n')
}

plot.sortino.omega(R, ' based on random data')

#*****************************************************************
# Load historical data
#*****************************************************************
data = env()
getSymbols('SPY', src = 'yahoo', from = '1970-01-01', env = data, auto.assign = T)

price = Ad(data$SPY)
ret = diff(log(price))

plot.sortino.omega(last(ret, 252), ' based on last year of SPY')

	#--------------------------------------------------------------------------
	# Create Efficient Frontier
	#--------------------------------------------------------------------------
	ia = aa.test.create.ia()
	n = ia$n		

	# 0 <= x.i <= 0.8 
	constraints = new.constraints(n, lb = 0, ub = 0.8)
	
	# SUM x.i = 1
	constraints = add.constraints(rep(1, n), 1, type = '=', constraints)		
	
	# Omega - http://en.wikipedia.org/wiki/Omega_ratio
		# do not set threshold too high
		ia$parameters.omega = 1/100 
		# convert annual to monthly
		ia$parameters.omega = ia$parameters.omega / 12

	max.omega.portfolio(ia, constraints)		

	# create efficient frontier(s)
	ef.risk = portopt(ia, constraints, 50, 'Risk')

	#plot.ef(ia, list(ef.risk), portfolio.risk, T, T)			
	
	
	# Plot Omega Efficient Frontiers and Transition Maps
	layout( matrix(1:4, nrow = 2, byrow=T) )
	
	# weights
	rownames(ef.risk$weight) = paste('Risk','weight',1:50,sep='_')
	plot.omega(ef.risk$weight[c(1,10,40,50), ], ia)
	
	# assets
	temp = diag(n)
	rownames(temp) = ia$symbols
	plot.omega(temp, ia)
		
	# portfolio
	plot.ef(ia, list(ef.risk), portfolio.omega, T, T)			

	#--------------------------------------------------------------------------
	# Create Efficient Frontier in Omega Ratio framework
	#--------------------------------------------------------------------------
	# Create maximum Omega Efficient Frontier
	ef.omega = portopt.omega(ia, constraints, 50)

	# Plot Omega Efficient Frontiers and Transition Maps
	layout( matrix(1:4, nrow = 2, byrow=T) )

	# weights
	plot.omega(ef.risk$weight[c(1,10,40,50), ], ia)	
	#plot.omega(ef.omega$weight[c(1,10,40,50), ], ia)
	

	# weights
	rownames(ef.omega$weight) = paste('Omega','weight',1:50,sep='_')	
	plot.omega(ef.omega$weight[c(1,10,40,50), ], ia)
		
	# portfolio
	plot.ef(ia, list(ef.omega, ef.risk), portfolio.omega, T, T)			

	#--------------------------------------------------------------------------
	# Compare Risk and Omega frontiers
	#--------------------------------------------------------------------------
	# Plot multiple Efficient Frontiers and Transition Maps
	layout( matrix(1:4, nrow = 2) )
	plot.ef(ia, list(ef.risk,ef.omega), portfolio.risk, F)			
	plot.ef(ia, list(ef.risk,ef.omega), portfolio.omega, F)			

	plot.transition.map(ef.risk)
	plot.transition.map(ef.omega)

#' @export 	
target.omega.portfolio <- function
(
	target.omega,
	type = 'mixed'
)
{
	target.omega = target.omega

	function
	(
		ia,			# input assumptions
		constraints	# constraints
	)
	{
		ia$parameters.omega = target.omega
		max.omega.portfolio(ia, constraints, type)
	}	
}
 
#*****************************************************************
# Load historical data
#*****************************************************************
library(SIT)
load.packages('quantmod')

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

tickers = '
LQD + VWESX
DBC + CRB
VTI +VTSMX # (or SPY)
ICF + VGSIX # (or IYR)
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)

	
#*****************************************************************
# Code Strategies
#****************************************************************** 					
obj = portfolio.allocation.helper(data$prices, 
		periodicity = 'months', lookback.len = 60, 
		min.risk.fns = list(
			EW=equal.weight.portfolio,
			RP=risk.parity.portfolio(),
			MD=max.div.portfolio,						
			
			MV=min.var.portfolio,
			MVE=min.var.excel.portfolio,
			MO=target.omega.portfolio(0/100, 'lp'),
			
			MC=min.corr.portfolio,
			MCE=min.corr.excel.portfolio,
			MC2=min.corr2.portfolio,
			
			MS=max.sharpe.portfolio()
		),
		silent=T
)
	
commission = list(cps = 0.01, fixed = 10.0, percentage = 0.0)

models = create.strategies(obj, data, dates='1997::', commission=commission, silent=T)$models
	
#*****************************************************************
# Create 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, perfromance.fn=engineering.returns.kpi))

data = bt.change.periodicity(data, periodicity = 'months')


#*****************************************************************
# Code Strategies
#****************************************************************** 					
obj = portfolio.allocation.helper(data$prices, 
		periodicity = 'months', lookback.len = 24, 
		min.risk.fns = list(
			EW=equal.weight.portfolio,
			RP=risk.parity.portfolio(),
			MD=max.div.portfolio,						
			
			MV=min.var.portfolio,
			MVE=min.var.excel.portfolio,
			MO=target.omega.portfolio(0/100, 'lp'),
			
			MC=min.corr.portfolio,
			MCE=min.corr.excel.portfolio,
			MC2=min.corr2.portfolio,
			
			MS=max.sharpe.portfolio()
		),
		silent=T
)
	
models = create.strategies(obj, data, dates='1997::', commission=commission, silent=T)$models
	
#*****************************************************************
# Create 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, perfromance.fn=engineering.returns.kpi))