Efficient Methods for Pricing American Options Across Volatility Surfaces

# Most efficient way to get American option prices from an implied volatility surface


# Most efficient way to get American option prices from an implied volatility surface












Using C++, I am generating end-of-day implied volatility surfaces with maturity and moneyness (strike/forward price) as the other 2 axes, via a 1D kernel (normalized gaussian) method and then smoothed with the same 1D kernel method, for american equity options for the Russell 3000 (3000 tickers) which is already a resource intensive process that takes a good amount of time. I then simply reprice the points in the surface but it takes a long time. We use finite differences method (BSM model) with discrete dividends to price and compute IV. Is there a more efficient way to get the corresponding option prices of an implied volatility surface instead of simply repricing each point in the surface?

## Answer by AKdemy (score 3)

https://quant.stackexchange.com/a/81628

Generally, it really shouldn't take too long. What language do you use. Is this all coded by yourself? There is a lot of difference between languages and the implementation itself. Usually C++ would be your number one choice.

What is most efficient will be depending on your language, but in general,

- try to avoid memory allocation

- global variables

- type instability

- use fast libraries (e.g. pythons numpy is mostly written in C, and uses neat things like SIMD instructions to boost performance)

Most efficient?

I cannot claim I know the fastest (most efficient) way but Voladynamics claims to price

> the whole US options universe on one box in a fraction of a second.

Commonly used pricing methods

- Solving the PDE should also be reasonably fast.

- What is frequently used is Cox Ross Rubenstein (CRR), which is a recombining binomial tree, introduced in 1979 and actually a particular case of an explicit FDM scheme (a special case of the FDM for the BS PDE). It's in essence just a discrete time approximation of the continuous process underlying the BS model.

- CRR has also been modified because the convergence of the binomial tree based value to the limit is not monotone but rather oscillatory. This observation was the basis for a method developed by Leisen and Reimer (LR) (1996) to compute accurate results with a "minimum" number of time steps.

- There are some papers when it comes to large scale pricing engines, see for example Andreasen, American Option Pricing in a Tick from Saxo Bank or Bank of America Merrill Lynch where the author claims to be able to improve the efficiency of American option pricing algorithms by at least 4 orders of magnitude.

Andreasen claims that

> The best we can do with FD pricing is probably O(1ms) CPU time per pricing.

On the other hand, their solution allows to

> price O(100,000) American options per second. simultaneously calibrate to dividends and implied volatilities of approx. 200 American calls and puts in O(0.01s).

Efficient languages

I'll copy paste something I did elsewhere about the impact of language and implementation to demonstrate the difference with simple examples. It was originally done to compare computation time between Julia and Python and explain why Julia is so much faster.

Many people use Python. However, the language is generally quite slow. Unlike Python, writing C, C++ or even Julia is significantly more complex when it comes to potential code implementation.

In a nutshell, everything that runs fast in Python is not Python, but C, Fortran and co. Julia can also turn slow when type inference fails or the JIT compiler has insufficient information to optimize effectively. If you notice no speedup to Python, it is either because you use examples (or benchmark that way) where compile time latency matters more than computation, or it is inefficient code.

Dynamic typing

In dynamic languages like Python, classes could be subclassed. This is also possible in some statically typed languages that are object oriented such as Java. For this reason, a single integer in `Python 3.x` actually contains four pieces:

- `ob_refcnt`, a reference count that helps Python silently handle memory allocation and deallocation

- `ob_type`, which encodes the type of the variable

- `ob_size`, which specifies the size of the following data members

- `ob_digit`, which contains the actual integer value that we expect the Python variable to represent.

This means that there is some overhead in storing an integer in Python as compared to an integer in say Julia or C. Python uses 28 bytes to store an integer, Julia only 8 bytes.

A great feature of Julia is that types are extremely lightweight, and user-defined types are exactly as performant as anything built-in: this is a big difference to languages like Python where constructing a new instance of a class is very expensive. The amount of data required to store a Struct (Class) called Foo is the same as the amount of data required to store a pointer to a string. At compilation, Julia (just like Fortran and C) doesn't have additional overhead contrary to how an object (in say Python) is represented in memory (includes type information, reference counter etc).

Julia:

Python:

You can see the massive overhead in Python. It's probably even better to say that Julia isn't slow (not that it is fast), meaning it doesn't get in the way of the hardware and lets the hardware do its job efficiently. Mathematically, the integer 5 is exactly identical to the floating point 5.0, but it has a different bitwise representation. This is why Types are important in computers. Slow programming languages have a very hard time to know in advance what types things are. So they test for all sorts of combinations to figure out what routine to use to call the correct CPU instructions, because at the lowest level, the microprocessor has to perform different instructions when dealing with for example integers or floats. A microprocessor does not know types. A CPU needs special instructions for all different types; e.g. signed and unsigned integers, floats with single precision or double precision, or combinations of floating points and integers, for which you will need conversion (promotion) of types. Usually there will be different special instructions for every possible type.

As long as your task is computational, Julia will always beat Python (provided both are optimized). Where it starts to matter is when you really work with lots of data and / or complex models.

Some examples

