Tuesday 24 April 2012

First Use of Rcpp Package

In preparation for the tests I want to conduct, I have finally started to use the Rcpp package. The reason for this decision is straightforward; I cannot begin to imagine how the code for the tests I want to conduct could be vectorised, and since loops in R are slow and not recommended there is no realistic alternative to biting the Rcpp bullet. Another reason is that some R packages, such as the ttrTests package, require an R function as part of the input and using Rcpp will enable me to transfer my .oct C++ function coding directly into R without having to rewrite in R code itself. A simple proof of concept programming exercise I set myself is shown below.
rm(list=ls()) # clear the entire workspace
library(xts)  # load the required library
library(zoo)  # load the required library
library(Rcpp) # load the required library
library(inline) # load the required library
library(compiler) # load the required library

# load the "indicator" file 
data <- read.csv(file="sind",head=FALSE,sep=,)
tick_size <- 0.025
tick_value <- 12.50

# extract other vectors of interest
#open <- data[,2]
market_mode <- data[,228]
kalman <- data[,283]

# source the Rcpp inline file for the C++ code for the test in question
# this file effectively creates a function that takes input and gives test 
# output in the desired form e.g. equity curve, position vector etc. for
# further statistical tests in the R environment, e.g the ttrTests package.

source("basic_market_mode_equity.r")
results <- basic_market_mode_equity(data[,2],market_mode,kalman,tick_size,tick_value)
This is "normal" R code which
  • loads the required libraries
  • loads the required file(s) from disk and takes other inputs for the test in question
  • extracts the required information from the loaded file(s)
  • sources and then calls the test function, named "basic_market_mode_equity"
This second code block contains the actual C++ test function code, which is contained in a file named "basic_market_mode_equity.r"
# This function takes as inputs vectors of opening prices, the "market mode," 
# the kalman filter of vwap, and single values for tick size and tick value.
# The values of "market_mode" are as follows:-
# 0 = cyclic
# 1 = up with retracement
# 2 = up with no retracement
# 3 = down with retracement
# 4 = down with no retracement
# 
# This first test is very simple - be long when market mode is 1 or 2, 
# short when 3 or 4, and when 0 be long or short depending on the direction
# of the kalman filter. This test is just practice in coding using Rcpp:- 
# don't expect any meaningful results. Equity curves for each market mode 
# are the function outputs

src <- '
Rcpp::NumericVector open(a) ;
Rcpp::NumericVector market_mode(b) ;
Rcpp::NumericVector kalman(c) ;
Rcpp::NumericVector tick_size(d) ;
Rcpp::NumericVector tick_value(e) ;
int n = open.size() ;
Rcpp::NumericVector cyc(n) ; // create output vector
Rcpp::NumericVector uwr(n) ; // create output vector
Rcpp::NumericVector unr(n) ; // create output vector
Rcpp::NumericVector dwr(n) ; // create output vector
Rcpp::NumericVector dnr(n) ; // create output vector

// fill the equity_curve with zeros for "burn in" period
for ( int ii = 0 ; ii < 100 ; ii++ ) {
    cyc[ii] = 0.0 ;
    uwr[ii] = 0.0 ;
    unr[ii] = 0.0 ;
    dwr[ii] = 0.0 ;
    dnr[ii] = 0.0 ; }

