Discrete dividens for american options

[mlambelho]

My goal is to compute american option prices with discrete dividends on a specific date. The example provided for American options does not indicate how could this be done.

# Define the coordinate grid
  s_min = 0.01
  s_max = 300.
  grid = pde.grids.uniform_grid(minimums=[s_min],
                                maximums=[s_max],
                                sizes=[number_grid_points],
                                dtype=dtype)

  # Define the values grid for the final condition
  s = grid[0]
  final_values_grid = tf.nn.relu(s - strike)

  # Define the PDE coefficient functions
  def second_order_coeff_fn(t, grid):
    del t
    s = grid[0]
    return [[volatility ** 2 * s ** 2 / 2]]

  def first_order_coeff_fn(t, grid):
    del t
    s = grid[0]
    return [risk_free_rate * s]

  def zeroth_order_coeff_fn(t, grid):
    del t, grid
    return -risk_free_rate

  # Define the boundary conditions
  @pde.boundary_conditions.dirichlet
  def lower_boundary_fn(t, grid):
    del t, grid
    return tf.constant(0.0, dtype=dtype)

  @pde.boundary_conditions.dirichlet
  def upper_boundary_fn(t, grid):
    del grid
    return tf.squeeze(s_max - strike * tf.exp(-risk_free_rate * (expiry - t)))

  # In order to price American option one needs to set option values to 
  # V(x) := max(V(x), max(x - strike, 0)) after each iteration
  def values_transform_fn(t, grid, values):
    del t
    s = grid[0]
    values_floor = tf.nn.relu(s - strike)
    return grid, tf.maximum(values, values_floor)

Given that discrete options happen at a specific date I thought it would be possible to change the value of spot (s) by the respective dividend after the ex_date occurs. Lets say that expiry = 0.8 and for example ex_date = 0.5, with a dividend of 1.5%, then the code would become:

# Define the coordinate grid
  s_min = 0.01
  s_max = 300.
  grid = pde.grids.uniform_grid(minimums=[s_min],
                                maximums=[s_max],
                                sizes=[number_grid_points],
                                dtype=dtype)

  # Define the values grid for the final condition
  s = grid[0]
  final_values_grid = tf.nn.relu(s - strike)

  # Define the PDE coefficient functions
  def second_order_coeff_fn(t, grid):
    s = grid[0]
    if t.numpy() > 0.5:
        s = s*(1 - 0.015)
    del t
    return [[volatility ** 2 * s ** 2 / 2]]

  def first_order_coeff_fn(t, grid):
    s = grid[0]
    if t.numpy() > 0.5:
        s = s*(1 - 0.015)
    del t
    return [risk_free_rate * s]

  def zeroth_order_coeff_fn(t, grid):
    del t, grid
    return -risk_free_rate

  # Define the boundary conditions
  @pde.boundary_conditions.dirichlet
  def lower_boundary_fn(t, grid):
    del t, grid
    return tf.constant(0.0, dtype=dtype)

  @pde.boundary_conditions.dirichlet
  def upper_boundary_fn(t, grid):
    del grid
    return tf.squeeze(s_max - strike * tf.exp(-risk_free_rate * (expiry - t)))

  # In order to price American option one needs to set option values to 
  # V(x) := max(V(x), max(x - strike, 0)) after each iteration
  def values_transform_fn(t, grid, values):
    s = grid[0]
    if t.numpy() > 0.5:
        s = s*(1 - 0.015)
    del t
    values_floor = tf.nn.relu(s - strike)
    return grid, tf.maximum(values, values_floor)

Unfortunately, this approach does not yield proper results when compared with Quantlib. Does the tff library allow for the computation of option prices using discrete dividends? Ideally the only thing that it is required is to shift all spot prices based on the dividend starting at ex date.

Thank you in advance

[mlambelho]

@cyrilchim If possible you could provide some help? This should be something really trivial, but I am honestly missing how to do it with the library

[cyrilchim]

Hi @mlambelho
I'd need to think a bit on this one. Could you please check, e.g., here. Seems you also need to update upper_boundary_fn

Could you please share the full code you are running?

