John Ehler's Sinewave Indicator Code

A reader recently inquired about my use of this indicator and so below I provide my Octave C++ .oct  version that I have been using for the past few years.

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. 

Denoising Autoencoder MATLAB/Octave Code

Following on from my last post I have been looking for Octave code for the denoising autoencoder to avoid reinventing the wheel and writing it myself from scratch, and luckily I have found two options. The first is a tutorial on autoencoders, by a Piotr Mirowski, which has a link to a Github page with code. The second is a toolbox for dimensionality reduction, by Laurens van der Maaten, which has autoencoders as one of its options.

As of now I'm not sure which one I'll use, and perhaps I might yet write my own code heavily reusing elements from both of the above. More in due course.

Recent Readings and New Directions.

Since my last post I have been doing a fair bit of online research and fortunately I have discovered the following papers, which mesh nicely with what I am trying to do with Conditional Restricted Boltzmann Machines to model time series:-

Deep Learning Architecture for Univariate Time Series Forecasting
Temporal Autoencoding Restricted Boltzmann Machine
Temporal Autoencoding Improves Generative Models of Time Series
Deep Modelling Complex Couplings Within Financial Markets
Predicting Time Series of Railway Speed Restrictions with Time Dependent Machine Learning Techniques

The big take away from these readings is to explicitly model the autoregressive components via a Denoising Autoencoder and secondly, not to model a univariate time series in isolation, but model a multivariate time series where the "other" time series are either informative measures taken from the univariate series itself (informative indicators?) and/or related time series e.g in forex one could use concepts similar to fractional product inefficiency or on all markets the concept of Intermarket analysis.

For the nearest future I have therefore set myself the task of adapting my CRBM code to include the denoising autoencoder and to investigate the multivariate time series approach.

Giving up on Runge-Kutta Methods (for now?)

Over the last few weeks I have been looking at using Runge-Kutta methods for the creation of features, but I have decided to give up on this for now simply because I think I have found a better way to accomplish what I want. I was alerted to this possible approach by this post over at http://dsp.stackexchange.com/ and following up on this I remembered that a few years ago I coded up John Ehler's Linear Prediction Filter (the white paper for which might still be available here) and my crudely transcribed code given below:

 DEFUN_DLD (linearpredict, args, , "Help String")
{
octave_value retval;

 ColumnVector a = args(0).column_vector_value ();                                   // The input price                          
 ColumnVector b(a);                                                                 // This will be the output column vector returned to Octave by "retval"
 const int n = args(1).int_value();
 int lngth = 10;                                                                 
 
 double g[30];
 int zz;
 for (zz = 0; zz < 30; zz++)
     g[zz] = 0.0;

 double sigPredict[30];
 for (zz = 0; zz < 30; zz++)
     sigPredict[zz] = 0.0;

 double sigPower = 0.0;
 double mu = 0.0;
 double xBar = 0.0;
 int jj = 0;                                                          
 for (octave_idx_type ii (10); ii < a.length="" 10="" a="" average="" factor="" for="" href="https://draft.blogger.com/null" if="" ii="" jj="" lngth="" loop="" normalization="" ompute="" onvergence="" power="" sigpower="" start="" the=""> 0)
         mu = 0.25 / (sigPower * 10);

      //Compute signal estimate
      xBar = 0;
      for (jj = 1; jj <= lngth; jj++)
          xBar = xBar + a(ii - jj) * g[jj];

     //Compute gain coefficients
     for (jj = 1; jj <= lngth; jj++)
         g[jj] = g[jj] + (mu * (a(ii) - xBar) * a(ii - jj));

     //Compute signal prediction waveform
     for (jj = 0; jj <= lngth; jj++)
         sigPredict[jj] = a(ii - (10 - jj));

     //Extend signal prediction into the future
     int kk = 0;
     for (jj = lngth + 1; jj <= lngth + 5; jj++)
     {
         sigPredict[jj] = 0;

         for (kk = 1; kk <= lngth; kk++)
	     sigPredict[jj] = sigPredict[jj] + sigPredict[jj - kk] * g[kk];
     }

     b(ii) = sigPredict[lngth + n];

     }
 
retval = b;                                                                         // Assign the output column vector to the return value

return retval;                                                                       // Return the output to Octave
}

which is very similar in concept to Burg's method. I think some application of these methods show more promise than concentrating on my naive implementation of Runge-Kutta. 

Runge-Kutta Example and Code

Following on from my last post I thought I would, as a first step, code up a "straightforward" Runge-Kutta function and show how to deal with the fact that there is no "magic mathematical formula" to calculate the slopes that are an integral part of Runge-Kutta.

My approach is to fit a quadratic function to the last n_bars of price and take the slope of this via my Savitzky-Golay filter convolution code, and in doing so the Runge-Kutta value k1 can easily be obtained. The extrapolation beyond the k1 point to the points k2 and k3 is, in effect, trying to fit to points that have a half bar lead over the last available price. To accommodate this I use the last n_bar - 1 points of a 2 bar simple moving average plus the "position" of points k2  and k3 to calculate the slopes at k2 and k3. A 2 bar simple moving average is used because this has a half bar lag and is effectively an interpolation of the known prices at the "half bar intervals," and therefore points k2 and k3 are one h step ahead of the last half bar interval. The k4 point is again simply calculated directly from prices plus the half bar interval projection from point k3. If all this seems confusing, hopefully the Octave code below will clear things up.

clear all

