Friday, 27 January 2012

Testing the Delta Phenomenon

The second book I ever bought about trading was Welles Wilder's "The Delta Phenomenon," which I bought as a result of a recommendation in the first book I ever bought. Subsequently I also bought the "The Adam Theory of Markets" and a few years later I bought the "Ocean Theory" book, so one could say I own the whole trilogy!

When I bought these I was doing my charting by hand on graph paper using prices from the Wall Street Journal, but in due course I got a computer and began using various software programs; Excel, Open Office Calc, QtStalker and finally ended up where I am today using OctaveR and Gnuplot. But however proficient I became at using these last three my programming skills weren't up to coding the Delta Phenomenon, until now that is. I had already quite easily coded the Adam Projection and the Natural Market Mirror, Natural Market River and Natural Moving Average from Ocean theory. Over the next few posts I am going to outline how I intend to test the Delta Phenomenon and show the eventual results of these tests, but before that I am going to present in this post the "breakthrough" piece of coding that finally allows me to do so. I think other users of R may find it useful.

The issue to be solved was overcoming the problem of missing days in the time series data, e.g. weekends, holidays, exchange closures, just plain missing data etc. because the Delta Phenomenon is predicated on counting days before or after specified dates, and of course any algorithmic counting over the data would be upset by such missing days. My forum query here and forum searches that turned up this led me to the R xts package, which finally resulted in me being able to write the following piece of R code
rm(list=ls(all=TRUE)) # remove all previous data from workspace
library(xts) # load the required library

# preparation for output to new file
sink(file="solution",append=FALSE,type=c("output"),split=FALSE)

# read in the raw data files
data_1 <- read.csv(file="sp",header=FALSE,sep=",")
itdmtd <- read.csv(file="itdmtd",header=FALSE,sep=",")

# create xts objects from above data:
x <- xts(data_1[c('V2','V3','V4','V5')],as.Date(data_1[,'V1']))
y <- xts(itdmtd[c('V2','V3','V4','V5','V6','V7','V8','V9')],as.Date(itdmtd[,'V1']))
z <- merge.xts(x,y) # merge x and y

# create a contiguous date vector to encompass date range of above data
d <- timeBasedSeq(paste(start(z),end(z),"d",sep="/"), retclass="Date")

# merge z with an "empty" xts object, xts(,d), filling with NA
prices <- merge(z,xts(,d),fill=NA)

# coerce prices xts object to a data frame object
prices_df <- data.frame(date=index(prices), coredata(prices))

# output to new file
write.table(prices_df,quote=FALSE,row.names=FALSE,col.names=FALSE,sep=",")
sink()
The code takes a csv file of the price time series, merges it with a csv file of the "Delta Solution," fills in any missing dates and then writes out the result to a final csv file that looks like this

1995-01-03,828.55,829.45,827.55,828.8,NA,NA,NA,NA,NA,NA,NA,NA
1995-01-04,830.2,831.25,827.6,831.1,NA,NA,NA,NA,NA,NA,NA,NA
1995-01-05,830.5,831.45,829.85,830.6,NA,NA,NA,NA,NA,NA,NA,NA
1995-01-06,830.75,833.05,829,830.35,NA,NA,NA,NA,NA,NA,NA,NA
1995-01-07,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA
1995-01-08,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA

This might not actually look like much, but using this as input to this Gnuplot script
reset
set title "Medium Term Delta Turning Points" textcolor rgb "#FFFFFF"
set object 1 rect from graph 0, 0, 0 to graph 1.1, 1.1, 0 behind lw 1.0 fc rgb "#000000" fillstyle solid 1.00
set datafile separator ","
set xdata time
set timefmt "%Y-%m-%d"
set format x
set y2range [0:1]

plot "solution" using 1:10 notitle with impulses linecolor rgb "#B0171F" axis x1y2, \
"solution" using 1:11 notitle with impulses linecolor rgb "#0000FF" axis x1y2, \
"solution" using 1:12 notitle with impulses linecolor rgb "#FFA500" axis x1y2, \
"solution" using 1:13 notitle with impulses linecolor rgb "#00EE00" axis x1y2, \
"solution" using 1:2:3:4:5 notitle with candlesticks linecolor rgb "#FFFFFF" axis x1y1, \
"solution" using 1:($$2>$$5?$$2:1/0):($$2>$$5?$$3:1/0):($$2>$$5?$$4:1/0):($$2>$$5?$$5:1/0) notitle with candlesticks lt 1 axis x1y1, \
"solution" using 1:($$2<$$5?$$2:1/0):($$2<$$5?$$3:1/0):($$2<$$5?$$4:1/0):($$2<$$5?$$5:1/0) notitle with candlesticks lt 3 axis x1y1
gives a nice plot thus

