"Trading is statistics and time series analysis." This blog details my progress in developing a systematic trading system for use on the futures and forex markets, with discussion of the various indicators and other inputs used in the creation of the system. Also discussed are some of the issues/problems encountered during this development process. Within the blog posts there are links to other web pages that are/have been useful to me.
Sunday, 26 February 2012
Testing Delta, Conclusion
One should be cautioned, however, that statistical significance is not necessarily of practical significance, and perhaps one should ask why the Delta Phenomenon should work at all? A lot of the criticism I have read on forums etc. obviously relates to the difficulty of actually using Delta in trading. Another theme of criticism relates to the supposed "astrological" aspect of Delta and therefore some rather dismissive comments relegating Delta to the realm of crystal-ball gazing, mumbo-jumbo. I think this second form of criticism is ill-founded and rather ignorant. If you think about it, the interaction of the Sun and Earth, the phases of the moon etc. are just another way of saying "seasonality" and I think it is pretty well accepted that there is seasonality in, for example, agricultural commodity futures, oil futures etc. This seasonality may well be reflected in the share prices/indices by the linkage of large commercial interests in these areas, and various knock-on effects in trade, currency movements etc. My take on Delta is that it can be regarded as a sophisticated form of seasonal analysis, and all the normal caveats that would apply to trading seasonal tendencies should also apply to Delta.
Accepting that there could be some reasonable fundamental justification for Delta and that it has been shown to be statistically significant in terms of its accuracy as described by the methodology, how does one actually trade it? Some guidelines might be gleaned from these sample reports from a Delta timing service. Reading these sample reports shows how a knowledge of Delta can be used to draw inferences from current market action and to perhaps formulate a suitable trading strategy. Of course, analysis such as this is not a 100% mechanical approach, but it may very well add real value to the bottom line.
Given this, I have decided to deploy all the other Delta solutions I have to create Delta Charts that look like this
where the solid vertical lines are the MTD solution and the dashed the ITD solution. Doing this across the 30 or so commodities I regularly track will probably take me a few weeks. Once done, I will think about how to combine Delta with the other indicators I have mentioned in previous posts and test a trading system based on them. More in due course.
Monday, 30 January 2012
Testing the Delta Phenomenon, Part 3
Firstly, I am going to use the statistical concept of Cross-validation, whereby a model is trained on one set of data and tested for accuracy on another set of data, which is actually quite common in the world of back testing. Since the "Delta solutions" I will be testing come from late 2006 (See Testing the Delta Phenomenon, Part 2) the 4 years from 2008 to 2011 inclusive can be considered to be the validation set for this test. The test(s) will assess the accuracy of the predicted highs and lows for this period, on each Delta time frame for which I have solutions, by counting the difference in days between actual highs and lows in the data and their predicted occurrences and then creating an average error for these differences. This will be the test statistic. Using R, a Null Hypothesis average error distribution for random predictions on the same data will be created, using the same number of predicted turning points as per the Delta solution being tested. The actual average error test statistic on real data will be compared with this Null Hypothesis distribution of the average error test statistic and the Null Hypothesis rejected or not, as the case may be. The Null Hypothesis may be stated as
- given a postulated number of turning points the accuracy of the Delta Phenomenon in correctly predicting when these turning points will occur, using the average error test statistic described above as the measure of accuracy, is no better than random guessing as to where the turning points will occur.
- given a postulated number of turning points the accuracy of the Delta Phenomenon in correctly predicting when these turning points will occur, using the average error test statistic described above as the measure of accuracy, is better than could be expected from random guessing as to where the turning points will occur.
All of this might be made clearer for readers by following the commented R code below.
# Assume a 365 trading day year, with 4 Delta turning points in this year
# First, create a Delta Turning points solution vector, the projected
# days on which the market will make a high or a low
proj_turns <- c(47,102,187,234) # day number of projected turning points
# now assume we apply the above Delta solution to future market data
# and identify, according to the principles of Delta, the actual turning
# points in the future unseen "real data"
real_turns <- c(42,109,193,226) # actual market turns occur on these days
# calculate the distance between the real_turns and the days on which
# the turn was predicted to occur and calculate the test statistic
# of interest, the avge_error_dist
avge_error_dist <- mean( abs(proj_turns - real_turns) )
print(avge_error_dist) # print for viewing
# calculate the theoretical probability of randomly picking 4 turning points
# in our 365 trading day year and getting an avge_error_dist that is equal
# to or better than the actual avge_error_dist calculated above.
# Taking the first projected turning point at 47 and the actual turning
# that occurs at 42, to get an error for this point that is as small as or
# smaller than that which actually occurs, we must randomly choose one of
# the following days: 42,43,44,45,46,47,48,49,50,51 or 52. The probability of
# randomly picking one of these numbers out of 1 to 365 inclusive is
a <- 11/365
# and similarly for the other 3 turning points
b <- 15/364 # turning point 2
c <- 13/363 # turning point 3
d <- 17/362 # turning point 4
# Note that the denominator decreases by 1 each time because we are
# sampling without replacement i.e. it is not possible to pick the same
# day more than once. Combining the 4 probabilities above, we get
rdn_prob_as_good <- (a*b*c*d)/100 # expressed as a %
print( rdn_prob_as_good ) # a very small % !!!
# but rather than rely on theoretical calculations, we are actually
# going to repeated, randomly choose 4 turning points and compare their
# accuracy with the "real accuracy", as measured by avge_error_dist
# Create our year vector to sample, consisting of 365 numbered days
year_vec <- 1:365
# predefine vector to hold results
result_vec <- numeric(100000) # because we are going to resample 100000 times
# count how many times a random selection of 4 turning points is as good
# as or better than our "real" results
as_good_as = 0
# do the random turning point guessing, resampling year_vec, in a loop
for(i in 1:100000) {
# randomly choose 4 days from year_vec as turning points
this_sample <- sample( year_vec , size=4 , replace=FALSE )
# sort this_sample so that it is in increasing order
sorted_sample <- sort( this_sample , decreasing=FALSE )
# calculate this_sample_avge_error_dist, our test statistic
this_sample_avge_error_dist <- mean( abs(proj_turns - sorted_sample) )
# if the test statistic is as good as or better that our real result
if( this_sample_avge_error_dist <= avge_error_dist ) {
as_good_as = as_good_as + 1 # increment as_good_as count
}
# assign this sample result to result_vec
result_vec[i] <- this_sample_avge_error_dist
}
# convert as_good_as to %
as_good_as_percent <- as_good_as/100000
# some summary statistics of result_vec
mean_of_result_vec <- mean( result_vec )
standard_dev_of_result_vec <- sd( result_vec )
real_result_from_mean <- ( mean_of_result_vec - avge_error_dist )/standard_dev_of_result_vec
print( as_good_as ) # print for viewing
print( as_good_as_percent ) # print for viewing
print( mean_of_result_vec ) # print for viewing
print( standard_dev_of_result_vec ) # print for viewing
print( real_result_from_mean ) # print for viewing
# plot histgram of the result_vec
hist( result_vec , freq=FALSE, col='yellow' )
abline( v=avge_error_dist , col='red' , lwd=3 )
Typical output of this code iswhich shows a histogram of the distribution of random prediction average errors in yellow, with the actual average error shown in red. This is for the illustrative hypothetical values used in the code box above. Terminal prompt output for this is
[1] 6.5
[1] 2.088655e-08
[1] 38
[1] 0.00038
[1] 63.78108
[1] 32.33727
[1] 1.771364
where
6.5 is actual average error in days
2.088655e-08 is "theoretical" probability of Delta being this accurate
38 is number of times a random prediction is as good as or better than 6.5
0.00038 is 38 expressed as a percentage of random predictions made
63.78108 is the mean of the random distribution histogram
32.33727 is the standard deviation of the random distribution histogram
1.771364 is the difference between 63.78108 and 6.5 expressed as a multiple of the 32.33727 standard deviation.
This would be an example of the Null Hypothesis being rejected due to the 0.00038 % figure for random prediction accuracy being better than actual accuracy; in statistical parlance - a low p-value. Note, however, gross the difference between this figure and the "theoretical" figure. Also note that despite the Null being rejected the actual average error falls well within 2 standard deviations from the mean of the random distribution. This of course is due to the extremely heavy right-tailedness of the distribution, which expands the standard deviation range.
This second plot
and
[1] 77.75
[1] 2.088655e-08
[1] 48207
[1] 0.48207
[1] 79.85934
[1] 27.60137
[1] 0.07642148
shows what a typical failure to reject the Null Hypothesis would look like - a 0.48 p-value - and an actual average error that is indistinguishable from random, typified by it being well within a nice looking bell curve distribution.
So there it is, the procedure I intend to follow to objectively test the accuracy of the Delta Phenomenon.
Friday, 27 January 2012
Testing the Delta Phenomenon, Part 2
The main "problem" with the Delta Phenomenon is that it can be very subjective. This assertion is borne out by this quote - "We have observed over time that market participants using Delta and applying it differently do not make buy & sell decisions en-masse at the same time. Also, while the Delta order is perfect, our interpretation of the order is never 100% accurate." which is taken from the FAQ section of the above mentioned message board. Also here are some random comments that I have cut and pasted from the above linked forum thread (views expressed are those of the original forum posters, not mine)
- ...as the delta phenomenom could be really nothing but a recurring coincidence that no one can actually use with a good accuracy rate.
- ...while there may be something in the theory, finding a practical application for it is impossible. To make sense of it everything has to be viewed in hindsight...
- I thought it was interesting, but just like drawing lines on a chart, I could find a turning point any sequence I decided to use. I think we humans are extraordinary at seeking out patterns in the world around us. Maybe not such a good thing though if there really is no pattern, just perception.
- Like most any indicator, it looks great in hindsight as you can apply each turning point to the nearest high/low. Problem is, of course, you never know this in real time, only after the fact.
- Trading with Delta is a lot like looking at an MT4 indicator that "repaints".
- Mind you, I'm not saying the concept behind Delta is not valid. Just that because of the latitude afforded on either side of the "turning point" day, plus the idea of inversions......it's just real tough to be sure until after the fact. Even the monthly newsletter that Wilder sent out with turning point calls was often subject to correction and change.....a lot of subjectivity with this method.
- Much is left to the trader's own subjective judgment.
Resampling
Monte Carlo Methods
Bootstrapping
Statistical Hypothesis Testing
Testing the Delta Phenomenon
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 Octave, R 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 this1995-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 thusMany 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.
Sunday, 5 February 2012
Testing the Delta Phenomenon, Part 4
The projected solution turning points are indicated with white, red and green lines whilst the actual turns in the data are indicated in yellow. For the purposes of counting, each yellow, actual turn line is paired with one of the projected solution lines, i.e. counting from the left the first yellow line is paired with the first white solution line, the second with the second, third with third etc. The red solution lines indicate the "inversion" window with the green line between them being the number 1 point of the solution.
The first thing that stands out is that the biggest move of the year occurred within the inversion time period window, and the Delta Phenomenon book states that big moves can be expected around the number 1 point (green line), so kudos to Delta so far. Another thing that stands out is that most turning points occur within one to two weeks of their projected turning points, well within the margin of error for this solution. Kudos plus 2 for Delta. However the litmus test is this
[1] 8.769231
[1] 278
[1] 0.0556
[1] 33.09061
[1] 14.54137
[1] 1.672565
the results of a 0.5 million random permutation test. With only 278 random permutations out of half a million matching or bettering the average error test statistic, giving a p-value of 0.0556%, the null hypothesis can be rejected at a statistical significance level of 0.1%. Delta has passed this first test.
Monday, 27 February 2012
Gold Delta Solution
It's interesting because, according to Delta, gold is now in its inversion window on the MTD time frame, shown on the above chart as the two solid red vertical lines. By my count the first MTD red line marks a low, the green MTD a high point 1, followed by a possible in-between-point (IBP), which if it holds means higher prices are indicated for the high MTD point 2 in about 2 to 4 weeks time. However, at that time the ITD time frame will be in its inversion window too. Delta theory says that big moves can be expected around point 1 (in the inversion time window) and here we have two inversion time windows coinciding! Gold may very well experience some big move or increased volatility in the next few weeks, maybe as a result of a flight to safety from the unusual recent run-up in the S&P (see the chart in my previous post) and its impending correction. I look forward to seeing how this interpretation based on Delta pans out.
Sunday, 19 February 2012
Testing the Delta Phenomenon, Part 6
5th March 2007 (4445) to 25th June 2007 (4557)
[1] 3.5
[1] 0.3668 %
3rd July 2007 (4565) to 22nd Oct 2007 (4676)
[1] 4.75
[1] 3.3826 %
30th Oct 2007 (4684) to 17th Feb 2008 (4794)
[1] 4.583333
[1] 2.8366 %
25th Feb 2008 (4802) to 14th June 2008 (4912)
[1] 3.916667
[1] 0.9496 %
17th June 2008 (4915) to 10th Oct 2008 (5030)
[1] 3.333333
[1] 0.2472 %
18th Oct 2008 (5038) to 5th Feb 2008 (5148)
[1] 3.333333
[1] 0.2524 %
9th Feb 2009 (5152) to 4th June 2009 (5267)
[1] 5.083333
[1] 4.7866 %
5th June 2009 (5268) to 2nd Oct 2009 (5387)
[1] 5.333333
[1] 5.477 %
9th Oct 2009 (5394) to 29th Jan 2010 (5506)
[1] 4.25
[1] 1.9232 %
2nd Feb 2010 (5510) to 4th June 2010 (5632)
[1] 4.230769
[1] 1.306 %
31st May 2010 (5628) to 21st Sept 2010 (5741)
[1] 3.583333
[1] 0.4506 %
23rd Sept 2010 (5743) to 28th Jan 2011 (5870)
[1] 3.538462
[1] 0.2902 %
31st Jan 2011 (5873) to 21st May 2011 (5983)
[1] 4.083333
[1] 1.5542 %
31st May 2011 (5993) to 8th Sept 2011 (6093)
[1] 4.272727
[1] 3.2184 %
16th Sept 2011 (6101) to 2nd Feb 2012 (6240)
[1] 3.928571
[1] 0.2982 %
As can be seen, if 5 % is taken as the level of statistical significance only one test fails to reject the null hypothesis, and 7 out of the remaining 14 are statistically significant at the 1 % level. I also repeated the test several times over the entire data period, encompassing a total of 182 separate turning points with a test statistic of 4.10989 for these 182 points. In these repeated tests (for a total of a few million permutations) not a single permutation was as good as or better than the given test statistic! These results are much better than I had anticipated and I therefore consider Delta to have passed these tests as well. For readers' interest a plot of the period in which the null hypothesis is not rejected is shown below,
where the yellow lines are my identification of the actual turning points and the white, red and green lines are the turning points projected by the solution.
Now that this series of tests is complete, and Delta has passed them, what does it all mean? This will be the subject of my next post.
Friday, 29 July 2011
Update on the Trend Vigor indicator
Firstly, I create my "ideal" market types by using an Octave .oct function, the code for which is given below. In this code can also be seen my implementation of the code for the Trend Vigor indicator, which I think is slightly different from Elher's. The code is commented, so no further description is required here.
.oct function code
// This function takes as arguments a single value input for a sinewave period and a single value input for a degrees increment value. A sinewave of the
// given period is created and then trends are added to the sinewave such that 5 hypothesised "ideal" market types are created: these are
// 1) a perfectly cyclic sideways market with no trend i.e. just the sinewave component
// 2) an uptrending market with cyclic retracements (uwr) such that retracements are down to the 50% Fibonacci ratio level
// 3) an uptrending market with no retracemnets (unr) i.e. the uptrend completely swamps the downward cyclic component counter to the trend
// 4) a downtrending market with cyclic retracements (dwr) such that retracements are up to the 50% Fibonacci ratio level
// 5) a downtrending market with no retracements (dnr) i.e. the downtrend completely swamps the upward cyclic component counter to the trend
// The cyclic component of these markets is then extracted using the bandpass indicator. These vector lengths are 500 in order to allow time for the
// bandpass calculations to settle down. The peak to peak amplitude is then recovered from the bandpass at the end of each market vector. The trend slope
// is also calculated at the end of each market vector. The trend vigor indicator value is then calculated thus
// trend_slope/p_to_p. The bandpass delta is set to 0.2.
// The idea is that
// for a sideways market the value should be about 0
// for a trending with retracement market the value should be > 0 && < 1 or < 0 && > -1, depending on whether an up trend or down trend
// for a trending with no retracement market the ratio should be > 1 or > -1, depending on whether an up trend or down trend
// The original sinewave is then repeatedly phase shifted by the degrees increment value, the above process repeated and new values are calculated.
// All the calculated values for each market type are the vector outputs of the function.
#include
#include
#include
#define PI 3.14159265
DEFUN_DLD (trend_vigor_dist, args, , "Inputs are period & degrees increment value, outputs are vectors of repeated median slope values")
{
octave_value_list retval_list;
if (args(0).length () < 1 | args(0).length () > 1)
{
error ("Invalid arguments. Inputs are single value period length and single value degrees increment value");
return retval_list;
}
if (args(1).length () < 1 | args(1).length () > 1)
{
error ("Invalid arguments. Inputs are single value period length and single value degrees increment value");
return retval_list;
}
if (error_state)
{
error ("Invalid arguments. Inputs are single value period length and single value degrees increment value");
return retval_list;
}
// read inputs
int period = args(0).int_value ();
double degrees_inc = args(1).double_value ();
double period_inc = 360.0 / double(period);
int length = period + 1; // the length of the "lookback". Is equal to period + 1 to encompass a full period
// vectors to hold created market values
ColumnVector sideways_vec(500); // vector to hold sideways market values
ColumnVector sideways_bandpass_vec(500); // vector to hold bandpass of sideways market values
ColumnVector uwr_vec(500); // vector to hold uwr market values
ColumnVector uwr_bandpass_vec(500); // vector to hold bandpass of uwr market values
ColumnVector unr_vec(500); // vector to hold unr market values
ColumnVector unr_bandpass_vec(500); // vector to hold bandpass of unr market values
ColumnVector dwr_vec(500); // vector to hold dwr market values
ColumnVector dwr_bandpass_vec(500); // vector to hold bandpass of dwr market values
ColumnVector dnr_vec(500); // vector to hold dnr market values
ColumnVector dnr_bandpass_vec(500); // vector to hold bandpass of dnr market values
// calculate the market trend_incs
double uwr_trend_inc = 12 / ( 5 * double(period) );
double unr_trend_inc = 4 / double(period);
double dwr_trend_inc = -( 12 / ( 5 * double(period) ) );
double dnr_trend_inc = -( 4 / double(period) );
// declare variables for bandpass and trend vigor calculations
double delta = 0.2;
double beta = cos( (360.0/period)*PI/180.0 );
double gamma = 1.0 / cos( (720.0*delta/period)*PI/180.0 );
double alpha = gamma - sqrt(gamma*gamma - 1.0);
double power_side;
double power_uwr;
double power_unr;
double power_dwr;
double power_dnr;
double rms;
double p_to_p;
// create output vectors
int output_vec_length = int ( 360 / degrees_inc );
ColumnVector sideways_dist ( output_vec_length ); // create output column of correct length for sideways market
ColumnVector uwr_dist ( output_vec_length ); // create output column of correct length for uwr market
ColumnVector unr_dist ( output_vec_length ); // create output column of correct length for unr market
ColumnVector dwr_dist ( output_vec_length ); // create output column of correct length for dwr market
ColumnVector dnr_dist ( output_vec_length ); // create output column of correct length for dnr market
for (octave_idx_type ii (0); ii < output_vec_length; ii++)
{
// Create the market types and their bandpasses for this ii iteration
for (octave_idx_type jj (0); jj < 500; jj++)
{
// First create the sideways market type
sideways_vec(jj) = sin( (degrees_inc*ii + period_inc*jj) * PI / 180 );
if ( jj < 2 )
{
sideways_bandpass_vec(jj) = 0;
}
else
{
sideways_bandpass_vec(jj) = 0.5*(1.0 - alpha)*(sideways_vec(jj) - sideways_vec(jj-2)) + beta*(1.0 + alpha)*sideways_bandpass_vec(jj-1) - alpha*sideways_bandpass_vec(jj-2);
}
// next, the uwr retracement market (uwr)
uwr_vec(jj) = sideways_vec(jj) + jj*uwr_trend_inc;
if ( jj < 2 )
{
uwr_bandpass_vec(jj) = 0;
}
else
{
uwr_bandpass_vec(jj) = 0.5*(1.0 - alpha)*(uwr_vec(jj) - uwr_vec(jj-2)) + beta*(1.0 + alpha)*uwr_bandpass_vec(jj-1) - alpha*uwr_bandpass_vec(jj-2);
}
// next, the unr retracement market (unr)
unr_vec(jj) = sideways_vec(jj) + jj*unr_trend_inc;
if ( jj < 2 )
{
unr_bandpass_vec(jj) = 0;
}
else
{
unr_bandpass_vec(jj) = 0.5*(1.0 - alpha)*(unr_vec(jj) - unr_vec(jj-2)) + beta*(1.0 + alpha)*unr_bandpass_vec(jj-1) - alpha*unr_bandpass_vec(jj-2);
}
// next, the dwr retracement market (dwr)
dwr_vec(jj) = sideways_vec(jj) + jj*dwr_trend_inc;
if ( jj < 2 )
{
dwr_bandpass_vec(jj) = 0;
}
else
{
dwr_bandpass_vec(jj) = 0.5*(1.0 - alpha)*(dwr_vec(jj) - dwr_vec(jj-2)) + beta*(1.0 + alpha)*dwr_bandpass_vec(jj-1) - alpha*dwr_bandpass_vec(jj-2);
}
// next, the dnr retracement market (dnr)
dnr_vec(jj) = sideways_vec(jj) + jj*dnr_trend_inc;
if ( jj < 2 )
{
dnr_bandpass_vec(jj) = 0;
}
else
{
dnr_bandpass_vec(jj) = 0.5*(1.0 - alpha)*(dnr_vec(jj) - dnr_vec(jj-2)) + beta*(1.0 + alpha)*dnr_bandpass_vec(jj-1) - alpha*dnr_bandpass_vec(jj-2);
}
} // end of jj loop to create the different markets and their bandpasses
// now loop over end of each market vector to create the distributions
power_side = 0.0;
power_uwr = 0.0;
power_unr = 0.0;
power_dwr = 0.0;
power_dnr = 0.0;
for (octave_idx_type jj (0); jj < length; jj++)
{
power_side = power_side + sideways_bandpass_vec(499-jj)*sideways_bandpass_vec(499-jj) + sideways_bandpass_vec(499-jj-int(period/4.0))*sideways_bandpass_vec(499-jj-int(period/4.0)) ;
power_uwr = power_uwr + uwr_bandpass_vec(499-jj)*uwr_bandpass_vec(499-jj) + uwr_bandpass_vec(499-jj-int(period/4.0))*uwr_bandpass_vec(499-jj-int(period/4.0)) ;
power_unr = power_unr + unr_bandpass_vec(499-jj)*unr_bandpass_vec(499-jj) + unr_bandpass_vec(499-jj-int(period/4.0))*unr_bandpass_vec(499-jj-int(period/4.0)) ;
power_dwr = power_dwr + dwr_bandpass_vec(499-jj)*dwr_bandpass_vec(499-jj) + dwr_bandpass_vec(499-jj-int(period/4.0))*dwr_bandpass_vec(499-jj-int(period/4.0)) ;
power_dnr = power_dnr + dnr_bandpass_vec(499-jj)*dnr_bandpass_vec(499-jj) + dnr_bandpass_vec(499-jj-int(period/4.0))*dnr_bandpass_vec(499-jj-int(period/4.0)) ;
}
// fill the distribution vectors
rms = sqrt( power_side / (period+1) ) ;
p_to_p = 2.0 * 1.414 * rms ;
sideways_dist(ii) = ( sideways_vec(499) - sideways_vec(499-period) ) / p_to_p ;
rms = sqrt( power_uwr / (period+1) ) ;
p_to_p = 2.0 * 1.414 * rms ;
uwr_dist(ii) = ( uwr_vec(499) - uwr_vec(499-period) ) / p_to_p ;
rms = sqrt( power_unr / (period+1) ) ;
p_to_p = 2.0 * 1.414 * rms ;
unr_dist(ii) = ( unr_vec(499) - unr_vec(499-period) ) / p_to_p ;
rms = sqrt( power_dwr / (period+1) ) ;
p_to_p = 2.0 * 1.414 * rms ;
dwr_dist(ii) = ( dwr_vec(499) - dwr_vec(499-period) ) / p_to_p ;
rms = sqrt( power_dnr / (period+1) ) ;
p_to_p = 2.0 * 1.414 * rms ;
dnr_dist(ii) = ( dnr_vec(499) - dnr_vec(499-period) ) / p_to_p ;
} // end of main ii loop
retval_list(4) = dnr_dist;
retval_list(3) = dwr_dist;
retval_list(2) = unr_dist;
retval_list(1) = uwr_dist;
retval_list(0) = sideways_dist;
return retval_list;
}
This function is called by this simple Octave scriptclear all
inc = input( "Enter phase increment: ");
for ii = 6:50
[sideways_dist,uwr_dist,unr_dist,dwr_dist,dnr_dist] = trend_vigor_dist(ii,inc);
A=[sideways_dist,uwr_dist,unr_dist,dwr_dist,dnr_dist];
file = strcat( int2str(ii),"_period_dist" );
dlmwrite(file,A)
endfor
which writes the output of the tests to named files which are to be used for further analysis in R.Firstly, using R, I wanted to see what the distribution of the results looks like, so this R script
rm(list=ls())
data <- as.matrix(read.csv(file="20_period_dist",head=FALSE,sep=,))
side <- density(data[,1])
uwr <- density(data[,2])
max_uwr_y <- max(uwr$y)
max_uwr_x <- max(uwr$x)
min_uwr_x <- min(uwr$x)
unr <- density(data[,3])
max_unr_y <- max(unr$y)
max_unr_x <- max(unr$x)
min_unr_x <- min(unr$x)
dwr <- density(data[,4])
max_dwr_y <- max(dwr$y)
max_dwr_x <- max(dwr$x)
min_dwr_x <- min(dwr$x)
dnr <- density(data[,5])
max_dnr_y <- max(dnr$y)
max_dnr_x <- max(dnr$x)
min_dnr_x <- min(dnr$x)
plot_max_y <- max(max_uwr_y,max_unr_y,max_dwr_y,max_dnr_y)
plot_max_x <- max(max_uwr_x,max_unr_x,max_dwr_x,max_dnr_x)
plot_min_x <- min(min_uwr_x,min_unr_x,min_dwr_x,min_dnr_x)
par(mfrow=c(2,1))
plot(uwr,xlim=c(plot_min_x,plot_max_x),ylim=c(0,plot_max_y),col="red")
par(new=TRUE)
plot(unr,xlim=c(plot_min_x,plot_max_x),ylim=c(0,plot_max_y),col="blue")
par(new=TRUE)
plot(dwr,xlim=c(plot_min_x,plot_max_x),ylim=c(0,plot_max_y),col="green")
par(new=TRUE)
plot(dnr,xlim=c(plot_min_x,plot_max_x),ylim=c(0,plot_max_y),col="black")
plot(side)
produces plots such as this outputwhere the top plot shows the distributions of the uwr, unr, dwr and dnr markets, and the lower plot the sideways market. In this particular case, the spread of each distribution is so narrow ( measured differences of the order of thousandths of a decimal place ) that I consider that for practical purposes the distributions can be treated as single values. This simple R boot strap script gets the average value of the distributions to be used as this single point value.
rm(list=ls())
data <- as.matrix(read.csv(file="trend_vigor_dist_results",head=FALSE,sep=,))
side <- data[,1]
uwr <- data[,2]
unr <- data[,3]
dwr <- data[,4]
dnr <- data[,5]
side_samples <- matrix(0,50000)
uwr_samples <- matrix(0,50000)
unr_samples <- matrix(0,50000)
dwr_samples <- matrix(0,50000)
dnr_samples <- matrix(0,50000)
for(ii in 1:50000) {
side_samples[ii] <- mean(sample(side, replace=T))
uwr_samples[ii] <- mean(sample(uwr, replace=T))
unr_samples[ii] <- mean(sample(unr, replace=T))
dwr_samples[ii] <- mean(sample(dwr, replace=T))
dnr_samples[ii] <- mean(sample(dnr, replace=T))
}
side_mean <- mean(side_samples)
uwr_mean <- mean(uwr_samples)
unr_mean <- mean(unr_samples)
dwr_mean <- mean(dwr_samples)
dnr_mean <- mean(dnr_samples)
For readers' interest, the actual values are 0.829, 1.329, -0.829 and -1.329 with 0 for the sideways market.This final plot is the same Natural Gas plot as in my previous post, but with the above values substituted for Ehler's default values of 1 and -1.
What I intend to do now is use these values as the means of normal distributions with varying standard deviations as inputs for my Naive Bayes classifier. Further Monte Carlo testing will be done such that values for the standard deviations are obtained that result in the classifier giving false classifications, when tested using the "ideal" markets code above, within acceptable limits, most probably a 5% classification error rate.
Saturday, 30 January 2021
Temporal Clustering Times on Forex Majors Pairs
In the following code box there are the results from the temporal clustering routine of my last few posts on the four forex majors pairs of EUR_USD, GBP_USD, USD_CHF and USD_JPY.
###### EUR_USD 10 minute bars #######
## In the following order
## Both Delta turning point filter and "normal" TPF combined ##
## Delta turning point filter only ##
## "Normal" turning point filter only
###################### Monday ##############################################
K_opt == 8, ix values == 13 38 63 89 112 135 162 186 ## averaged over all 15 n_bars 1 to 15 inclusive
00 4:10 8:20 12:40 16:30 20:20 00:50 4:50
K_opt == 8, ix values == 13 39 64 89 112 135 161 186 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt == 5, ix_values == 21 60 97 134 175 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
K == 6, ix values == 21 59 94 125 158 184
K_opt == 11, ix values == 9 26 43 60 78 95 113 132 151 169 185 ## averaged over all 15 n_bars 1 to 15 inclusive
K_opt == 8, ix values == 13 36 61 86 111 136 161 186 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt == 8, ix values == 13 34 61 87 110 137 164 187 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
K_opt == 8, ix values == 13 38 63 88 112 137 162 186 ## averaged over all 15 n_bars 1 to 15 inclusive
K_opt == 10, ix values == 10 31 52 72 91 112 131 150 169 188 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt == 8, ix values == 12 35 62 88 112 137 164 187 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
###################### Tuesday #############################################
K_opt == 6, ix values == 131 169 206 244 283 322 ## averaged over all 15 n_bars 1 to 15 inclusive
19:40 02:00 8:10 14:30 21:00 03:30
K_opt == 6, ix values == 131 170 207 245 284 323 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt == 7, ix values == 131 168 206 243 274 305 330 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
K_opt == 11, ix values == 124 143 164 184 205 226 247 268 289 310 331 ## averaged over all 15 n_bars 1 to 15 inclusive
K_opt == 11, ix values == 124 144 164 185 204 225 246 267 288 309 332 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt == 7, ix values = 133 169 206 241 273 304 329 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
K_opt == 9, ix values == 127 152 175 202 228 253 278 305 330 ## averaged over all 15 n_bars 1 to 15 inclusive
K_opt == 9, ix values == 127 152 177 202 228 253 278 304 329 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt == 7, ix values == 132 168 205 242 273 304 329 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
###################### Wednesday ###########################################
K_opt == 6, ix values == 275 312 351 389 426 465 ## averaged over all 15 n_bars 1 to 15 inclusive
19:40 01:50 08:20 14:40 20:50 03:20
K_opt == 6, ix values == 275 313 352 391 428 466 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt == 6, ix values == 274 312 350 389 424 463 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
K_opt == 9, ix values == 272 299 322 347 372 397 422 449 474 ## averaged over all 15 n_bars 1 to 15 inclusive
K_opt == 11, ix values == 268 288 308 329 348 369 390 411 432 453 476 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt == 6, ix values == 275 312 351 388 424 463 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
K_opt == 9, ix values == 272 297 322 348 373 398 423 449 474 ## averaged over all 15 n_bars 1 to 15 inclusive
K_opt == 9, ix values == 271 297 322 348 373 398 423 448 473 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt == 6, ix values == 276 311 350 389 426 465 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
####################### Thursday ###########################################
K_opt == 6, ix values == 420 457 495 532 570 609 ## averaged over all 15 n_bars 1 to 15 inclusive
19:50 02:00 08:20 14:30 20:50 03:20
K_opt == 6, ix values == 420 457 494 531 570 610 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt == 6, ix values == 420 457 495 532 568 607 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
K_opt == 9, ix values == 416 443 466 492 518 543 568 593 618 ## averaged over all 15 n_bars 1 to 15 inclusive
K_opt == 10, ix values == 414 437 460 483 506 527 550 573 596 619 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt == 9, ix values == 416 443 466 493 520 543 568 595 618 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
K_opt == 9, ix values == 415 440 465 492 518 543 568 593 618 ## averaged over all 15 n_bars 1 to 15 inclusive
K_opt == 9, ix values == 415 440 465 492 518 543 568 593 618 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt == 7, ix values == 420 457 494 529 561 592 617 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
####################### Friday #############################################
K_opt == 5, ix values == 564 599 635 670 703 ## averaged over all 15 n_bars 1 to 15 inclusive
19:50 01:40 07:40 13:30 19:00
K_opt == 6, ix values == 563 596 627 654 680 707 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K == 5, ix values == 564 599 635 668 703
K_opt == 5, ix values == 564 601 639 674 705 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
K_opt == 9, ix values == 556 575 595 614 633 652 672 691 711 ## averaged over all 15 n_bars 1 to 15 inclusive
K_opt == 11, ix values == 554 570 587 602 619 634 651 667 682 698 713 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt == 9, ix values == 556 575 595 614 633 652 671 691 711 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt == 9, ix values == 556 575 596 613 634 652 672 691 711 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
K_opt == 9, ix values == 556 575 594 613 633 652 672 691 710 ## averaged over all 15 n_bars 1 to 15 inclusive
K_opt == 9, ix values == 556 575 594 613 634 653 672 691 710 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt == 5, ix values == 564 600 637 674 705 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
############################################################################
###### GBP_USD 10 minute bars #######
## In the following order
## Both Delta turning point filter and "normal" TPF combined ##
###################### Monday ##############################################
K_opt = 8, ix_values = 13 36 61 86 111 136 162 186 ## averaged over all 15 n_bars 1 to 15 inclusive
0:00 3:50 8:00 12:10 16:20 20:30 0:50 4:50
K_opt = 9, ix_values = 12 34 56 78 99 120 141 164 187 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt = 8, ix_values = 12 35 61 86 110 136 163 186 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
###################### Tuesday #############################################
K_opt = 12, ix_values = 124 143 162 180 199 216 235 254 274 293 312 332 ## averaged over all 15 n_bars 1 to 15 inclusive
18:30 21:40 0:50 3:50 7:00 9:50 13:00 16:10 19:30 22:40 1:50 5:10
K_opt = 11, ix_values = 124 143 164 185 206 227 248 269 290 311 332 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt = 9, ix_values = 128 154 177 205 230 254 279 307 330 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
###################### Wednesday ###########################################
K_opt = 11, ix_values = 269 290 311 331 352 373 394 415 434 455 476 ## averaged over all 15 n_bars 1 to 15 inclusive
18:40 22:10 1:40 5:00 8:30 12:00 15:30 19:00 22:10 1:40 5:10
K_opt = 11, ix_values = 269 289 310 330 351 372 393 413 434 455 476 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt = 8, ix_values = 275 310 341 367 394 422 451 475 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
###################### Thursday ############################################
K_opt = 9, ix_values = 415 440 465 492 517 542 568 594 618 ## averaged over all 15 n_bars 1 to 15 inclusive
19:00 23:10 3:20 7:50 12:00 16:10 20:30 0:50 4:50
K_opt = 9, ix_values = 415 440 465 491 517 542 568 593 618 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt = 9, ix_values = 416 441 464 492 519 542 569 596 619 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
###################### Friday ##############################################
K_opt = 9, ix_values = 557 576 595 614 633 652 671 690 711 ## averaged over all 15 n_bars 1 to 15 inclusive
18:40 21:50 1:00 4:10 7:20 10:30 13:40 16:50 20:20
K_opt = 9, ix_values = 557 576 595 614 633 652 671 691 711 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt = 8, ix_values = 557 576 599 621 642 665 686 709 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
############################################################################
###### USD_CHF 10 minute bars #######
## In the following order
## Both Delta turning point filter and "normal" TPF combined ##
###################### Monday ##############################################
K_opt = 11, ix_values = 8 25 42 61 79 96 113 131 150 169 188 ## averaged over all 15 n_bars 1 to 15 inclusive
23:10 2:00 4:50 8:00 11:00 13:50 16:40 19:40 22:50 2:00 5:10
K_opt = 11, ix_values = 9 26 43 60 79 96 114 133 151 170 189 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt = 7, ix_values = 13 38 66 99 127 157 184 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
###################### Tuesday #############################################
K_opt = 9, ix_values = 127 152 177 202 228 253 279 306 330 ## averaged over all 15 n_bars 1 to 15 inclusive
19:00 23:10 3:20 7:30 11:50 16:00 20:20 0:50 4:50
K_opt = 11, ix_values = 124 144 165 185 204 225 246 267 288 309 331 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt = 7, ix_values = 133 170 205 240 270 301 328 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
###################### Wednesday ###########################################
K_opt = 10, ix_values = 270 293 316 342 365 388 411 432 454 475 ## averaged over all 15 n_bars 1 to 15 inclusive
18:50 22:40 2:30 6:50 10:40 14:30 18:20 21:50 1:30 5:00
K_opt = 12, ix_values = 268 287 308 327 346 365 384 401 420 439 458 477 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt = 7, ix_values = 276 313 349 383 414 444 471 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
###################### Thursday ############################################
K_opt = 11, ix_values = 413 432 452 471 491 512 533 554 575 598 619 ## averaged over all 15 n_bars 1 to 15 inclusive
18:40 21:50 1:10 4:20 7:40 11:10 14:40 18:10 21:40 1:30 5:00
K_opt = 12, ix_values = 412 431 450 469 488 507 526 545 563 582 601 621 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt = 9, ix_values = 415 440 463 491 518 543 570 597 619 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
###################### Friday ##############################################
K_opt = 9, ix_values = 557 576 596 615 634 653 672 691 710 ## averaged over all 15 n_bars 1 to 15 inclusive
18:40 21:50 1:10 4:20 7:30 10:40 13:50 17:00 20:10
K_opt = 9, ix_values = 556 575 595 614 633 652 671 690 710 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt = 7, ix_values = 558 579 602 629 652 677 705 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
############################################################################
###### USD_JPY 10 minute bars #######
## In the following order
## Both Delta turning point filter and "normal" TPF combined ##
###################### Monday ##############################################
K_opt = 12, ix_values = 8 24 41 58 73 90 107 124 141 158 173 190 ## averaged over all 15 n_bars 1 to 15 inclusive
23:10 1:50 4:40 7:30 10:00 12:50 15:40 18:30 21:20 0:10 2:40 5:30
K_opt = 12, ix_values = 8 24 41 56 73 90 107 124 141 158 173 190 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt = 5, ix_values = 20 60 99 136 175 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
###################### Tuesday #############################################
K_opt = 9, ix_values = 128 154 179 204 229 254 279 306 331 ## averaged over all 15 n_bars 1 to 15 inclusive
19:10 23:30 3:40 7:50 12:00 16:10 20:20 0:50 5:00
K_opt = 9, ix_values = 128 153 178 203 228 254 279 305 330 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt = 7, ix_values = 133 168 205 240 271 302 329 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
###################### Wednesday ###########################################
K_opt = 11, ix_values = 269 289 310 331 352 373 394 414 433 454 476 ## averaged over all 15 n_bars 1 to 15 inclusive
18:40 22:00 1:30 5:00 8:30 12:00 15:30 18:50 22:00 1:30 5:10
K_opt = 9, ix_values = 272 297 322 348 374 399 424 449 474 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt = 10, ix_values = 269 288 309 331 352 376 398 423 450 475 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
###################### Thursday ############################################
K_opt = 9, ix_values = 416 442 467 492 518 543 568 593 618 ## averaged over all 15 n_bars 1 to 15 inclusive
19:10 23:30 3:40 7:50 12:10 16:20 20:30 0:40 4:50
K_opt = 12, ix_values = 412 431 450 469 488 507 526 545 564 583 602 621 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt = 7, ix_values = 420 455 492 527 560 591 618 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
###################### Friday ##############################################
K_opt = 7, 8 or 9
ix_values 7 = 561 588 613 638 663 686 709 ## averaged over all 15 n_bars 1 to 15 inclusive
ix_values 8 = 557 578 599 622 643 666 687 710
ix_values 9 = 557 576 596 616 635 653 672 691 711 ## timings are for this bottom row
18:40 21:50 1:10 4:30 7:40 10:40 13:50 17:00 20:20
K_opt = 8, ix_values = 558 579 600 621 644 665 687 709 ## averaged over n_bars 1 to 6 inclusive ( upto and include 1 hour )
K_opt = 6, ix_values = 563 594 621 646 676 705 ## averaged over n_bars 7 to 15 inclusive ( over 1 hour )
############################################################################This is based on 10 minute bars over the last year or so. Readers should read my last few previous posts for background.
The first set of results, EUR_USD, are what the charts of my previous posts were based on and include combined results of my "Delta Turning Point Filter" and "Normal Turning Point Filter" and the results for each filter separately. Since there doesn't appear to be significant differences between these, the other three pairs' results are the combined filter results only.
The K_opt variable is the optimal number of clusters (see my temporal-clustering-part-3 post for how "optimal" is decided) and the ix_values are also described in this post. For convenience the first set of ix_values per day have the relevant times anotated underneath and therefore it is a simple matter to count forwards/backwards in 10 minute increments to place times to the other ix_values. The variable n_bars is an input to the turning point filter functions and essentially indicates the lookback/lookforward period (n_bar == 2 would mean 2 x 10 minute periods) used for determining a local high/low according to each function's logic.
As to how to interpret this, a typical sequence of times per day might look like this:
18:40 22:00 1:30 5:00 8:30 12:00 15:30 18:50 22:00 1:30 5:10
where the highlighted times represent the BST times for the period covering the London session open to the New York session close for one day. The preceding and following times are the two "book-ending" Asian sessions.
Close inspection of these results reveals some surprising regularities. In even just the above single example (an actual copy and paste of a code box example) there appear to be definite times per day at which a local high/low occurs. I hopefully will be able to incorporate this into some type of chart for a nice visual presentation of the data.
More in due course. Enjoy.
Monday, 28 May 2012
Update on Gold Delta Solution
Below is an updated chart of Gold showing the next few turning points following on from where the previous chart ended.
My read of this chart is that MTD 1 (solid green line) is a high which came in early March, followed by a low MTD 2 which came in a week late (the second solid red line). The market then struggled to move up to a high MTD 3, indicating a very weak market, and now the market is moving down to a low MTD 4 (the second solid yellow line). Since the most recent MTD 3 high is lower than the previous MTD 1 high and it looks like the MTD 4 low will come in lower that the MTD 2 low, I'd say that on the MTD time frame Gold is now a bear market.
On an unrelated note, it has been more than a month since my last post. During this time I have been busy working through the free, online version of Andrew Ng's Machine Learning course. More on this in a future post.
Wednesday, 15 February 2012
Delta Training Starts Tomorrow!
...there is a Delta Phenomenon multi-week training session via IRC by a group of DP experts!.. The classes are Monday, Wednesday and Friday 2-5pm EST starting on February 16th. irc server irc.forex.com #training
I shall try to attend these training sessions, and look forward to them.
Wednesday, 2 October 2013
Review of Quantpedia.com
The thing that strikes me is the wide and eclectic range of studies available, which can be easily seen by clicking on the keywords dropdown box on the screener tab. There is something for almost all trading styles and I would be surprised if a visitor to the premium section found nothing of value. As always the devil is in the details, and of course I haven't read anywhere near all the information that is available in the various studies, but a brief visual overview of the various studies' performance taken from Quantpedia's screener page is provided in the scatterchart below:
(n.b. Those studies that lie exactly on the y-axis (volatility = 0%) do not have no volatility, but no % figure for volatility was given in the screener.)
As can be seen there are some impressive performers, which are fully disclosed in the premium section. A few studies, irrespective of performance, did catch my immediate attention. Firstly, there are a couple of studies that look at lunar effects in stock markets and precious metals. Long time readers of this blog will know that some time back I did a series of tests on the Delta Phenomenon, which is basically a lunar effect model of price turning points. The conclusion of the tests I conducted was that the Delta Phenomenon did have statistically significant predictive ability. The above mentioned studies come to the same conclusion with regard to lunar effects, although via a different testing methodology. It is comforting to have one's own research conclusions confirmed by independent researchers, and it's a valuable resource to be able to see how other researchers approach the testing of similar market hypotheses.
Secondly, there is a study on using Principal Components Analysis to characterise the current state of the market. This is very much in tune with what I'm working on at the moment with my neural net market classifier, and the idea of using PCA as an input is a new one to me and one that I shall almost certainly look into in more detail.
This second point I think neatly sums up the main value of the studies on the Quantpedia site - they can give you new insights as to how one might develop one's own trading system(s), with the added benefit that the idea is not a dead end because it has already been tested by the original paper's authors and the Quantpedia site people. You could also theoretically just take one of the studies as a stand alone system and tweak it to suit your own needs, or add it as a form of diversification to an existing set of trading systems. Given the wide range of studies available, this would be a much more robust form of diversification than merely adjusting a look back length, parameter value or some other such cosmetic adjustment.
In conclusion, since I can appreciate the value in the Quantpedia site, I would like to thank Martin Nizny of Quantpedia for extending me the opportunity to review the premium section of the site.
Tuesday, 24 November 2020
Temporal Clustering on Real Prices
Having now had time to run the code shown in my previous post, Temporal Clustering, part 3, in this post I want to show the results on real prices.
Firstly, I have written two functions in Octave to identify market turning points and each function takes as input an n_bar argument which determines the lookback/lookforward length along price series to determine local relative highs and lows. I ran both these for n_bar values of 1 to 15 inclusive on EUR_USD forex 10 minute bars from July 2012 upto and including last week's set of 10 minute bars. I created 3 sets of turning point data per function by averaging the function outputs over n_bar 1 - 15, 1 - 6 and 7 - 15, and also averaged the outputs over the average of the 2 functions over the same ranges. In total this gives 9 slightly different sets of turning point data.
I then ran the optimal K clustering code, shown in previous posts, over each set of data to get the "solutions" per set of data. Six of the sets had an optimal K value of 8 and a combined plot of these is shown below.
For each "solution" turning point ix (ix ranges from 1 to 198) a turning point value of 1 is added to get a sort of spike train plot through time. The ix = 1 value is 22:00 BST on Sunday and ix = 198 is 06:50 BST on Tuesday. I chose this range so that there would be a buffer at each end of the time range I am really interested in: 7:00 BST to 22:00 BST, which covers the time from the London open to the New York close. The vertical blue lines are plotted for clarity to help identify the the turns and are plotted as 3 consecutive lines 10 minutes apart. The added text shows the time of occurence of the first bar of each triplet of lines, the time being London BST. The following second plot is the same as above but with the other 3 "solutions" of K = 5, 10 and 11 added.For those readers who are familiar with the Delta Phenomenon the main vertical blue lines could conceptually be thought of as MTD lines with the other lines being lower timeframe ITD lines, but on an intraday scale. However, it is important to bear in mind that this is NOT a Delta solution and therefore rules about numbering, alternating highs and lows and inversions etc. do not apply. It is more helpful to think in terms of probability and see the various spikes/lines as indicating times of the day at which there is a higher probability of price making a local high or low. The size of a move after such a high or low is not indicated, and the timings are only approximate or alternatively represent the centre of a window in which the high or low might occur.The proof of the pudding is in the eating, however, and the following plots are yesterday's (23 November 2020) out of sample EUR_USD forex pair price action with the lines of the above "solution" overlaid. The first plot is just the K = 8 solution plot
whilst this second plot has all lines shown.Given the above caveats about caution with regards to the lines only being probabilities, it seems uncanny how accurately the major highs and lows of the day are picked out. I only wish I had done this analysis sooner as then yesterday could have been one of my best trading days ever!More soon.
Sunday, 19 February 2012
I've Missed the Delta Training :'-(
Perhaps the anonymous reader who alerted me to the course in the first place has more information? I would be grateful to hear from you again.
Monday, 13 February 2012
Testing the Delta Phenomenon, Part 5
Years 2007 to 2008 - 16th Aug 2007 to 6th Jan 2009 incl.
[1] 9.722222 - test statistic, the average error of identified turns, in days
[1] 103 - number of random permutations as good as or better than test stat
[1] 0.0206 - above expressed as % of permutations, the p-value
[1] 37.6425 - average of permutation error distribution, in days
[1] 16.32295 - standard deviation of permutation error distribution, in days
[1] 1.710493 - test statistic distance from average of permutation error distribution, expressed as a multiple of standard deviation of permutation error distribution
Year 2009 - 20th Dec 2008 to 16th Jan 2010 incl.
[1] 10.42857 - test statistic, the average error of identified turns, in days
[1] 809 - number of random permutations as good as or better than test stat
[1] 0.1618 - above expressed as % of permutations, the p-value
[1] 33.43553 - average of permutation error distribution, in days
[1] 14.18181 - standard deviation of permutation error distribution, in days
[1] 1.622286 - test statistic distance from average of permutation error distribution, expressed as a multiple of standard deviation of permutation error distribution
Year 2010 - 27th Nov 2009 to 28th Jan 2011 incl.
[1] 11.35714 - test statistic, the average error of identified turns, in days
[1] 1706 - number of random permutations as good as or better than test stat
[1] 0.3412 - above expressed as % of permutations, the p-value
[1] 34.92687 - average of permutation error distribution, in days
[1] 15.36711 - standard deviation of permutation error distribution, in days
[1] 1.533777 - test statistic distance from average of permutation error distribution, expressed as a multiple of standard deviation of permutation error distribution
All out of sample data - 16th Aug 2007 to 28th Jan 2012
[1] 9.8 - test statistic, the average error of identified turns, in days
[1] 0 - number of random permutations as good as or better than test stat
[1] 0 - above expressed as % of permutations, the p-value
[1] 67.47298 - average of permutation error distribution, in days
[1] 29.38267 - standard deviation of permutation error distribution, in days
[1] 1.962823 - test statistic distance from average of permutation error distribution, expressed as a multiple of standard deviation of permutation error distribution
The histogram plot of the final, all out of sample data
I think the results are unambiguous - given the consistently low p-values the null hypothesis can be rejected and the alternative hypothesis accepted i.e. the S & P 500 Medium term solution is NOT random and therefore has some predictive value.
My next round of tests will be of the Intermediate term solution for the S & P 500. Before conducting these tests however, I would like to state my expectation that the results will not be as conclusive as those given above. This is due to the fact that the Intermediate term solution also has 12 points, but these 12 points occur within a time period that is approximately 120 days long, so the actual average error on this time frame will have to be 2 or 3 bars/days or better to match the above results.
Tuesday, 20 October 2020
A Temporal Clustering Function
Recently a reader contacted me with a view to collaborating on some work regarding the Delta phenomenon but after a brief exchange of e-mails this seems to have petered out. However, for my part, the work I have done has opened a few new avenues of investigation and this post is about one of them.
One of the problems I set out to solve was clustering in the time domain, or temporal clustering as I call it. Take a time series and record the time of occurance of an event by setting to 1, in an otherwise zero filled 1-dimensional vector the same length as the original time series, the value of the vector at time index tx and repeat for all occurances of the event. In my case the event(s) I am interested in are local highs and lows in the time series. This vector is then "chopped" into segments representing distinct periods of time, e.g. 1 day, 1 week etc. and stacked into a matrix where each row is one complete period and the columns represent the same time in each period, e.g. the first column is the first hour of the trading week, the second column the second hour etc. Sum down the columns to get a final 1-dimensional vector of counts of the timing of events happening within each period over the entire time series data record. A chart of such is shown below.
The coloured vertical lines show the opening and closing times of the London and New York sessions (7am to 5pm in their respective local times) for one complete week at a 10 minute bar time scale, in this case for the GBP_USD forex pair. This is what I want to cluster.The solution I have come up with is the Octave function in the code box below
## Copyright (C) 2020 dekalog
##
## This program is free software: you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program. If not, see
## .
## -*- texinfo -*-
## @deftypefn {} {@var{centre_ix} =} blurred_maxshift_1d_linear (@var{train_vec}, @var{bandwidth})
##
## Clusters the 1 dimensional vector TRAIN_VEC using a "centred" sliding window of length 2 * BANDWIDTH + 1.
##
## Based on the idea of the Blurred Meanshift Algorithm.
##
## The "centre ix" value of the odd length sliding window is assigned to the
## maximum value ix of the sliding window. The centre_ix, if it is not the
## maximum value, is then set to zero. A pass through the whole length of
## TRAIN_VEC is completed before any assignments are made.
##
## @seealso{}
## @end deftypefn
## Author: dekalog
## Created: 2020-10-17
function new_train_vec = blurred_maxshift_1d_linear ( train_vec , bandwidth )
if ( nargin < 2 )
bandwidth = 1 ;
endif
if ( numel( train_vec ) < 2 * bandwidth + 1 )
error( 'Bandwidth too wide for length of train_vec.' ) ;
endif
length_train_vec = numel( train_vec ) ;
assigned_cluster_centre_ix = zeros( size( train_vec ) ) ;
## initialise the while condition variable
has_converged = 0 ;
while ( has_converged < 1 )
new_train_vec = zeros( size( train_vec ) ) ;
## do the beginning and end of train_vec first
[ ~ , ix ] = max( train_vec( 1 : 2 * bandwidth + 1 ) ) ;
new_train_vec( ix ) = sum( train_vec( 1 : bandwidth ) ) ;
[ ~ , ix ] = max( train_vec( end - 2 * bandwidth : end ) ) ;
new_train_vec( end - 2 * bandwidth - 1 + ix ) = sum( train_vec( end - bandwidth + 1 : end ) ) ;
for ii = 2 * bandwidth + 1 : numel( train_vec ) - bandwidth
[ ~ , ix ] = max( train_vec( ii - bandwidth : ii + bandwidth ) ) ;
new_train_vec( ii - bandwidth - 1 + ix ) += train_vec( ii ) ;
endfor
if ( sum( ( train_vec == new_train_vec ) ) == length_train_vec )
has_converged = 1 ;
else
train_vec = new_train_vec ;
endif
endwhile
endfunction
I have named the function "blurred_maxshift_1d_linear" as it is inspired by the "blurred" version of the Mean shift algorithm, operates on a 1-dimensional vector and is "linear" in that there is no periodic wrapping of the data within the function code. The two function inputs are the above type of data, obviously, and an integer parameter "bandwidth" which controls the size of a moving window over the data in which the data is shifted according to a maximum value, hence maxshift rather than meanshift. I won't discuss the code further as it is pretty straightforward.
A chart of a typical clustering solution is (bandwidth setting == 2)
where the black line is the original count data and red the clustering solution. The bandwidth setting in this case is approximately equivalent to clustering with a 50 minute moving window.The following heatmap chart is a stacked version of the above where the bandwidth parameter is varied from 1 to 10 inclusive upwards, with the original data being at the lowest level per pane.
The intensity reflects the counts at each time tx index per bandwidth setting. The difference between the panes is that in the upper pane the raw data is the function input per bandwidth setting, whilst the lower pane shows hierarchical clustering whereby the output of the function is used as the input to the next function call with the next higher bandwidth parameter setting.More in due course.
Monday, 21 December 2015
John Ehler's Sinewave Indicator Code
DEFUN_DLD ( sinewave_indicator, args, nargout )
{
octave_value_list retval_list ;
int nargin = args.length () ;
int vec_length = args(0).length () ;
// check the input argument
if ( nargin != 1 )
{
error ("Invalid argument. Input is a single price vector.") ;
return retval_list ;
}
if ( vec_length < 50 )
{
error ("Invalid argument. Input is a single price vector.") ;
return retval_list ;
}
if ( error_state )
{
error ("Invalid argument. Input is a single price vector.") ;
return retval_list ;
}
// end of input checking
// inputs
ColumnVector price = args(0).column_vector_value () ;
// outputs
ColumnVector sinewave( vec_length ) ;
ColumnVector sinewave_lead_1( vec_length ) ;
ColumnVector smoothperiod_out( vec_length ) ;
ColumnVector dcphase_vec( vec_length ) ;
ColumnVector sumperiod( vec_length ) ;
ColumnVector sum_period( vec_length ) ;
ColumnVector deltaphase( vec_length ) ;
// Declarations for calculations of period, phase & sine wave measurements
ColumnVector smooth( vec_length ) ;
ColumnVector period( vec_length ) ;
ColumnVector smoothperiod( vec_length ) ;
ColumnVector detrender( vec_length ) ;
ColumnVector Q1( vec_length ) ;
ColumnVector I1( vec_length ) ;
ColumnVector jI( vec_length ) ;
ColumnVector jQ( vec_length ) ;
ColumnVector I2( vec_length ) ;
ColumnVector Q2( vec_length ) ;
ColumnVector sI2( vec_length ) ;
ColumnVector sQ2( vec_length ) ;
ColumnVector Re( vec_length ) ;
ColumnVector Im( vec_length ) ;
ColumnVector sRe( vec_length ) ;
ColumnVector sIm( vec_length ) ;
int dcperiod ;
double realpart ;
double imagpart ;
double dcphase ;
double sum_deltaphase ;
int count ;
// unrolled loop to fill the first 5 elements of above calculation vectors ( unrolled for speed optimisation )
sinewave(0) = 0.0 ; sinewave(1) = 0.0 ; sinewave(2) = 0.0 ; sinewave(3) = 0.0 ; sinewave(4) = 0.0 ;
sinewave_lead_1(0) = 0.0 ; sinewave_lead_1(1) = 0.0 ; sinewave_lead_1(2) = 0.0 ; sinewave_lead_1(3) = 0.0 ; sinewave_lead_1(4) = 0.0 ;
smoothperiod_out(0) = 0.0 ; smoothperiod_out(1) = 0.0 ; smoothperiod_out(2) = 0.0 ; smoothperiod_out(3) = 0.0 ; smoothperiod_out(4) = 0.0 ;
dcphase_vec(0) = 0.0 ; dcphase_vec(1) = 0.0 ; dcphase_vec(2) = 0.0 ; dcphase_vec(3) = 0.0 ; dcphase_vec(4) = 0.0 ;
smooth(0) = 0.0 ; smooth(1) = 0.0 ; smooth(2) = 0.0 ; smooth(3) = 0.0 ; smooth(4) = 0.0 ;
period(0) = 0.0 ; period(1) = 0.0 ; period(2) = 0.0 ; period(3) = 0.0 ; period(4) = 0.0 ;
smoothperiod(0) = 0.0 ; smoothperiod(1) = 0.0 ; smoothperiod(2) = 0.0 ; smoothperiod(3) = 0.0 ; smoothperiod(4) = 0.0 ;
detrender(0) = 0.0 ; detrender(1) = 0.0 ; detrender(2) = 0.0 ; detrender(3) = 0.0 ; detrender(4) = 0.0 ;
Q1(0) = 0.0 ; Q1(1) = 0.0 ; Q1(2) = 0.0 ; Q1(3) = 0.0 ; Q1(4) = 0.0 ;
I1(0) = 0.0 ; I1(1) = 0.0 ; I1(2) = 0.0 ; I1(3) = 0.0 ; I1(4) = 0.0 ;
jI(0) = 0.0 ; jI(1) = 0.0 ; jI(2) = 0.0 ; jI(3) = 0.0 ; jI(4) = 0.0 ;
jQ(0) = 0.0 ; jQ(1) = 0.0 ; jQ(2) = 0.0 ; jQ(3) = 0.0 ; jQ(4) = 0.0 ;
I2(0) = 0.0 ; I2(1) = 0.0 ; I2(2) = 0.0 ; I2(3) = 0.0 ; I2(4) = 0.0 ;
Q2(0) = 0.0 ; Q2(1) = 0.0 ; Q2(2) = 0.0 ; Q2(3) = 0.0 ; Q2(4) = 0.0 ;
sI2(0) = 0.0 ; sI2(1) = 0.0 ; sI2(2) = 0.0 ; sI2(3) = 0.0 ; sI2(4) = 0.0 ;
sQ2(0) = 0.0 ; sQ2(1) = 0.0 ; sQ2(2) = 0.0 ; sQ2(3) = 0.0 ; sQ2(4) = 0.0 ;
Re(0) = 0.0 ; Re(1) = 0.0 ; Re(2) = 0.0 ; Re(3) = 0.0 ; Re(4) = 0.0 ;
Im(0) = 0.0 ; Im(1) = 0.0 ; Im(2) = 0.0 ; Im(3) = 0.0 ; Im(4) = 0.0 ;
sRe(0) = 0.0 ; sRe(1) = 0.0 ; sRe(2) = 0.0 ; sRe(3) = 0.0 ; sRe(4) = 0.0 ;
sIm(0) = 0.0 ; sIm(1) = 0.0 ; sIm(2) = 0.0 ; sIm(3) = 0.0 ; sIm(4) = 0.0 ;
for ( octave_idx_type ii (5) ; ii < vec_length ; ii++ ) // Start the main loop
{
// smooth the price for hilbert calculations
smooth(ii) = (4.0 * price(ii) + 3.0 * price(ii-1) + 2.0 * price(ii-2) + price(ii-3) ) / 10.0 ;
// Detrend the input
detrender(ii) = (0.0962 * smooth(ii) + 0.5769 * smooth(ii-2) - 0.5769 * smooth(ii-4) - 0.0962 * smooth(ii-6)) * (0.075 * period(ii-1) + 0.54) ;
// Compute InPhase and Quadrature components
Q1(ii) = (0.0962 * detrender(ii) + 0.5769 * detrender(ii-2) - 0.5769 * detrender(ii-4) - 0.0962 * detrender(ii-6)) * (0.075 * period(ii-1) + 0.54) ;
I1(ii) = detrender(ii-3) ;
// Advance the phase of I1 and Q1 by 90 degrees
jI(ii) = (0.0962 * I1(ii) + 0.5769 * I1(ii-2) - 0.5769 * I1(ii-4) - 0.0962 * I1(ii-6)) * (0.075 * period(ii-1) + 0.54) ;
jQ(ii) = (0.0962 * Q1(ii) + 0.5769 * Q1(ii-2) - 0.5769 * Q1(ii-4) - 0.0962 * Q1(ii-6)) * (0.075 * period(ii-1) + 0.54) ;
// Phasor addition for 3 bar averaging
I2(ii) = I1(ii) - jQ(ii) ;
Q2(ii) = Q1(ii) + jI(ii) ;
// Smooth the I and Q components before applying the discriminator
sI2(ii) = 0.2 * I2(ii) + 0.8 * sI2(ii-1) ;
sQ2(ii) = 0.2 * Q2(ii) + 0.8 * sQ2(ii-1) ;
// Homodyne Discriminator
Re(ii) = sI2(ii) * sI2(ii-1) + sQ2(ii) * sQ2(ii-1) ;
Im(ii) = sI2(ii) * sQ2(ii-1) - sQ2(ii) * sI2(ii-1) ;
sRe(ii) = 0.2 * Re(ii) + 0.8 * sRe(ii-1) ;
sIm(ii) = 0.2 * Im(ii) + 0.8 * sIm(ii-1) ;
if ( (sIm(ii) > 0.0 || sIm(ii) < 0.0) && (sRe(ii) > 0.0 || sRe(ii) < 0.0) )
{
period(ii) = 360.0 / ( ((atan(sIm(ii) / sRe(ii))) * 180.0) / PI ) ;
}
else
{
period(ii) = period(ii-1) ;
}
if ( period(ii) > 1.5 * period(ii-1) )
{
period(ii) = 1.5 * period(ii-1) ;
}
if ( period(ii) < 0.67 * period(ii-1) )
{
period(ii) = 0.67 * period(ii-1) ;
}
if ( period(ii) < 6.0 )
{
period(ii) = 6.0 ;
}
if ( period(ii) > 50.0 )
{
period(ii) = 50.0 ;
}
period(ii) = 0.2 * period(ii) + 0.8 * period(ii-1) ;
smoothperiod(ii) = 0.33 * period(ii) + 0.67 * smoothperiod(ii-1) ;
// Compute Dominant Cycle
dcperiod = int ( smoothperiod(ii) + 0.5 ) ;
realpart = 0.0 ;
imagpart = 0.0 ;
dcphase = 0.0 ;
for ( octave_idx_type jj (0) ; jj <= ( dcperiod - 1 ) ; jj++ )
{
realpart += sin( PI/180.0 * 360.0 * jj / dcperiod ) * ( smooth(ii-jj) ) ;
imagpart += cos( PI/180.0 * 360.0 * jj / dcperiod ) * ( smooth(ii-jj) ) ;
}
if ( fabs( imagpart ) > 0.0 )
{
dcphase = atan( realpart / imagpart ) * 180.0 / PI ;
}
else if ( fabs( imagpart ) < 0.001 )
{
if ( realpart < 0.0 )
{
dcphase -= 90.0 ;
}
else if ( realpart > 0.0 )
{
dcphase += 90.0 ;
}
}
dcphase += 90.0 ;
// Compensate for one bar lag of the 4 bar weighted moving average
dcphase += 360.0 / smoothperiod(ii) ;
if ( imagpart < 0.0 )
dcphase += 180.0 ;
if ( dcphase > 315.0 )
dcphase -= 360.0 ;
// phase output
dcphase_vec(ii) = dcphase ;
//Now compute a differential phase, resolve phase wraparound, and limit delta phase errors
deltaphase(ii) = dcphase_vec(ii) - dcphase_vec(ii-1) ;
if ( dcphase_vec(ii-1) > 270.0 && dcphase_vec(ii) < 90.0 )
{
deltaphase(ii) = 360.0 - dcphase_vec(ii-1) + dcphase_vec(ii) ;
}
if ( deltaphase(ii) < 1.0 )
{
deltaphase(ii) = 1.0 ;
}
if ( deltaphase(ii) > 60.0 )
{
deltaphase(ii) = 60.0 ;
}
// Sum Deltaphases to reach 360 degrees. The sum is the instantaneous period.
sum_period(ii) = 0.0 ;
sum_deltaphase = 0.0 ;
count = 0 ;
while ( sum_deltaphase < 360.0 )
{
sum_deltaphase += deltaphase(ii-count) ;
count ++ ;
sum_period(ii) = count ;
}
// Resolve Instantaneous Period errors and smooth
if ( sum_period(ii) == 0.0 )
{
sum_period(ii) = sum_period(ii-1) ;
}
sumperiod(ii) = 0.25 * sum_period(ii) + 0.75 * sum_period(ii-1) ;
// sinewave output
sinewave(ii) = sin( dcphase * PI / 180.0 ) ;
// one bar leading function
sinewave_lead_1(ii) = sin( ( dcphase + 360.0 / smoothperiod(ii) ) * PI / 180.0 ) ;
// period output
smoothperiod_out(ii) = floor ( smoothperiod(ii) + 0.5 ) ;
} // end of main ii loop
retval_list(3) = dcphase_vec ;
retval_list(2) = smoothperiod_out ;
retval_list(1) = sinewave_lead_1 ;
retval_list(0) = sinewave ;
return retval_list ;
} // end of function
This is a straightforward conversion of the code available from here. A nice intro to how it can be used is here and Ehler's own website can be found here.



















