Estimating American SOFR Option Volatility with a Consistent Tree

# Calculate implied volatility of american option on interest rate futures


# Calculate implied volatility of american option on interest rate futures












I want to calculate implied volatility of american option of a short term interest rate future.

Let's take for example a put option for a SOFR future with $K=95, price=0.105, T=0.750685, underlying=95.505$

I currently use as a first approximation the implied vol by using finding implied vol using the Bachelier model (used to price European options from my understanding). The undiscounted price for a put option is

$$ P(K) = (K - F_0)N(-d) + \sigma\sqrt{T}n(d) $$

where $K, F_0, \sigma, T$ is the strike, forward price, volatility and time to expiry, respectively. $N(d)$ is the CDF of a standard normal distribution, $n(d)$ is the PDF and $d=(F_0 - K)/\sigma\sqrt{T}$.

Obviously calculating the implied vol using above formula will give implied vol for european option which the SOFR option is not.

I want to now improve on this and calculate implied vol for american option. I now attempt to price the american option using a binomial tree and calculate implied vol like that.

The following is my attempt in python but im running into some problems;

- Im unsure of my assumption of the probability being 1/2 or if my binomial pricing is correct. The alternative I have seen is $p = (1 - d)/(u - d)$ where $u,d$ are the up and down price moves.

- i cant seem to find root for the implied vol for the binomial pricing even though i can for the analytical formula

Some guidance on where I'm going wrong is appreciated

```
import numpy as np
from scipy.stats import norm

def bachelier_binomial_tree_vectorized(F0, K, T, sigma, N, opt_type='C'):
    
    dt = T / N
    u = sigma * np.sqrt(dt)
    d = -u
    p = 0.5  
    
    # initialise asset prices at maturity - Time step N
    V = F0 * d ** (np.arange(N,-1,-1)) * u ** (np.arange(0,N+1,1))
    
    # Initialize option value tree at maturity 
    if opt_type == 'C':
        V = np.maximum( V - K , np.zeros(N+1) )
    elif opt_type=='P':
        V = np.maximum( K - V , np.zeros(N+1) )
    else:
        raise NotImplementedError(f'Unexpected type {opttype}')
    
    # Backward induction

    for i in np.arange(N,0,-1):
        V = ( p * V[1:i+1] + (1-p) * V[0:i] )
    
    return V[0]

def bachelier_price(F0, K, T, sigma, opt_type):
    d = (F0 - K) / (sigma * np.sqrt(T))
    if opt_type == 'C':
        price = ((F0 - K) * norm.cdf(d) + sigma * np.sqrt(T) * norm.pdf(d))
    elif opt_type == 'P':
        price = ((K - F0) * norm.cdf(-d) + sigma * np.sqrt(T) * norm.pdf(d))
    else:
        raise NotImplementedError(f'Unexpected type {opttype}')
    return price

def solve_for_sigma(F0, K, T, price, N, opt_type, model):

    def premium_error(sigma):
        if model == 'b':
            model_price = bachelier_binomial_tree_vectorized(F0, K, T, sigma, N, opt_type)
        elif model == 'a':
            model_price = bachelier_price(F0, K, T, sigma, opt_type)
        return model_price - price
        
    res = brentq(premium_error, a=0.000001, b=1000, full_output=True)

    if res[1].converged:
        return res[0]
    else:
        return np.nan
        

# Example usage
F0 = 95.505  # Initial forward price
K = 95   # Strike price
T = 0.750685     # Time to maturity 
N = 10  # Number of time steps
opt_type = 'P'
price = 0.105

sigma_analytical = solve_for_sigma(F0, K, T, price, N, opt_type, model='a')
print(f"sigma analytical: {sigma_analytical}")

option_price = bachelier_binomial_tree_vectorized(F0, K, T, sigma_analytical, N, opt_type)
print(f"Bachelier Binomial Tree Option Price using analytical sigma: {option_price:.4f}")

analytical_price = bachelier_price(F0, K, T, sigma_analytical, opt_type)
print(f"Bachelier Analytical Price: {analytical_price:.4f}")

sigma_binomial = solve_for_sigma(F0, K, T, price, N, opt_type, model='b')
print(f"sigma binomial: {sigma_binomial}")
```

The output i get is

