Conditional Restricted Boltzmann Machine

I have recently been looking at using a Conditional Restricted Boltzmann Machine, and in particular Graham Taylor's thesis paper, Composable Distributed-State Models for High-Dimensional Time Series and the associated code available from here. My adaptation of this is very much a work in progress, but I think it holds great promise.

Instead of using the "label" of my MFE/MAE indicator as determined by the matches found by my Cauchy Schwarz matching algorithm, I'm thinking of using the above CRBM to directly model the price action that immediately follows said matches, use the generative ability of the CRBM to create a set of modelled price bars and then apply the MFE/MAE indicator to these modelled price bars. In this way, I should be able to generate a sort of distribution of future MFE/MAE values as input for making trading decisions.

In my adaptation of the above linked code I have replaced the somewhat complicated motion capture input data with a very simple set of 5 features of the log changes in price between consecutive bars, and also made the code easier to read by spacing it out, using underscores to separate compound word named variables, added some extra comments and renamed some variables to be more relevant to the intended use. I have also split the code into separately callable functions, with an eye on future coding into .oct C++ functions. However, in essence, the code is much the same as is downloadable from the above code link. The code blocks below show my code adaptations so far:-
first, the basic training file

clear all

% load price file of interest
filename = input( 'Enter filename for prices, e.g. es or esmatrix: ' , 's' ) ;
data = load( "-ascii" , filename ) ;

% get tick size
switch filename

case { "cc" }
tick = 1 ;

case { "gc" "lb" "pl" "sm" "sp" }
tick = 0.1 ;

case { "ausyen" "bo" "cl" "ct" "dx" "euryen" "gbpyen" "sb" "usdyen" }
tick = 0.01 ;

case { "ng" }
tick = 0.001 ;

case { "auscad" "aususd" "euraus" "eurcad" "eurchf" "eurgbp" "eurusd" "gbpchf" "gbpusd" "ho" "rb" "usdcad" "usdchf" }
tick = 0.0001 ;

case { "c" "o" "s" "es" "nd" "w" }
tick = 0.25 ;

case { "fc" "lc" "lh" "pb" }
tick = 0.025 ;

case { "ed" }
tick = 0.0025 ;

case { "si" }
tick = 0.5 ;

case { "hg" "kc" "oj" "pa" }
tick = 0.05 ;

case { "ty" "us" }
tick = 0.015625 ;

case { "ccmatrix" }
tick = 1 ;

case { "gcmatrix" "lbmatrix" "plmatrix" "smmatrix" "spmatrix" }
tick = 0.1 ;

case { "ausyenmatrix" "bomatrix" "clmatrix" "ctmatrix" "dxmatrix" "euryenmatrix" "gbpyenmatrix" "sbmatrix" "usdyenmatrix" }
tick = 0.01 ;

case { "ngmatrix" }
tick = 0.001 ;

case { "auscadmatrix" "aususdmatrix" "eurausmatrix" "eurcadmatrix" "eurchfmatrix" "eurgbpmatrix" "eurusdmatrix" "gbpchfmatrix" "gbpusdmatrix" "homatrix" "rbmatrix" "usdcadmatrix" "usdchfmatrix" }
tick = 0.0001 ;

case { "cmatrix" "omatrix" "smatrix" "esmatrix" "ndmatrix" "wmatrix" }
tick = 0.25 ;

case { "fcmatrix" "lcmatrix" "lhmatrix" "pbmatrix" }
tick = 0.025 ;

case { "edmatrix" }
tick = 0.0025 ;

case { "simatrix" }
tick = 0.5 ;

case { "hgmatrix" "kcmatrix" "ojmatrix" "pamatrix" }
tick = 0.05 ;

case { "tymatrix" "usmatrix" }
tick = 0.015625 ;

endswitch

%  get data
open = data( : , 4 ) ;
high = data( : , 5 ) ;
low = data( : , 6 ) ;
close = data( : , 7 ) ;

clear data

%  create log series for crbm training
vwap = vwap_all_period_detrend( open , high , low , close , tick ) ;
open_from_vwap = log( open ./ shift( vwap , 1 ) ) ; open_from_vwap(1) = 0 ;
open_from_previous_close = log( open ./ shift( close , 1 ) ) ; open_from_previous_close(1) = 0 ;
high_from_vwap = log( high ./ vwap ) ; high_from_vwap(1) = 0 ;
low_from_vwap = log( low ./ vwap ) ; low_from_vwap(1) = 0 ;
close_from_vwap = log( close ./ vwap ) ; close_from_vwap(1) = 0 ;

%  assemble into a features matrix
features = [ open_from_vwap open_from_previous_close high_from_vwap low_from_vwap close_from_vwap ] ;

%  now prepare data for the NN training
%  visually inspect chart to choose a training target
pkg load financial
clf () ;
figure(1) ; highlow( high , low , close , open ) ;
training_point_index = input( 'Choose a training point for NN training by entering an index from figure 1: ' ) ;
bar_number = 0 ;

%  load the pre-computed top 500 Cauchy Schwarz Cross validation index matches
load eurusd_top_500_cv_matches

%  select how many of the above matches to use - an optimisable parameter via pso?
number_of_training_examples_to_use = 450 ;

for ii = 0 : 2
%  bar_number = input( 'Enter a 1, 2 or 3 to identify the bar sequence order' ) ;
bar_number = bar_number + 1 ;

% using training_point_index get the top number_of_training_examples_to_use matches
features_index = eurusd_top_500_cv_matches( training_point_index+ii , 1 : number_of_training_examples_to_use )' ;

%  use the above features_index to create batchdata for training, assembled into form suitable for crbm training
%  how-many timesteps do we look back for directed connections - this is what we call the "order" of the model 
n1 = 3 ; % first layer
n2 = 3 ; % second layer

%  create batchdataindex to index bars that include bars at timesteps t-2, t-1, t and t+1 where t is the bar at features_index
batchdataindex = [ features_index ; features_index.- 2 ; features_index.-1 ; features_index.+1 ] ;
%  sorted_batchdata_index = sort( batchdataindex ) ;

%  shuffle training set by shuffling the index values
batchdataindex = batchdataindex( randperm( size( batchdataindex , 1 ) ) , : ) ;

%  normalise the features matrix by the mean and std of just the batchdataindex values in features matrix
data_mean = mean( features( batchdataindex ) , 1 ) ;
data_std = std( features( batchdataindex ) ) ;
batchdata = ( features - repmat( data_mean , size( features , 1 ) , 1 ) ) ./ repmat( data_std , size( features , 1 ) , 1 ) ;

if ( bar_number == 1 )
data_mean_1 = data_mean ; data_std_1 = data_std ;
save bar_1_mean_std.mat data_mean_1 data_std_1 ;
elseif ( bar_number == 2 )
data_mean_2 = data_mean ; data_std_2 = data_std ;
save bar_2_mean_std.mat data_mean_2 data_std_2 ;
elseif ( bar_number == 3 )
data_mean_3 = data_mean ; data_std_3 = data_std ;
save bar_3_mean_std.mat data_mean_3 data_std_3 ;
else
fprintf( 'Have entered wrong value for bar number.\n' ) ;
end

%  now create minibatches from the batchdataindex in the form of a cell array
batchsize = 100 ;
minibatchindex = reshape( batchdataindex , size( batchdataindex , 1 ) / batchsize , batchsize ) ;
minibatch = num2cell( minibatchindex , 2 ) ;

%  Set network properties
numdims = size( features , 2 ) ;
%  numhid1 = 150 ; numhid2 = 150 ; numepochs = 2000 ; % original
numhid1 = 100 ; numhid2 = 100 ; numepochs = 2000 ;
gsd = 1 ;          % fixed standard deviation for Gaussian units
nt = n1 ;          % crbm "order"
numhid = numhid1 ;
%  fprintf( 1 , 'Training Layer 1 CRBM, order %d: %d-%d \n', nt , numdims , numhid ) ;
restart = 1 ;      % initialize weights
tic() ;
[ w1 , bj1 , bi1 , A1 , B1 ] = gaussian_crbm( batchdata , minibatch , nt , gsd , numepochs , numhid , restart ) ;
gaussian_crbm_training_time = toc()

if ( bar_number == 1 )
w11 = w1 ; bj11 = bj1 ; bi11 = bi1 ; A11 = A1 ; B11 = B1 ;
save layer1.mat w11 bj11 bi11 A11 B11 ;
elseif ( bar_number == 2 )
w12 = w1 ; bj12 = bj1 ; bi12 = bi1 ; A12 = A1 ; B12 = B1 ;
save layer1_2.mat w12 bj12 bi12 A12 B12 ;
elseif ( bar_number == 3 )
w13 = w1 ; bj13 = bj1 ; bi13 = bi1 ; A13 = A1 ; B13 = B1 ;
save layer1_3.mat w13 bj13 bi13 A13 B13 ;
else
fprintf( 'Have entered wrong value for bar number.\n' ) ;
end