Many readers might say "So what! And...?" but to those readers who know what the Delta Phenomenon is, the coloured lines will be significant. Furthermore, using free, open source software I have created a free, as in gratis, alternative to the software that the Delta Society sells on its website. Most importantly of all, of course, is that I am now in a position to do computerised testing of the Delta Phenomenon. More on these tests in upcoming posts.

Tuesday, 10 January 2012

Formula for Approximating VWAP from OHLC Data

An answer to another question on the Quantitative Finance forum here gives a formula that approximates the volume weighted average price (VWAP) from OHLC data only, with an exhortation in a comment to do one's own research before using it. This blog post is my quick attempt at such research.

My simple test uses the relationship between the VWAP and the Pivot point calculated as (H+L+C)/3. This wiki entry on pivot points states "Trading above or below the pivot point indicates the overall market sentiment" and my test of being above or below is the VWAP being above or below: i.e. if today's VWAP is above today's pivot point (calculated from yesterday's H,L & C) then this is s bullish sign and therefore we can expect tomorrow's VWAP to be higher than today's. The converse is true for bearish sentiment. In Octave it is trivially simple to write a script to test whether this hypothesis is true and the terminal output of the results is shown in the box below.

The first line identifies the market being tested (a forex pair or commodity), the second "sign_correct" is the percentage of time the hypothesis appears to be true: i.e. for the first AUSCAD entry the figure 67.228 means that 67.228 % of the time the following day's VWAP is actually higher/lower than today's according to the "prediction" based on the categorisation of today's bar being bullish or bearish. The third "result" is the p-value of a 5000 repetition Monte Carlo permutation test to test the statistical significance of the given % accuracy. The permutation test was the "sign" test described in Evidence Based Technical Analysis and in the PDF available from its companion website.

Given that all the p-values are zero it can be stated that the VWAP being above or below the pivot point has a non random predictive accuracy for the next day's VWAP being higher or lower to a degree of accuracy that is statistically significant.
auscad 
sign_correct =  65.695
result = 0
------------------ 
aususd 
sign_correct =  67.228
result = 0
------------------ 
ausyen 
sign_correct =  67.100
result = 0
------------------ 
bo 
sign_correct =  63.989
result = 0
------------------ 
c 
sign_correct =  63.530
result = 0
------------------ 
cc 
sign_correct =  60.565
result = 0
------------------ 
cl 
sign_correct =  63.318
result = 0
------------------ 
ct 
sign_correct =  60.273
results = 0
------------------ 
dx 
sign_correct =  63.811
result = 0
------------------ 
ed 
sign_correct =  63.945
results = 0
------------------ 
euraus 
sign_correct =  68.024
results = 0
------------------ 
eurcad 
sign_correct =  66.854
result = 0
------------------ 
eurchf 
sign_correct =  65.707
result = 0
------------------ 
eurgbp 
sign_correct =  66.760
result = 0
------------------ 
eurusd 
sign_correct =  65.544
result = 0
------------------ 
euryen 
sign_correct =  66.444
result = 0
------------------ 
fc 
sign_correct =  61.905
result = 0
------------------ 
gbpchf 
sign_correct =  67.497
result = 0
------------------ 
gbpusd 
sign_correct =  66.936
result = 0
------------------ 
gbpyen 
sign_correct =  66.936
result = 0
------------------ 
gc 
sign_correct =  60.667
result = 0
------------------ 
hg 
sign_correct =  58.554
result = 0
------------------ 
ho 
sign_correct =  62.685
result = 0
------------------ 
kc 
sign_correct =  61.732
result = 0
------------------ 
lb 
sign_correct =  61.765
result = 0
------------------ 
lc 
sign_correct =  62.372
result = 0
------------------ 
lh 
sign_correct =  61.601
result = 0
------------------ 
ng 
sign_correct =  62.356
result = 0
------------------ 
o 
sign_correct =  60.705
result = 0
------------------ 
oj 
sign_correct =  61.848
result = 0
------------------ 
pa 
sign_correct =  62.497
result = 0
------------------ 
pb 
sign_correct =  59.116
result = 0
------------------ 
pl 
sign_correct =  60.737
result = 0
------------------ 
rb 
sign_correct =  63.107
result = 0
------------------ 
s 
sign_correct =  64.091
result = 0
------------------ 
sb 
sign_correct =  61.106
result = 0
------------------ 
si 
sign_correct =  59.563
result = 0
------------------ 
sm 
sign_correct =  63.810
result = 0
------------------ 
sp 
sign_correct =  66.954
result = 0
------------------ 
es 
sign_correct =  66.744
result = 0
------------------ 
nd 
sign_correct =  66.221
result = 0
------------------ 
ty 
sign_correct =  65.260
result = 0
------------------ 
us 
sign_correct =  65.893
result = 0
------------------ 
usdcad 
sign_correct =  67.357
result = 0
------------------ 
usdchf 
sign_correct =  67.088
result = 0
------------------ 
usdyen 
sign_correct =  66.947
result = 0
------------------ 
w 
sign_correct =  65.118
result = 0
At the moment I use the value given by the pivot point calculation as my "price input" for the day but, based on these test results and the future results of some other tests I can think of, I may change to the VWAP approximation as my price input. More in a future post.

Tuesday, 3 January 2012

Creating Synthetic Data Using the Fast Fourier Transform

This question was recently asked on the Quantitative Finance forum and one answerer's link to here piqued my interest. My comments to said answer were
"Interesting link! I theorise that one way it could be done is to apply a Fast Fourier Transform (e.g. FFT function in Octave/MATLAB) to the original series, then divide the resulting vector into segments and within each separate segment randomly permute the values, and finally reconstruct the time series by applying the Inverse Fast Fourier Transform (IFFT function). By randomly permuting the values you will not alter the amplitudes of the composite sine waves but will slightly alter their periodicity - this will naturally tend to restrict the range of the new series to that of the original."
and
"This is such an intriguing approach that I may "knock up" some Octave code on my blog and post a link to it from here."
This blog post is my attempt to "knock up" the Octave code.

The "guts" of the code are a .oct function that randomly permutes the FFT output vector that is given to it as its only input argument and the C++ code for this .oct function is given below.
Note: due to formatting issues the headers are not correctly displayed; the octave headers should be enclosed between a "<" and ">" after the #include.
#include octave/oct.h
#include octave/dColVector.h
#include octave/CNDArray.h
#include "MersenneTwister.h"
     
DEFUN_DLD (permute_vector, args, , "Input is a vector that is to be permuted once")
{
octave_value_list retval;

int nargin = args.length () ;
int vec_length = args(0).length () ;

// check the input arguments
    if ( nargin > 1 )
    {
        error ("Invalid arguments. Input is a vector that is to be permuted once") ;
        return retval ;
    }
    
    if (vec_length < 2 )
    {
        error ("Invalid arguments. Input is a vector that is to be permuted once") ;
        return retval ;
    }
    
    if (error_state)
    {
        error ("Invalid arguments. Input is a vector that is to be permuted once") ;
        return retval ;
    }
// end of input checking

 ComplexNDArray input_vector = args(0).complex_array_value () ;
 ComplexNDArray output_vector = args(0).complex_array_value () ;
 int k1;
 int k2;
 
 MTRand mtrand1; // Declare the Mersenne Twister Class - will seed from system time

       k1 = vec_length - 1; // initialise prior to shuffling the vector

       while (k1 > 0) // While at least 2 left to shuffle
             {          
             k2 = mtrand1.randInt( k1 ); // Pick an int from 0 through k1 

               if (k2 > k1) // check if random vector index no. k2 is > than max vector index - should never happen 
                  {
                  k2 = k1 - 1; // But this is cheap insurance against disaster if it does happen
                  }

             output_vector(k1) = input_vector(k2) ; // allocate random pick k2 from input_vector to the k1 vector index of output_vector
             input_vector(k2) = input_vector(k1) ; // replace random pick k2 content of input_vector with content k1 of input_vector
             k1 = k1 - 1; // count down 
             } // Shuffling is complete when this while loop exits

 retval(0) = output_vector ; 
 
return retval; // Return the output to Octave
} 
The Octave script that calls the above function is
clear all

contract = input( "Enter contract symbol e.g. sp: ","s") ;
data = load("-ascii",contract) ;
n = rows(data)
index_begin = input( "Enter index_begin: ") ;
index_end = input( "Enter index_end, value not greater than n: ") ;

close = data(index_begin:index_end,7) ;

% detrend the close vector prior to applying the fft
slope = ( close(end) - close(1) ) / ( length(close) - 1 ) ;
v1 = (0:1:length(close)-1)' ;
detrended_close = close .- ( v1 .* slope ) ;
close_index_begin = close(1)
detrended_close_index_begin = detrended_close(1)
detrended_close_index_end = detrended_close(end)

% create zero padded vector for fft
L2 = 2^nextpow2( length(close) ) ; half_L2 = L2/2 ;
y2 = zeros( 1,L2 ) ; y2( 1:length(close) ) = detrended_close ;

% apply the fft
transform = fft( y2 ) ;

% permute the first half of the transform vector in "chunks" of 10
max_ii_value = floor( half_L2 / 10 ) ;
for ii = 1:max_ii_value
transform( (ii*10):(ii*10)+9 ) = permute_vector( transform( (ii*10):(ii*10)+9 ) ) ;
endfor

% take the inverse fft
ifft_vec = real( ifft( transform ) ) ;

% retrend the ifft_permuted_vec
retrended_ifft_vec = ifft_vec( 1:length(close) )' .+ ( v1 .* slope ) ;

% statistics
correl = corrcoef (close, retrended_ifft_vec)
spear = spearman (close, retrended_ifft_vec)
tau = kendall (close, retrended_ifft_vec)

plot( close,'b',retrended_ifft_vec,'r' ) ; legend( 'close','permute' ) ; 
This script could undoubtedly be made more efficient by putting the for loop on lines 28 to 30 within the .oct function itself, but this is just a proof of concept. The script allows one to extract any contiguous section of data by use of the index_begin and index_end inputs or, of course, to select the whole data range. Some simple statistics are calculated using in-built Octave functions and can be seen in the terminal windows in the screen shots. Finally the script plots the original input data along with the FFT permuted output data. Typical screen shots are shown below, with the original data in blue and the FFT permuted data in red.

All S&P data

The last few years of the above

The same last few years indexed using the index_begin and index_end method to apply the function just to these years and not the whole data set
The close_index_begin, detrended_close_index_begin and detrended_close_index_end values also seen in the terminal windows are a simple sanity check to see that the original data has been detrended correctly prior to performing the FFT.

Wednesday, 14 December 2011

Sunday, 11 December 2011

The Perfect Oscillator on S&P E-Mini Data

Following on from my previous posts the chart below is a snapshot of the last 150 days of the S&P E-mini continuous futures contract, up to and including the 9th December.


The top pane is obviously a plot of price showing triangles to indicate cyclic highs and lows when the Perfect Oscillator (middle pane) crosses its SMA (see previous post for explanation). As explained earlier the algorithm used in calculating the PO makes use of projected/predicted prices, and since these "forward prices" are available from the .oct function the bars are coloured according to the 5 and 10 day "forward prices." Based on a very simple if/or statement in the code essentially we should be long when the bars are blue and short when red. The few whites bars that are visible are bars where the conditions for being blue or red are not met.

Superimposed over the middle pane plot of the PO are the full period EMA and plus and minus 1 exponential standard deviation lines, a sort of Bollinger Bands over the indicator. The bottom plot is of the bandpass indicator which is used in the "forward prices" algorithm.

Based on the logic and theory of the PO the rules for use are essentially:-
1) go long/short according to the PO SMA crossings
2) confirm the above by the bar colours
3) exit/close up stops when the PO crosses its exponential standard deviation lines

