High Resolution Tools for Spectral Analysis - Update

Following on from my initial enthusiasm for the code on the High Resolution Tools for Spectral Analysis page, I have say that I have been unable to get the code performing as I would like it for my intended application to price time series.

My original intent was to use the zero crossing period estimation function, the subject of my last few posts, to get a rough idea of the dominant cycle period and then use the most recent data in a rolling window of this length as input to the high resolution code. This approach, however, ran into problems.

Firstly, windows of just the dominant cycle length (approximately 10 to 30 data points only) would lead to all sorts of errors being thrown from the toolkit functions as well as core Octave functions, such as divide by zero warnings and cryptic error messages that even now I don't understand. My best guess here is that the amount of data available in such short windows is simply insufficient for the algorithm to work, in much the same way as the Fast Fourier transform may fail to work if given too little data that is not a power of 2 in length. It might be possible to substantially rewrite the relevant functions, but my understanding of the algorithm and the inner workings of Octave means this is well beyond my pay grade.

My second approach was to simply extend the amount of data available by using the Octave repmat function on the windowed data so that all the above errors disappeared. This had very hit and miss results - sometimes they were very accurate, other times just fair to middling, and on occasion just way off target. I suspect here that the problem is the introduction of signal artifacts via the use of the repmat function, which results in Aliasing of the underlying signal.

As a result I shall not continue with investigating this toolbox for now. Although I only investigated the use of the me.m function (there are other functions available) I feel that my time at the moment can be used more productively.

Inputs For Zero Crossing Function

Continuing with the work on the zero crossing function, for now I have decided to use a "roofing filter" to detrend and smooth the raw price data and then apply some simple cycle extractor code

DEFUN_DLD ( cycle_extractor, args, nargout,
"-*- texinfo -*-\n\
@deftypefn {Function File} {} cycle_extractor (@var{price})\n\
This function takes a single price vector input and outputs\n\
a vector of the cycle extracted from the input.\n\
@end deftypefn" )

{
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 () ;
ColumnVector cycle = args(0).column_vector_value () ; cycle(0) = 0.0 ; cycle(1) = 0.0 ;
double alpha = 0.07 ;

 for ( octave_idx_type ii (2) ; ii < vec_length ; ii ++ ) // Start the main loop
     {
     cycle(ii) = ( 1.0 - 0.5 * alpha ) * ( 1.0 - 0.5 * alpha ) * ( price(ii) - 2.0 * price(ii-1) + price(ii-2) ) + 2.0 * ( 1.0 - alpha ) * cycle(ii-1) - ( 1.0 - alpha ) * ( 1.0 - alpha ) * cycle(ii-2) ;
     }
                                                                 
 retval_list(0) = cycle ;

return retval_list ; 
                                                                       
} // end of function

to create the inphase input for the zero cross function. For the imaginary or quadrature input I have decided to use the simple trigonometric identity of the derivative of a sine wave being the cosine (i.e. 90 degree phase lead), easily implemented using the bar to bar difference of the sine wave, or in our case the above simple cycle. I might yet change this, but for now it seems to work. The upper pane in the screen shot below shows the raw price in blue, the extracted cycle in black and the cycle derivative in red.

It can be seen that most of the time the cycle (inphase) either leads or is approximately in sync with the price, whilst the quadrature is nicely spaced between the zero crossings of the cycle.

The lower pane shows the measured periods as in the previous posts. The zero crossing measured periods are now a bit more erratic than before, but some simple tests show that, on average, the zero crossing measurement is consistently closer to the real period than the sine wave indicator period; however, this improvement cannot said to be statistically significant at any p-value.

Now I would like to introduce some other work I've been doing recently. For many years I have wanted to measure the instantaneous period using Maximum entropy spectral estimation, which has been popularised over the years by John Ehlers. Unfortunately I had never found any code in the public domain which I could use or adapt, until now this is. My discovery is High Resolution Tools For Spectral Analysis. This might not actually be the same as Ehlers' MESA, but it certainly covers the same general area and, joy of joys, it has freely downloadable MATLAB code along with an accessible description of the theoretical background along with some academic papers.

As is usual in situations like this, I have had to refactor some of the MATLAB code so that it can run in Octave without error messages; specifically this non Octave function

function [msg,A,B,C,D] = abcdchk(A,B,C,D)
%ABCDCHK Checks dimensional consistency of A,B,C,D matrices.
%   ERROR(ABCDCHK(A,B,C,D)) checks that the dimensions of A,B,C,D
%   are consistent for a linear, time-invariant system model.
%   An error occurs if the nonzero dimensions are not consistent.
%
%   [MSG,A,B,C,D] = ABCDCHK(A,B,C,D) also alters the dimensions
%   any 0-by-0 empty matrices to make them consistent with the others.

%   Copyright 1984-2001 The MathWorks, Inc.
%   $Revision: 1.21 $  $Date: 2001/04/15 11:59:11 $

if nargin < 4, D = []; end
if nargin < 3, C = []; end
if nargin < 2, B = []; end
if nargin < 1, A = []; end

[ma,na] = size(A);
[mb,nb] = size(B);
[mc,nc] = size(C);
[md,nd] = size(D);

if mc==0 & nc==0 & (md==0 | na==0)
   mc = md; nc = na; C = zeros(mc,nc);
end
if mb==0 & nb==0 & (ma==0 | nd==0)
   mb = ma; nb = nd; B = zeros(mb,nb);
end
if md==0 & nd==0 & (mc==0 | nb==0)
   md = mc; nd = nb; D = zeros(md,nd);
end
if ma==0 & na==0 & (mb==0 | nc==0)
   ma = mb; na = nc; A = zeros(ma,na);
end

if ma~=na & nargin>=1
   msg = 'The A matrix must be square';
elseif ma~=mb & nargin>=2
   msg = 'The A and B matrices must have the same number of rows.';
elseif na~=nc & nargin>=3
   msg = 'The A and C matrices must have the same number of columns.';
elseif md~=mc & nargin>=4
   msg = 'The C and D matrices must have the same number of rows.';
elseif nd~=nb & nargin>=4
   msg = 'The B and D matrices must have the same number of columns.';
else
   msg = '';
end

should be in the loadpath, and in the function file dlsim.m, on line 71, the call to the internal MATLAB function

x = ltitr( a , b , u , x0 ) ;

should be replaced by

x = lsim ( ss  (a , b , [ ] , [ ] , -1 ) , u , [ ] , x0 ) ;
(thanks to Lukas for this)

Having made these changes, this much simplified demo.m script

%% SIGNAL = 2 x sinusoids + noise
%  Setting up the signal parameters and time history
N = 100 ;

mag0 = 1.8 ; mag1 = 1.5 ; o1 = 1.3 ; mag2 = 2 ; o2 = 1.35 ;

t = 0 : N - 1 ; t = t(:) ;

y0 = mag0 * randn(N,1) ;

y1 = mag1 * exp( 1i * ( o1 * t + 2 * pi * rand) ) ;

y2 = mag2 * exp( 1i * ( o2 * t + 2 * pi * rand )) ;

y = real( y0 .+ y1 .+ y2 ) ;

NN = 2048 ; th = linspace( 0, 2 * pi, NN ) ;

%% setting up filter parameters and the svd of the input-to-state response
thetamid = 1.325 ; 

[ A , B ] = rjordan( 5, 0.88 * exp( thetamid * 1i ) ) ;

%% obtaining state statistics
R = dlsim_real( A, B, y' ) ;

%% maximum entropy
spectrum = me( R, A, B, th ) ;

spectrum = spectrum / max( spectrum ) * 1.2 ;

figure(1) ; 
subplot(211) ; plot( real(y1), 'r', real(y2), 'm', y, 'b' ) ; legend( 'Underlying Sinewave 1', 'Underlying Sinewave 2', '2 Sinewaves Plus Noise' ) ; title( 'Noisy Price series' ) ;
subplot(212) ; plot( th(400:475), spectrum(400:475) ) ; title( 'The Spectrum - Peaks should be centred at 1.3 and 1.35' ) ;

produces this plot

The upper pane shows two sine waves (red and magenta), very close to each other in period and amplitude, combined with random noise to create price (in blue). Looking at just the blue price, it would seem to be almost impossible that there are in fact two sine waves hidden within noise in this short time series, yet the algorithm clearly picks these out as shown by the two peaks in the spectral plot in the lower pane. Powerful stuff!

It is my intention over the coming days to investigate using the zero crossing function to select the data length prior to using this spectral analysis to determine the instantaneous period of price.

Zero Crossing Frequency/Period Measuring Function

Following on from my last post I have now finished coding the first two versions of this function, which are shown in the code boxes below.

DEFUN_DLD ( iq_zc_period_indicator_v1, args, nargout,
"-*- texinfo -*-\n\
@deftypefn {Function File} {} iq_zc_period_indicator_v1 (@var{Realpart,Quadrature})\n\
This function takes two input vectors of the real part and quadrature of price,\n\
the analytic signal, and outputs the period calculated from the zero crossings of\n\
the two input vectors.\n\
@end deftypefn" ) {
octave_value_list retval_list ;
int nargin = args.length () ;
int vec_length = args(0).length () ;

// check the input arguments
if ( nargin != 2 )
   {
   error ( "Invalid input argument(s). Inputs are two vectors, the realpart and quadrature." ) ;
   return retval_list ;
   }

if ( vec_length < 50 )
   {
   error ( "Invalid input argument length(s). Input vectors must be > 50 in length." ) ;
   return retval_list ;
   }
   
if ( args(1).length () != vec_length )
   {
   error ( "Invalid input argument length(s). Input vectors differ in length." ) ;
   return retval_list ;
   }   

if ( error_state )
   {
   error ( "Invalid input argument(s). Inputs are two vectors, the realpart and quadrature." ) ;
   return retval_list ;
   }
// end of input checking

// inputs
ColumnVector inphase_1 = args(0).column_vector_value () ;
ColumnVector quadrature_1 = args(1).column_vector_value () ;
ColumnVector f_hat ( vec_length ) ;
ColumnVector inphase_1_state ( vec_length ) ;
ColumnVector quadrature_1_state ( vec_length ) ;

for ( octave_idx_type ii (0) ; ii < 50 ; ii ++ ) // Start the initialising loop
{
f_hat(ii) = 20.0 ; // a default period length
inphase_1_state(ii) = 0.0 ; // assuming exactly on zero cross line
quadrature_1_state(ii) = 0.0 ; // assuming exactly on zero cross line
} // end of initialising ii loop

double latest_i_1_zc_time = 49.0 ; // as the "current" bar is bar no. 49
double previous_i_1_zc_time = 49.0 ; // as the "current" bar is bar no. 49
double latest_q_1_zc_time = 49.0 ; // as the "current" bar is bar no. 49
double previous_q_1_zc_time = 49.0 ; // as the "current" bar is bar no. 49

// get the initial states of the inphase_1 and quadrature_1 vectors. A value of
// 1 is assigned if > 0, -1 for < 0 & 0.0 if == 0
inphase_1_state(49) = inphase_1(49) > 0.0 ? 1.0 : ( inphase_1(49) < 0.0 ? -1.0 : 0.0 ) ;
quadrature_1_state(49) = quadrature_1(49) > 0.0 ? 1.0 : ( quadrature_1(49) < 0.0 ? -1.0 : 0.0 ) ;

for ( octave_idx_type ii (50) ; ii < vec_length ; ii ++ ) // Start the main loop
{
  
// first, assign current states for the inphase_1 and quadrature_1 vectors on this
// new (ii) bar
inphase_1_state(ii) = inphase_1(ii) >= 0.0 ? 1.0 : ( inphase_1(ii) < 0.0 ? -1.0 : 0.0 ) ;
quadrature_1_state(ii) = quadrature_1(ii) >= 0.0 ? 1.0 : ( quadrature_1(ii) < 0.0 ? -1.0 : 0.0 ) ;

  // check for a change in state, which indicates a zero crossing
  if ( inphase_1_state(ii) == 1.0 && inphase_1_state(ii-1) == -1.0 )
  { // inphase_1 has crossed up over the zero line
    latest_i_1_zc_time = double(ii-1) - inphase_1(ii-1) / ( inphase_1(ii) - inphase_1(ii-1) ) ; // calculate time of zero cross
    f_hat(ii) = ( latest_i_1_zc_time - previous_i_1_zc_time ) * 2.0 ; // calculate the period from times of zero crosses
    previous_i_1_zc_time = latest_i_1_zc_time ; // update the previous_zc_time to equal the latest_zc_time
  }
  else if ( inphase_1_state(ii) == -1.0 && inphase_1_state(ii-1) == 1.0 )
  { // inphase_1 has crossed down over the zero line
    latest_i_1_zc_time = double(ii-1) + inphase_1(ii-1) / ( inphase_1(ii-1) - inphase_1(ii) ) ; // calculate time of zero cross
    f_hat(ii) = ( latest_i_1_zc_time - previous_i_1_zc_time ) * 2.0 ; // calculate the period from times of zero crosses
    previous_i_1_zc_time = latest_i_1_zc_time ; // update the previous_zc_time to equal the latest_zc_time    
  }
  else if ( quadrature_1_state(ii) == 1.0 && quadrature_1_state(ii-1) == -1.0 )
  { // quadrature_1 has crossed up over the zero line
    latest_q_1_zc_time = double(ii-1) - quadrature_1(ii-1) / ( quadrature_1(ii) - quadrature_1(ii-1) ) ; // calculate time of zero cross
    f_hat(ii) = ( latest_q_1_zc_time - previous_q_1_zc_time ) * 2.0 ; // calculate the period from times of zero crosses
    previous_q_1_zc_time = latest_q_1_zc_time ; // update the previous_zc_time to equal the latest_zc_time
  }
  else if ( quadrature_1_state(ii) == -1.0 && quadrature_1_state(ii-1) == 1.0 )
  { // quadrature_1 has crossed down over the zero line
    latest_q_1_zc_time = double(ii-1) + quadrature_1(ii-1) / ( quadrature_1(ii-1) - quadrature_1(ii) ) ; // calculate time of zero cross
    f_hat(ii) = ( latest_q_1_zc_time - previous_q_1_zc_time ) * 2.0 ; // calculate the period from times of zero crosses
    previous_q_1_zc_time = latest_q_1_zc_time ; // update the previous_zc_time to equal the latest_zc_time 
  }
  else
  { // neither the inphase_1 or quadrature_1 have crossed the zeroline on this bar
    f_hat(ii) = f_hat(ii-1) ; // set to most recently calcualted zero crossing period length
  }
  
} // end of main ii loop
 
retval_list(0) = f_hat ;

return retval_list ; 
                                                                       
} // end of function

and, almost identical to the above

DEFUN_DLD ( iq_zc_period_indicator_v1_1, args, nargout,
"-*- texinfo -*-\n\
@deftypefn {Function File} {} iq_zc_period_indicator_v1_1 (@var{Realpart,Quadrature})\n\
This function takes two input vectors of the real part and quadrature of price,\n\
the analytic signal, and outputs the period calculated from the zero crossings of\n\
the two input vectors.\n\
@end deftypefn" )

{
octave_value_list retval_list ;
int nargin = args.length () ;
int vec_length = args(0).length () ;

// check the input arguments
if ( nargin != 2 )
   {
   error ( "Invalid input argument(s). Inputs are two vectors, the realpart and quadrature." ) ;
   return retval_list ;
   }

if ( vec_length < 50 )
   {
   error ( "Invalid input argument length(s). Input vectors must be > 50 in length." ) ;
   return retval_list ;
   }
   
if ( args(1).length () != vec_length )
   {
   error ( "Invalid input argument length(s). Input vectors differ in length." ) ;
   return retval_list ;
   }   

if ( error_state )
   {
   error ( "Invalid input argument(s). Inputs are two vectors, the realpart and quadrature." ) ;
   return retval_list ;
   }
// end of input checking

// inputs
ColumnVector inphase_1 = args(0).column_vector_value () ;
ColumnVector quadrature_1 = args(1).column_vector_value () ;
ColumnVector f_hat ( vec_length ) ;
ColumnVector inphase_1_state ( vec_length ) ;
ColumnVector quadrature_1_state ( vec_length ) ;

for ( octave_idx_type ii (0) ; ii < 50 ; ii ++ ) // Start the initialising loop
{
f_hat(ii) = 20.0 ; // a default period length
inphase_1_state(ii) = 0.0 ; // assuming exactly on zero cross line
quadrature_1_state(ii) = 0.0 ; // assuming exactly on zero cross line
} // end of initialising ii loop

double latest_iq_1_zc_time = 49.0 ; // as the "current" bar is bar no. 49
double previous_iq_1_zc_time = 49.0 ; // as the "current" bar is bar no. 49

// get the initial states of the inphase_1 and quadrature_1 vectors. A value of
// 1 is assigned if > 0, -1 for < 0 & 0.0 if == 0
inphase_1_state(49) = inphase_1(49) > 0.0 ? 1.0 : ( inphase_1(49) < 0.0 ? -1.0 : 0.0 ) ;
quadrature_1_state(49) = quadrature_1(49) > 0.0 ? 1.0 : ( quadrature_1(49) < 0.0 ? -1.0 : 0.0 ) ;

for ( octave_idx_type ii (50) ; ii < vec_length ; ii ++ ) // Start the main loop
{
  
// first, assign current states for the inphase_1 and quadrature_1 vectors on this
// new (ii) bar
inphase_1_state(ii) = inphase_1(ii) >= 0.0 ? 1.0 : ( inphase_1(ii) < 0.0 ? -1.0 : 0.0 ) ;
quadrature_1_state(ii) = quadrature_1(ii) >= 0.0 ? 1.0 : ( quadrature_1(ii) < 0.0 ? -1.0 : 0.0 ) ;

  // check for a change in state, which indicates a zero crossing
  if ( inphase_1_state(ii) == 1.0 && inphase_1_state(ii-1) == -1.0 )
  { // inphase_1 has crossed up over the zero line
    latest_iq_1_zc_time = double(ii-1) - inphase_1(ii-1) / ( inphase_1(ii) - inphase_1(ii-1) ) ; // calculate time of zero cross
    f_hat(ii) = ( latest_iq_1_zc_time - previous_iq_1_zc_time ) * 4.0 ; // calculate the period from times of zero crosses
    previous_iq_1_zc_time = latest_iq_1_zc_time ; // update the previous_zc_time to equal the latest_zc_time
  }
  else if ( inphase_1_state(ii) == -1.0 && inphase_1_state(ii-1) == 1.0 )
  { // inphase_1 has crossed down over the zero line
    latest_iq_1_zc_time = double(ii-1) + inphase_1(ii-1) / ( inphase_1(ii-1) - inphase_1(ii) ) ; // calculate time of zero cross
    f_hat(ii) = ( latest_iq_1_zc_time - previous_iq_1_zc_time ) * 4.0 ; // calculate the period from times of zero crosses
    previous_iq_1_zc_time = latest_iq_1_zc_time ; // update the previous_zc_time to equal the latest_zc_time    
  }
  else if ( quadrature_1_state(ii) == 1.0 && quadrature_1_state(ii-1) == -1.0 )
  { // quadrature_1 has crossed up over the zero line
    latest_iq_1_zc_time = double(ii-1) - quadrature_1(ii-1) / ( quadrature_1(ii) - quadrature_1(ii-1) ) ; // calculate time of zero cross
    f_hat(ii) = ( latest_iq_1_zc_time - previous_iq_1_zc_time ) * 4.0 ; // calculate the period from times of zero crosses
    previous_iq_1_zc_time = latest_iq_1_zc_time ; // update the previous_zc_time to equal the latest_zc_time
  }
  else if ( quadrature_1_state(ii) == -1.0 && quadrature_1_state(ii-1) == 1.0 )
  { // quadrature_1 has crossed down over the zero line
    latest_iq_1_zc_time = double(ii-1) + quadrature_1(ii-1) / ( quadrature_1(ii-1) - quadrature_1(ii) ) ; // calculate time of zero cross
    f_hat(ii) = ( latest_iq_1_zc_time - previous_iq_1_zc_time ) * 4.0 ; // calculate the period from times of zero crosses
    previous_iq_1_zc_time = latest_iq_1_zc_time ; // update the previous_zc_time to equal the latest_zc_time 
  }
  else
  { // neither the inphase_1 or quadrature_1 have crossed the zeroline on this bar
    f_hat(ii) = f_hat(ii-1) ; // set to most recently calcualted zero crossing period length
  }
  
} // end of main ii loop
 
retval_list(0) = f_hat ;

return retval_list ; 
                                                                       
} // end of function

The first measures the period every quarter of a cycle, based on a half period measure of zero crossings, and the second every quarter of a cycle based on quarter cycle measurements. The screen shot below shows typical performance on ideal sine and cosine inputs (top pane) with the period measurements shown in the lower pane.

It can be seen that they both perform admirably within their respective theoretical constraints on this ideal input. The subject of my next post will concentrate on the inputs that can be realistically extracted from live price data.

Update on Zero Crossing Frequency Measurement

Having now returned from a long summer hiatus I've got back to work on the subject matter of my last post and I think I've made some significant progress. Below is a screen shot of the fruits of my work so far.

The top pane shows, in blue, a series of randomly selected sine waves stitched together to form a sideways "price" channel of various periods of 10 to 30 inclusive, whilst the bottom pane shows the period measurements of the "price." The blue period line is the actual period of the blue "price" in the upper pane whilst the red line is the measured period of "price" using my current implementation of a zero crossing period measurement indicator using the actual phase of "price" as input.

Of course in reality one would not know the actual phase but would have to estimate it, which I'm currently doing via a Hilbert Sine Wave indicator and which is shown in green in the upper pane above. It can be seen that this is not exactly in sync with the blue price, and the period of price as extracted via the sine wave indicator calculations is shown in green in the the lower pane above.The phase as extracted by the Sine Wave indicator is compared with the actual phase in the screen shot below, where blue is the actual and red the extracted phase.

Obviously there is a difference, which manifests itself as an extremely erratic zero crossing measurement when this extracted phase is used as input to the zero crossing indicator, shown in magenta in the screen shot below.

It is very encouraging that the zero crossing indicator performs in an almost theoretically perfect manner when supplied with good inputs, but it is disappointing that such inputs might not be able to be supplied by the Sine Wave indicator. I surmise that the delay induced by the calculations of the Hilbert transform might be partly to blame, but there might also be some improvements I can make to my implementation of the Sine Wave indicator (elimination of coding bugs or perhaps some coding kluge?). I will focus on this over the coming days.