tic() ;
[ filtering_distance ] = get_filtering_distance( batchdata , unique( batchdataindex ) , numhid , w1 , bj1 , n1 , n2 , B1 , gsd ) ;
get_filtering_distance_timing = toc()

%  now use filtering_distance as new features for training of second layer
%  first create a new batchdata set
batchdata_2 = zeros( size( features , 1 ) , size( filtering_distance , 2 ) ) ;
%  and fill with the relevant, indexed values from filtering_distance
batchdata_2( unique( batchdataindex ) , : ) = filtering_distance ;

%  create new minibatch to index batchdata_2
batchsize = 25 ;
minibatchindex = reshape( features_index .+ 1 , size( features_index , 1 ) / batchsize , batchsize ) ;
minibatch = num2cell( minibatchindex , 2 ) ;

%  set network properties for the 2nd layer training
numhid = numhid2 ; nt = n2 ;
numdims = size( batchdata_2 , 2 ) ; % data (visible) dimension
%  fprintf( 1 , 'Training Layer 2 CRBM, order %d: %d-%d \n' , nt , numdims , numhid ) ;
restart = 1 ; % initialize weights
tic() ;
[ w2 , bj2 , bi2 , A2 , B2 ] = binary_crbm( batchdata_2 , minibatch , nt , numepochs , numhid , restart ) ;
binary_crbm_training_time = toc()

if ( bar_number == 1 )
w21 = w2 ; bj21 = bj2 ; bi21 = bi2 ; A21 = A2 ; B21 = B2 ;
save layer2.mat w21 bj21 bi21 A21 B21 ;
elseif ( bar_number == 2 )
w22 = w2 ; bj22 = bj2 ; bi22 = bi2 ; A22 = A2 ; B22 = B2 ;
save layer2_2.mat w22 bj22 bi22 A22 B22 ;
elseif ( bar_number == 3 )
w23 = w2 ; bj23 = bj2 ; bi23 = bi2 ; A23 = A2 ; B23 = B2 ;
save layer2_3.mat w23 bj23 bi23 A23 B23 ;
else
fprintf( 'Have entered wrong value for bar number.\n' ) ;
end

end % end of training_point_index loop

which in turn calls gaussian_crbm.m

function [ w , bj , bi , A , B ] = gaussian_crbm( batchdata , minibatch , nt , gsd , num_epochs , num_hid , restart )

%  Version 1.01 
%  %
%  Code provided by Graham Taylor, Geoff Hinton and Sam Roweis 
%  %
%  For more information, see:
%     http://www.cs.toronto.edu/~gwtaylor/publications/nips2006mhmublv
%  %
%  Permission is granted for anyone to copy, use, modify, or distribute this
%  program and accompanying programs and documents for any purpose, provided
%  this copyright notice is retained and prominently displayed, along with
%  a note saying that the original programs are available from our
%  web page.
%  The programs and documents are distributed without any warranty, express or
%  implied.  As the programs were written for research purposes only, they have
%  not been tested to the degree that would be advisable in any important
%  application.  All use of these programs is entirely at the user's own risk.
%  
%  This program trains a Conditional Restricted Boltzmann Machine in which
%  visible, Gaussian-distributed inputs are connected to
%  hidden, binary, stochastic feature detectors using symmetrically
%  weighted connections. Learning is done with 1-step Contrastive Divergence.
%  Directed connections are present, from the past nt configurations of the
%  visible units to the current visible units (A), and the past nt
%  configurations of the visible units to the current hidden units (B)
%  
%  The program assumes that the following variables are set externally:
%  nt        -- order of the model
%  gsd       -- fixed standard deviation of Gaussian visible units
%  num_epochs -- maximum number of epochs
%  num_hid    --  number of hidden units 
%  batchdata --  a matrix of data (num_cases,num_dims) 
%  minibatch -- a cell array of dimension batchsize, indexing the valid
%  frames in batchdata
%  restart   -- set to 1 if learning starts from beginning 

%  batchdata is a big matrix of all the frames
%  we index it with "minibatch", a cell array of mini-batch indices
num_batches = length( minibatch ) ; 

num_dims = size( batchdata , 2 ) ; % visible dimension

%  Setting learning rates
epsilon_w = 1e-3 ; % undirected
epsilon_bi = 1e-3 ; % visibles
epsilon_bj = 1e-3 ; % hidden units
epsilon_A = 1e-3 ; % autoregressive
epsilon_B = 1e-3 ; % prev visibles to hidden

w_decay = 0.0002 ; % currently we use the same weight decay for w, A, B
mom = 0.9 ;       % momentum used only after 5 epochs of training

if ( restart == 1 )  
    restart = 0 ;
    epoch = 1 ;

    % Randomly initialize weights
    w = 0.01 * randn( num_hid , num_dims ) ;
    bi = 0.01 * randn( num_dims , 1 ) ;
    bj = -1 + 0.01 * randn( num_hid , 1 ) ; % set to favour units being "off"

    % The autoregressive weights; A(:,:,j) is the weight from t-j to the vis
    A = 0.01 * randn( num_dims ,num_dims ,nt ) ;

    % The weights from previous time-steps to the hiddens; B(:,:,j) is the
    % weight from t-j to the hidden layer
    B = 0.01 * randn( num_hid , num_dims , nt ) ;

    clear w_grad bi_grad bj_grad A_grad B_grad
    clear neg_w_grad neg_bi_grad neg_bj_grad neg_A_grad neg_B_grad

    % keep previous updates around for momentum
    w_update = zeros( size( w ) ) ;
    bi_update = zeros( size( bi ) ) ;
    bj_update = zeros( size( bj ) ) ;
    A_update = zeros( size( A ) ) ;
    B_update = zeros( size( B ) ) ;
    
end % end of restart if

% Main loop
for epoch = epoch : num_epochs

  errsum = 0 ; % keep a running total of the difference between data and recon
  
  for batch = 1 : num_batches     
   