[mlambelho]

@cyrilchim here it is. I included the quantlib code for comparison.

https://colab.research.google.com/drive/1WiS1Lz259vv1UZZTYK5TBafG1MFir3WW?usp=sharing

[cyrilchim]

I looked in the colab. As I understand, you indeed need to adjust the grid by the dividend at the dividend date. You also need to adjust boundary conditions, values_transform_fn and the final_values_grid.

@tf.function(autograph=True)
def american_option(number_grid_points,
                    strike,
                    volatility,
                    risk_free_rate,
                    expiry,
                    dividend_time,
                    dtype=tf.float64):
  """ Computes American Call options prices."""
  # Define the coordinate grid
  s_min = 250
  s_max = 700
  grid = pde.grids.uniform_grid(minimums=[s_min],
                                maximums=[s_max],
                                sizes=[number_grid_points],
                                dtype=dtype)

  # Define the values grid for the final condition
  s = grid[0]
  final_values_grid = tf.nn.relu(
      s - 1.57 * tf.exp(-risk_free_rate * (dividend_time - expiry)) - strike)

  # Define the PDE coefficient functions
  def second_order_coeff_fn(t, grid):
    #del t
    #s = grid[0]
    s = grid[0]
    if t > dividend_time:
      s = s - 1.57
    del t
    return [[volatility ** 2 * s ** 2 / 2]]

  def first_order_coeff_fn(t, grid):
    #del t
    #s = grid[0]
    s = grid[0]
    if t > dividend_time:
      s = s  - 1.57
    del t
    return [risk_free_rate * s]

  def zeroth_order_coeff_fn(t, grid):
    del t, grid
    return -risk_free_rate

  # Define the boundary conditions
  @pde.boundary_conditions.dirichlet
  def lower_boundary_fn(t, grid):
    del t, grid
    return tf.constant(0.0, dtype=dtype)

  @pde.boundary_conditions.dirichlet
  def upper_boundary_fn(t, grid):
    #del grid
    if t > dividend_time:
      res = tf.squeeze(s_max - strike * tf.exp(-risk_free_rate * (expiry - t)) 
              - (1.57 * tf.exp(-risk_free_rate * (dividend_time - t))))
    else:
       res = tf.squeeze(s_max - strike * tf.exp(-risk_free_rate * (expiry - t)))
    return res


  # In order to price American option one needs to set option values to 
  # V(x) := max(V(x), max(x - strike, 0)) after each iteration
  def values_transform_fn(t, grid, values):
    s = grid[0]
    # Ensure shape of t is the same on both branches
    if t > dividend_time:
      s = s - 1.57 * tf.exp(-risk_free_rate * (dividend_time - t))
    else:
      s = tf.broadcast_to(s, tf.shape(values))
    values_floor = tf.nn.relu(s - strike)

    return grid, tf.maximum(values, values_floor)


  # Solve
  estimate_values, estimate_grid, _, _ = \
    pde.fd_solvers.solve_backward(
      start_time=expiry_dte,
      end_time=0,
      values_transform_fn=values_transform_fn,
      coord_grid=grid,
      values_grid=final_values_grid,
      num_steps=500,
      boundary_conditions=[(lower_boundary_fn, upper_boundary_fn)],
      second_order_coeff_fn=second_order_coeff_fn,
      first_order_coeff_fn=first_order_coeff_fn,
      zeroth_order_coeff_fn=zeroth_order_coeff_fn,
      dtype=dtype
    )

  return estimate_values, estimate_grid[0]

I estimate prices as

prices = (estimate[:, spot_index + 1] + estimate[:, spot_index - 1]) / 2

and relative absolute error of 1e-3 with QuantLib (I do 500 time steps on a 1024 point grid). I am sure QuantLib does a few more things like adjusting the grid on the fly which you could probably do as well. You'll need to check QuantLib source code for details. You could also ensure that the dividend time lies on the time grid.

Also, in your code you can try speeding up both pricing and Brent search with XLA compilation (I get 4x speed up for Brent and 6x for option pricing):

