Heston Model Reference Prices

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 \mbox{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:

StrikeType (P=Put, C=Call)Option Price
80P7.958878113256768285213263077598987193482161301733
80C26.774758743998854221382195325726949201687074848341
90P12.017966707346304987709573290236471654992071308187
90C20.933349000596710388139445766564068085476194042256
100P17.055270961270109413522653999411000974895436309183
100C16.070154917028834278213466703938231827658768230714
110P23.017825898442800538908781834822560777763225722188
110C12.132211516709844867860534767549426052805766831181
120P29.811026202682471843340682293165857439167301370697
120C9.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:

StrikeType (P=Put, C=Call)Option Price
90P4.5183603586861772614990106188215872180542*10^-8
90C9.989001595065276544935948045293485530832966049263
95P0.000461954855653851579672612557018857858641926937
95C4.989963479738160122154264702582719627807098780529
100P0.477781171629504680023239655436072890669645669297
100C0.467782671512844263098248405184095087949465507760
105P5.009501052563650299130635110520904481889436667608
105C2.527447823194706060519991248106500619490942*10^-6
110P10.008998550115123724684210555728039829315964456261
110C1.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.

Introducing XGBM — a new solvable stochastic volatility model having a stationary distribution

I have posted a new research paper to the arXiv titled “Exact Solutions for a GBM-type Stochastic Volatility Model having a Stationary Distribution”. The article may be found here.

Let me explain the title and motivation. The most popular stochastic volatility model is certainly the Heston ’93 model. Its virtues are that it is exactly solvable and it has some very nice features: a mean-reverting, stochastic volatility process  that can be negatively correlated with equity returns. This process has the well-known “square-root model” form

$\displaystyle \mbox{Heston:} \quad d v_t = (\omega - \theta v_t) \, dt + \xi \,\sqrt{ v_t} \, dW_t,$

where  $v_t$ is the instantaneous variance rate of the asset in question. We could also write this evolution in terms of the instantaneous volatility $\sigma_t = \sqrt{v_t}$. However, time series analysis of the volatility of broad-based indices, such as the S&P500 index, suggests that a better specification would have $d \sigma_t \sim \xi \sigma_t dW_t$, which is a geometric Brownian motion (GBM)-type volatility. Taking $d \sigma_t = \xi \sigma_t dW_t$ is the well-known SABR model. The problem with that one is that the volatility is not mean-reverting and does not admit a stationary density $\psi(\sigma)$. To get a stationary density, you need a mean-reverting drift term. Also, it would be very nice to also have a model that has exact solutions. The XGBM model may be the first to combine all these properties. It is a standard bivariate model for the pair $(S_t,\sigma_t)$, but here I’m just going to write the volatility evolution, which is

$\displaystyle \mbox{XGBM:} \quad d \sigma_t = \sigma_t (\omega - \theta \sigma_t) \, dt + \xi \, \sigma_t \, dW_t.$

To see how this model is closer to the “real-world” than the Heston ’93 model, take a look at the figures at the end, which show Maximum Likelihood fits to a proxy series for $\{\sigma_t\}$. The proxy is the daily (annualized) volatility for the S&P500 taken from the Oxford-Man Institutes “realized library”. They maintain a number of estimators — I am using a basic one (rv5) which is simply the daily (intraday) volatility using 5-minute log-return observations. I am using all the data available at the time of this study, which is January 3, 2000 – September 28, 2018 (4705 volatility observations). The first figure shows this volatility time series. You can see there is a maximum of around 140% which is the annualized volatility at the height of the 2007-2008 Financial Crisis.

The next figures show the stationary volatility fits. For these histograms, I am using the annualized volatility in decimal (not percent), so the time series maximum is a histogram entry at $\sigma \approx 1.4.$ That small bump is not really visible, but you can see it is accounted for by the axis extension out to that value.

As one can see, the visual fit is better for XGBM vs. Heston, with corresponding log-likelihoods:  5356 (Heston) vs 5920 (XGBM). In fact, the fit is even better (LL = 6055) for the figure labeled “GARCH -3/2 model”, which is the stationary density corresponding to both the GARCH diffusion model and the 3/2-model. In any event, the point of the exercise is to motivate my new paper and interested readers may find it at the link above.

Update. This research paper is now published. The cite is “Lewis, Alan L. (2019), “Exact solutions for a GBM-type stochastic volatility model having a stationary distribution”, Wilmott mag, Vol. 2019, Issue 101, May, 20-41. The editors were kind enough to make it a
cover story. As part of that publication, I wrote an additional Introduction, which is available here: XGBMIntro.WMag.Final.

Calibration of the GARCH Diffusion Model

The GARCH diffusion model  is one of the running examples of bivariate stochastic volatility models in my first book. (Others include the well-known Heston ’93 model and the so-called 3/2-model). As with most finance models, it comes in two versions: a real-world (aka P-measure) version and a risk-neutral (aka Q-measure) version. The latter is used to value options.

For equity applications, the stochastic process pair is $\{S_t,v_t\}$. Here $S_t \ge 0$ is a stock price and $v_t = \sigma_t^2 > 0$ is the associated stochastic variance rate. Then, the risk-neutral version, as an SDE (stochastic differential equation) system, reads

$\displaystyle \left \{ \begin{array}{l}dS_t = (r-q) S_t dt + \sigma_t S_t dB_t, \\d v_t = \kappa (\bar{v} - v_t) + \xi \, v_t \, dW_t, \quad dB_t \, dW_t = \rho \, dt. \end{array}\right. \quad\quad\quad\quad\quad (1)$

It’s a very simple and natural model which, unfortunately, lacksÂ  exact solutions for option values or transition densities. To determine the unknown parameters in (1), one needs to calculate option prices (numerically) and fit the model to option chains, a procedure generally called “calibration”. To do this efficiently and accurately for this model (and many others) requires a PDE approach. (There have been some earlier calibrations of this model via Monte Carlo).

The lack of an efficient — or, apparently, any — PDE calibration for this model prompted Yiannis Papadopoulos and me to perform one. Our methods and first results were recently posted on the arXiv:

A First Option Calibration of the GARCH Diffusion Model by a PDE Method

Yiannis’ own announcement may be found in his blog here. You will also find at that link a free downloadable GARCH diffusion calibrator demo that Yiannis has developed. (Windows executable). You can run it on a sample option chain that he supplies, and see a calibration in well under a minute (11 seconds on my desktop). Or you can run it on your own data, simply by imitating the provided data file format.