Hello opengamma developers,
You’ll note I put a question mark in the subject because I am not entirely sure of my proper use of the model and of my interpretation of results. However, what I think I see is that the class 1) calculates incorrect values; 2) fails on sigma of zero (a valid input); and 3) which causes the implied volatility function to crash in some cases.
- Here are actual market values for an IBM July call, with fabricated rate and volatility to test with:
Strike = 185.0
Spot = 196.25
Rate = .03
Volatility = .16
Time = .25924
Test code:
public static void main(String[] args) {
BjerksundStenslandModel model = new BjerksundStenslandModel();
Double tv = model.price(196.25, 185, .03, 0, .25924, 0.16, true);
System.out.println("TV = " + tv);
Double iv = model.impliedVolatility(14.3, 196.25, 185, .03, .0, .25924, true);
System.out.println("IV = " + iv);
}
The resulting theoretical value is 13.289.
Here is why I suspect this value. In my search for Java option value and IV code that I can trust (unsuccessful so far!), I have tested (besides opengamma) javaquant.net, jquantlib, and maygard. (Maygard is merely a transcription of the undocumented code from the book “Java Methods for Financial Engineering (Barker)”). Here are the results:
opengamma (BS) 13.289
maygard (BS) 13.270
jquantlib (BS) 14.464
binomial (javaquant & various online including CBOE): 14.465 to 14.475
This is a huge discrepancy (~10%) that can’t be attributed to roundoff or inherent inaccuracies in the underlying approximation. The reasons I suspect that both opengamma and maygard are delivering an incorrect calculation are:
- Both are far from the binomial and javaquant BS calculations, which agree.
- The javaquant test class verifies its calculations with those from a textbook, while opengamma’s tester seems to check only for internal consistency, without checking for correctness of result.
- I made corresponding tests of the Barone-Adesi & Whaley (BAW) estimation. The same pattern of results obtain: opengamma and maygard agree at 13.295, while javaquant and jquantlib yield 14.612 and 14.464.
What appears to be happening is that opengamma and maygard use the same calculations for BS and BAW. The results from jquantlib, javaquant, and online calculators are consistent with each other across the BS, BAW, and binomial methods, as well as in one case verified against textbook examples. Therefore I am inclined to conclude that opengamma’s and maygard’s calculations are incorrect.
A question for the developers of this class: did you derive the method from the original Bjerksund and Stensland paper, or take it from maygard or the Barker book? Frankly, I suspect that the method and code from the book may be wrong. I myself do not have the skills to verify the methods against the original paper, especially since not one of the four source codes contains even a word of documentation explaining their translation from formulas to code. But I’ve heard that Java is self-documenting (;/.
I’d appreciate hearing back your thoughts on this analysis. Or if you think the model is correct, please show me the flaw in my programming or thinking.
The second problem with BjerksundStenslandModel is that the price function fails with an uncaught exception when sigma==0. Simply edit the code line from the example above:
Double tv = model.price(196.25, 185, .03, 0, .25924, 0.0, true);
The problem is division by zero in method getCallPrice:
final double y = 0.5 - b / sigmaSq;
A volatility of zero is entirely valid: it’s the intrinsic value of the option, 11.25 in this example. The model correctly approaches 11.25 as sigma->0.0, just try .0001 for sigma. The issue is that the method blows up at 0.0.
The problem with the above bug is when you use the GenericImpliedVolatiltySolver with price. The root search correctly sends volatility 0.0 to the pricing function, causing it to fail.
In this case, I substituted price for getPriceAndVega, in an attempt to gain some efficiency. But the original pricing function, getPriceAndVega, also blows up at sigma==0, most likely for the same mathematical reasons that price does. Then the solver falls back to a bisection method which also fails. Try it (code above).
BTW, there is a valid IV of 0.193, which can be found by a manual root search. I did try altering the price method so that when sigma==0, it is set to .001. This converts price to a continuous function over all its proper domain, and allows the root to be found. However, that’s an egregious hack that I would never offer as a solution. The mathematical formula cited in the class also blows up at sigma==0, so I think you will need a correct mathematical solution for that case. It’s beyond my skills to derive one.
Thanks for taking the time to look at these claims, and looking forward to comments or solutions.
Doug