for ( int ii = 100 ; ii < n ; ii++ ) {
    
    // uwr market type
    if ( market_mode[ii-2] == 1 ) { 
    cyc[ii] = cyc[ii-1] ;
    uwr[ii] = uwr[ii-1] + tick_value[0] * ( (open[ii]-open[ii-1])/tick_size[0] ) ;
    unr[ii] = unr[ii-1] ;
    dwr[ii] = dwr[ii-1] ;
    dnr[ii] = dnr[ii-1] ; }

    // unr market type
    if ( market_mode[ii-2] == 2 ) { 
    cyc[ii] = cyc[ii-1] ;
    uwr[ii] = uwr[ii-1] ; 
    unr[ii] = unr[ii-1] + tick_value[0] * ( (open[ii]-open[ii-1])/tick_size[0] ) ;
    dwr[ii] = dwr[ii-1] ;
    dnr[ii] = dnr[ii-1] ; }

    // dwr
    if ( market_mode[ii-2] == 3 ) { 
    cyc[ii] = cyc[ii-1] ;
    uwr[ii] = uwr[ii-1] ;
    unr[ii] = unr[ii-1] ;
    dwr[ii] = dwr[ii-1] + tick_value[0] * ( (open[ii-1]-open[ii])/tick_size[0] ) ; 
    dnr[ii] = dnr[ii-1] ; }

    // dnr
    if ( market_mode[ii-2] == 4 ) { 
    cyc[ii] = cyc[ii-1] ;
    uwr[ii] = uwr[ii-1] ;
    unr[ii] = unr[ii-1] ;
    dwr[ii] = dwr[ii-1] ; 
    dnr[ii] = dnr[ii-1] + tick_value[0] * ( (open[ii-1]-open[ii])/tick_size[0] ) ; }

    // cyc long
    if ( market_mode[ii-2] == 0 && kalman[ii-2] > kalman[ii-3] ) { 
    cyc[ii] = cyc[ii-1] + tick_value[0] * ( (open[ii]-open[ii-1])/tick_size[0] ) ;
    uwr[ii] = uwr[ii-1] ;
    unr[ii] = unr[ii-1] ;
    dwr[ii] = dwr[ii-1] ; 
    dnr[ii] = dnr[ii-1] ; }

    // cyc short
    if ( market_mode[ii-2] == 0 && kalman[ii-2] < kalman[ii-3] ) { 
    cyc[ii] = cyc[ii-1] + tick_value[0] * ( (open[ii-1]-open[ii])/tick_size[0] ) ;
    uwr[ii] = uwr[ii-1] ;
    unr[ii] = unr[ii-1] ;
    dwr[ii] = dwr[ii-1] ; 
    dnr[ii] = dnr[ii-1] ; }

} // end of main for loop

return List::create(
  _["cyc"] = cyc ,
  _["uwr"] = uwr ,
  _["unr"] = unr ,
  _["dwr"] = dwr ,
  _["dnr"] = dnr ) ; '

basic_market_mode_equity <- cxxfunction(signature(a = "numeric", b = "numeric", c = "numeric",
                                d = "numeric", e = "numeric"), body=src, 
                                plugin = "Rcpp")
which basically just outputs five equity curves corresponding to the identified market mode (see the comments in the code for more details).

This final code block is again "normal" R code
# coerce the above results list object to a data frame object
results_df <- data.frame( results )
df_max <- max(results_df) # for scaling of results plot
df_min <- min(results_df) # for scaling of results plot

# and now create an xts object for plotting
results_xts <- xts(results_df,as.Date(data[,'V1']))

# a nice plot of the results_xts object
plot(results_xts[,'cyc'],ylim=c(df_min,df_max),type="l")
par(new=TRUE,col="cyan")
plot(results_xts[,'uwr'],ylim=c(df_min,df_max),type="l")
par(new=TRUE,col="blue")
plot(results_xts[,'unr'],ylim=c(df_min,df_max),type="l")
par(new=TRUE,col="green")
plot(results_xts[,'dwr'],ylim=c(df_min,df_max),type="l")
par(new=TRUE,col="red")
plot(results_xts[,'dnr'],ylim=c(df_min,df_max),type="l")
which produces this typical plot output.
This test is just a toy test and the results are not really important. The important thing is that I can now easily import the output of my Octave .oct C++ functions into R and conduct a whole range of tests utilising R packages from CRAN or simply write my own test routines in C++ (my intention) to run in R and not have to struggle with vectorising code for speed optimisation purposes.

No comments: