Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
36,364 characters · 7 sections · 30 citation commands
A Bayesian take on option pricing with Gaussian processes
The local volatility model is widely used in finance within option pricing. One reason is its capability of reproducing any given surface of observed prices from European options. The flexibility is due to a latent function---the local volatility function---which uniquely defines the model's pricing operator. The crux of the matter is that it can only be unambiguously retrieved in the limit of infinite data. Numerous calibration approaches have therefore been suggested for its (re)construction. Deterministic approaches prevail, typically parametric methods that involve steps of interpolation and extrapolation.
We look at calibration from a probabilistic angle with an approach based on Gaussian processes. This gives a way of encoding prior believes about local volatility with a nonparametric model which is flexible, yet prone to overfitting. In addition to providing point-estimate(s), we draw posterior inference about local volatility to better understand uncertainty attached with the calibrates. Our approach also gives means for probabilist prediction of local volatility as well as for deriving quantities from the model which takes the posterior information into account.
In mathematical finance, option pricing is often based on It\^{o}-process models of the financial market. We consider the local volatility model derman1994,dupire1994pricing where the risk-neutral dynamics of a stock price $S$ is
(see oksendal2003stochastic for an introduction to stochastic differential equations). The latent {local volatility function} $\sigma:\mathbb{R}_+^2 \rightarrow\mathbb{R}_+$ is the essential component of (ref) and our main target for modelling. $S$ is driven by a Wiener process $W$ and the constant $r$ is the interest rate of a bank account with value $dB_t = rB_tdt$. Together, $(S,B)$ constitute the underlying financial market on which options are derivative.
A European call option ({call} for short) with strike price $K$ and maturity time $T$ gives its holder a payoff $\max(S_T-K,0)$ at maturity, where the terminal stock price $S_T$ is unknown at $t<T$. Thus, one need a method for computing the time-$t$ fair price $C(t,S_t,T,K)$ of the option. No-arbitrage pricing theory (e.g., bjork2009arbitrage) yields that
where $E[\cdot|S_t]$ denotes conditional expectation given the current stock price $S_t$. Dupire dupire1994pricing shows how to represent (ref) as the solution to a partial differential equation in the variables $(T,K)$
with boundary condition $C(T=0,K) = \max(S_t-K,0)$.
In practice, an option-pricing model like (ref) may be used for generating prices of options not quoted on the market, deriving hedging strategies and evaluating risk measures. Before employing a model for such purposes, it has to be fitted to data with a {calibration} approach modelCalibration. Here, data is market prices of liquid vanilla options--- European puts and calls---due to their role as benchmark instruments and since they are frequently (and publicly) traded. Figure (ref) shows the data we use for this paper.
A theoretical appeal with the local volatility model is that it can perfectly fit any given set of consistent call (and/or put) prices. This means that there are no internal arbitrage opportunities in the set carr2005note. In practice, however, market prices are quoted only up to their bid-ask spread (left Figure (ref)) from which it is common to use mid-prices for calibration. This is often formulated as a least-squares problem
where model prices $C_i(\sigma)$ correspond to mid-market prices $\hat{c}_i$. When optimising (ref) one needs to evaluate $C_i(\sigma)$, $\forall i$, for each candidate local volatility function. This can only be done with numerical methods. Dupire's equation speeds up this computation---thereby the overall calibration procedure---since a numerical solution of (ref) gives $C_{i}$ for all $i$ in a single sweep. Following gupta2014robust,cont2004recovering,jackson1998computation among others, we use a Crank–Nicolson scheme hirsa2012computational for this purpose, i.e., for the mapping $\sigma\mapsto \{C_i(\sigma)\}_{i}$.
We place local volatility in a Bayesian framework by a probabilistic generative model of data. To this end, we take the space of strikes and maturities, $x:=(T,K)\in\mathcal{X} := \mathbb{R}^2_+$, to be the input space of a Gaussian process prior. The Gaussian process defines a distribution over functions from $\mathcal{X}$ to $\mathbb{R}$. It is uniquely specified by a covariance function $k_{\kappa}: \mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ parametrised by $\kappa$, and a mean function which we take to be zero rasmussen2006gaussian. A distribution over local volatility functions is then defined through a positive transformation
The prior model (ref) provides means of encoding prior knowledge about local volatility. The stock price diffusion (ref) and Dupire's equation (ref) pertain such knowledge (to guarantee a solution of (ref) with continuous transition density) but it is only minimal: $\sigma:\mathcal{X}\rightarrow\mathbb{R}$ being continuous and positive suffices. Stronger regularity assumptions are often motivated from a practical viewpoint. A smooth volatility surface is desired when using the model for hedging and risk management, and to obtain a more robust calibration method. This is conveniently encoded by the squared exponential covariance
since (ref)-(ref) has sample paths with continuous derivatives of all orders lindgren2012stationary.
While the Gaussian process has the title role in our approach, we accommodate it with a hyperprior distribution over $\kappa=(\sigma_f,l_T,l_K)$. These parameters specify the general shape of the functional output: $l_T,l_K$ control variability along input dimensions and $\sigma_f$ controls output variance. Specifying sensible values for $\sigma_f,l_T,l_K$ is less straightforward and hyperparameter inference is therefore integral in the learning process. We assume independent sigmoid Gaussians
where $\phi$ is the standard normal density. A limited support of the parameter space make the model prone to identifiability issues and it helps improving sampling efficiency aki. We will also leverage that (ref) defines a transformed Gaussian in the design of our sampling strategy.
For the final component of the generative procedure we take the least-squares calibration (ref) as basis for a likelihood. Market data $\mathcal{D}=\{\hat{c}_i,\hat{x}_i\}_{i=1}^{n}$ (the \{mid-price, maturity, strike\} of quoted options) is assumed to be fair prices conditional on local volatility, and observed with additive Gaussian noise. In the simple case of independent noise with common standard deviation\footnote{We use boldface to denote a vector, or more generally, a set $\textbf{x}=\{x_i\}_{i}$.}
In the likelihood (ref) we explicitly write out the dependency of the model price $C$ on a strike-maturity point $\hat{x}_i$ (corresponding to $\hat{c}_i$) and a finite set of $N\geq n$ volatility/functional values, $\textbf{f}=\log\boldsymbol{\sigma}$, taken at an input set $\textbf{x}$. For the finite difference scheme, the input set is a Cartesian product $\textbf{x} = \{(T_i,K_j)\}_{i,j}$ and it is constructed to include all market inputs of $\mathcal{D}$, i.e., $\hat{\textbf{x}}\subset\textbf{x}$. The reader can visualise a lattice in the $x$-$y$ plane of Figure (ref) with nodes $\textbf{x}$ covering all red dots $\hat{\textbf{x}}$.
We are now in a position to stage the play and apply Bayes' theorem for a posterior distribution over all unknowns. First we make a simple augmentation of the likelihood: $p(\hat{\textbf{c}}|\textbf{f},\mu_f,\sigma_{\epsilon})\equiv p(\hat{\textbf{c}}|\textbf{f}+\mu_f,\sigma_{\epsilon})$ with $\mu_f$ allowing for a constant (non-zero) mean of the Gaussian process. The joint posterior over hyperparameters and function values is then
where $Z$ is a normalising constant. The prior over $\textbf{f}$ is induced by the Gaussian process, $p(\textbf{f}|\kappa) = \mathcal{N}(\textbf{f}|\textbf{0},\boldsymbol{\Sigma}_{\kappa})$ with a $N\times N$ covariance matrix with elements $[\boldsymbol{\Sigma}_{\kappa}]_{i,j}=k_{\kappa}(x_i,x_j)$ for all input pairs in $\textbf{x}$. For the hyperprior over $(\mu_f,\sigma_{\epsilon})$ we also assume independent sigmoid Gaussians.
Calibration in our probabilistic model amounts to posterior inference about $\textbf{f}$. Since the posterior distribution is analytically intractable, we perform approximate inference by generating samples from (ref) (Section (ref)). From such a sample, we can extract the mean, mode or {maximum a posteriori} (MAP) surface for stand-in replacements of “classical” deterministic point-estimates. What more is, we obtain a fully probabilistic representation of local volatility. This is beneficial for several reasons: (i) posterior credible intervals provide a probabilistic notion of uncertainty for our estimates, (ii) the posterior uncertainty may be consistently incorporated in derivative quantities of the model, $g(\sigma)$, through Bayes' estimator $\int g(\textbf{f})p(\textbf{f}|\mathcal{D})d\textbf{f}$ and, (iii) it provides consistent probabilistic prediction of local volatility at unseen inputs ${x}_{\star}\in\mathcal{X}$ by
The conditional prior in (ref) follows from the predictive equations of Gaussian processes, i.e., $\mathcal{N}(f({x}_{\star})|\textbf{k}_{\kappa}^{\top}\boldsymbol{\Sigma}_{\kappa}^{-1}\textbf{f},k_{\kappa}({x}_{\star},{x}_{\star})-\textbf{k}_{\kappa}^{\top}\boldsymbol{\Sigma}_{\kappa}^{-1}\textbf{k}_{\kappa})$ with $\textbf{k}_{\kappa} = \{k_{\kappa}({x}_{\star},x)\}_{x\in\textbf{x}}$ a column vector rasmussen2006gaussian. The likelihood ratio does not cancel out since (ref) is non-factorisable over the local volatility surface: a model price $C( \{\textbf{f},{f}_{\star}\},\hat{x})$ depends on all elements of $\{\textbf{f},{f}_{\star}\}$ even if the dependency is strongest on values taken {locally} around $\hat{x}$. To first order, we may approximate it to be constant, in particular when we predict at a point ${x}_{\star}$ lying outside of $\hat{\textbf{x}}$. It is then possible to generate predictions with a Gaussian mixture approximating (ref), or by direct simulation.
Models with latent variables like (ref) are intractable in general. We use Markov chain Monte Carlo to represent the posterior with samples. These are asymptotically exact methods in that they provide samples from a distribution which is arbitrarily close to the target distribution in the long run robert2004monte.
We generate states of $(\textbf{f},\kappa,\mu_f,\sigma_{\epsilon})$ by Gibbs sampling of three blocks. First, we update $\textbf{f}$ by targeting its conditional posterior
We use elliptical slice sampling for this purpose murray2010elliptical. Second, we update $\kappa$ with surrogate data slice sampling murray2010slice. This means targeting the conditional
where $\textbf{g}$ augments the model with a noisy version of $\textbf{f}$, and $\textbf{f}(\kappa,\boldsymbol{\eta},\textbf{g})$ is a representation of $\textbf{f}|\textbf{g}$ with a standard normal variable $\boldsymbol{\eta}\sim\mathcal{N}(\textbf{0},\mathbb{I})$---see murray2010slice for details. Since (ref) defines a transformed Gaussian, $\kappa_i = \kappa_{i,\text{max}}(1+\exp(-\textbf{z}_i))^{-1}$ with $\textbf{z}\sim\mathcal{N}(\textbf{0},\mathbb{I})$, we target (ref) with elliptical slice sampling in $\textbf{z}$-space. Third, we target the conditional of $(\mu_f,\sigma_{\epsilon})$
and again we use elliptical slice sampling for updating the corresponding bivariate Gaussian.
In terms of computational complexity, the second update of covariance parameters is the most expensive with a cost $\mathcal{O}(N^3)$ for (standard) inversion of the covariance matrix each time a new $\kappa$ is proposed. Proposing a new $\textbf{f}$ also costs $\mathcal{O}(N^3)$ for the Cholesky decomposition of $\boldsymbol{\Sigma}_{\kappa}$, although the same decomposition can be recycled for consecutive proposals. The third update is the cheapest since it only requires computing the likelihood. One such evaluation involves solving $I$ systems of $J\times J$ tridiagonal matrices for the Crank-Nicolson scheme of the call price, where $I$ ($J$) is the number of maturities (strikes) of the input set $\textbf{x}$ i.e., a cost $\mathcal{O}(IJ)=\mathcal{O}(N)$. In fact, we can readily exploit that $\textbf{x}$ is a Cartesian product and use Kronecker methods when computing matrix decompositions and inversions saatcci2012scalable. This reduces the cost to $\mathcal{O}(N^{3/2})$ for both the second (for each $\kappa$ proposal) and the first update (recycled throughout a $\textbf{f}$-update with elliptical slice sampling).
To the best of our knowledge, there is no previous work that treats local volatility calibration with a full Bayesian approach based on functional priors. In fact, we try to achieve what Hamida and Cont cont2004recovering deemed too complicated just over a decade ago “The computational complexity of the Bayesian approach has prevented it from being applied beyond a simple case such as Black-Scholes.” While a Bayesian analysis of Black-Scholes model is given in jacquier2000bayesian, the authors propose an evolutionary algorithm to generate a population of local volatility. All surfaces which, on average, reproduces benchmark prices within an average bid-ask spread are treated as equal calibrates, and used to span price-intervals of derivatives. They use cubic splines to represent the surfaces, a parametric model of local volatility surfaces adopted before by jackson1998computation.
The spline model is also used by gupta2014robust in a Bayesian approach to estimate its parameters. Prior beliefs are here encoded by penalising local volatility deviations from the earliest at-the-money implied volatility with respect to the Sobolev norm of first order. Hence, only first order differentiability is considered as a smoothness requirement. By defining the negative log-prior proportional to this distance they obtain a Gaussian over the parameters of the spline model. They take the least-squares error to define a Gaussian likelihood with (noise) standard deviation fixed at the average bid-ask spread. Similar to cont2004recovering they assign zero likelihood to surfaces which do not respect the bid-ask spread criterion.
On a conceptual level, our approach is similar to gupta2014robust. The major difference is that we use a nonparametric functional prior and perform full Bayesian analysis. While the Gaussian process defines a random function with properties specified by the covariance kernel, cubic splines are fundamentally deterministic interpolation polynomials, even if a distribution over the nodes induces a distribution over functions. Interpretability of the models is therefore different: encoding prior knowledge in a Gaussian process is straightforward by choosing an appropriate covariance (this may even be part of the learning process, see e.g. wilson2015kernel). In the spline model, the choice of norm and attached Gaussian is ad hoc such that prior beliefs about the polynomial function are somewhat indirect. Moreover, the interpretation of hyperparameters is clear in the Gaussian process model. We include these in the inference process---and stress their importance---for a full Bayesian analysis. In contrast, pre-specified constants control the prior behaviour of the spline model, with an ex post robustness analysis of their impact on model outputs---not on the posterior over local volatility. Similarly, the impact of the fixed likelihood variance is analysed in terms of output robustness with respect to the data-noise (not the variance parameter itself). Here we emphasise that noise variance do indeed affect the posterior over local volatility to a great extent---so do all hyperparameters. If we fix a corresponding $\sigma_{\epsilon}$ from bid-ask spreads and exclude it from inference, we obtain much wider confidence regions compared to Figure (ref).
Finally we just mention that there is a large body of literature on calibration through regularised problems of the form (ref), e.g., chiarella2000calibration, crepey2003calibration,egger2005tikhonov, jackson1998computation. Albeit deterministic, these are related to our approach: unconstrained minimisation of (ref) corresponds to maximising (ref) (for fixed $\sigma_{\epsilon}$), i.e., the maximum likelihood estimator---or a completely uninformative prior. Similarly, minimising (ref) with a quadratic regularisation term corresponds to maximising (ref) (for fixed hyperparameters), i.e., the MAP estimator.
We demonstrate our approach with market data from call options on the S&P 500 index, shown in Figure (ref). The left figure shows the (best) bid and ask price of 473 options quoted on 4 January 2010. The mid-market prices are converted to implied Black-Scholes volatilities, shown in the right figure. {These two representations are equivalent as the Black-Scholes price is monotonically increasing in implied volatility.} The latter expresses option prices in a standardised “currency”, often preferred by option traders. It can be though of as representing a market sentiment in that it is the market's expectation of future fluctuations in the underlying stock-price. A high level of implied volatility attached to a strike-maturity yields an expensive option. It anticipates a volatile stock-price (high risk), such that trading in options to protect from (adverse) price movements will be expensive.
\paragraph{Calibration} From the market data, we take a subset $\mathcal{D}$ of $n=117$ options for calibration with strike-maturities $\hat{\textbf{x}}$ distributed over 20 strikes and 8 maturities. These are marked as blue strike-maturity points in right Figure (ref). The input set $\textbf{x}$ thus forms a $20\times8$ grid covering these points. This is to have a reasonably sized input of $N=160$ points for computational efficiency in the MCMC algorithm (compared with $103\times14=1442$ for the full data). We still cover the original ranges of strikes and maturities; we just take a less dense input set evenly distributed over the original quotes.
To set-up an operational model, we scale $\textbf{x}$ to lie in the unit square and set $\kappa_{i,\max}=1$, $\forall i$, for the covariance hyperprior, and $\sigma_{\epsilon,\max}=0.5$ and $ \mu_{f,\max} = \log 0.5$ for the likelihood . These values stem from heuristics and some initial runs of the MCMC algorithm. Length-scales around one (maximum) give slowly varying volatility surfaces much smoother than the implied volatility in right Figure (ref): prior $\kappa$-ranges $(0,1)$ provide enough flexibility. For the likelihood, the maximum noise standard-deviation 0.5 is in the same range as the average bid-ask spread of 2.5. Note, however, that the noise models the deviance of the mid-market price from the (unobserved) fair price, not the (half) bid-ask spread.
We run 50,000 iterations of the three-step algorithm described in Section (ref). We discard the first 10,000 states as burn-in, and subsample the reminder to obtain a thinned posterior sample of 1000 local volatility surfaces. The result is represented in Figure (ref), where we show four cross-sections of the volatility surfaces taken along the strike dimension, for maturities of 0.13, 0.48, 1.5 and 3 years. In each panel, we show a slice from the MAP surface and a credible interval of $\pm2$ standard deviations around the mean (estimated from the posterior sample point-wise at each $x$-input).
The shapes of the credible intervals are symptomatic of two main factors: non-linearity of the pricing/likelihood function and data availability. For early maturities (0.13 and 0.48 years), the price is not very sensitive to model parameters when looking at strikes far from the current asset price $S_t = \$1130$. For $K>>S_t$, there is not much chance of the option reaching at/in-the-money before maturity, i.e., $K\leq S_T$. This keeps the price close to zero (since the option is likely to expire worthless; see left Figure (ref)). Similarly, for $K<<S_t$ there is little chance of going out-of-money and the option's price is close to its intrinsic value, $S_t-K$ (its time value is small). In effect, the price insensitivity makes the likelihood relatively uninformative about local volatility over these points. For the same reasons there is little trading in the options. Market quotes are concentrated around at-the money, $K=S_t$, which also contribute to increased uncertainty outside this region as there is lack of data. The effects can be seen as a high degree of uncertainty for large and small strikes in the top panes of Figure (ref).
On the other hand, posterior confidence is high at strikes close to at-the-money where the price is most sensitive to local volatility, and where there is a high concentration of quotes. For late maturities (3 years) the price is also sensitive to model parameters for small and large strikes. With time, there is an increasing probability that the option will leave its in- or out-of-money region. Considering the dynamics (ref), local volatility is a parameter of diffuseness of the stock-price. For long maturities, the volatility level therefore has higher impact on option prices over all strikes. In effect, the credible region of bottom right Figure (ref) shows a rather constant uncertainty across strikes. Compared to earlier maturites, the band is wider around at-the-money. Another reason for this is data availability with a larger gap between maturities of 2 and 3 years (see Figure (ref)).
The posterior sample of local volatility surfaces gives a probabilistic picture of the local volatility function. This is interesting in its own right as it provides information on what can be learned (and not learned) from data about the latent function, i.e., to what degree of certainty one can pin down $\sigma$-values over different strike-maturity points. The result---point-estimate(s) and associated uncertainty---can then be used when employing the model for other purposes. For example, if pricing a derivative which is sensitive to $\sigma$ over points with high calibration uncertainty, one better takes this uncertainty into account for a robust price-estimate.
\paragraph{Re-pricing market data} To see how well the calibrated model explains market data we transform the posterior sample over local volatility to samples over call prices and implied volatilities. The result is shown in Figure (ref) for implied volatilities; corresponding call prices (not shown for brevity) show confidence regions which are completely tight (prices are \$0--500 with standard deviations $\sim$ \$0.25), with all market prices indistinguishable from the MAP surface. In terms of implied volatilities, Figure (ref) show that the MAP surface is close to the data used for calibration (blue dots) and to the market data excluded from this set (red asterisks). The only exemptions are some implied volatilities at the shortest maturity, as shown in the top left pane of Figure (ref). This confirms the fact that the local volatility model is often inconsistent with market prices at short maturities gatheral2011volatility. An explanation is that the underlying stock price is modelled as a continuous process, while the market anticipate jump-like behaviour which contribute to more pronounced implied-volatility skews (as observed in top-left Figure (ref)).
Implied volatility MAP-to-market errors are within $0.0003\pm0.02$ and within $0\pm0.01$ if excluding the shortest maturity; both indicate a strong re-pricing performance of our approach (recall that the likelihood is based on prices, not implied volatilities) and competitive to what is reported in the literature: cont2004recovering report $\sim0\pm0.04$ for a set of DAX options, luo2010local a few percentage points.
Finally we note that the uncertainty of the local volatility estimate transfers to modest credible intervals for implied volatility, especially over the in-the-money region. The reason is, again, model price/implied volatility (in)sensitivity to the local volatility parameter---even if there is a high level of local-volatility uncertainty over some strike-maturity regions, it is not propagated back to uncertainty over prices/implied volatilities. This property of call options should, however, not be taken for granted for other derivatives, especially if they are designed to be heavily dependent on volatility.
\paragraph{Prediction} As mentioned, it is common to employ a calibrated model to price off-exchange products, i.e., possibly non-standard derivatives not (publicly) traded. These may dependent on unquoted strike-maturity points which require interpolation and extrapolation of the calibrated local volatility surface. We readily have a machinery for this purpose that provides a predictive posterior over local volatility at any strike-maturity point of our choice. As an illustrative example we consider long dated options with maturities 4--10 years. We generate a prediction sample of 1000 local volatility surfaces by simulation: each of 100 randomly selected states $(\textbf{f},\kappa,\mu_f)$ from the posterior sample are used to calculate the moments of the conditional prior in (ref). We then generate 10 surfaces from every such set of moments. The result is shown in left Figure (ref) as a point-wise credible envelope of $\pm2$ standard deviation around the mean. Comparing with the calibration (Figure (ref)) the prediction uncertainty is substantial over small/large strikes and increasing in maturity up to a point $\sim7$ years, from which the surfaces level off to their prior confidence $\sim 0.2\pm0.3$. The sample is used to a predict call prices (middle Figure (ref)) which in turn are transformed to a sample over implied volatilities (right Figure (ref)). Both illustrates how prediction uncertainty of local volatility propagate: long-dated options can only be consistently priced subject to an increasing uncertainty-range which is substantial for the latest maturities.
We have presented an approach for Bayesian inference in the local volatility model based on Gaussian process priors. To the best of our knowledge, this is the first full Bayesian analysis with a nonparametric functional prior. We discussed several benefits of our approach such as straightforward encoding of prior beliefs and interpretability, the gain of a probabilistic notion of calibration uncertainty, posterior predictions consistent with data, and the evaluation of derivate model-quantities with posterior information taken into account. We also proposed an MCMC algorithm for inference. For real market data, we demonstrated encouraging calibration performance and discussed how the pricing model and benchmark data interact to generate the resulting posterior, with focus on uncertainty. Finally, we predicted local volatility over unseen, long maturities and used the result to generate a predictive distribution for long-dated call options.
As immediate extensions of this work, one might consider alternations of the setup and model. A prior Matérn covariance would probably achieve higher likelihood values (smaller re-pricing errors) to the cost of less smooth volatility surfaces. In terms of the likelihood, the noise model (ref) is the simplest choice. Market features such as liquidity and transaction costs could motivate more intricate models with structural and also non-Gaussian noise. Such could readily be adopted in the likelihood and their parameters added to the inference process. We also note the freedom to choose other calibration benchmarks---bid and ask prices for example---as long they permit a pricing operator which can be evaluated (Monte Carlo is a go-to method, although involving SDE simulation in the the likelihood would be terrible in terms of computational effort). Several extensions exist for the local volatility model itself---such as stochastic interest rates and discrete dividends--- which could be investigated. A more fundamental and interesting extension is to consider a Bayesian Gaussian-process approach to the Brunick-Shreve model brunick2013mimicking. The latent diffusion function of this model is extended to depend on a third state, the running maximum of $S$. Since the payoff of barrier options depends on the running maximum, barriers take a similar role for this model as calls have for local volatility. With a pricing equation corresponding to (ref) derived in hambly2016forward, we could thus mimic our prevailing approach with an extended input space of the Gaussian process. We leave this for future work.