american_option_xla = tf.function(american_option, jit_compile=True)
brent_search = tf.function(tff.math.root_search.brentq, jit_compile=True)

Hope this helps.

[mlambelho]

Thank you for your help. If possible could you share the colab that you are running? jit_compile is much slower for me every time that I test it.

Could this be expandable for more than one dividend at one point in time? So two or three dividends?

On another note, is it possible to compute Greeks here? Things like theta, vega, gamma, delta?

Thank you in advance

[cyrilchim]

Edited your colab. Please check it out (see the last section there). Delta and Gamma should be estimated with finite difference as you already have price estimate on the grid. I've added that too. You can try computing gradients directly though it will be quite slow.

I've added theta and vega calculation using standard AD.

I have changed how dividend date is incorporated. This should make it clear how to add multiple dividend dates.

Hope this helps!

[mlambelho]

Hi @cyrilchim thank you very very much for your reply, it is super helpful. A couple of questions though:

1. If I increase the expiry date of the option to for example 2022-11-18 (so 5 months after the original one), the prices that the model shows are quite different than the quantlib ones (attach picture for comparison). As the DTE seems to increase, the mismatch between quantlib and tff seems to increase as well.
   
   Could it be that there is something missing in the formula? I was not able to find an error though.
   american_div.pdf

2. Is it possible to price both calls and put options at the same time? If so, how exactly? Both require different final values grids, boundaries and also values_transform_fn.

3. Is it possible to price different options at the same time as well? If an option requires a grid from 0 to 100 dollars and other option requires a grid from 1000 to 4000 dollars how can we incorporate pricing of many different options at the same time on such an approach?

Apologize for the questions, thank you very much for your help.

[cyrilchim]

Hi @mlambelho

1. Your grid boundaries are too narrow. Set something like s_min=10, s_max=900. The boundaries should be selected so that to cover most of the possible terminal values S(T). For Black-Scholes model S(T) has a normal distribution which can be used to build an appropriate interval. You can check how it is done in QuantLib (link)

2. Yes, you can do so. You'd need to use tf.where, something like done here

3. You can have a batch of grids like shown here. All grids should have the same number of points though.

Please let me know if this helps!

[mlambelho]

Thank you @cyrilchim, solutions 1 and 2 worked perfectly! Thank you! Regarding solution 3, I am sure your suggestion will also prove successful, I just can't seem to make it work. So, my understanding is that if we want to run in batches 20 options for example (10 puts and 10 calls), we need to provide a list with size 20 for s_min, s_max and sizes (where the number of grid points is always the same), correct?
So, in other words, imagine that we want to price 100 options and 50 of those options will have the same grid, but the other 50 should have a different grid other than the first one. The s_min and s_max lists length should be 100 with the first 50 elements being the same in all of them and the last 50 elements also being the same, right?

https://colab.research.google.com/drive/10Vu37T5f45E9yb1-LZNODFuAcTofmoy-?usp=sharing

I provide here what I was able to accomplish for other people possibly watching this thread as well. So the pricing of the put/call options at the same time seems to work pretty well. I used both tf.where and tf.multiply to be able to make it work.

However, in order to test the batch of grids, I converted the grid into three lists of 20 s_min, s_max and sizes (of course this does not make much sense since all batches have the same s_min and s_max, but that is not relevant for testing), but I am seeing other problems:

Shape must be at most rank 1 but is rank 2 for '{{node solve_backward/while/composite_scheme_step/cond/BroadcastTo_6}} = BroadcastTo[T=DT_DOUBLE, Tidx=DT_INT32](solve_backward/while/composite_scheme_step/cond/Squeeze, solve_backward/while/composite_scheme_step/cond/BroadcastTo_6/shape)' with input shapes: [20,20], [1].

I am not exactly sure of what could be causing the problem. I tested changing multiple other things, but nothing seems to be working for me. Thank you in advance again for your help. Much appreciated!

[cyrilchim]

Sorry for the late reply. Typically shape error means there is some shape mismatch somewhere. The best way to debug that is to remove tf.function and step through the code inspecting places where things might go wrong.

In your case, I simply printed out shapes for each of the input functions to the pde solver (second_order_coeff_fn, etc)

