A few years back, over at the Wilmott forum, somebody requested some high precision reference prices for the Heston ’93 stochastic volatility model, so as to compare to an implementation. I responded to that with some prices computed in Mathematica to high precision. A number of people found that thread quite useful, and other cases were added over time. Unfortunately, the thread became corrupted after a forum revision. I have gotten a recent request for the same data, so I reproduce some of those results here as perhaps a more permanent record.
In SDE notation, the Heston model is the system
\displaystyle \left \{ \begin{array}{l}dS_t = (r-q) \, S_t \, dt + \sqrt{V_t} S_t \, dB_t, \\dV_t = (\omega - \theta V_t) \, dt + \xi \, \sqrt{V_t} \, dW_t, \quad dB_t \, dW_t = \rho \, dt. \end{array}\right. \quad\quad\quad (1)I use the parameter set:
\displaystyle r = \frac{1}{100},\quad q = \frac{2}{100}, \quad S_0 = 100, \quad T = 1, \quad V_0 = \frac{4}{100}, \omega = 1, \quad \theta = 4, \quad \xi = 1, \quad and \quad \rho = -\frac{1}{2}.With WorkingPrecision=50, I asked Mathematica for 20 good digits, which is probably a conservative estimate of how many good digits there are in the result. (The results have been confirmed by others to at least 15-16 good digits). With that, the results for Put and Call option values were:
| Strike | Type (P=Put, C=Call) | Option Price |
|---|---|---|
| 80 | P | 7.958878113256768285213263077598987193482161301733 |
| 80 | C | 26.774758743998854221382195325726949201687074848341 |
| 90 | P | 12.017966707346304987709573290236471654992071308187 |
| 90 | C | 20.933349000596710388139445766564068085476194042256 |
| 100 | P | 17.055270961270109413522653999411000974895436309183 |
| 100 | C | 16.070154917028834278213466703938231827658768230714 |
| 110 | P | 23.017825898442800538908781834822560777763225722188 |
| 110 | C | 12.132211516709844867860534767549426052805766831181 |
| 120 | P | 29.811026202682471843340682293165857439167301370697 |
| 120 | C | 9.024913483457835636553375454092357136489051667150 |
Subsequently, Mark Joshi, a well-known researcher who has since passed away, asked for more extreme cases with small T and small volatility, specifically T = V_0 = 0.01. With those changes, but otherwise the same parameters, I found:
| Strike | Type (P=Put, C=Call) | Option Price |
|---|---|---|
| 90 | P | 4.5183603586861772614990106188215872180542*10^-8 |
| 90 | C | 9.989001595065276544935948045293485530832966049263 |
| 95 | P | 0.000461954855653851579672612557018857858641926937 |
| 95 | C | 4.989963479738160122154264702582719627807098780529 |
| 100 | P | 0.477781171629504680023239655436072890669645669297 |
| 100 | C | 0.467782671512844263098248405184095087949465507760 |
| 105 | P | 5.009501052563650299130635110520904481889436667608 |
| 105 | C | 2.527447823194706060519991248106500619490942*10^-6 |
| 110 | P | 10.008998550115123724684210555728039829315964456261 |
| 110 | C | 1.29932760052624920704881258510264466*10^-13 |
Mark found agreement with his results to 6 digits.
Update (Feb 21, 2019). Recently, Wilmott forum member zukimaten did additional checks. He confirmed the accuracy of the first panel of results to machine precision (say 15 digits). For the more difficult second panel, he confirmed the call values with strikes (90,95,100) to 15 leading digits, and the strikes (105, 110) to 13 leading digits. Of course, by put-call parity, similar results can be inferred for the corresponding puts.