Finite Difference Approximation of Derivatives

A while ago, someone asked me to reference him in a paper of mine because I used formulas of a finite difference approximation of a derivative on a non uniform grid. I was shocked as those formula are very widespread (in countless papers, courses and books) and not far off elementary mathematics.

There are however some interesting old papers on the technique. Usually people approximate the first derivative by the central approximation of second order:

$$ f'(x) = \frac{f(x_{i+1})-f(x_{i-1})}{x_{i+1} - x_{i-1}} $$
However there are some other possibilities. For example one can find a formula directly out of the Taylor expansions of f(x(i+1)) and f(x(i-1)). This paper and that one seems to indicate it is more precise, especially when the grid does not vary smoothly (a typical example is uniform by parts).

This can make a big difference in practice, here is the example of a Bond priced under the Cox-Ingersoll-Ross model by finite differences. EULER is the classic central approximation, EULER1 uses the more refined approximation based on Taylor expansion, EULER2 uses Taylor expansion approximation as well as a higher order boundary condition. I used the same parameters as in the Tavella-Randall book example and a uniform grid between [0,0.2] except that I have added 2 points at the far end at 0.5 and 1.0. So the only difference between EULER and EULER1 lies in the computation of derivatives at the 3 last points.

I also computed the backward 2nd order first derivative on a non uniform grid (for the refined boundary). I was surprised not to find this easily on the web, so here it is:

$$ f'(x_i) = \left(\frac{1}{h_i}+\frac{1}{h_i+h_{i-1}}\right) f(x_i)- \left(\frac{1}{h_{i-1}}+\frac{1}{h_i}\right) f(x_{i-1})+ \left(\frac{1}{h_{i-1}} - \frac{1}{h_i+h_{i-1}} \right) f(x_{i-2}) + ...$$

Incidently while writing this post I found out it was a pain to write Math in HTML (I initially used a picture). MathML seems a bit crazy, I wonder why they couldn't just use the LaTeX standard. Update January 3rd 2013 - I now use Mathjax. It's not very good solution as I think this should typically be handled by the browser directly instead of huge javascript library, but it looks a bit better

Finite Difference Approximation of Derivatives

A while ago, someone asked me to reference him in a paper of mine because I used formulas of a finite difference approximation of a derivative on a non uniform grid. I was shocked as those formula are very widespread (in countless papers, courses and books) and not far off elementary mathematics.

There are however some interesting old papers on the technique. Usually people approximate the first derivative by the central approximation of second order:

$$ f'(x) = \frac{f(x_{i+1})-f(x_{i-1})}{x_{i+1} - x_{i-1}} $$
However there are some other possibilities. For example one can find a formula directly out of the Taylor expansions of f(x(i+1)) and f(x(i-1)). This paper and that one seems to indicate it is more precise, especially when the grid does not vary smoothly (a typical example is uniform by parts).

This can make a big difference in practice, here is the example of a Bond priced under the Cox-Ingersoll-Ross model by finite differences. EULER is the classic central approximation, EULER1 uses the more refined approximation based on Taylor expansion, EULER2 uses Taylor expansion approximation as well as a higher order boundary condition. I used the same parameters as in the Tavella-Randall book example and a uniform grid between [0,0.2] except that I have added 2 points at the far end at 0.5 and 1.0. So the only difference between EULER and EULER1 lies in the computation of derivatives at the 3 last points.

I also computed the backward 2nd order first derivative on a non uniform grid (for the refined boundary). I was surprised not to find this easily on the web, so here it is:

$$ f'(x_i) = \left(\frac{1}{h_i}+\frac{1}{h_i+h_{i-1}}\right) f(x_i)- \left(\frac{1}{h_{i-1}}+\frac{1}{h_i}\right) f(x_{i-1})+ \left(\frac{1}{h_{i-1}} - \frac{1}{h_i+h_{i-1}} \right) f(x_{i-2}) + ...$$

Incidently while writing this post I found out it was a pain to write Math in HTML (I initially used a picture). MathML seems a bit crazy, I wonder why they couldn't just use the LaTeX standard. Update January 3rd 2013 - I now use Mathjax. It's not very good solution as I think this should typically be handled by the browser directly instead of huge javascript library, but it looks a bit better

A Discontinuity

I am comparing various finite difference schemes on simple problems and am currently stumbling upon a strange discontinuity at the boundary for some of the schemes (Crank-Nicolson, Rannacher, and TR-BDF2) when I plot an American Put Option Gamma using a log grid. It actually is more pronounced with some values of the strike, not all. The amplitude oscillates with the strike. And it does not happen on a European Put, so it's not a boundary approximation error in the code. It might well be due to the nature of the scheme as schemes based on implicit Euler work (maybe monotonicity preservation is important). This appears on this graph around S=350.