For example look at the question Why is this task faster in Python than Julia? on stackoverflow. It is almost always code written in non-performant ways if you get slow Julia code. Without any optimization, just getting rid of compile time latency, resulted in Julia taking 8s vs 30s in Python. Regarding the stackoverflow question, Quantinsti has some good tests for Multiple operations on large datasets. For 100 groups of ~10,000,000 rows, Python (pandas package) and R (dplyr package) resulted in an internal error and out of memory error, respectively while Julia took 2.4 seconds the first time and 1.8 seconds the second time.

What makes a language fast?

Julia is fast because of its design decisions that allow a compiler to make efficient code. This is achieved with type-stability through specialization via multiple-dispatch as a core design decision. To illustrate this, you can try to compute `2^3`.

In Julia, you can run `@code_llvm` as well as `@code_native` macros to print the LLVM bitcodes as well as the native assembly instructions generated for running the method matching the given generic function and type.

If you use `@code_llvm exponential(2,3)` as well as `@code_llvm exponential(2,3)` and compared it to the 2^3 counterpart, you will be surprised how messy the “manual” exponential solution is compared to Julia’s built in.

Manual (ignoring almost half the output):

Built-in (entire output):

The built in solution is something you can decipher fairly quickly with some basic assembly.

- The lines beginning with `;` are comments, and explain which section of the code the following instructions come from. Calling it with `debuginfo=:none` removes the comments in the assembly code

- `pushq` , `movq`, `subq` and `addq` are to save stack frames, move content to or from registers (the suffix `q` is short for `quad`, which means 64-bit integer), allocate to registers and dealocate. Neither is important for our purpose. For short, a CPU operates only on data present in registers. These are small, fixed size slots (e.g. 8 bytes in size) inside the CPU itself. A register is meant to hold one single piece of data, like an integer or a floating point number.

- `movabsq` is the so called Exponentiation by squaring. For example if you compute $x^{15}$, it will look like this

```
x^15 = (x^7)*(x^7)*x 
x^7 = (x^3)*(x^3)*x 
x^3 = x*x*x
```

Or only the last line in our simple example above, which is $x^3$ only.

The reason this works is that Julia's built in function made some clever design decisions. For instance, `x^-n` works for an integer `x` and literal `n` but not for expressions because Julia developers defined a special syntactic treatment of literals.

- The benefit is quick and reliable exponentiation for something like $x^-3$. The advantage is that literal powers automatically become type-stable, because they turn into a fixed number of * calls. This is why Julia now lowers `x^-3` to `inv(x)^3`.

- The downside is that $p = -3$ and $x^p$ does not work and throws an error which seems potentially quite confusing. However, the `x^literal` has a different meaning than `x^expression`. In essence, referential transparency was sacrificed, and type stability "extended": That is why `^` to a literal integer power is different than raising to a variable with the same integer value.

Note: The crux here is that type-stability doesn’t mean that the type returned by a function is the same as the input. It means that the type can be inferred at every step of the way. If inv of an integer is always a floating point number, then that’s type stable.

The problem however, is that most casual Python users will expect x^p to work for p = -3.

Similarly, scope matters a lot. You can look at this simple function.

```
b = 1.0
function g(a)
    global b
    tmp = a
    for i in 1:1_000_000
        tmp = tmp + a + b
    end
    return tmp
end
```

b is global here, which is like poison for performance because the type of a global variable is never certain.

If you instead write

```
function g(a, b)
    tmp = a
    for i in 1:1_000_000
        tmp = tmp + a + b
    end
    return tmp
end
```

you are eliminating the global variable which

• reduces the number of allocations to zero (from 3000001 allocations: 45.78 MiB) • speeds up execution and • produces clean and fairly simple machine code.

For short, with some care, Julia can be on par with the fastest languages.

These allocations matter massively. On my private laptop for example, finding a single file (including opening and closing the file) takes about 1.4 milliseconds. Accessing 1,000,000 integers from memory takes 75.7 milliseconds. So my RAM is almost 20,000 times faster than my disk. I ran this with a, at least back in the time, reasonably fast SSD but even the newer Optane technology disks are usually thousands of times slower than RAM.

Simple speed comparison between C, Julia and Python

I ran a speed test between C, Python and Julia, which is based on a notebook from an introductory MIT math course.

It is a speed comparisons between C, Python and Julia by implementing a sum function `sum(a)`, which computes

$$ \mathrm{sum}(a) = \sum_{i=1}^n a_i $$

for an array `a` with `n` elements. I compared built in `sum` functions in Julia and Python along with hand-coded implementations in C, Python, and Julia. If you are a Windows user and want to replicate this on your own, you will need to install the MinGW (GCC/G++) compiler.

The C code is:

```
C_code = """
#include <stddef.h>
double c_sum(size_t n, double *X) {
    double s = 0.0;
    for (size_t i = 0; i < n; ++i) {
        s += X[i];
    }
    return s;
}
"""

const Clib = tempname()   # make a temporary file

# compile to a shared library by piping C_code to gcc
open(`gcc -fPIC -O3 -msse3 -xc -shared -o $(Clib * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code) 
end
```

Julia's

```
function my_sum1(x)
    result = zero(eltype(x))
    for element in x
        result += element
    end
    return result
