A time series {x_k} can be represented as a discrete real-time domain function f, x = f(t). In the 18^th century, Jean Baptiste Joseph Fourier demonstrated that any continuous periodic function f can be represented as a linear combination of sine and cosine functions. The discrete Fourier transform (DFT) is a linear transformation that converts a time series into a list of coefficients of a finite combination of complex or real trigonometric functions, ordered by their frequencies.

The frequency ω of each trigonometric function defines one of the harmonics of the signal. The space that represents the signal amplitude versus frequency of the signal is known as the frequency domain. The generic DFT transforms a time series into a sequence of frequencies defined as complex numbers a + j.φ (j² = -1), where a is the amplitude of the frequency and φ is the phase.

This section is dedicated to the real DFT that converts a time series into an ordered sequence of frequencies with real values.

Note

Real discrete Fourier transform

M5: A periodic function f can be represented as an infinite combination of sine and cosine functions:

M6: The Fourier cosine transform of a function f is defined as:

M7: The discrete real cosine series of a function f(-x) = f(x) is defined as:

M8: The Fourier sine transform of a function is defined as:

M9: The discrete real sine series of a function f(-x) = f(x) is defined as:

The computation of the Fourier trigonometric series is time consuming with an asymptotic time complexity of O(n²). Scientists and mathematicians have been working to make the computation as effective as possible. The most common numerical algorithm used to compute the Fourier series is the Fast Fourier Transform (FFT) created by J.W. Cooley and J. Tukey [3:7].

The algorithm called Radix-2 version recursively breaks down the Fourier transform for a time series of N data points into any combination of N₁ and N₂ sized segments such as N = N₁ N₂. Ultimately, the discrete Fourier transform is applied to the deeper-nested segments.

The Radix-2 implementation requires that the number of data points is N=2ⁿ for even functions (sine) and N=2ⁿ+1 for cosine. There are two approaches to meet this constraint:

Padding the time series is the preferred option because it does not affect the original set of observations.

Let's define a DTransform trait for any variant of the discrete Fourier transform. The first step is to wrap the default configuration parameters used in the Apache Commons Math library into a Config singleton:

trait DTransform {
  object Config {
     final val FORWARD = TransformType.FORWARD
     final val INVERSE = TransformType.INVERSE
     final val SINE = DstNormalization.STANDARD_DST_I
     final val COSINE = DctNormalization.STANDARD_DCT_I
   }
   …
}

The main purpose of the DTransform trait is to pad the vec time series with zero values:

def pad(vec: DblVector, 
    even: Boolean = true)(implicit f: T =>Double): DblArray = {
  val newSize = padSize(vec.size, even)   //1
  val arr: DblVector = vec.map(_.toDouble)
  if( newSize > 0) arr ++ Array.fill(newSize)(0.0) else arr //2
}

def padSize(xtSz: Int, even: Boolean= true): Int = {
  val sz = if( even ) xtSz else xtSz-1  //3
  if( (sz & (sz-1)) == 0) 0
  else {
    var bitPos = 0
    do { bitPos += 1 } while( (sz >> bitPos) > 0) //4
    (if(even) (1<<bitPos) else (1<<bitPos)+1) - xtSz
  }
}

The pad method computes the optimal size of the frequency vector as 2^N by invoking the padSize method (line 1). It then concatenates the padding with the original time series or vector of observations (line 2). The padSize method adjusts the size of the data depending on whether the time series has initially an even or odd number of observations (line 3). It relies on bit operations to find the next radix, N (line 4).

The fast implementation of the padding method, pad, consists of detecting the number of N observations as a power of 2 (the next highest radix). The method evaluates if N and (N-1) are zero after it shifts the number of bits in the value, N. The code illustrates the effective use of implicit conversion to make the code readable in the pad method:

 val arr: DblVector = vec.map(_.toDouble)

The next step is to write the DFT class for the real sine and cosine discrete transforms by subclassing DTransform. The class relies on the padding mechanism implemented in DTransform whenever necessary:

class DFT[@specialized(Double) T <: AnyVal](
    eps: Double)(implicit f: T => Double)
    extends ETransform[Double](eps) with DTransform { //5
  type U = XSeries[T]   //6
  type V = DblVector
  
  override def |> : PartialFunction[U, Try[V]] = { //7
    case xv: U if(xv.size >= 2) => fwrd(xv).map(_._2.toVector) 
  }
}

We treat the discrete Fourier transform as a transformation on the time series using an explicit ETransform configuration (line 5). The U data type of the input and the V type of the output have to be defined (line 6). The |> transformation function delegates the computation to the fwrd method (line 7):

def fwrd(xv: U): Try[(RealTransformer, DblArray)] = {
  val rdt = if(Math.abs(xv.head) < config)  //8
      new FastSineTransformer(SINE)  //9
  else  new FastCosineTransformer(COSINE)  //10
  
  val padded = pad(xv.map(_.toDouble), xv.head == 0.0).toArray
  Try( (rdt, rdt.transform(padded, FORWARD)) )
}

