This is more of an academic question. The results are SO close, I think they are ACTUALLY THE SAME FORMULAS. So someone published a paper with a "new" method to adjust Kirk's formula to make it more accurate, and the reality is, I think it's just a rearranged version of another paper's formula. Or that's what it appears to be... check my tests and judge for yourself.
I read up on "Adjusted Kirk's formula" for spread options and decided why not code it up in Python, as it's supposed to fix the errors that come with the original Kirk's formula and tie out almost exactly to a MC simulation. Paper by Elisa Alòs & Jorge A. León, 2013. "On the closed-form approximation of short-time random strike options," pages 24-26: https://www.tandfonline.com/doi/full/10.1080/14697688.2015.1013499, with code adapted from: https://arxiv.org/ftp/arxiv/papers/1812/1812.04272.pdf. Now after doing so, I had already coded up Bjerksund-Stensland's 2008 spread option model https://papers.ssrn.com/sol3/papers.cfm?abstract_id=1145206 and started to compare the results... and I found, surprising enough, they were always within machine precision of each other. I tried over numerous strikes and underlyings, correlations, and call/puts. I tried to recreate the comparisons in the paper I read with a range of strikes and correlations, and added some varying underlying prices on top of it.
Try for yourself. Just finding this quite amusing... I added weights to both and also figured out the put formula for the Adjusted Kirk's, although it wasn't in the paper. I'm just shocked that this paper got published with all the math geniuses out there - why did no one else compare the 2 models? I'm certainly not a mathematician, and don't claim to be. Or did the paper I copied get the formulas wrong?
""" Bjerksund-Stensland spread option model callput has to be 1 call or -1 put modifications made to input a leg2 weight (positive) for example: heat rate option spreads are assumed to be +leg1 -leg2 """ from math import sqrt, erf, log, exp def SpreadBJS(F1, F2, vol1, vol2, K, T, IR, corr, leg2_weight_ratio, callput): a = (F2*leg2_weight_ratio) + K b = (F2*leg2_weight_ratio) / (F2*leg2_weight_ratio + K) spreadvol = sqrt(vol1*vol1 - 2*b * corr * vol1 * vol2 + b*b * vol2*vol2) d1 = (log(F1 / a) + (0.5 * vol1*vol1 - b* corr * vol1 * vol2 + 0.5 * b*b * vol2*vol2 ) * T) / (spreadvol * sqrt(T)) d2 = (log(F1 / a) + (- 0.5 * vol1*vol1 + corr * vol1 * vol2 + 0.5 * b*b * vol2*vol2 - b * vol2*vol2) * T) / (spreadvol * sqrt(T)) d3 = (log(F1 / a) + (-0.5 * vol1*vol1 + 0.5 * b*b * vol2*vol2) * T) / (spreadvol * sqrt(T)) return exp(-IR * T) * (callput * F1 * normcdf(callput * d1) - callput * (F2*leg2_weight_ratio) * normcdf(callput * d2) - callput * K * normcdf(callput * d3)) """ Modified Kirk's supposed to be more accurate than Kirk's https://arxiv.org/ftp/arxiv/papers/1812/1812.04272.pdf MOST RESULTS MATCH Bjerksund-Stensland within 1e-7, guessing it's really the same formula... """ def SpreadKirkMod(F1, F2, vol1, vol2, K, T, IR, corr, leg2_weight_ratio, callput): F2 = F2*leg2_weight_ratio a_kirk = sqrt(vol1**2-2*corr*vol1*vol2*(F2/(F2+K))+vol2**2*((F2/(F2+K))**2)) xt = log(F1) xstart = log(F2+K) ihatt = sqrt(a_kirk**2) + 0.5 *(((vol2 * F2/(F2+K)) - corr*vol1)**2 ) * (1/((sqrt(a_kirk**2))**3)) * (vol2**2) * ((F2*K)/((F2+K)**2)) * (xt - xstart) S = F1/(F2+K) d1 = (log(S)+ 0.5 *(ihatt**2)*T) / (ihatt * (sqrt(T))) d2 = d1 - ihatt * sqrt(T) return exp(-IR*T) * ((callput * F1 * normcdf(callput * d1)) - callput * ((F2+K) * normcdf(callput * d2))) # fast norm cdf function def normcdf(x): return (1+erf(x/sqrt(2)))/2
Give it a shot with this test script:
import numpy as np F1 = np.arange(10,100,5) F2 = np.arange(10,46,2) vol2 = np.arange(0.95,0.05,-0.05) vol1 = np.arange(0.05,0.95,0.05) corr = np.arange(-0.9,0.9,0.1) for i in range(F1.shape): print(SpreadBJS(F1[i],F2[i],vol1[i],vol2[i],0,0.05,0.5,corr[i],1,1)-SpreadKirkMod(F1[i],F2[i],vol1[i],vol2[i],0,0.05,0.5,corr[i],1,1)) for i in range(F1.shape): print(SpreadBJS(F1[i],F2[i],vol1[i],vol2[i],0,0.05,0.5,corr[i],1,-1)-SpreadKirkMod(F1[i],F2[i],vol1[i],vol2[i],0,0.05,0.5,corr[i],1,-1))
Literally no difference:
0.0 1.7763568394002505e-15 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.7763568394002505e-15 4.440892098500626e-16 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
And I covered a bunch of option values:
0.8641109341058706 3.157008938949211 5.922236376518944 8.801436936706379 11.71190938228034 14.632399936782699 17.556431842920915 20.48174418431956 23.407495315452923 26.33338023040011 29.259300092551516 32.185227795685435 35.111157105827544 38.0370867785629 40.96301665970048 43.88894723101272 46.81488185033222 49.74084103409464 0.8641109341058706 0.23107920286421468 0.0703769043489477 0.02364772845138316 0.00819043794034862 0.0027512563577136934 0.0008534264109282218 0.00023603172457514483 5.742677293366057e-05 1.2605635129521401e-05 2.7317015360618936e-06 6.987504636992223e-07 2.7280756718442266e-07 2.0945792974138515e-07 3.545105183808255e-07 1.1897377637357556e-06 6.0729722538421966e-06 3.552064967292224e-05