```
sigma analytical: 0.8397460469518556
Bachelier Binomial Tree Option Price using analytical sigma: 95.0000
Bachelier Analytical Price: 0.1050
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[62], line 73
     70 analytical_price = bachelier_price(F0, K, T, sigma_analytical, opt_type)
     71 print(f"Bachelier Analytical Price: {analytical_price:.4f}")
---> 73 sigma_binomial = solve_for_sigma(F0, K, T, price, N, opt_type, model='b')
     74 print(f"sigma binomial: {sigma_binomial}")

Cell In[62], line 48, in solve_for_sigma(F0, K, T, price, N, opt_type, model)
     45         model_price = bachelier_price(F0, K, T, sigma, opt_type)
     46     return model_price - price
---> 48 res = brentq(premium_error, a=0.000001, b=1000, full_output=True)
     50 if res[1].converged:
     51     return res[0]

File ~\AppData\Local\miniconda3\envs\analytics\Lib\site-packages\scipy\optimize\_zeros_py.py:806, in brentq(f, a, b, args, xtol, rtol, maxiter, full_output, disp)
    804     raise ValueError(f"rtol too small ({rtol:g} < {_rtol:g})")
    805 f = _wrap_nan_raise(f)
--> 806 r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp)
    807 return results_c(full_output, r, "brentq")

ValueError: f(a) and f(b) must have different signs
```
```

## Answer by Rojolithos (score 0)

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

You can’t use 0.5 for the risk-neutral probability as it would be negating the effects of the risk-free rate. The equation below is correct, $$ p = \frac{1 - d}{u - d} $$

However, you also need to take into account early exercise, as the value of an American option at each node should be the maximum of the intrinsic value and the continuation value. Below is the fixed code

```
import numpy as np
from scipy.stats import norm
from scipy.optimize import brentq

def bachelier_binomial_tree_vectorized(F0, K, T, sigma, N, opt_type='C'):
    dt = T / N
    u = np.exp(sigma * np.sqrt(dt))
    d = 1 / u
    p = (1 - d) / (u - d)  # Risk-neutral probability

    # Initialize asset prices at maturity
    asset_prices = F0 * d ** np.arange(N, -1, -1) * u ** np.arange(0, N + 1, 1)

    # Initialize option values at maturity
    if opt_type == 'C':
        option_values = np.maximum(asset_prices - K, 0)
    elif opt_type == 'P':
        option_values = np.maximum(K - asset_prices, 0)
    else:
        raise NotImplementedError(f'Unexpected type {opt_type}')

    # Backward induction
    for i in range(N, 0, -1):
        asset_prices = asset_prices[1:] * u  # Adjust asset prices
        option_values = (p * option_values[1:] + (1 - p) * option_values[:-1]) * np.exp(-0.0 * dt)  # Discount factor can be added if needed
        if opt_type == 'C':
            option_values = np.maximum(option_values, asset_prices - K)
        elif opt_type == 'P':
            option_values = np.maximum(option_values, K - asset_prices)

    return option_values[0]

def bachelier_price(F0, K, T, sigma, opt_type):
    d = (F0 - K) / (sigma * np.sqrt(T))
    if opt_type == 'C':
        price = ((F0 - K) * norm.cdf(d) + sigma * np.sqrt(T) * norm.pdf(d))
    elif opt_type == 'P':
        price = ((K - F0) * norm.cdf(-d) + sigma * np.sqrt(T) * norm.pdf(d))
    else:
        raise NotImplementedError(f'Unexpected type {opt_type}')
    return price

def solve_for_sigma(F0, K, T, price, N, opt_type, model):

    def premium_error(sigma):
        if model == 'b':
            model_price = bachelier_binomial_tree_vectorized(F0, K, T, sigma, N, opt_type)
        elif model == 'a':
            model_price = bachelier_price(F0, K, T, sigma, opt_type)
        return model_price - price

    res = brentq(premium_error, a=0.000001, b=1.0, full_output=True)  # Adjusted bounds for better convergence

    if res.converged:
        return res.root
    else:
        return np.nan

# Example usage
F0 = 95.505  # Initial forward price
K = 95   # Strike price
T = 0.750685  # Time to maturity
N = 10  # Number of time steps
opt_type = 'P'
price = 0.105