The fwrd method selects the discrete Fourier sine series if the first value of the time series is 0.0, otherwise it selects the discrete cosine series. This implementation automates the selection of the appropriate series by evaluating xt.head (line 8). The transformation invokes the FastSineTransformer (line 9) and FastCosineTransformer (line 10) classes of the Apache Commons Math library [3:8] introduced in the first chapter.

This example uses the standard formulation of the cosine and sine transformations, defined by the COSINE argument. The orthogonal normalization that normalizes the frequency by a factor of 1/sqrt(2(N-1)), where N is the size of the time series, generates a cleaner frequency spectrum for a higher computation cost.

In order to illustrate the different concepts behind DFTs, let's consider the case of a time series generated by a h sequence of sinusoidal functions:

val F = Array[Double](2.0, 5.0, 15.0)
val A = Array[Double](2.0, 1.0, 0.33)

def harmonic(x: Double, n: Int): Double =  
      A(n)*Math.cos(Math.PI*F(n)*x)
val h = (x: Double) => 
    Range(0, A.size).aggregate(0.0)((s, i) => 
          s + harmonic(x, i), _ + _)

As the signal is synthetically created, we can select the size of the time series to avoid padding. The first value in the time series is not null, so the number of observations is 2ⁿ+1. The data generated by the h function is plotted as follows:

An example of sinusoidal time series

Let's extract the frequencies' spectrum for the time series generated by the h function. The data points are created by tabulating the h function. The frequencies spectrum is computed with a simple invocation of the explicit |> data transformation of the DFT class:

val OUTPUT1 = "output/chap3/simulated.csv"
val OUTPUT2 = "output/chap3/smoothed.csv"
val FREQ_SIZE = 1025; val INV_FREQ = 1.0/FREQ_SIZE

val pfnDFT = DFT[Double] |> //11
for {
  values <- Try(Vector.tabulate(FREQ_SIZE)
               (n => h(n*INV_FREQ))) //12
  output1 <- DataSink[Double](OUTPUT1).write(values)
  spectrum <- pfnDFT(values)
  output2 <- DataSink[Double](OUTPUT2).write(spectrum) //13
} yield {
  val results = format(spectrum.take(DISPLAY_SIZE), 
"x/1025", SHORT)
  show(s"$DISPLAY_SIZE frequencies: ${results}")
}

The execution of the data simulator follows these steps:

The spectrum of frequencies of the time series, plotted for the first 32 points, clearly shows three frequencies at k = 2, 5, and 15. This result is expected because the original signal is composed of three sinusoidal functions. The amplitude of these frequencies are 1024/1, 1024/2, and 1024/6, respectively. The following plot represents the first 32 harmonics for the time series:

The frequency spectrum for a three-frequency sinusoidal

The next step is to use the frequencies spectrum to create a low-pass filter using DFT. There are many algorithms available to implement a low or pass band filter in the time domain, from autoregressive models to the Butterworth algorithm. However, the discrete Fourier transform is still a very popular technique to smooth signals and identify trends.

The purpose of this section is to introduce, describe, and implement a noise filtering mechanism that leverages the discrete Fourier transform. The idea is quite simple: the forward and inverse Fourier series are used sequentially to convert the raw data from the time domain to the frequency domain and back. The only input you need to supply is a function g that modifies the sequence of frequencies. This operation is known as the convolution of the filter g and the frequencies' spectrum.

A convolution is similar to an inner product of two time series in the frequencies domain. Mathematically, the convolution is defined as follows:

Note

Convolution

M10: The convolution of two functions f and g is defined as:

M11: The convolution F of a time series x = (x_i} with a frequency spectrum ω^x and a filter f in frequency domain ω^f is defined as:

Let's apply the convolution to our filtering problem. The filtering algorithm using the discrete Fourier transform consists of five steps:

The most commonly used low-pass filter functions are known as the sinc and sinc2 functions, which are defined as a rectangular function and triangular function, respectively. These functions are partially applied functions that are derived from a generic convol method. The simplest sinc function returns 1 for frequencies below a cutoff frequency, fC, and 0 if the frequency is higher:

val convol = (n: Int, f: Double, fC: Double) => 
     if( Math.pow(f, n) < fC) 1.0 else 0.0
val sinc = convol(1, _: Double, _:Double)
val sinc2 = convol(2, _: Double, _:Double)
val sinc4 = convol(4, _: Double, _:Double)

Note

Partially applied functions versus partial functions

Partial functions and partially applied functions are not actually related.

A partial function f' is a function that is applied to a subset X' of the input space X. It does not execute all possible input values:

A partially applied function f" is a function value for which the user supplies the value for one or more arguments. The projection reduces the dimension of the input space (X, Z):

The DFTFilter class inherits from the DFT class in order to reuse the fwrd forward transform function. The g frequency domain function is an attribute of the filter. The g function takes the fC frequency cutoff value as the second argument (line 14). The two sinc and sinc2 filters defined in the previous section are examples of filtering functions:

class DFTFilter[@specialized(Double) T <: AnyVal]( 
    fC: Double,
    eps: Double)
    (g: (Double, Double) =>Double)(implicit f: T => Double)
  extends DFT[T](eps) { //14

  override def |> : PartialFunction[U, Try[V]] = {
    case xt: U if( xt.size >= 2 ) => {
      fwrd(xt).map{ case(trf, freq) => {  //15
        val cutOff = fC*freq.size
        val filtered = freq.zipWithIndex
                     .map{ case(x, n) => x*g(n, cutOff) } //16
        trf.transform(filtered, INVERSE).toVector }) //17
    }
  }
}

The filtering process follows three steps:

Let's evaluate the impact of the cutoff values on the filtered data. The implementation of the test program consists of loading the data from the file (line 19) and then invoking the DFTFilter of the pfnDFTfilter partial function (line 19):

import YahooFinancials._

val inputFile = s"$RESOURCE_PATH$symbol.csv"
val src = DataSource(input, false, true, 1)
val CUTOFF = 0.005
val pfnDFTfilter = DFTFilter[Double](CUTOFF)(sinc) |>
for {
  price <- src.get(adjClose)  //18
  filtered <- pfnDFTfilter(price)  //19
} 
yield { /* ... */ }

Filtering out the noise is accomplished by selecting the cutoff value between any of the three harmonics with the respective frequencies of 2, 5, and 15. The original and the two filtered time series are plotted on the following graph:

Plotting of the discrete Fourier filter-based smoothing

As you would expect, the low-pass filter with a cutoff value of 12 eliminates the noise with the highest frequencies. The filter with the cutoff value 4 cancels out the second harmonic (low-frequency noise), leaving out only the main trend cycle.

Using the discrete Fourier transform to generate the frequencies spectrum of a periodical time series is easy. However, what about real-world signals such as the time series that represent the historical price of a stock?

The purpose of the next exercise is to detect, if any, the long term cycle(s) of the overall stock market by applying the discrete Fourier transform to the quote of the S&P 500 index between January 1, 2009 and December 31, 2013, as illustrated in the following graph:

Historical S&P 500 index prices

The first step is to apply the DFT to extract a frequencies spectrum for the S&P 500 historical prices, as shown in the following graph, with the first 32 harmonics:

Frequencies spectrum for the historical S&P index

The frequency domain chart highlights some interesting characteristics regarding the S&P 500 historical prices:

A low-pass filter is limited to reduce or cancel out the noise in the raw data. In this case, a band-pass filter using a range or window of frequencies is appropriate to isolate the frequency or the group of frequencies that characterize a specific cycle. The sinc function, which was introduced in the previous section to implement a low-pass filter, is modified to enforce the band-pass filter within a window, [w₁, w₂], as follows:

def sinc(f: Double, w: (Double, Double)): Double = 
    if(f > w._1 && f < w._2) 1.0 else 0.0

Let's define a DFT-based band-pass filter with a window of width 4, w=(i, i+4), with i ranging between 2 and 20. Applying the window [4, 8] isolates the impact of the second harmonic on the price. As we eliminate the main upward trend with frequencies less than 4, all the filtered data varies within a short range relative to the main trend. The following graph shows the output of this filter:

The output of a band-pass DFT filter range 4-8 on the historical S&P index

In this case, we filter the S&P 500 index around the third group of harmonics with frequencies ranging from 18 to 22; the signal is converted into a familiar sinusoidal function, as shown here:

The output of a band-pass DFT filter range 18-22 on the historical S&P index

There is a possible rational explanation for the shape of the S&P 500 data filtered by a band-pass filter with a frequency of 20, as illustrated in the previous graph. The S&P 500 historical data plot shows that the frequency of the fluctuation in the middle of the uptrend (trading sessions 620 to 770) increases significantly.

This phenomenon can be explained by the fact that the S&P 500 index reaches a resistance level around the trading session 545 when the existing uptrend breaks. A tug of war starts between the bulls, betting the market nudges higher, and the bears, who are expecting a correction. The back and forth between the traders ends when the S&P 500 index breaks through its resistance and resumes a strong uptrend characterized by a high amplitude low frequency, as shown in the following graph:

An illustration of support and resistance levels for the historical S&P 500 index prices

One of the limitations of using the discrete Fourier-based filters to clean up data is that it requires the data scientist to extract the frequencies spectrum and modify the filter on a regular basis, as he or she is never sure that the most recent batch of data does not introduce noise with a different frequency. The Kalman filter addresses this limitation.