Wednesday, 4 September 2013

FFT Smoother Function Test

Following on from my last post I now want to show the FFT smoother function results on my usual idealised market data. The Octave script code for this test is given in the code box below.
clear all

% create the "price" vector
period = input( 'Enter period between 10 & 50 for underlying sine: ' ) ;
underlying_sine = sinewave( 500 , period ) ; % create underlying sinewave function with amplitude 2
trend_mult = input( 'Enter value for trend (a number between -4 & 4 incl.) to 4 decimal places: ' ) ;
trend = ( trend_mult / period ) .* (0:1:499) ;
noise_mult = input( 'Enter value for noise multiplier (a number between 0 & 1 incl.) to 2 decimal places: ' ) ;
noise = noise_mult * randn(500,1)' ;
smooth_price = underlying_sine .+ trend ;
price = underlying_sine .+ trend .+ noise ;
super_smooth = super_smoother( price ) ;

smooth = price ;
mean_value = price ;
mean_smooth = price ;

for ii = period+1 : length( price )

mean_value(ii) = mean( price( ii-period : ii ) ) ;
fft_input = price( ii-(period-1) : ii ) .- mean_value(ii) ;
N = length( fft_input ) ;

data_in_freq_domain = fft( fft_input ) ;

if mod( N , 2 ) == 0 % N even
% The DC component is unique and should not be altered
% adjust value of element N / 2 to compensate for it being the sum of frequencies
data_in_freq_domain( N / 2 ) = abs( data_in_freq_domain( N / 2 ) ) / 2 ;

[ max1 , ix1 ] = max( abs( data_in_freq_domain( 2 : N / 2 ) ) ) ; % get the index of the first max value
ix1 = ix1 + 1 ;

[ max2 , ix2 ] = max( abs( data_in_freq_domain( N / 2 + 1 : end ) ) ) ; % get index of the second max value
ix2 = N / 2 + ix2 ;

non_zero_values_ix = [ 1 ix1 ix2 ] ; % indices for DC component, +ve freq max & -ve freq max

else % N odd
% The DC component is unique and should not be altered

[ max1 , ix1 ] = max( abs( data_in_freq_domain( 2 : ( N + 1 ) / 2 ) ) ) ; % get the index of the first max value
ix1 = ix1 + 1 ;

[ max2 , ix2 ] = max( abs( data_in_freq_domain( ( N + 1 ) / 2 + 1 : end ) ) ) ; % get index of the second max value
ix2 = ( N + 1 ) / 2 + ix2 ;

non_zero_values_ix = [ 1 ix1 ix2 ] ; % indices for DC component, +ve freq max & -ve freq max

end

% extract the constant term and frequencies of interest
non_zero_values = data_in_freq_domain( non_zero_values_ix ) ;

% now set all values to zero
data_in_freq_domain( 1 : end ) = 0.0 ;

% replace the non_zero_values
data_in_freq_domain( non_zero_values_ix ) = non_zero_values ;

inv = ifft( data_in_freq_domain ) ; % reverse the FFT

recovered_sig = real( inv ) ; % recover the real component of ifft for plotting

smooth( ii ) = recovered_sig( end ) + mean_value( ii ) ;

mean_smooth( ii ) = mean( smooth( ii-period : ii ) ) ;

end

smooth = smooth .+ ( mean_value .- mean_smooth ) ;

plot( smooth_price , 'r' ,price , 'c' , smooth , 'm' , super_smooth , 'g' )
legend( 'smooth price' , 'price' , 'fft smooth' , 'super smooth' )
This script prompts for terminal input for the values for an underlying sine wave period, a trend multiplier and a noise multiplier component to construct the "price" time series. Here is a plot of output with the settings
• period = 20
• trend multiplier = 1.3333
• noise multiplier = 0.0
The cyan line is the "price," the magenta is the FFT smoother and the green is the "Super Smoother" filter shown for comparative purposes. As can be seen, after the "burn in" period for the FFT smoother calculations to settle down, the smoother is a perfect recovery of the price time series. However, this is without noise. With a noise multiplier value of 0.25 a typical plot is
where the underlying smooth "true" signal can now be seen in red. I think the performance of the magenta FFT smoother in recovering this "true" signal speaks for itself. My next post will show its performance on real price data.