sigma_analytical = solve_for_sigma(F0, K, T, price, N, opt_type, model='a')
print(f"sigma analytical: {sigma_analytical}")

option_price = bachelier_binomial_tree_vectorized(F0, K, T, sigma_analytical, N, opt_type)
print(f"Bachelier Binomial Tree Option Price using analytical sigma: {option_price:.4f}")

analytical_price = bachelier_price(F0, K, T, sigma_analytical, opt_type)
print(f"Bachelier Analytical Price: {analytical_price:.4f}")

sigma_binomial = solve_for_sigma(F0, K, T, price, N, opt_type, model='b')
print(f"sigma binomial: {sigma_binomial}")
```

## Answer by carry_and_pray (score 0)

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

From what I see, the root-finding is not the real problem but instead that your Bachelier tree is not actually a Bachelier tree.

In your code you set $u = \sigma \sqrt{\Delta t}$ and $d = -u$ but then you build terminal nodes at $F_N(j) = F_0 \cdot d^{N - j}u^j$ which is a multiplicative tree formula while Bachelier / normal dynamics are additive. You're mixing a normal-model formula with a CRR-style lattice which is why you see your tree price blow up to something like 95 instead of something near 0.105 and then `brentq` cannot bracket a root.

For a Bachelier tree, the futures price should evolve additively like $F_{n + 1} = F_n \pm \sigma \sqrt{\Delta t}$ so the terminal nodes are $F_N(j) = F_0 + (2j - N) \sigma \sqrt{\Delta t}$ for $j = 0, \dots, N$.

If you want an American option, you then do backward induction with early exercise $$ V_n(j) = \max(\text{intrinsic}, D_{n, n+1} (pV_{n+1}(j+1) + (1-p)V_{n+1}(j))). $$ If you are working in the same undiscounted setup as your Bachelier formula, then $D_{n, n+1} = 1$. If you want PV instead, discount consistently in both the analytic formula and the tree.

Also, do NOT replace $p = \frac{1}{2}$ by $p = \frac{1-d}{u-d}$ unless you also switch the whole tree to a multiplicative Black / CRR tree. That probability formula belongs to a log-normal lattice and not a normal one. In a centered additive Bachelier tree for a futures price with zero drift, $p = \frac{1}{2}$ is the consistent choice.

There is also a contract-convention point, that is, CME says that the Three-Month SOFR options are American-style and their SOFR options were designed to mirror Eurodollar options in nearly all respects. OpenGamma's CME Eurodollar note also makes the standard STIR convention clear which is that a call on the quoted IMM price is a put on the rate because $\text{rate} = 100 - \text{price}$.

So you have to essentially decide whether you want a normal/Bachelier model or a lognormal/Black-CRR model, then build the tree to match that model and add early exercise with `max(intrinsic, continuation)` and only then solve for implied vol.

```
def bachelier_price(F0, K, T, sigma, opt_type='P'):
    d = (F0 - K) / (sigma * np.sqrt(T))
    if opt_type == 'C':
        return (F0 - K) * norm.cdf(d) + sigma * np.sqrt(T) * norm.pdf(d)
    else:
        return (K - F0) * norm.cdf(-d) + sigma * np.sqrt(T) * norm.pdf(d)

def bachelier_american_tree(F0, K, T, sigma, N, opt_type='P', discount=1.0):
    dt = T / N
    step = sigma * np.sqrt(dt)
    p = 0.5

    F = F0 + (2 * np.arange(N + 1) - N) * step

    if opt_type == 'C':
        V = np.maximum(F - K, 0.0)
    else:
        V = np.maximum(K - F, 0.0)

    for i in range(N - 1, -1, -1):
        F = F0 + (2 * np.arange(i + 1) - i) * step
        cont = discount * (p * V[1:i+2] + (1 - p) * V[0:i+1])

        if opt_type == 'C':
            exer = np.maximum(F - K, 0.0)
        else:
            exer = np.maximum(K - F, 0.0)

        V = np.maximum(cont, exer)

    return V[0]

def implied_vol_tree(F0, K, T, price, N, opt_type='P'):
    f = lambda sigma: bachelier_american_tree(F0, K, T, sigma, N, opt_type) - price
    return brentq(f, 1e-8, 5.0)
```