end
```

Not only is the Julia code significantly easier, but also, unlike C, it works in Julia for any type of array (any iterable that provides an `eltype` method).

Python (within a Julia Kernel using PyCall)

```
py"""
def mysum(a):
    s = 0.0
    for x in a:
        s = s + x
    return s
"""
mysum_py = py"mysum"
```

The results are

There are a few interesting aspects.

- Starting with C, testing the function on an array of $10^7$ random numbers in [0,1] gives a result of ≈ 10ms for summing these numbers, which means that we have about 1 billion additions per second. This sounds a lot, but is well below the clock-speed of 2.5GHZ my processor has. This is because each Floating-Point addition needs to perform several additional calculations to load the next element of the array from memory, plus the time needed to access memory itself.

- Probably most interesting is that Python's built in sum is fairly slow, even though the Python sum function is written in C (Julia's sum is built using only Julia). It takes almost 4x longer than the equivalent C code and allocates memory. Python code pays a price for being generic and being able to handle arbitrary iterable data structures. See above for memory allocation differences between Julia and Python. Therefore, there are not only the computations involving addition, but also the overhead from fetching each item from memory.

- Julia equivalent of a Python list is a Vector{Any}: internally, this is an array of pointers to "boxes" that can hold any type (Any). This makes things much slower: each + computation on an Any value must dynamically look up the type of object, figure out what + function to call, and allocate a new "box" to store the result. If you run this logic in Julia, it is even slower than Python. This is the same problem as with the manual Numpy array implementation. Python is faster (better optimized for dealing) with untyped `Any` values. This is because Julia expects you to use "concretely" typed arrays in all performance-critical cases and why you should always avoid abstract types in Julia.

- On the other hand, Numpy arrays can take advantage of the fact that all of the elements are of the same type. This is actually faster than C! The reason is that NumPy gets and extra turbo boost from exploiting SIMD instructions.

You can run SIMD for C as well as Julia as well.

Julia

```
function mysum(A)
    s = zero(eltype(A)) 
    @simd for a in A
        s += a
    end
    return s
end
```

C

```
const Clib_fastmath = tempname()   # make a temporary file

# The same as above but with a -ffast-math flag added
open(`gcc -fPIC -O3 -msse3 -xc -shared -ffast-math -o $(Clib_fastmath * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code) 
end

# define a Julia function that calls the C function:
c_sum_fastmath(X::Array{Float64}) = ccall(("c_sum", Clib_fastmath), Float64, (Csize_t, Ptr{Float64}), length(X), X)
```

However, once you work with SIMD, you really need to know what you do when you implement and use this. Incorrect use of the @simd macro may cause unexpected results.

While this may sound like esoteric computer science stuff, you can look at the following question to see what can happen with floating point math results in a reasonable normal situation. The most intriguing example demonstrating the dangers of floating point math I have seen so far comes from Stefan Karpinksi. You can use this to illustrate the difference between left-to-right summation and left-to-right summation with a @simd annotation.

Essentially, if you use this code snippet in Julia,

```
function sumsto(x::Float64)
    0 <= x < exp2(970) || throw(ArgumentError("sum must be in [0,2^970)"))
    n, p₀ = Base.decompose(x) # integers such that `n*exp2(p₀) == x`
    [floatmax(); [exp2(p) for p in reverse(-1074:969) if iseven(n >> (p-p₀))]

    -floatmax(); [exp2(p) for p in reverse(-1074:969) if isodd(n >> (p-p₀))]]
end    

s = rand()*10^6
v = sumsto(s)
```

you order the random vector (10^6) from highest to lowest value. None of the standard algorithms in Julia or Python provide the true sum, although interestingly, in this example SIMD comes closest in this case. The reason is that in the case of summation, strict left-to-right summation is actually one of the least accurate algorithms, so putting @simd on the loop makes it faster and more accurate at computing the true sum because SIMD largely allows for reassociation of mathematically associative operations which are not actually associative in floating point. This includes `+` and `*`. Generally, `(a + b) + c` can be evaluated as `a + (b + c)` instead, but this does not compute the same thing for floating-point values.

The correct answer can be obtained with the Kahan-Babuska-Neumaier (KBN) summation algorithm for computing sums (and cumulative sums). This comes at a speed difference of about 20x, which is why none of the conventional implementations uses it.

Bottom line is, while Julia comes with a great mix of speed and usability, the main problem is that it is not trivial to reach C (C++ or Fortran speed) and you need to pay a lot of attention to detail to not slow the language down.

If you require real speed, firms will prefer C++, also because it's significantly easier go get C++ devslopers compared to Julia (at least in my field, quant finance, C++ is the by far most frequently used language, with some niche products like OCAML used by Lexifi or Bloomberg's DLIB, which focus on complex derivs).

## Answer by marie albertini (score 2)

https://quant.stackexchange.com/a/81629

rather than use finite difference you could have a look at the ODgrid method, which is flexible in terms of dividend and rates and fairly quick and accurate.

it is a mix between trinomial trees and pde solving

you'll fins more at https://arxiv.org/pdf/2411.05425