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.505K=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−F0)N(−d)+σTn(d)P(K) = (K - F_0)N(-d) + \sigma\sqrt{T}n(d) where K,F0,σ,TK, F_0, \sigma, T is the strike, forward price, volatility and time to expiry, respectively. N(d)N(d) is the CDF of a standard normal distribution, n(d)n(d) is the PDF and d=(F0−K)/σTd=(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)p = (1 - d)/(u - d) where u,du,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 ```