Calculate implied volatility of american option on interest rate futures

External question — Quantitative Finance Stack Exchange Author: Naim Hussain Original post: https://quant.stackexchange.com/questions/80125 License: CC BY-SA 4.0 — https://creativecommons.org/licenses/by-sa/4.0/ Adaptation: HTML converted to plain text; contact email addresses removed. 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 ```

External answer — Quantitative Finance Stack Exchange Author: Rojolithos Original post: https://quant.stackexchange.com/a/80148 License: CC BY-SA 4.0 — https://creativecommons.org/licenses/by-sa/4.0/ Adaptation: HTML converted to plain text; contact email addresses removed. 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}")

External answer — Quantitative Finance Stack Exchange Author: carry_and_pray Original post: https://quant.stackexchange.com/a/85611 License: CC BY-SA 4.0 — https://creativecommons.org/licenses/by-sa/4.0/ Adaptation: HTML converted to plain text; contact email addresses removed. 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)