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[0]):
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[0]):
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
