# Vhayu/REvolution Computing High-Performance Computing with R
# Supplementary Material

# This script requires the following R packages:
install.packages('quantmod')
require ('quantmod')
install.packages('PerformanceAnalytics')
require ('PerformanceAnalytics')
require ('TTR')
install.packages('xts')
require(xts)

# We define a few convenience functions to convert between Vhayu SQL
# bar results and the xts time series.

# VHTime2ct: convert Vhayu time to POSIXct value
VHTime2ct <- function (v)
{
  t <- as.POSIXct (v/10000000, origin="1600-12-31")
  return (t)
}

# vbar2xts: convert Vhayu bar to xts bar
vbar2xts <- function (v)
{
  x <- xts (v[,2:ncol(v)], order.by=VHTime2ct(v[,1]))
  return (x)
}

# 1. Interacting with Vhayu Velocity Engine and basic data manipulation
# -----------------------------------------------------------------------
# The Vhayu databse supports SQL queries like the one shown below (note
# the order of of the requested fields--they match the xts bar order).
# conn <- odbcDriverConnect("DSN=Vhayu")
# AAPL <- sqlQuery (conn, "SELECT TICKTIME, Open, High, Low, Close, Volume FROM BarDaily WHERE SYMBOL = 'AAPL' AND TICKTIME BETWEEN VT('20051001') AND  VT()")

# The sample 'Vhayu.Rdata' file includes a small set of hourly data for a
# few equities, just as they would have been extracted from Vhayu with the above
# SQL commands.

load ('Vhayu.Rdata')

# Convert the example Vhayu SQL bar data to xts:
AAPL       <- (vbar2xts(AAPL_1h))
AAPL_clean <- (vbar2xts(AAPL_1h_clean))

# Elementary data manipulation:
head (AAPL)
AAPL[,]       # all rows, columns selected (cf Matlab (:,:) notation)

# We can extract rows by index:
AAPL[1:5,]
# or, by date and time index (note the out of range index selection):
# (see in a minute below...)

# We can extract columns by index:
AAPL[1:5, 2:4]
# or, by name:
AAPL[1:5, 'Close']
# or, by a list of names in any order:
AAPL[1:5, c('Close','Low')]


# The Vhayu Velocity engine includes very high performance data
# processing capabilities, including intervalization and data cleaning
# functions. We find that using Velocity for data pre-processing tasks
# and R for analysis and backtesting computation is one reasonable way
# to balance the computational abilities of each math engine.

# Illustration of Vhayu  data pre-processing (outlier removal)
# Raw, intervalied data:
chartSeries(AAPL['2005-10-03::2005-10-12'],TA=NULL)
# With outliers removed by Velocity:
chartSeries(AAPL_clean['2005-10-03::2005-10-12'],TA=NULL)

# Illustration of fancy quantmod charts!
addBBands()
addMACD(fast=12, slow=26, signal=9)

# The xts time-series class is flexible and fast. The following
# example illustrates a join on time-series defined over different 
# periods:
x <- AAPL['2005-10-01::2005-10-12', 'Close']
y <- AAPL_clean['2005-10-03::2005-10-08', 'Close']
z <- cbind (x,y)
colnames (z) <- c('AAPL','A')
z

# NAs were introduced by the join operation. We can easily expunge them:
na.omit (z)

# 2. Technical trading rules
# -----------------------------------------------------------------------
# For this part of the demo, we consider the INTC and IEF symbols on daily
# intervals over a longer time period. 

# First, let's chart INTC:
png(file="intc-macd.png",width=700, height=800)
chartSeries(INTC)
addMACD(fast=12, slow=26, signal=9)
dev.off()

# In order to emulate trading-like
# strategies, we use IEF (the iShares 10-year treasury instrument) as a
# benchmark.  We consider only the processed close price data, and 
# over a consistent time period:
z <- cbind (INTC$Close, IEF$Close)
colnames(z) <- c("INTC", "IEF")
z <- na.omit (z)