%%%%%%%%% START POSITIVE PHASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%      
    num_cases = length( minibatch{ batch } ) ;
    mb = minibatch{ batch } ; % caches the indices   

    % data is a nt+1-d array with current and delayed data
    % corresponding to this mini-batch
    data = zeros( num_cases , num_dims , nt + 1 ) ;
    data( : , : , 1 ) = batchdata( mb , : ) ;
    for hh = 1 : nt
        data( : , : , hh + 1 ) = batchdata( mb - hh , : ) ;
    end
     
    % calculate contributions from directed autoregressive connections
    bi_star = zeros( num_dims , num_cases ) ; 
    for hh = 1 : nt       
        bi_star = bi_star +  A( : , : , hh ) * data( : , : , hh + 1 )' ;
    end   
    
    % Calculate contributions from directed visible-to-hidden connections
    bj_star = zeros( num_hid , num_cases ) ;
    for hh = 1:nt
        bj_star = bj_star + B( : , : , hh ) * data( : , : , hh + 1 )' ;
    end
      
    % Calculate "posterior" probability -- hidden state being on 
    % N ote that it isn't a true posterior   
    eta = w * ( data( : , : , 1 ) ./ gsd )' + ... % bottom-up connections
          repmat( bj , 1 , num_cases ) + ...      % static biases on unit
          bj_star ;                               % dynamic biases
    
    h_posteriors = 1 ./ ( 1 + exp( -eta ) ) ;    % logistic
    
    % Activate the hidden units    
    hid_states = double( h_posteriors' > rand( num_cases , num_hid ) ) ; 
    
    % Calculate positive gradients (note w.r.t. neg energy)
    w_grad = hid_states' * ( data( : , : , 1 ) ./ gsd ) ;
    bi_grad = sum( data( : , : , 1 )' - repmat( bi , 1 , num_cases ) - bi_star , 2 ) ./ gsd^2 ;
    bj_grad = sum( hid_states , 1 )' ;
           
    for hh = 1 : nt      
        A_grad( : , : , hh ) = ( data( : , : , 1 )' - repmat( bi , 1 , num_cases ) - bi_star ) ./ gsd^2 * data( : , : , hh + 1 ) ;
        B_grad( : , : , hh ) = hid_states' * data( : , : , hh + 1 ) ;      
    end
    
%%%%%%%%% END OF POSITIVE PHASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
    % Activate the visible units
    % Find the mean of the Gaussian
    topdown = gsd .* ( hid_states * w ) ;

    % This is the mean of the Gaussian 
    % Instead of properly sampling, negdata is just the mean
    % If we want to sample from the Gaussian, we would add in
    % gsd .* randn( num_cases , num_dims ) ;
    negdata = topdown + ...                       % top down connections
              repmat( bi' , num_cases , 1 ) + ... % static biases
              bi_star' ;                          % dynamic biases

    %N ow conditional on negdata, calculate "posterior" probability
    % for hiddens
    eta = w * ( negdata ./ gsd )' + ...      % bottom-up connections
          repmat( bj , 1 , num_cases ) + ... % static biases on unit (no change)
          bj_star ;                          % dynamic biases (no change)

    h_posteriors = 1 ./ ( 1 + exp( -eta ) ) ; % logistic
    
    % Calculate negative gradients
    neg_w_grad = h_posteriors * ( negdata ./ gsd ) ; % not using activations
    neg_bi_grad = sum( negdata' - repmat( bi , 1 , num_cases ) - bi_star , 2 ) ./ gsd^2 ;
    neg_bj_grad = sum( h_posteriors , 2 ) ;
    
    for hh = 1 : nt
        neg_A_grad( : , : , hh ) = ( negdata' - repmat( bi , 1 , num_cases ) - bi_star ) ./ gsd^2 * data( : , : , hh + 1 ) ;
        neg_B_grad( : , : , hh ) = h_posteriors * data( : , : , hh + 1 ) ;        
    end

%%%%%%%%% END NEGATIVE PHASE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    
    err = sum( sum( ( data( : , : , 1 ) - negdata ) .^2 ) ) ;
    errsum = err + errsum ;
    
    if ( epoch > 5 ) % use momentum
        momentum = mom ;
    else % no momentum
        momentum = 0 ;
    end
    
%%%%%%%%% UPDATE WEIGHTS AND BIASES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
    
    w_update = momentum * w_update + epsilon_w * ( ( w_grad - neg_w_grad ) / num_cases - w_decay * w ) ;
    bi_update = momentum * bi_update + ( epsilon_bi / num_cases ) * ( bi_grad - neg_bi_grad ) ;
    bj_update = momentum * bj_update + ( epsilon_bj / num_cases ) * ( bj_grad - neg_bj_grad ) ;

    for hh = 1 : nt
        A_update( : , : , hh ) = momentum * A_update( : , : , hh ) + epsilon_A * ( ( A_grad( : , : , hh ) - neg_A_grad( : , : , hh ) ) / num_cases - w_decay * A( : , : , hh ) ) ;
        B_update( : , : , hh ) = momentum * B_update( : , : , hh ) + epsilon_B * ( ( B_grad( : , : , hh ) - neg_B_grad( : , : , hh ) ) / num_cases - w_decay * B( : , : , hh ) ) ;
    end

    w = w + w_update ;
    bi = bi + bi_update ;
    bj = bj + bj_update ;

    for hh = 1 : nt
        A( : , : , hh ) = A( : , : , hh ) + A_update( : , : , hh ) ;
        B( : , : , hh ) = B( : , : , hh ) + B_update( : , : , hh ) ;
    end
    
%%%%%%%%%%%%%%%% END OF UPDATES  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
  
  end
    
%    % every 10 epochs, show output
%    if ( mod( epoch , 10 ) == 0 )
%       fprintf( 1 , 'epoch %4i error %6.1f  \n', epoch , errsum ) ; 
%       % Could see a plot of the weights every 10 epochs
%       % figure(3); weightreport
%       % drawnow;
%    end
    
end % end of epoch for loop
    
end % end of function

get_filtering_distance.m

function [ filtering_distance ] = get_filtering_distance( batchdata , batchdataindex , num_hid , w , bj , n1 , n2 , B , gsd )

% Version 1.01 
%
% Code provided by Graham Taylor, Geoff Hinton and Sam Roweis 
%
% For more information, see:
%     http://www.cs.toronto.edu/~gwtaylor/publications/nips2006mhmublv
%
% Permission is granted for anyone to copy, use, modify, or distribute this
% program and accompanying programs and documents for any purpose, provided
% this copyright notice is retained and prominently displayed, along with
% a note saying that the original programs are available from our
% web page.
% The programs and documents are distributed without any warranty, express or
% implied.  As the programs were written for research purposes only, they have
% not been tested to the degree that would be advisable in any important
% application.  All use of these programs is entirely at the user's own risk.

% This program runs the entire training set on the trained 1-hidden-layer
% network and saves the filtering distribution vectors in mini-batch format
% This is done before training a second CRBM on top of the first

% The program assumes that the following variables are set externally:
% n2    -- order of the next layer CRBM

%  batchsize = 100 ; % size of mini-batches

%  take all valid examples ( indexed by batchdataindex )
num_cases = length( batchdataindex ) ;

%  Calculate contributions from directed visible-to-hidden connections
bj_star = zeros( num_hid , num_cases ) ;

for hh = 1 : n1

  bj_star = bj_star + B( : , : , hh ) * batchdata( batchdataindex - hh , : )' ;
  
end % end for loop

%  Calculate "posterior" probability -- hidden state being on
%  Note that it isn't a true posterior
bottom_up = w * ( batchdata( batchdataindex , : ) ./ gsd )' ;

eta = bottom_up + ...                    % bottom-up connections
      repmat( bj , 1 , num_cases ) + ... % static biases on unit
      bj_star ;                          % dynamic biases

filtering_distance = 1 ./ ( 1 + exp( -eta' ) ) ; % logistic

end % end of function

binary_crbm.m

function [ w , bj , bi , A , B ] = binary_crbm( input_data , minibatch , nt , num_epochs , num_hid , restart )

% Version 1.01 
%
% Code provided by Graham Taylor, Geoff Hinton and Sam Roweis 
%
% For more information, see:
%     http://www.cs.toronto.edu/~gwtaylor/publications/nips2006mhmublv
%
% Permission is granted for anyone to copy, use, modify, or distribute this
% program and accompanying programs and documents for any purpose, provided
% this copyright notice is retained and prominently displayed, along with
% a note saying that the original programs are available from our
% web page.
% The programs and documents are distributed without any warranty, express or
% implied.  As the programs were written for research purposes only, they have
% not been tested to the degree that would be advisable in any important
% application.  All use of these programs is entirely at the user's own risk.

% This program trains a Conditional Restricted Boltzmann Machine in which
% visible, binary, stochastic inputs are connected to
% hidden, binary, stochastic feature detectors using symmetrically
% weighted connections. Learning is done with 1-step Contrastive Divergence.
% Directed connections are present, from the past nt configurations of the
% visible units to the current visible units (A), and the past nt
% configurations of the visible units to the current hidden units (B)

% The program assumes that the following variables are set externally:
% nt        -- order of the model
% num_epochs -- maximum number of epochs
% num_hid    --  number of hidden units 
% input_data --  a matrix of data (num_cases,num_dims) 
% minibatch -- a cell array of dimension batchsize, indexing the valid
% frames in input_data
% restart   -- set to 1 if learning starts from beginning 

% input_data is a big matrix of all the frames
% we index it with "minibatch", a cell array of mini-batch indices
num_batches = length( minibatch ) ; 

num_dims = size( input_data , 2 ) ; % visible dimension

%  Setting learning rates
epsilon_w = 1e-3 ; % undirected
epsilon_bi = 1e-3 ; % visibles
epsilon_bj = 1e-3 ; % hidden units
epsilon_A = 1e-3 ; % autoregressive
epsilon_B = 1e-3 ; % prev visibles to hidden

w_decay = 0.0002 ; % currently we use the same weight decay for w, A, B
mom = 0.9 ;        % momentum used only after 5 epochs of training

if ( restart == 1 )  
   restart = 0 ;
   epoch = 1 ;
  
  % Randomly initialize weights
  w = 0.01 * randn( num_hid , num_dims ) ;
  bi = 0.01 * randn( num_dims , 1 ) ;
  bj = -1 + 0.01 * randn( num_hid , 1 ) ; % set to favour units being "off"
  
  % The autoregressive weights; A(:,:,j) is the weight from t-j to the vis
  A = 0.01 * randn( num_dims , num_dims , nt ) ;
 
  % The weights from previous time-steps to the hiddens; B(:,:,j) is the
  % weight from t-j to the hidden layer
  B = 0.01 * randn( num_hid , num_dims , nt ) ;
  
  clear w_grad bi_grad bj_grad A_grad B_grad
  clear neg_w_grad neg_bi_grad neg_bj_grad neg_A_grad neg_B_grad
  
  % keep previous updates around for momentum
  w_update = zeros( size( w ) ) ;
  bi_update = zeros( size( bi ) ) ;
  bj_update = zeros( size( bj ) ) ;
  A_update = zeros( size( A ) ) ;
  B_update = zeros( size( B ) ) ;
    
end % end of restart if

% Main loop
for epoch = epoch : num_epochs

  errsum = 0 ; % keep a running total of the difference between data and recon
  
  for batch = 1 : num_batches     

%%%%%%%%% START POSITIVE PHASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    num_cases = length( minibatch{ batch } ) ;
    mb = minibatch{ batch } ; % caches the indices   

    % data is a nt+1-d array with current and delayed data
    % corresponding to this mini-batch
    data = zeros( num_cases , num_dims , nt + 1 ) ;
    data( : , : , 1 ) = input_data( mb , : ) ;
    for hh = 1 : nt
        data( : , : , hh + 1 ) = input_data( mb - hh , : ) ;
    end

    % Calculate contributions from directed autoregressive connections
    bi_star = zeros( num_dims , num_cases ) ; 
    for hh = 1 : nt       
        bi_star = bi_star +  A( : , : , hh ) * data( : , : , hh + 1 )' ;
    end   
    
    % Calculate contributions from directed visible-to-hidden connections
    bj_star = zeros( num_hid , num_cases ) ;
    for hh = 1 : nt
        bj_star = bj_star + B( : , : , hh ) * data( : , : , hh + 1 )' ;
    end
     
    bottom_up = w * data( : , : , 1 )' ;
    
    % Calculate "posterior" probability -- hidden state being on 
    % Note that it isn't a true posterior   
    eta = bottom_up + ...                    % bottom-up connections
          repmat( bj , 1 , num_cases ) + ... % static biases on unit
          bj_star ;                          % dynamic biases
    
    h_posteriors = 1 ./ ( 1 + exp( -eta ) ) ; % logistic
       
    % Activate the hidden units    
    hid_states = h_posteriors' > rand( num_cases , num_hid ) ; 
    
    % Calculate positive gradients (note w.r.t. neg energy)
    w_grad = hid_states' * data( : , : , 1 ) ;
    bi_grad = sum( data( : , : , 1 )' - repmat( bi , 1 , num_cases ) - bi_star , 2 ) ;
    bj_grad = sum( hid_states , 1 )' ;
           
    for hh = 1 : nt      
        A_grad( : , : , hh ) = ( data( : , : , 1 )' - repmat( bi , 1 , num_cases ) - bi_star ) * data( : , : , hh + 1 ) ;
        B_grad( : , : , hh ) = hid_states' * data( : , : , hh + 1 ) ;      
    end
    
%%%%%%%%% END OF POSITIVE PHASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%% START NEGATIVE PHASE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    
    % Activate the visible units
    topdown = hid_states * w ;

    eta = topdown + ...                       % top down connections
          repmat( bi' , num_cases , 1 ) + ... % static biases
          bi_star' ;                          % dynamic biases
    
    negdata = 1 ./ ( 1 + exp( -eta ) ) ; % logistic
    
    % Now conditional on negdata, calculate "posterior" probability
    % for hiddens
    bottom_up = w * negdata' ;
    eta = bottom_up + ...                    % bottom-up connections
          repmat( bj , 1 , num_cases ) + ... % static biases on unit (no change)
          bj_star ;                          % dynamic biases (no change)

    h_posteriors = 1 ./ ( 1 + exp( -eta ) ) ; %logistic
    
    % Calculate negative gradients
    neg_w_grad = h_posteriors * negdata ; % not using activations
    neg_bi_grad = sum( negdata' - repmat( bi , 1 , num_cases ) - bi_star , 2 ) ;
    neg_bj_grad = sum( h_posteriors , 2 ) ;
    
    for hh = 1 : nt
        neg_A_grad( : , : , hh ) = ( negdata' - repmat( bi , 1 , num_cases ) - bi_star ) * data( : , : , hh + 1 ) ;
        neg_B_grad( : , : , hh ) = h_posteriors * data( : , : , hh + 1 ) ;        
    end
   
%%%%%%%%% END NEGATIVE PHASE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    err = sum( sum( ( data( : , : , 1 ) - negdata ) .^2 ) ) ;
    errsum = err + errsum ;
    
    if ( epoch > 5 ) % use momentum
        momentum = mom ;
    else % no momentum
        momentum = 0 ;
    end
    
%%%%%%%%% UPDATE WEIGHTS AND BIASES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    w_update = momentum * w_update + epsilon_w * ( ( w_grad - neg_w_grad ) / num_cases - w_decay * w ) ;
    bi_update = momentum * bi_update + ( epsilon_bi / num_cases ) * ( bi_grad - neg_bi_grad ) ;
    bj_update = momentum * bj_update + ( epsilon_bj / num_cases ) * ( bj_grad - neg_bj_grad ) ;

    for hh = 1 : nt
        A_update( : , : , hh ) = momentum * A_update( : , : , hh ) + epsilon_A * ( ( A_grad( : , : , hh ) - neg_A_grad( : , : , hh ) ) / num_cases - w_decay * A( : , : , hh ) ) ;
        B_update( : , : , hh ) = momentum * B_update( : , : , hh ) + epsilon_B * ( ( B_grad( : , : , hh ) - neg_B_grad( : , : , hh ) ) / num_cases - w_decay * B( : , : , hh ) ) ;
    end

    w = w +  w_update ;
    bi = bi + bi_update ;
    bj = bj + bj_update ;

    for hh = 1 : nt
        A( : , : , hh ) = A( : , : , hh ) + A_update( : , : , hh ) ;
        B( : , : , hh ) = B( : , : , hh ) + B_update( : , : , hh ) ;
    end
    
%%%%%%%%%%%%%%%% END OF UPDATES  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
  end
    
%    % every 10 epochs, show output
%    if ( mod( epoch , 10 ) == 0 )
%        fprintf( 1 , 'epoch %4i error %6.1f  \n' , epoch , errsum ) ; 
%        % Could see a plot of the weights every 10 epochs
%        % figure(3); weightreport
%        % drawnow;
%    end
    
end % end main loop
    
end % end of function

generates new data with generate_visible_data.m

function [ visible_data  ] = generate_visible_data( init_data , hidden_1 , num_hid_1 , num_hid_2 , w1 , w2 , bj1 , bi1 , bj2 , bi2 , A1 , A2 , B1 , B2 , n1 , n2 , gsd )

% Version 1.01 
%
% Code provided by Graham Taylor, Geoff Hinton and Sam Roweis 
%
% For more information, see:
%     http://www.cs.toronto.edu/~gwtaylor/publications/nips2006mhmublv
%
% Permission is granted for anyone to copy, use, modify, or distribute this
% program and accompanying programs and documents for any purpose, provided
% this copyright notice is retained and prominently displayed, along with
% a note saying that the original programs are available from our
% web page.
% The programs and documents are distributed without any warranty, express or
% implied.  As the programs were written for research purposes only, they have
% not been tested to the degree that would be advisable in any important
% application.  All use of these programs is entirely at the user's own risk.
%
% This program uses the 2-level CRBM to generate data
% More efficient version than the original

num_frames = 1 ;
max_clamped = n1 + n2 ;
A2_flat = reshape( A2 , num_hid_1 , n2 * num_hid_1 ) ;
B2_flat = reshape( B2 , num_hid_2 , n2 * num_hid_1 ) ;

num_Gibbs = 30 ; % number of alternating Gibbs iterations

%  initialize second layer ( first n1 + n2 frames padded )
hidden_2 = ones( num_frames , num_hid_2 ) ;

%  keep the recent past in vector form
%  for input to directed links
%  it's slightly faster to pre-compute this vector and update it ( by
%  shifting ) at each time frame, rather than to fully recompute each time
%  frame
past = reshape( hidden_1( max_clamped : -1 : max_clamped + 1 - n2 , : )' , num_hid_1 * n2 , 1 ) ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  First generate a hidden sequence (top layer)
%  Then go down through the first CRBM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf( 'Generating hidden states\n' ) ;

%  for tt = max_clamped + 1 : num_frames
  tt = size( hidden_1 , 1 ) ;

  % initialize using the last frame
  % noise is not added for binary units
  hidden_1( tt , : ) = hidden_1( tt - 1 , : ) ;
  
  % Dynamic biases aren't re-calculated during Alternating Gibbs
  bi_star = A2_flat * past ;
  bj_star = B2_flat * past ;

  % Gibbs sampling
  for gg = 1 : num_Gibbs
  
    % Calculate posterior probability -- hidden state being on (estimate)
    % add in bias
    bottomup = w2 * hidden_1( tt , : )' ;
               eta = bottomup + ...       % bottom-up connections
               bj2 + ...                  % static biases on unit
               bj_star ;                  % dynamic biases
    
    hposteriors = 1 ./ ( 1 + exp( -eta ) ) ; % logistic
    
    hidden_2( tt , : ) = double( hposteriors' > rand( 1 , num_hid_2 ) ) ; % Activating hiddens
    
    % Downward pass; visibles are binary logistic units     
    topdown = hidden_2( tt , : ) * w2 ;
        
    eta = topdown + ... % top down connections
          bi2' + ...    % static biases
          bi_star';     % dynamic biases   
    
    hidden_1( tt , : ) = 1 ./ ( 1 + exp( -eta ) ) ; % logistic 
    
  end % end gibbs sampling

  % If we are done Gibbs sampling, then do a mean-field sample
  
  topdown = hposteriors' * w2 ; % Very noisy if we don't do mean-field here

  eta = topdown + ... % top down connections
        bi2' + ...    % static biases
        bi_star';     % dynamic biases
        hidden_1( tt , : ) = 1 ./ ( 1 + exp( -eta ) ) ;

  % update the past
%    past( num_hid_1 + 1 : end ) = past( 1 : end - num_hid_1 ) ; % shift older history down
%    past( 1 : num_hid_1 ) = hidden_1( tt , : ) ;                % place most recent frame at top
  
  
%    if ( mod( tt , 10 ) == 0 )
%       fprintf( 'Finished frame %d\n' , tt ) ;
%    end
  
%  end % end tt loop

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Now that we've decided on the "filtering distribution", generate visible
%  data through CRBM 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf( 'Generating visible data\n' ) ;

%  for tt = max_clamped + 1 : num_frames

    % Add contributions from autoregressive connections
%      bi_star = zeros( num_dims , 1 ) ; % original
    bi_star = zeros( size( init_data , 2 ) , 1 ) ;
    for hh = 1 : n1
%          bi_star = bi_star + A1( : , : , hh ) * visible( tt - hh , : )' ; % original
        bi_star = bi_star + A1( : , : , hh ) * init_data( tt - hh , : )' ;
    end

    % Mean-field approx; visible units are Gaussian
    % ( filtering distribution is the data we just generated )
    topdown = gsd .* ( hidden_1( tt , : ) * w1 ) ;
%      visible( tt , : ) = topdown + ... % top down connections ( original )
    visible_data = topdown + ...
                   bi1' + ...    % static biases
                   bi_star' ;    % dynamic biases
                        
%  end % end tt loop

end % end of function

and then finally plots this new data as price bars, crbm_visual_test.m

clear all

% load price file of interest
filename = input( 'Enter filename for prices, e.g. es or esmatrix: ' , 's' ) ;
data = load( "-ascii" , filename ) ;

% get tick size
switch filename

case { "cc" }
tick = 1 ;

case { "gc" "lb" "pl" "sm" "sp" }
tick = 0.1 ;

case { "ausyen" "bo" "cl" "ct" "dx" "euryen" "gbpyen" "sb" "usdyen" }
tick = 0.01 ;

case { "ng" }
tick = 0.001 ;

case { "auscad" "aususd" "euraus" "eurcad" "eurchf" "eurgbp" "eurusd" "gbpchf" "gbpusd" "ho" "rb" "usdcad" "usdchf" }
tick = 0.0001 ;

case { "c" "o" "s" "es" "nd" "w" }
tick = 0.25 ;

case { "fc" "lc" "lh" "pb" }
tick = 0.025 ;

case { "ed" }
tick = 0.0025 ;

case { "si" }
tick = 0.5 ;

case { "hg" "kc" "oj" "pa" }
tick = 0.05 ;

case { "ty" "us" }
tick = 0.015625 ;

case { "ccmatrix" }
tick = 1 ;

case { "gcmatrix" "lbmatrix" "plmatrix" "smmatrix" "spmatrix" }
tick = 0.1 ;

case { "ausyenmatrix" "bomatrix" "clmatrix" "ctmatrix" "dxmatrix" "euryenmatrix" "gbpyenmatrix" "sbmatrix" "usdyenmatrix" }
tick = 0.01 ;

case { "ngmatrix" }
tick = 0.001 ;

case { "auscadmatrix" "aususdmatrix" "eurausmatrix" "eurcadmatrix" "eurchfmatrix" "eurgbpmatrix" "eurusdmatrix" "gbpchfmatrix" "gbpusdmatrix" "homatrix" "rbmatrix" "usdcadmatrix" "usdchfmatrix" }
tick = 0.0001 ;

case { "cmatrix" "omatrix" "smatrix" "esmatrix" "ndmatrix" "wmatrix" }
tick = 0.25 ;

case { "fcmatrix" "lcmatrix" "lhmatrix" "pbmatrix" }
tick = 0.025 ;

case { "edmatrix" }
tick = 0.0025 ;

case { "simatrix" }
tick = 0.5 ;

case { "hgmatrix" "kcmatrix" "ojmatrix" "pamatrix" }
tick = 0.05 ;

case { "tymatrix" "usmatrix" }
tick = 0.015625 ;

endswitch

%  get data
open = data( : , 4 ) ;
high = data( : , 5 ) ;
low = data( : , 6 ) ;
close = data( : , 7 ) ;

%  clear data

%  create log series for crbm training
vwap = vwap_all_period_detrend( open , high , low , close , tick ) ;
open_from_vwap = log( open ./ vwap ) ; open_from_vwap(1) = 0 ;
open_from_previous_close = log( open ./ shift( close , 1 ) ) ; open_from_previous_close(1) = 0 ;
high_from_vwap = log( high ./ vwap ) ; high_from_vwap(1) = 0 ;
low_from_vwap = log( low ./ vwap ) ; low_from_vwap(1) = 0 ;
close_from_vwap = log( close ./ vwap ) ; close_from_vwap(1) = 0 ;

%  assemble into a features matrix
features = [ open_from_vwap open_from_previous_close high_from_vwap low_from_vwap close_from_vwap ] ;

training_point_index = input( 'The training point index for NN training from figure 1: ' ) ;

%  load the weights etc.
load layer1.mat ;
load layer1_2.mat ;
load layer1_3.mat ;
load layer2.mat ;
load layer2_2.mat ;
load layer2_3.mat ;
load bar_1_mean_std ;
load bar_2_mean_std ;
load bar_3_mean_std ;

%  how-many timesteps do we look back for directed connections - this is what we call the "order" of the model 
n1 = 3 ; % first layer
n2 = 3 ; % second layer

%  Set network properties
numhid1 = 100 ; numhid2 = 100 ; % numepochs = 2000 ;
gsd = 1 ;          % fixed standard deviation for Gaussian units

%  the first bar
init_data = features( training_point_index - 5 : training_point_index , : ) ;

%  which must be normalised using the same data_mean and data_std used to normalise the training data
init_data = ( init_data - repmat( data_mean_1 , size( init_data , 1 ) , 1 ) ) ./ repmat( data_std_1 , size( init_data , 1 ) , 1 ) ; 

%  and we also need the 3 fixed hidden units of the crbm's first hidden layer
[ hidden_1 ] = get_fixed_hidden_layer( init_data , numhid1 , w11 , bj11 , B11 , n1 , n2 , gsd ) ;

%  do Gibbs sampling and generate the visible layer from the above, clamped, init_data and hidden_1 layers
[ visible_data ] = generate_visible_data( init_data , hidden_1 , numhid1 , numhid2 , w11 , w21 , bj11 , bi11 , bj21 , bi21 , A11 , A21 , B11 , B21 , n1 , n2 , gsd ) ;

%  now post process the visible_data - scale back to original using data_mean & data_std from above
visible_data_1 = data_std_1 .* visible_data .+ data_mean_1 ;

%  repeat for 2nd bar
init_data = [ features( training_point_index - 4 : training_point_index , : ) ; visible_data_1 ] ;

%  which must be normalised using the same data_mean and data_std used to normalise the training data
init_data = ( init_data - repmat( data_mean_2 , size( init_data , 1 ) , 1 ) ) ./ repmat( data_std_2 , size( init_data , 1 ) , 1 ) ; 

%  and we also need the 3 fixed hidden units of the crbm's first hidden layer
[ hidden_1 ] = get_fixed_hidden_layer( init_data , numhid1 , w12 , bj12 , B12 , n1 , n2 , gsd ) ;

%  do Gibbs sampling and generate the visible layer from the above, clamped, init_data and hidden_1 layers
[ visible_data ] = generate_visible_data( init_data , hidden_1 , numhid1 , numhid2 , w12 , w22 , bj12 , bi12 , bj22 , bi22 , A12 , A22 , B12 , B22 , n1 , n2 , gsd ) ;

%  now post process the visible_data - scale back to original using data_mean & data_std from above
visible_data_2 = data_std_2 .* visible_data .+ data_mean_2 ;

%  repeat for 3rd bar
init_data = [ features( training_point_index - 3 : training_point_index , : ) ; visible_data_1 ; visible_data_2 ] ;

%  which must be normalised using the same data_mean and data_std used to normalise the training data
init_data = ( init_data - repmat( data_mean_3 , size( init_data , 1 ) , 1 ) ) ./ repmat( data_std_3 , size( init_data , 1 ) , 1 ) ; 

%  and we also need the 3 fixed hidden units of the crbm's first hidden layer
[ hidden_1 ] = get_fixed_hidden_layer( init_data , numhid1 , w13 , bj13 , B13 , n1 , n2 , gsd ) ;

%  do Gibbs sampling and generate the visible layer from the above, clamped, init_data and hidden_1 layers
[ visible_data ] = generate_visible_data( init_data , hidden_1 , numhid1 , numhid2 , w13 , w23 , bj13 , bi13 , bj23 , bi23 , A13 , A23 , B13 , B23 , n1 , n2 , gsd ) ;

%  now post process the visible_data - scale back to original using data_mean & data_std from above
visible_data_3 = data_std_3 .* visible_data .+ data_mean_3 ;

%  create new bars for plotting
new_opens = zeros( 3 , 1 ) ; new_highs = zeros( 3 , 1 ) ; new_lows = zeros( 3 , 1 ) ; new_closes = zeros( 3 , 1 ) ;

visible_data_1 = exp( visible_data_1 ) ; visible_data_2 = exp( visible_data_2 ) ; visible_data_3 = exp( visible_data_3 ) ;

new_opens( 1 , 1 ) = ( visible_data_1( 1 , 1 ) * vwap(training_point_index) + visible_data_1( 1 , 2 ) * close(training_point_index) ) / 2 ;
new_highs( 1 , 1 ) = visible_data_1( 1 , 3 ) * vwap(training_point_index) ;
new_lows( 1 , 1 ) = visible_data_1( 1 , 4 ) * vwap(training_point_index) ;
new_closes( 1 , 1 ) = visible_data_1( 1 , 5 ) * vwap(training_point_index) ;
vwap_1 = ( new_opens(1,1) + new_closes(1,1) + ( new_highs(1,1) + new_lows(1,1) )/2 ) / 3 ;

new_opens( 2 , 1 ) = ( visible_data_2( 1 , 1 ) * vwap_1 + visible_data_2( 1 , 2 ) * new_closes(1,1) ) / 2 ;
new_highs( 2 , 1 ) = visible_data_2( 1 , 3 ) * vwap_1 ;
new_lows( 2 , 1 ) = visible_data_2( 1 , 4 ) * vwap_1 ;
new_closes( 2 , 1 ) = visible_data_2( 1 , 5 ) * vwap_1 ;
vwap_2 = ( new_opens(2,1) + new_closes(2,1) + ( new_highs(2,1) + new_lows(2,1) )/2 ) / 3 ;

new_opens( 3 , 1 ) = ( visible_data_3( 1 , 1 ) * vwap_2 + visible_data_3( 1 , 2 ) * new_closes(2,1) ) / 2 ;
new_highs( 3 , 1 ) = visible_data_3( 1 , 3 ) * vwap_2 ;
new_lows( 3 , 1 ) = visible_data_3( 1 , 4 ) * vwap_2 ;
new_closes( 3 , 1 ) = visible_data_3( 1 , 5 ) * vwap_2 ;

figure(1) ;
plot_opens = open( training_point_index-20:training_point_index+3 , : ) ; plot_highs = high( training_point_index-20:training_point_index+3 , : ) ; plot_lows = low( training_point_index-20:training_point_index+3 , : ) ; plot_closes = close( training_point_index-20:training_point_index+3 , : ) ;
pkg load financial ;
clf() ;
highlow( plot_highs , plot_lows , plot_closes , plot_opens )

figure(2) ;
plot_opens = [ open( training_point_index-20:training_point_index , : ) ; new_opens ] ; plot_highs = [ high( training_point_index-20:training_point_index , : ) ; new_highs ] ; plot_lows = [ low( training_point_index-20:training_point_index , : ) ; new_lows ] ; plot_closes = [ close( training_point_index-20:training_point_index , : ) ; new_closes ] ; ;
clf() ;
highlow( plot_highs , plot_lows , plot_closes , plot_opens )

which produces this as output

and

where the last 3 bars on the right of the second chart are the modelled bars of the 3 rightmost bars in the upper chart. At the moment the modelling isn't very good because I'm using simplistic features - this is just proof of concept coding. Future work to do will include finding much better input features, as well as speeding up the running time of the code by compiling into .oct files, and then optimising parameters using my particle swarm optimisation code.

Particle Swarm Optimisation, Part 2.

Following on from my last post, here is an Octave .oct file implementation of the one dimensional Particle swarm optimisation routine, with one slight twist: instead of using a for loop I've implemented it within a while loop with a stopping condition that the algorithm should cease once there has been no improvement in the global_best value for 25 iterations.

DEFUN_DLD ( pso_conversion_code, args, nargout,
"-*- texinfo -*-\n\
@deftypefn {Function File} {} pso_conversion_code (@var{target_val,max_lambda})\n\
The output of this test function should be half the input target_val.\n\
@end deftypefn" )

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

// check the number of input arguments
if ( nargin != 2 )
{
    error ( "Invalid number of arguments." ) ;
    return retval_list ;
}

if ( args(0).length () != 1 )
{
error ( "Invalid target_val. Should be a single value." ) ;
return retval_list ;
}

if ( args(1).length () != 1 )
{
error ( "Invalid max_lambda value. Should be a single value for the maximum 'guess'." ) ;
return retval_list ;
}

if (error_state)
{
    error ( "Invalid arguments." ) ;
    return retval_list ;
}
// end of input checking

double target_val = args(0).double_value() ;
double max_lambda = args(1).double_value() ;

double loocv_value ;

// the the pso algorithm
int no_iterations_until_cease = 25 ;
int no_of_particles = 100 ;

double global_best = std::numeric_limits::infinity() ;
double global_best_lambda = 0.0 ; // initially set to unregularised

ColumnVector local_best( no_of_particles , 1 ) ; local_best.fill( global_best ) ;
ColumnVector local( no_of_particles , 1 ) ;
ColumnVector local_best_so_far( no_of_particles , 1 ) ;

ColumnVector velocity( no_of_particles , 1 ) ; velocity.fill( 0.0 ) ; // particle velocity vector

// A Mersenne Twister random number generator can be declared with a simple from the #include "MersenneTwister.h" 
MTRand mtrand1 ;

// values for the random updating process
double r1 ; 
double r2 ;

// an inertial constant. Good values are usually slightly less than 1. Or it could be randomly initialized for each particle.
double w_ic = 0.9 ;

// c1 and c2 are constants that say how much the particle is directed towards good positions. They represent a "cognitive" and a "social" component, respectively, 
// in that they affect how much the particle's personal best and the global best (respectively) influence its movement. Usually we take c_1, c_2  approx = 2. 
// Or they could be randomly initialized for each particle.
double c1 = 2.0 ;
double c2 = 2.0 ; 

// fill the local vector with initial random values < max_lambda, temporarily using r1
for ( octave_idx_type ii (0) ; ii < no_of_particles ; ii++ )
    {
    r1 = mtrand1.randDblExc () ;
    local(ii) = r1 * max_lambda ;
    }

int while_counter = 0 ;
while ( while_counter < no_iterations_until_cease )
{
  
  // loop once over local_best and local vectors
  for ( octave_idx_type jj (0) ; jj < no_of_particles ; jj++ )
      {
 
          // Replace this code box as necessary 
 
          //*************************************************************//
   //                                                             // 
          // fitness function evaluation                                 //
   loocv_value = local(jj) * local(jj) - target_val * local(jj) ; //
   //                                                             //
   //*************************************************************//

   // check if the local_best has improved
   if ( loocv_value < local_best(jj) )
      {
      // update local_best and local_best_so_far vector if it has
      local_best(jj) = loocv_value ;
      local_best_so_far(jj) = local(jj) ;
      }
     
   // check if the above local_best has also improved the global_best  
   if ( local_best(jj) < global_best )
      {
      // update global_best and global_best_lambda if it has
      global_best = local_best(jj) ;
      global_best_lambda = local_best_so_far(jj) ;
      while_counter = 0 ;
      }
    
      } // end of loop once over local_best and local vectors
      
  // now update the particle velocity and position vectors
  for ( octave_idx_type jj (0) ; jj < no_of_particles ; jj++ )
      {
      r1 = mtrand1.randDblExc () ;
      r2 = mtrand1.randDblExc () ;
      velocity(jj) = w_ic * velocity(jj) + c1 * r1 * ( local_best_so_far(jj) - local(jj) ) + c2 * r2 * ( global_best - local(jj) ) ;
      local(jj) = local(jj) + velocity(jj) ;
      } // end of particle velocity and position vectors updates loop
 
while_counter += 1 ;
 
} // end of main while loop 

retval_list(1) = global_best ;
retval_list(0) = global_best_lambda ;

return retval_list ;
  
} // end of function

In the "commented" function section there is the same test function as in the vectorised code in my previous post. In real life application, of course, this code would be replaced - the above is just a test of my conversion of the pso algorithm from Octave to C++ code.

I've been looking at pso because I think I can easily use it as a simple Hyperparameter optimisation tool to tune the regularisation of the weights in my proposed neural net trading system. The reason I chose the fitness function I did in the above code is simply that it has a global minimum, which is what neural net training is all about - minimisation of an error function. 

Particle Swarm Optimisation

Having decided that I'm going to use my mfe_mae indicator as a target for my neural net, over the last couple of months I've been doing some research on what would make good features for this target. In the course of this I've also decided that Particle swarm optimization would be a useful tool to use.

Thanks to the pseudo-code on this wiki and the code on in this stackoverflow thread I've been able to write some Octave code to perform pso over one dimension, which is shown in the code box below:

clear all

function [out] = evaluate_fitness( R , val )
    out = R .* R .- val .* R ;
end

val = input( 'Enter test val: ' ) ;

% Initialize the particle positions and their velocities
n_particles = 100 ;
upper_limit = val ;
n_iterations = 50 ;

w = 0.9 .* ones( n_particles , 1 ) ; % is an inertial constant. Good values are usually slightly less than 1. Or it could be randomly initialized for each particle.

% c1 and c2 are constants that say how much the particle is directed towards good positions. They represent a "cognitive" and a "social" component, respectively, 
% in that they affect how much the particle's personal best and the global best (respectively) influence its movement. Usually we take c_1, c_2  approx = 2. 
% Or they could be randomly initialized for each particle.
C1 = 2.0 ;
C2 = 2.0 ;

X = upper_limit * rand( n_particles , 1 ) ; % == particle position vector containing lambda values
V = zeros( size(X) ) ; % particle velocity vector
X_lbest = X ; % == particle position vector containing lambda values for local best
 
% Initialize the global and local fitness to the worst possible. Fitness == LOOCV "press" statistic
fitness_gbest = Inf ; % _gbest == global best
fitness_lbest = fitness_gbest * ones( n_particles , 1 ) ; % _lbest == local best
 
% Loop until convergence, in this example a finite number of iterations chosen
for ii = 1 : n_iterations

    % evaluate the fitness of each particle, i.e. do the linear regression and get
    % the LOOCV statistic
    fitness_X = evaluate_fitness( X , val ) ;
 
    % Update the local bests and their fitness 
    ix = find( fitness_X < fitness_lbest ) ; % if local LOOCV "press" statistic improves
    fitness_lbest( ix ) = fitness_X( ix ) ;  % record this better LOOCV statistic value 
    X_lbest( ix ) = X( ix ) ;                % and the lambda value that gave rise to it    
 
    % Update the global best and its fitness 
    [ min_fitness min_fitness_index ] = min( fitness_X ) ;
    
    if ( min_fitness < fitness_gbest )       % if global LOOCV "press" statistic improves
        fitness_gbest = min_fitness ;        % record this better LOOCV statistic value
        X_gbest = X( min_fitness_index ) ;   % and the lambda value that gave rise to it
    end % end if    
 
    % Update the particle velocity and position vectors
    R1 = rand( n_particles , 1 ) ;
    R2 = rand( n_particles , 1 ) ;
    V = w .* V + C1 .* R1 .* ( X_lbest .- X ) .+ C2 .* R2 .* ( X_gbest .* ones( n_particles , 1 ) .- X ) ;
    X = X .+ V ;
    
end % end main ii loop

and which is vectorised as much as I can make it at the moment. The evaluate_fitness function I've chosen to use is a Quadratic function of the form 

f(x) = x^2 - bx 

which, when a positive value for "val" is input as the test value, ensures the function curve goes through the origin and has a minimum y-axis value at a point on the x-axis that is half the input test value. This make it easy to quickly verify that the pso code is performing as expected, with the global minimum x_axis value found by the algorithm being given by the variable X_gbest. My reasons for choosing a test function of this form, and for looking at pso in general, will be given in my next post.  

MFE/MAE Indicator Test Results

Following on from the previous post, the test I outlined in that post wasn't very satisfactory, which I put down to the fact that the Sigmoid transformation of the raw MFE/MAE indicator values is not amenable to the application of standard deviation as a meaningful measure. Instead, I changed the test to one based on the standard error of the mean, an example screen shot of which is shown below:-

The top pane shows the the long version of the indicator and the bottom pane the short version. In each there are upper and lower limits of the sample standard error of the mean above and below the population mean (mean of all values of the indicator) along with the cumulative mean value of the top N matches as shown on the x-axis. In this particular example it can be seen that around the 170-180 samples mark the cumulative mean moves inside the standard error limits, never to leave them again. The meaning I ascribe to this is that there is no value to be gained from using more than approximately 180 samples for machine learning purposes, for this example, as to use more samples would be akin to training on all available data, which makes the use of my Cauchy Schwarz matching algo superfluous. I repeated the above on all instances of sigmoid transformed and untransformed MFE/MAE indicator values to get an average of 325 samples for transformed, and an average of 446 samples for the untransformed indicator values across the 4 major forex pairs. Based on this, I have decided to use the top 450 Cauchy Schwarz matches for training purposes, which has ramifications for model complexity will be discussed shortly.

Returning to the above screen shot, the figure 2 inset shows the price bars that immediately follow the price bar for which the main screen shows the top N matches. Looking at the extreme left of the main screen it can be seen that the lower pane, short indicator has an almost maximum reading of 1 whilst the upper pane, long indicator shows a value of approx. 2.7, which is not much above the global minimum for this indicator and well below the 0.5 neutral level. This strongly suggests a short position, and looking at the inset figure it can be seen that over the 3 days following the extreme left matched bar a short position was indeed the best position to hold. This is a pattern that seems to frequently present itself during visual inspection of charts, although I am unable to quantify this in any way.

On the matter of model complexity alluded to above, I found the Learning From Data course I have recently completed on the edX platform to be very enlightening, particularly the concept of the VC dimension, which is nicely explained in the Learning From Data Video library. I'll leave it to interested readers to follow the links, but the big take away for me is that using 450 samples as described above implies that my final machine learning model must have an upper bound of approximately 45 on the VC dimension, which in turn implies a maximum of 45 weights in the neural net. This is a design constraint that I will discuss in a future post. 



  

Test of the MFE/MAE Indicator

Continuing from my last post, wherein I stated I was going to conduct a more pertinent statistical test of the returns of the bars(s) immediately following the N best, Cauchy Schwarz matching algorithm matched bars in the price history, readers may recall that the basic premise behind this algorithm is that by matching current price action to the N best matches, the price action after these matches can be used to infer what will occur after the current price action. However, rather than test the price action directly I have decided to apply the test to the MFE/MAE indicator. There are several reasons for this, which are enumerated below.

1. I intend to use the indicator as the target function for future Neural net training
2. the indicator represents a reward to risk ratio, which indirectly reflects price action itself, but without the noise of said action
3. this reward to risk ratio is of much more direct concern, from a trading perspective, than accurately predicting price
4. since the indicator is now included as a feature in the matching algorithm, testing the indicator is, very indirectly, a test of the matching algorithm too

The test I have in mind is a hybrid of a hypothesis test, Cross validation and the application of Statistical process control.  Consider the chart below:

This shows two sampling distributions of the mean for Long MFE/MAE indicator values > 0.5, the upper pane for sample sizes of 20 and the lower pane for 75. For simplicity I shall only discuss the Long > 0.5 version of the indicator, but everything that follows applies equally to the Short version. As expected the upper pane shows greater variance, and for the envisioned test a whole series of these sampling distributions will be produced for different sampling rates. The way I intend it to work is as follows:

• take a single bar in the history and see what the value of the MFE/MAE indicator value is 3 bars later (assume > 0.5 for this exposition, so we compare to long sampling distributions only)
• get the top 20 matched bars for the above selected bar and the corresponding 20 indicator values for 3 bars later and take the mean of these 20 indicator values
• check if this mean falls within the sampling distribution of the mean of 20, as shown in the upper pane above by the vertical black line at 0.8 on the x axis. If it does fall with the sampling distribution, we accept the null hypothesis that the 20 best matches in history future indicator values and the value of the indicator after the bar to be matched come from the same distribution
• repeat the immediately preceding step for means of 21, 22, ... etc until such time as the null hypothesis can be rejected, shown in the lower pane above. At this point, we then then declare an upper bound on the historical number of matches for the bar to be predicted

For any single bar to be predicted we can then produce the following chart, which is completely artificial and just for illustrative purposes:

where the cyan and red lines are the +/- 2 standard deviations above/below a notional mean value for the whole distribution of approximately 0.85, and the chart can be considered to be a type of control chart. The upper and lower control lines converge towards the right, reflecting the decreasing variance of increasingly large N sample means, as shown in the first chart above. The green line represents the cumulative N sample mean of the best N historical matches' future values. I have shown it as decreasing as it is to be expected that as more N matches are included, the greater the chance that incorrect matches, unexpected price reversals etc. will be caught up in this mean calculation, resulting in the mean value moving into the left tail of the sampling distribution. This effect combines with the shrinking variance to reach a critical point (rejection of the null hypothesis) at which the green line exits below the lower control line.

The purpose of all the above is provide a principled manner to choose the number N matches from the Cauchy-Schwarz matching algorithm to supply instances of training data to the envisioned neural net training. An incidental benefit of this approach is that it is indirectly a hypothesis test of the fundamental assumption underlying the matching algorithm; namely that past price action has predictive ability for future price action, and furthermore, it is a test of the MFE/MAE indicator. Discussion of the results of these tests in a future post.       

First Use for the MFE/MAE Indicator

This first use is as an input to my Cauchy-Schwarz matching algorithm, previous posts about which can be read here, here and here. The screen shot below shows what I would characterise as a "good" set of matches: 

The top left pane shows the original section of the price series to be matched, and the panes labelled #1, #5, etc. are the best match, 5th best match and so on respectively. The last 3 rightmost bars in each pane are "future" price bars, i.e. the 4th bar in from the right is the target bar that is being matched, matched over all the bars to the left or in the past of this target bar.

I consider the above to be a set of "good" matches because, for the #1 through #25 matches for "future" bars:

• if one considers the logic of the mfe/mae indicator each pane gives indicator readings of "long," which all agree with the original "future" bars
• similarly the mae (maximum adverse excursion) occurs on the day immediately following the matched day
• the mfe (maximum favourable excursion) occurs on the 3rd "future" bar, with the slight exception of pane #10
• the marked to market returns of an entry at the open of the 1st "future" bar to the close of the 3rd "future" bar all show a profit, as does the original pane

However, it can be seen that the above noted "goodness" breaks down for panes #25 and #30, which leads me to postulate that there is an upper bound on the number of matches for which there is predictive ability for "future" returns.

In the above linked posts the test statistic used to judge the predictive efficacy of the matching algorithm was effect size. However, I think a more pertinent test statistic to use would be the average bar return over the bars immediately following a matched bar, and a discussion of this will be the subject of my next post.  

A New MFE/MAE Indicator.

After stopping my investigation of tools for spectral analysis over the last few weeks I have been doing another mooc, this time Learning from Data, and also working on the idea of one of my earlier posts.

In the above linked post there is a video showing the idea as a "paint bar" study. However, I thought it would be a good idea to render it as an indicator, the C++ Octave .oct code for which is shown in the code box below.

DEFUN_DLD ( adjustable_mfe_mae_from_open_indicator, args, nargout,
"-*- texinfo -*-\n\
@deftypefn {Function File} {} adjustable_mfe_mae_from_open_indicator (@var{open,high,low,close,lookback_length})\n\
This function takes four input series for the OHLC and a value for lookback length. The main outputs are\n\
two indicators, long and short, that show the ratio of the MFE over the MAE from the open of the specified\n\
lookback in the past. The indicators are normalised to the range 0 to 1 by a sigmoid function and a MFE/MAE\n\
ratio of 1:1 is shifted in the sigmoid function to give a 'neutral' indicator reading of 0.5. A third output\n\
is the max high - min low range over the lookback_length normalised by the range of the daily support and\n\
resistance levels S1 and R1 calculated for the first bar of the lookback period. This is also normalised to\n\
give a reading of 0.5 in the sigmoid function if the ratio is 1:1. The point of this third output is to give\n\
some relative scale to the unitless MFE/MAE ratio and to act as a measure of strength or importance of the\n\
MFE/MAE ratio.\n\
@end deftypefn" )

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

// check the input arguments
if ( nargin != 5 )
   {
   error ( "Invalid arguments. Arguments are price series for open, high, low and close and value for lookback length." ) ;
   return retval_list ;
   }

if ( args(4).length () != 1 )
   {
   error ( "Invalid argument. Argument 5 is a scalar value for the lookback length." ) ;
   return retval_list ;
   }
   
   int lookback_length = args(4).int_value() ;

if ( args(0).length () < lookback_length )
   {
   error ( "Invalid argument lengths. Argument lengths for open, high, low and close vectors should be >= lookback length." ) ;
   return retval_list ;
   }
   
if ( args(1).length () != args(0).length () )
   {
   error ( "Invalid argument lengths. Argument lengths for open, high, low and close vectors should be equal." ) ;
   return retval_list ;
   }
   
if ( args(2).length () != args(0).length () )
   {
   error ( "Invalid argument lengths. Argument lengths for open, high, low and close vectors should be equal." ) ;
   return retval_list ;
   }
   
if ( args(3).length () != args(0).length () )
   {
   error ( "Invalid argument lengths. Argument lengths for open, high, low and close vectors should be equal." ) ;
   return retval_list ;
   }   

if (error_state)
   {
   error ( "Invalid arguments. Arguments are price series for open, high, low and close and value for lookback length." ) ;
   return retval_list ;
   }
// end of input checking
  
// inputs
ColumnVector open = args(0).column_vector_value () ;
ColumnVector high = args(1).column_vector_value () ;
ColumnVector low = args(2).column_vector_value () ;
ColumnVector close = args(3).column_vector_value () ;
// outputs
ColumnVector long_mfe_mae = args(0).column_vector_value () ;
ColumnVector short_mfe_mae = args(0).column_vector_value () ;
ColumnVector range = args(0).column_vector_value () ;

// variables
double max_high = *std::max_element( &high(0), &high( lookback_length ) ) ;
double min_low = *std::min_element( &low(0), &low( lookback_length ) ) ;
double pivot_point = ( high(0) + low(0) + close(0) ) / 3.0 ;
double s1 = 2.0 * pivot_point - high(0) ;
double r1 = 2.0 * pivot_point - low(0) ;

for ( octave_idx_type ii (0) ; ii < lookback_length ; ii++ ) // initial ii loop
    {

      // long_mfe_mae
      if ( open(0) > min_low ) // the "normal" situation
      {
 long_mfe_mae(ii) = 1.0 / ( 1.0 + exp( -( ( max_high - open(0) ) / ( open(0) - min_low ) - 1.0 ) ) ) ;
      }
      else if ( open(0) == min_low )
      {
 long_mfe_mae(ii) = 1.0 ;
      }
      else
      {
 long_mfe_mae(ii) = 0.5 ;
      }
           
      // short_mfe_mae
      if ( open(0) < max_high ) // the "normal" situation
      {
 short_mfe_mae(ii) = 1.0 / ( 1.0 + exp( -( ( open(0) - min_low ) / ( max_high - open(0) ) - 1.0 ) ) ) ;
      }
      else if ( open(0) == max_high )
      {
 short_mfe_mae(ii) = 1.0 ;
      }
      else
      {
 short_mfe_mae(ii) = 0.5 ;
      }
      
      range(ii) = 1.0 / ( 1.0 + exp( -( ( max_high - min_low ) / ( r1 - s1 ) - 1.0 ) ) ) ;
      
    } // end of initial ii loop

for ( octave_idx_type ii ( lookback_length ) ; ii < args(0).length() ; ii++ ) // main ii loop
    { 
    // assign variable values  
    max_high = *std::max_element( &high( ii - lookback_length + 1 ), &high( ii + 1 ) ) ;
    min_low = *std::min_element( &low( ii - lookback_length + 1 ), &low( ii + 1 ) ) ;
    pivot_point = ( high(ii-lookback_length) + low(ii-lookback_length) + close(ii-lookback_length) ) / 3.0 ;
    s1 = 2.0 * pivot_point - high(ii-lookback_length) ;
    r1 = 2.0 * pivot_point - low(ii-lookback_length) ;

      // long_mfe_mae
      if ( open( ii - lookback_length + 1 ) > min_low && open( ii - lookback_length + 1 ) < max_high ) // the "normal" situation
      {
 long_mfe_mae(ii) = 1.0 / ( 1.0 + exp( -( ( max_high - open( ii - lookback_length + 1 ) ) / ( open( ii - lookback_length + 1 ) - min_low ) - 1.0 ) ) ) ;
      }
      else if ( open( ii - lookback_length + 1 ) == min_low )
      {
 long_mfe_mae(ii) = 1.0 ;
      }
      else
      {
 long_mfe_mae(ii) = 0.0 ;
      }
   
      // short_mfe_mae
      if ( open( ii - lookback_length + 1 ) > min_low && open( ii - lookback_length + 1 ) < max_high ) // the "normal" situation
      {
 short_mfe_mae(ii) = 1.0 / ( 1.0 + exp( -( ( open( ii - lookback_length + 1 ) - min_low ) / ( max_high - open( ii - lookback_length + 1 ) ) - 1.0 ) ) ) ;
      }
      else if ( open( ii - lookback_length + 1 ) == max_high )
      {
 short_mfe_mae(ii) = 1.0 ;
      }
      else
      {
 short_mfe_mae(ii) = 0.0 ;
      }
      
      range(ii) = 1.0 / ( 1.0 + exp( -( ( max_high - min_low ) / ( r1 - s1 ) - 1.0 ) ) ) ;

    } // end of main ii loop

retval_list(2) = range ;    
retval_list(1) = short_mfe_mae ;
retval_list(0) = long_mfe_mae ;

return retval_list ;
  
} // end of function

The way to interpret this is as follows:

• if the "long" indicator reading is above 0.5, go long
• if the "short" is above 0.5, go short
• if both are below 0.5, go flat

An alternative, if the indicator reading is flat, is to maintain any previous non flat position. I won't show a chart of the indicator itself as it just looks like a very noisy oscillator, but the equity curve(s) of it, without the benefit of foresight, on the EURUSD forex pair are shown below.

The yellow equity curve is the cumulative, close to close, tick returns of a buy and hold strategy, the blue is the return going flat when indicated, and the red maintaining the previous position when flat is indicated. Not much to write home about. However, this second chart shows the return when one has the benefit of the "peek into the future" as discussed in my earlier post.

The colour of the curves are as before except for the addition of the green equity curve, which is the cumulative, vwap value to vwap value tick returns, a simple representation of what an equity curve with realistic slippage might look like. This second set of equity curves shows the promise of what could be achievable if a neural net to accurately predict future values of the above indicator can be trained. More in an upcoming post.