Update December 13, 2012: after a close look at what was happening. It was after all a boundary issue. It's more visible on the American because the Gamma is more spread out. But I reproduced it on a European as well.

A Discontinuity

I am comparing various finite difference schemes on simple problems and am currently stumbling upon a strange discontinuity at the boundary for some of the schemes (Crank-Nicolson, Rannacher, and TR-BDF2) when I plot an American Put Option Gamma using a log grid. It actually is more pronounced with some values of the strike, not all. The amplitude oscillates with the strike. And it does not happen on a European Put, so it's not a boundary approximation error in the code. It might well be due to the nature of the scheme as schemes based on implicit Euler work (maybe monotonicity preservation is important). This appears on this graph around S=350.

Update December 13, 2012: after a close look at what was happening. It was after all a boundary issue. It's more visible on the American because the Gamma is more spread out. But I reproduced it on a European as well.

Scala is Mad

I spent quick a bit of time to figure out why something that is usually simple to do in Java did not work in Scala: Arrays and ArrayLists with generics.

For some technical reason (type erasure at the JVM level), Array sometimes need a parameter with a ClassManifest !?! a generic type like [T :< Point : ClassManifest] need to be declared instead of simply [T :< Point].

And then the quickSort method somehow does not work if invoked on a generic... like quickSort(points) where points: Array[T]. I could not figure out yet how to do this one, I just casted to points.asInstanceOf[Array[Point]], quite ugly.

In contrast I did not even have to think much to write the Java equivalent. Generics in Scala, while having a nice syntax, are just crazy. This is something that goes beyond generics. Some of the Scala library and syntax is nice, but overall, the IDE integration is still very buggy, and productivity is not higher.

Update Dec 12 2012: here is the actual code (this is kept close to the Java equivalent on purpose):

object Point {
  def sortAndRemoveIdenticalPoints[T <: Point : ClassManifest](points : Array[T]) : Array[T] = {
      Sorting.quickSort(points.asInstanceOf[Array[Point]])
      val l = new ArrayBuffer[T](points.length)
      var previous = points(0)
      l += points(0)
      for (i <- 1 until points.length) {
        if(math.abs(points(i).value - previous.value)< Epsilon.MACHINE_EPSILON_SQRT) {
          l += points(i)
        }
      }
      return l.toArray
    }
    return points
  }
}

class Point(val value: Double, val isMiddle: Boolean) extends Ordered[Point] {
  def compare(that: Point): Int = {
    return math.signum(this.value - that.value).toInt
  }
}

In Java one can just use Arrays.sort(points) if points is a T[]. And the method can work with a subclass of Point.

Scala is Mad

I spent quick a bit of time to figure out why something that is usually simple to do in Java did not work in Scala: Arrays and ArrayLists with generics.

For some technical reason (type erasure at the JVM level), Array sometimes need a parameter with a ClassManifest !?! a generic type like [T :< Point : ClassManifest] need to be declared instead of simply [T :< Point].

And then the quickSort method somehow does not work if invoked on a generic... like quickSort(points) where points: Array[T]. I could not figure out yet how to do this one, I just casted to points.asInstanceOf[Array[Point]], quite ugly.

In contrast I did not even have to think much to write the Java equivalent. Generics in Scala, while having a nice syntax, are just crazy. This is something that goes beyond generics. Some of the Scala library and syntax is nice, but overall, the IDE integration is still very buggy, and productivity is not higher.

Update Dec 12 2012: here is the actual code (this is kept close to the Java equivalent on purpose):

object Point {
  def sortAndRemoveIdenticalPoints[T <: Point : ClassManifest](points : Array[T]) : Array[T] = {
      Sorting.quickSort(points.asInstanceOf[Array[Point]])
      val l = new ArrayBuffer[T](points.length)
      var previous = points(0)
      l += points(0)
      for (i <- 1 until points.length) {
        if(math.abs(points(i).value - previous.value)< Epsilon.MACHINE_EPSILON_SQRT) {
          l += points(i)
        }
      }
      return l.toArray
    }
    return points
  }
}

class Point(val value: Double, val isMiddle: Boolean) extends Ordered[Point] {
  def compare(that: Point): Int = {
    return math.signum(this.value - that.value).toInt
  }
}

In Java one can just use Arrays.sort(points) if points is a T[]. And the method can work with a subclass of Point.