I noticed that upper_boundary_fn outputs a tensor of shape (20, 20) but it should be (20, 1). The issue is with

    res = tf.squeeze(tf.multiply((grid[0][..., -1] - strike * tf.exp(-risk_free_rate * (expiry - t))
                     + extra_term), is_call_options))

on line 104. There should be grid[0][..., :-1] and not grid[0][..., -1] .

    res = tf.squeeze(tf.multiply((grid[0][..., :-1] - strike * tf.exp(-risk_free_rate * (expiry - t))
                     + extra_term), is_call_options))

This should be it.

Also, where you get spot prices, you need a slice on the grid

spot_price = grid_locations[:, spot_index]

Please let me know if this works (it seems to work for me!)

[cyrilchim]

Also, I think you need to adjust the boundary conditions for puts. The lower boundary for put is

tf.squeeze((strike * tf.exp(-risk_free_rate * (expiry - t)) - grid[0][..., -1] 
                     + extra_term) * is_call_options)

and the upper boundary is zero.

[mlambelho]

I was able to make the adjustments you recommended, you can see the changes below:

https://colab.research.google.com/drive/10Vu37T5f45E9yb1-LZNODFuAcTofmoy-?usp=sharing

I think the pricing of the put options looks ok now.

If possible could you elaborate on how should I safely specify a grid? For example if a stock has spot 100, what should be the safest grid to specify that is neither too big nor too small? In my mind it would be possible to specify a grid just by knowing the current spot price.

Also I am seeing an issue in the notebook above where the brent optimizer with jit compile seems to be way slower vs without jit compile. Finding this quite strange.

Lastly, is it possible within this framework to price options at the same time with different dividends? I think it should be, the issue here is that I am not sure how is it possible to make each row of the tensor instead of using the same tf.where statement with the extra term, use a different tf.where statement per row (since each row would require different dividends).

Again, thank you very much for all the help.

[cyrilchim]

You could use similar approach to QuantLib (ql). Under Black-Scholes the distribution of the terminal values S(T) is known (it is Geometric normal). That means we can put some upper and lower boundaries such that S(T) is very unlikely to cross these. For eaxample, we can follow ql and put a boundary such that there is only 0.0001 chance crossing the boundary. We can compute the boundaries as

# shape [num_options]
vol = ... # volatilities
# shape [num_options]
r = ... # risk free rate
T = expiry_dte  # expiry

# shape [num_options]
spot = ... # desired spot values

# probability threshhold
prob = tf.constant(0.0001, dtype=spot.dtype)

# quantile for the normal distribution of log(S(T))
quantile = (r - vol**2 / 2) * T + np.sqrt(2) * vol * tf.sqrt(T) * tf.math.ndtri(1 - eps)

# Get the boundaries
min_val = spot * tf.math.exp(-quantile)
max_val = spot * tf.math.exp(quantile)

The grid can then be constructed as