# Let's think about returns instead of prices...
# Ra is the log return for a buy-and-hold strategy, Rb the 'benchmark.'
# We will extensively use these log-return series in the sequel.
Ra <- Return.calculate(z$INTC)
Rb <- Return.calculate(z$IEF)
png(file="intc-ief-returns.png",width=700, height=400)
chart.CumReturns (cbind (Ra,Rb), main='Returns', legend.loc='topright')
dev.off()

# Define a very simple MACD-based example trading rule.  This rule 
# takes a buy or sell position on the 'long' instrument,
# converting into the 'benchmark' instrument on the sell signal.
# Note the use of named function arguments, optionally supplied with
# default values.
simpleRule <- function (z, fast=12, slow=26, signal=9, long, benchmark)
{
  x <- MACD(z, nFast=fast, nSlow=slow, nSig=signal, maType="EMA")
  position <- sign(x[,1]-x[,2])
  s <- xts(position,order.by=index(z))
  return (long*(s>0) + benchmark*(s<=0))
}

# Now let's evaluate the return from this rule and compute
# a Sharpe-like information ratio versus the benchmark.
# R is the log return from our simple rule, Dt the differential log return
R <- simpleRule (z$INTC, fast=12, slow=26, signal=9, long=Ra, benchmark=Rb)
png(file="default-returns.png",width=500, height=350)
chart.CumReturns(R, main="nFast=12 nSlow=26 nSig=9")
dev.off()

Dt <- na.omit(R - Rb)
sharpe <- mean(Dt)/sd(Dt)
print (paste("Ratio = ",sharpe))

# Compare against a buy-and-hold strategy over the time period:
SharpeRatio (na.omit(Ra))

# We see that the MACD trading rule with default parameter values outperforms
# just holding the INTC stock over this period. Can we find values of the
# parameters that do even better?

# One approach is to simply try a brute-force optimization of the 
# parameters:
# (Warning: This will take a long time if M is large!!!)
M <- 100
S <- matrix(0,M,M)
system.time(
for (j in 5:(M-1)) {
  for (k in min((j+2),M):M) {
    R <- simpleRule(z$INTC,j,k,9, Ra, Rb)
    Dt <- na.omit(R - Rb)
    S[j,k] <- mean(Dt)/sd(Dt)
  }
}
)

# The best Sharpe-like ratio and corresponding parameters
print (paste("Best ratio is ",max(S)))
which (S==max(S), arr.ind=TRUE)

# OK, the optimal parameters as expected perform a lot better. Typically,
# we would then like to test those parameters on other date ranges.



# 3. Parallel R Example: Parallel for loop
# ----------------------------------------------------------------------
# Unfortunately, brute force optimization and validation across dates
# can be very time consuming, even with simple computations as illustrated
# above. The bottle-neck in this kind of operation often consists of
# loops--in particular, nested for loops. In many cases, one of the
# for loops can simply be replaced by 'foreach.'
#
# The 'foreach' and '%do%'/'%dopar%' operators provide a looping
# construct that can be viewed as a hybrid of the standard 'for'
# loop and 'lapply' function. It looks similar to the 'for' loop,
# and it evaluates an expression, rather than a function (as in
# 'lapply'), but it ultimately returns a list as it's result, rather
# than working by side-effects. 

# The code below can run on one machine or across clusters, on Windows
# GNU/Linux and/or Mac OS X platforms. 
install.packages('doNWS')
require ('doNWS') # load the ParallelR packages

# For distributed computation on a cluster:
# s <- sleigh (workerCount=2, launch=hpcserv)

# For parallel processing on a dual-core workstation 
s <- sleigh(workerCount=2)

setSleigh(s)
registerDoNWS(s) # use NetWorkSpaces as the back-end engine
system.time (
SS <- foreach (j=5:(M-1), .combine=rbind, .packages=c('xts','TTR')) %dopar% {
    x <- rep(0,M)
    for (k in min((j+2),M):M) {
      R <- simpleRule(z$INTC,j,k,9,Ra,Rb)
      Dt <- na.omit(R - Rb)
      x[k] <- mean(Dt)/sd(Dt)
    }
    return(x)
  }
)

stopSleigh(s) # Shutdown parallel workers