and based on these simple rules a short entry is indicated for 12th December, particularly since the bandpass is also indicating a short trade with its leading function cross. Will be interesting to see how this indicated trade works out.

Friday, 18 November 2011

Update #2 on Perfect Moving Average and Oscillator

Below are screenshots of the PMA and PO based upon the improved price projection algorithm discussed in the previous post, shown on an idealised time series of a sinewave (plus trend) in cyan in all screenshots.

The first 3 are plots of the 10 bar PMA in red with a normal 10 bar EMA in green shown for comparative purposes.
As the theory suggests it can be seen that the PMA has the equivalent smoothing of the 10 bar EMA but without any lag.

The next series of screenshots are of the PO, where again the "price" is cyan, the various red lines are 2 bar, 3 bar etc. POs up to an 11 bar PO, and the green PO is a simple average of all these POs.
As I had anticipated these POs lead price by a quarter period which means the cycle highs and lows are indicated by the green "average PO" crossing its centre line, which in turn is a full cycle period length simple moving average of the average PO. My leading oscillator indicator should give timely warnings of these upcoming turns in price. It is also interesting to note that the individual, longer length PO crossings (greater amplitude) also give advance warning of the upcoming turn and the fact that the shorter length POs NOT crossing the centre line and being entirely above or below said centre line seem to indicate the strength of any trend.

Wednesday, 16 November 2011

Update on the Perfect Moving Average and Oscillator

It has been some time since my last post and the reason for this was, in large part, due to my changing over from an Ubuntu to a Kubuntu OS and the resultant teething problems etc. However, all that is over so now I return to the Perfect Moving Average (PMA) and the Perfect Oscillator (PO).

In my previous post I stated that I would show these on real world data, but before doing this there was one thing I wanted to investigate - the accuracy of the price projection algorithm upon which the average and oscillator calculations are based. Below are four sample screenshots of these investigations.

Sideways Market

Up Trending Market

Down Trending Market

The charts show a simulated "ideal" time series (cyan) with two projected price series of 10 points each (red and green) plotted such that the red and green projections are plotted concurrently with the actual series points they are projections, or predictions, of. Both projections/predictions start at the point at which the red "emerges" from the actual series. As can be seen the green projection/prediction is much more accurate than the red, and the green is the improved projection algorithm I have recently been working on.

As a result of this improvement and due to the nature of the internal calculations I expect that the PMA will follow prices much closer and that the PO will lead prices by a quarter cycle compared with my initial price projection algorithm. More on this in the next post.