# `number_grid_points // 2` points below spot value and number_grid_points // 2 above spot value
grid = tf.transpose(
    tf.concat([tf.linspace(min_val, spot, number_grid_points // 2),
               tf.linspace(spot, max_val, number_grid_points // 2)[1:]], axis=0))[0]

That is about it!

You will also need to update how prices and gradients are computed (see, e.g., here)

prices = estimate[:, spot_index]# (estimate[:, spot_index + 1] + estimate[:, spot_index - 1]) / 2
prices_plus = estimate[:, spot_index + 1]
prices_minus = estimate[:, spot_index - 1]

# Estimated deltas
d_plus = grid_locations[spot_index + 1] - grid_locations[spot_index]
d_minus = grid_locations[spot_index] - grid_locations[spot_index - 1]
c_plus = d_minus / ((d_minus + d_plus) * d_plus)
c_minus = -d_plus / ((d_minus + d_plus) * d_minus)

grid_delta = (grid_locations[spot_index + 1] - grid_locations[spot_index - 1])
D = 1 / grid_delta**2

deltas = (c_plus * prices_plus + c_minus * prices_minus
          - prices * (c_plus + c_minus)) 
# Estimated gammas
D_plus = 2 / ((d_minus + d_plus) * d_plus)
D_minus = 2 / ((d_minus + d_plus) * d_minus)
gammas = (D_plus * prices_plus + D_minus * prices_minus
          - prices * (D_plus + D_minus))

Let me know if this works (I get 9e-5 relative error compared to QL).

[mlambelho]

@cyrilchim I am coming back to you this late because I was trying for so long to make the grid work as you suggested, however I am not able to get it to work. I think something really small must be going on... But can't figure out what is. from the error message

https://colab.research.google.com/drive/1TsbOkm_LH_SgRRICaBUfh7nigEVyp6A8?usp=sharing

Basically there are 3 things that I would like to do that I still wasn't able:
1 - make a different grid for each option as per quantlib as a function of the spot value. So if there are 20 different options, ideally we would have 20 different grids
2 - apply different dividends per option on the same run. I can do it fine if it is just one. But more than one it is getting trickier to me. Well, I actually can make more than one dividend work... the problem for me is when one of the options has only one dividend, but others have 2 or 3
3 - make a fast "dummy" run to warm up jti so that when I want to run it right away, it becomes fast from the get go.

I tried to prevent asking for help again, but really not figuring out what is happening. Many apologies

[cyrilchim]

1. I've added a few edits to the colab. I made a mistake earlier. The shape of boundary function results should be [num_options] (that is 20 in your case). Please check that the results are correct
2. Dividends in your case can take the same shape as spot so each option can have its own dividend (see edit in the colab). If you have options of various number of dividends, the easiest way is to make all options have the same number dividends and set to zero the values that are not needed.
3. Function calibration_fn_xla needs to accept tensor values, not numpy arrays. You need to run calibration_fn_xla with some random inputs to 'warm' it up but with the shapes you are going to use in practice (in your case for 20 options and 1023 grid points) . Check this thread.

Please let me know if that makes sense

[cyrilchim]

I have added a section on CPU utilization optimization.

[mlambelho]

This was extremely helpful @cyrilchim. Thank you very much.

If possible, I have one additional question regarding how to calculate the deltas where the stocks have different grids.

  spot_index = int(number_grid_points/2)
  prices = estimate[:, spot_index]# (estimate[:, spot_index + 1] + estimate[:, spot_index - 1]) / 2
  prices_plus = estimate[:, spot_index + 1]
  prices_minus = estimate[:, spot_index - 1]

  # Estimated deltas
  d_plus = grid_locations[:, spot_index + 1] - grid_locations[:, spot_index]
  d_minus = grid_locations[:, spot_index] - grid_locations[:, spot_index - 1]
  c_plus = d_minus / ((d_minus + d_plus) * d_plus)
  c_minus = -d_plus / ((d_minus + d_plus) * d_minus)

  grid_delta = (grid_locations[:, spot_index + 1] - grid_locations[:, spot_index - 1])
  D = 1 / grid_delta**2

  deltas = (c_plus * prices_plus + c_minus * prices_minus
            - prices * (c_plus + c_minus)) 
  # Estimated gammas
  D_plus = 2 / ((d_minus + d_plus) * d_plus)
  D_minus = 2 / ((d_minus + d_plus) * d_minus)
  gammas = (D_plus * prices_plus + D_minus * prices_minus
            - prices * (D_plus + D_minus))

On this example, we are assuming that the spot_index is exactly at the center of the grid. What happens if that is not the case? How can deltas be computed if that is not the case?

The formula to calculate the grid relies on volatilities values. Meaning that even on the same stock/same expiry, the grids will vary. In such a case how can deltas be computed with AD?

quantile = (risk_free_rate - volatility**2 / 2) * expiry + np.sqrt(2) * volatility * tf.sqrt(expiry) * tf.math.ndtri(1 - tf.constant(0.0001, dtype=tf.float64))

Uploading here my link where I try to calculate deltas for 20 options which although they are the same stock, they have slightly different grids due to the different volatility term in the quantile formula above.

https://colab.research.google.com/drive/1WkFza_cfH1fC-B9jMT7RRQvi3vIaPZP9?usp=sharing

The 'correct' deltas were computed in a url above, leaving it here as well:

https://colab.research.google.com/drive/10Vu37T5f45E9yb1-LZNODFuAcTofmoy-?usp=sharing

Deltas that I am observing now:
<tf.Tensor: shape=(20,), dtype=float64, numpy=
array([ 0.63432159, 0.23296208, 0.25946409, 0.55911974, 0.56797372,
0.55901103, 0.51600366, 0.42287884, 0.41631046, 0.50195433,
-0.36743661, -0.77497464, -0.74827531, -0.44312381, -0.43424029,
-0.44358379, -0.48676025, -0.58105902, -0.58871937, -0.50144817])>

I suspect the reason for this discrepancy is related with the spot_index that is not aligned in this new approach of calculating the grid. If you could help me how can I calculate the delta taking into account that the spot will assume a different position in the grid (i.e spot_index is not int(number_grid_points/2) ) it would be very appreciated. many many apologies again. thank you

[cyrilchim]

Sorry for the late reply. I think you just need to adjust spot_index to a correct value. Is there an index which corresponds to the spot?

[mlambelho]

Hi @cyrilchim, my problem was indeed that spot_index was not in the correct index. Here it is the correct final solution.

https://colab.research.google.com/drive/1tfWjWQzP5ZS92x3zoMV7GhBbTlKUYeLk?usp=sharing

Everything seems to be working as expected, however it seems that I am only able to price more than 1 option. If I try to run it with a single option the code is failing.

ValueError: slice index 1 of dimension 0 out of bounds. for '{{node solve_backward/while/composite_scheme_step/strided_slice_6}} = StridedSlice[Index=DT_INT32, T=DT_DOUBLE, begin_mask=0, ellipsis_mask=0, end_mask=0, new_axis_mask=0, shrink_axis_mask=1](solve_backward/while/composite_scheme_step/sub, solve_backward/while/composite_scheme_step/strided_slice_6/stack, solve_backward/while/composite_scheme_step/strided_slice_6/stack_1, solve_backward/while/composite_scheme_step/strided_slice_6/stack_2)' with input shapes: [1,1022], [1], [1], [1] and with computed input tensors: input[1] = <1>, input[2] = <2>, input[3] = <1>.

Do you think you could help? I have tried debugging it a lot, but not sure where the error is coming from. I bet is just an indexing in one of the statements but I can't seem to be able to find it.

[cyrilchim]

Looking into it. There is probably some bug. Will report back once investigate

[cyrilchim]

Ok, should be fixed now! Please check it out. Thank you for spotting this one

[mlambelho]

Hi @cyrilchim . Thank you very much. It is working perfectly. The only issue that I am seeing is that some runs take a lot longer than other runs. I leave here a snippet as an example:

https://colab.research.google.com/drive/1L4YXephA0tg3XkSVbX4GD5o8klAxXCa2?usp=sharing

Is there any possibility to use dynamic shapes with the current structure of the code? Perhaps dynamic shapes might also make the first run go faster if the first run only has 1 option? I tried to use all of your suggestions across the code

Sorry again to bother again, if you could provide a bit of help would be really appreciated. Thank you very much in advance

[cyrilchim]

Sorry for keeping you waiting. You need to add input_signatures (check this) thread. This way your performance will not be impacted by input tensor shapes.

Could you please give more details on the use case? Seems like what you need is to wrap the pricer in tf.Module calss (see here).

In your case this will be something like

class AmericanPricer(tf.Module):

@tf.function(input_signature=[...])
def pricer(...):
  # Plug in you code here
  ...

model = AmericanPricer()

Save the model with tf.saved_mode.save(model, local_path, signatures={'serving_default': model.pricer})).

Once that is done, use TF serving to start local server for your model. The link contains copy-pastable code to do that.

Please let me know if that makes sense?

[magic-man]

Hi @cyrilchim is there any news on the American examples being formalized and made available? If not, what needs to be done as I am currently trying to piecemeal together a working example from this thread to get the IV, price and greeks for an american option with discrete dividends.