%  create the raw price series
period = input( 'Period? ' ) ;
sideways = zeros( 1 , 2*period ) ;
uwr = ( 1 : 1 : 2*period ) .* 1.333 / period ; uwr_end = uwr(end) ;
unr = ( 1 : 1 : 2*period ) .* 4 / period ; unr_end = unr(end) ;
dwr = ( 1 : 1 : 2*period ) .* -1.333 / period ; dwr_end = dwr(end) ;
dnr = ( 1 : 1 : 2*period ) .* -4 / period ; dnr_end = dnr(end) ;

trends = [ sideways , uwr , unr.+uwr_end , sideways.+uwr_end.+unr_end , dnr.+uwr_end.+unr_end , dwr.+uwr_end.+unr_end.+dnr_end , sideways ] .+ 2 ;
noise = randn( 1 , length(trends) ) .* 0.0 ;
price = sinewave( length( trends ) , period ) .+ trends .+ noise ;
ma_2 = sma( price , 2 ) ;
 
%  regress over 'n_bar' bars
n_bar = 9 ;
%  and a 'p' order fit
p = 2 ;

%  get the relevant coefficients
slope_coeffs = generalised_sgolay_filter_coeffs( n_bar , p , 1 ) ; 

%  container for 1 bar ahead projection
projection_1_bar = zeros( 1 , length( price ) ) ;

for ii = n_bar : length( price )

%  calculate k1 value i.e. values at price(ii), the most recent price
k1 = price( ii-(n_bar-1) : ii ) * slope_coeffs( : , end ) ;
projection_of_point_k2 = price(ii) + k1 / 2 ;

%  calculate k2 value
k2 = [ ma_2( ii-(n_bar-2) : ii ) ; projection_of_point_k2 ]' * slope_coeffs( : , end ) ;
projection_of_point_k3 = price(ii) + k2 / 2 ;

%  calculate k3 value
k3 = [ ma_2( ii-(n_bar-2) : ii ) ; projection_of_point_k3 ]' * slope_coeffs( : , end ) ;
projection_of_point_k4 = price(ii) + k3 / 2 ;

%  calculate k4 value
k4 = [ price( ii-(n_bar-2) : ii ) , projection_of_point_k4 ] * slope_coeffs( : , end ) ;

%  the runge-kutta weighted moving average
projection_1_bar(ii) = price(ii) + ( k1 + 2 * ( k2 + k3 ) + k4 ) / 6 ;

end

%  shift for plotting
projection_1_bar = shift( projection_1_bar , 1 ) ;
projection_1_bar( : , 1:n_bar ) = price( : , 1:n_bar ) ;

plot( price , 'c' , projection_1_bar , 'r' ) ;

This code produces a plot like this, without noise,

and this with noise ( line 12 of the code ).

The cyan line is the underlying price and the red line is the Runge-Kutta 1 bar ahead projection. As can be seen, when the price is moving in rather straight lines the projection is quite accurate, however, at turnings there is some overshoot, which is to be expected. I'm not unduly concerned about this overshoot as my intention is simply to get the various k values as features, but this overshoot does have some implications which I will discuss in a future post.

Runge-Kutta Methods

As stated in my previous post I have been focusing on getting some meaningful features as possible inputs to my machine learning based trading system, and one of the possible ideas that has caught my attention is using Runge-Kutta methods to project ( otherwise known as "guessing" ) future price evolution. I have used this sort of approach before in the construction of my perfect oscillator ( links here, here, here and here ) but I deem this approach unsuitable for the type of "guessing" I want to do for my mfe-mae indicator as I need to "guess" rolling maximum highs and minimum lows, which quite often are actually flat for consecutive price bars. In fact, what I want to do is predict/guess upper and lower bounds for future price evolution, which might actually turn out to be a more tractable problem than predicting price itself.

The above wiki link to Runge-Kutta methods is a pretty dense mathematical read and readers may be wondering how approximation of solutions to ordinary differential equations can possibly relate to my stated aim, however the following links visualise Runge-Kutta in an accessible way:

• https://jmckennonmth212s09.wordpress.com/2009/02/03/runge-kutta-method-vs-euler-method/
• https://wigglewave.wordpress.com/2012/11/09/intuitive-description-of-runge-kutta-integration/
• https://www.youtube.com/watch?v=0iRZqmzs9dk
• http://buttersblog.com/runge-kutta/

Readers should hopefully see that what Runge-Kutta basically does is compute a "future" position given a starting value and a function to calculate slopes. When it comes to prices our starting value can be taken to be the most recent price/OHLC/bid/offer/tick in the time series, and whilst I don't possess a mathematical function to exactly calculate future evolution of price slopes, I can use my Savitzky Golay filter convolution code to approximate these slopes. The final icing on this cake is the weighting scheme for the Runge-Kutta method, as shown in the buttersblog link above, i.e.

which is just linear regression on the k slope values with the most recent price set as the intercept term! Hopefully this will be a useful way to generate features for my conditional restricted boltzmann machine, and if I use regularized linear regression I might finally be able to use my particle swarm optimisation code.

More in due course.

Recent reading

In my last post I mentioned that I was going away for the summer, but now I'm back. During the summer I didn't get to do as much reading etc. as I had hoped, but I did manage to play around with the Rssa package for singular spectrum analysis and this is still an ongoing investigation. I also briefly looked at independent component analysis and the FastICA package.

The common theme of the above is the extraction of meaningful time series features, and this general area is what I will be looking into for my next set of posts.