Algorithmic market making for options

TL;DR

This paper develops a tractable approach for optimal market making in options by reducing a high-dimensional control problem to a low-dimensional one using vega approximation, enabling efficient numerical solutions.

Contribution

It introduces a novel method to simplify the complex stochastic control problem of options market making into a low-dimensional equation using vega approximation, applicable to models like Heston.

Findings

The method simplifies high-dimensional control problems to low-dimensional equations.

Numerical solutions are feasible using Euler schemes and interpolation.

The approach is demonstrated with numerical examples.

Abstract

In this article, we tackle the problem of a market maker in charge of a book of options on a single liquid underlying asset. By using an approximation of the portfolio in terms of its vega, we show that the seemingly high-dimensional stochastic optimal control problem of an option market maker is in fact tractable. More precisely, when volatility is modeled using a classical stochastic volatility model -- e.g. the Heston model -- the problem faced by an option market maker is characterized by a low-dimensional functional equation that can be solved numerically using a Euler scheme along with interpolation techniques, even for large portfolios. In order to illustrate our findings, numerical examples are provided.

Figures

Equations

$$\displaystyle\left\{\begin{array}[]{ll}dS_{t}=\mu S_{t}dt+\sqrt{\nu_{t}}S_{t}dW_{t}^{S}\\ d\nu_{t}=a_{\mathbb{P}}(t,\nu_{t})dt+\xi\sqrt{\nu_{t}}dW_{t}^{\nu},\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{ll}dS_{t}=\sqrt{\nu_{t}}S_{t}d\widehat{W}_{t}^{S}\\ d\nu_{t}=a_{\mathbb{Q}}(t,\nu_{t})dt+\xi\sqrt{\nu_{t}}d\widehat{W}_{t}^{\nu},\end{array}\right.$$

$$0$$

$$\displaystyle dq_{t}^{i}:=\int_{\mathbb{R}_{+}^{*}}z\big{(}N^{i,b}(dt,dz)-N^{i,a}(dt,dz)\big{)},\forall i\in\{1,\ldots,N\},$$

$${\lambda }_t^{i,b}(dz):={\Lambda }^{i,b}({\delta }_t^{i,b}(z))\text{mathds}{1}_{\{q_{t-}+z{e}^i\in \mathcal{Q}\}}{\mu }^{i,b}(dz){\lambda }_t^{i,a}(dz):={\Lambda }^{i,a}({\delta }_t^{i,a}(z))\text{mathds}{1}_{\{q_{t-}-z{e}^i\in \mathcal{Q}\}}{\mu }^{i,a}(dz)$$

$${\Delta }_t:=\underset{i=1}{\sum ^N}{\partial }_S{O}^i(t,S_t,{\nu }_t)q_t^i\text{ for all }t\in [0,T].$$

$$\displaystyle dX_{t}:=\underset{i=1}{\overset{N}{\sum}}\bigg{(}\int_{\mathbb{R}_{+}^{*}}z\Big{(}\delta_{t}^{i,b}(z)N^{i,b}(dt,dz)+\delta_{t}^{i,a}(z)N^{i,a}(dt,dz)\Big{)}-\mathcal{O}_{t}^{i}dq_{t}^{i}\bigg{)}+S_{t}d\Delta_{t}+d\big{\langle}\Delta,S\big{\rangle}_{t}.$$

$$V_t:=X_t-{\Delta }_t{S}_t+\underset{i=1}{\sum ^N}{q}_t^i{\mathcal{O}}_t^i.$$

$$d{V}_t$$

$${\mathcal{V}}_t^i:={\partial }_{\sqrt{\nu }}{O}^i(t,S_t,{\nu }_t)=2\sqrt{{\nu }_t}{\partial }_{\nu }{O}^i(t,S_t,{\nu }_t)\text{ for all }t\in [0,T].$$

$$\displaystyle dV_{t}=\underset{i=1}{\overset{N}{\sum}}\left(\int_{\mathbb{R}_{+}^{*}}z\Big{(}\delta_{t}^{i,b}(z)N_{t}^{i,b}(dt,dz)+\delta_{t}^{i,a}(z)N_{t}^{i,a}(dt,dz)\Big{)}+q_{t}^{i}\mathcal{V}_{t}^{i}\frac{a_{\mathbb{P}}(t,\nu_{t})-a_{\mathbb{Q}}(t,\nu_{t})}{2\sqrt{\nu_{t}}}dt+\frac{\xi}{2}q_{t}^{i}\mathcal{V}_{t}^{i}dW_{t}^{\nu}\right).$$

$$\sup_{\delta\in\mathcal{A}}\mathbb{E}\bigg{[}-\exp\big{(}-\gamma V_{T}\big{)}\bigg{]},$$

$$\underset{\delta \in \mathcal{A}}{\sup }\mathbb{E}\left[V_T-\frac{\gamma }{2}{\int }_0^T{\left(\underset{i=1}{\sum ^N}\frac{\xi }{2}{q}_t^i{\mathcal{V}}_t^i\right)}^{2}dt\right].$$

$$\underset{\delta \in \mathcal{A}}{\sup }\mathbb{E}$$

$$+q_t^i{\mathcal{V}}_t^i\frac{a_{\mathbb{P}}(t,{\nu }_t)-a_{\mathbb{Q}}(t,{\nu }_t)}{2\sqrt{{\nu }_t}}))dt-\frac{\gamma {\xi }^2}{8}{\int }_0^T{\left(\underset{i=1}{\sum ^N}{q}_t^i{\mathcal{V}}_t^i\right)}^{2}dt],$$

$$\displaystyle\psi(j):=\left\{\begin{array}[]{ll}+1\text{ if }j=a\\ -1\text{ if }j=b.\end{array}\right.$$

$$u:\big{(}t,S,\nu,q\big{)}\in[0,T]\times\mathbb{R}_{+}^{2}\times\mathcal{Q}\mapsto u\big{(}t,S,\nu,q\big{)}$$

$$\displaystyle u\big{(}t,S,\nu,q\big{)}$$

$${\mathcal{V}}_t^i={\mathcal{V}}_0^i=:{\mathcal{V}}^i\in \mathbb{R},\text{ for all }i\in \{1,\dots ,N\}.$$

$$\mathcal{Q}=\bigg{\{}q\in\mathbb{R}^{N}\left|\underset{i=1}{\overset{N}{\sum}}q^{i}\mathcal{V}^{i}\in\left[-\overline{\mathcal{V}},\overline{\mathcal{V}}\right]\right.\bigg{\}},$$

$$d\mathcal{V}^{\pi}_{t}=\sum_{i=1}^{N}\int_{\mathbb{R}_{+}^{*}}z\mathcal{V}^{i}\big{(}N^{i,b}(dt,dz)-N^{i,a}(dt,dz)\big{)}.$$

$$\forall(t,S,\nu,q\big{)}\in[0,T]\times\mathbb{R}_{+}^{2}\times\mathcal{Q},u(t,S,\nu,q\big{)}=v\left(t,\nu,\sum_{i=1}^{N}q^{i}\mathcal{V}^{i}\right),$$

$$\begin{array}{rl}v\left(t,\nu ,{\mathcal{V}}^{\pi }\right)=\underset{({\delta }_s{)}_{s\in [t,T]}\in {\mathcal{A}}_t}{\sup }{\mathbb{E}}_{(t,\nu ,{\mathcal{V}}^{\pi })}[{\int }_t^T& (\left(\underset{i=1}{\sum ^N}\sum _{j=a,b}{\int }_{{\mathbb{R}}_{+}^{\ast }}z{\delta }_s^{i,j}(z){\Lambda }^{i,j}({\delta }_s^{i,j}(z))\text{mathds}{1}_{\left\{|{\mathcal{V}}_s^{\pi }-\psi (j)z{\mathcal{V}}^i|\le \overline{\mathcal{V}}\right\}}{\mu }^{i,j}(dz)\right)\\ & +{\mathcal{V}}_s^{\pi }\frac{a_{\mathbb{P}}(s,{\nu }_s)-a_{\mathbb{Q}}(s,{\nu }_s)}{2\sqrt{{\nu }_s}}-\frac{\gamma {\xi }^2}{8}{{\mathcal{V}}_s^{\pi }}^2)ds].\end{array}$$

$$\begin{split}0=\ &\partial_{t}v(t,\nu,\mathcal{V}^{\pi})+a_{\mathbb{P}}(t,\nu)\partial_{\nu}v(t,\nu,\mathcal{V}^{\pi})+\frac{1}{2}\nu\xi^{2}\partial^{2}_{\nu\nu}v(t,\nu,\mathcal{V}^{\pi})+\mathcal{V}^{\pi}\frac{a_{\mathbb{P}}(t,\nu)-a_{\mathbb{Q}}(t,\nu)}{2\sqrt{\nu}}-\frac{\gamma\xi^{2}}{8}{\mathcal{V}^{\pi}}^{2}\\ &+\underset{i=1}{\overset{N}{\sum}}\underset{j=a,b}{\sum}\int_{\mathbb{R}_{+}^{*}}z\mathds{1}_{\left\{|\mathcal{V}^{\pi}-\psi(j)z\mathcal{V}^{i}|\leq\overline{\mathcal{V}}\right\}}H^{i,j}\bigg{(}\frac{v\big{(}t,\nu,\mathcal{V}^{\pi}\big{)}-v\big{(}t,\nu,\mathcal{V}^{\pi}-\psi(j)z\mathcal{V}^{i}\big{)}}{z}\bigg{)}\mu^{i,j}(dz),\end{split}$$

$$H^{i,j}(p):=\underset{{\delta }^{i,j}\ge {\delta }_{\infty }}{\sup }{\Lambda }^{i,j}({\delta }^{i,j})({\delta }^{i,j}-p),i\in \{1,\dots ,N\},j=a,b.$$

$$\delta^{i,j*}_{t}(z)=\max\left(\delta_{\infty},\left(\Lambda^{i,j}\right)^{-1}\left(-H^{i,j^{\prime}}\left(\frac{v\big{(}t,\nu_{t},\mathcal{V}_{t-}^{\pi}\big{)}-v\big{(}t,\nu_{t},\mathcal{V}_{t-}^{\pi}-\psi(j)z\mathcal{V}^{i}\big{)}}{z}\right)\right)\right),\text{ }i\in\{1,\ldots,N\},\text{ }j=a,b.$$

$$0=\partial_{t}w(t,\mathcal{V}^{\pi})-\frac{\gamma\xi^{2}}{8}{\mathcal{V}^{\pi}}^{2}+\underset{i=1}{\overset{N}{\sum}}\underset{j=a,b}{\sum}\int_{\mathbb{R}_{+}^{*}}z\mathds{1}_{\left\{|\mathcal{V}^{\pi}-\psi(j)z\mathcal{V}^{i}|\leq\overline{\mathcal{V}}\right\}}H^{i,j}\bigg{(}\frac{w\big{(}t,\mathcal{V}^{\pi}\big{)}-w\big{(}t,\mathcal{V}^{\pi}-\psi(j)z\mathcal{V}^{i}\big{)}}{z}\bigg{)}\mu^{i,j}(dz),$$

$$\mathcal{K}=\{8\text{€},9\text{€},10\text{€},11\text{€},12\text{€}\}\text{and}\mathcal{T}=\{1\text{ year},1.5\text{ years},2\text{ years},3\text{ years}\}.$$

$${\Lambda }^{i,j}(\delta )=\frac{{\lambda }^i}{1+e^{\alpha +\frac{\beta }{{\mathcal{V}}^i}\delta }},i\in \{1,\dots ,N\},j=a,b,$$

$$dX_{t}=\underset{i=1}{\overset{N}{\sum}}\bigg{(}\int_{\mathbb{R}_{+}^{*}}z\Big{(}\delta_{t}^{i,b}(z)N_{t}^{i,b}(dt,dz)+\delta_{t}^{i,a}(z)N_{t}^{i,a}(dt,dz)\Big{)}-\mathcal{O}_{t}^{i}dq_{t}^{i}\bigg{)}-S_{t}dq^{S}_{t}-d\big{\langle}q^{S},S\big{\rangle}_{t}.$$

$$V_t=X_t+q_t^S{S}_t+\underset{i=1}{\sum ^N}{q}_t^i{\mathcal{O}}_t^i$$

Full text

Key words: Market making, Algorithmic trading, Options, Stochastic optimal control.

1 Introduction

The electronification of financial markets started in the seventies with stock exchanges and now affects each and every asset class. For asset classes that are usually traded in a centralized way (stocks, futures, etc.), exchanges and other all-to-all trading platforms – based or not based on limit order books – are now fully automated. For assets that are still traded over the counter (OTC), electronification occurs through the introduction of new platforms, for instance single- and multi-dealer-to-client platforms.

This electronification is associated with a trend towards the automation of the trading process for many players of the financial industry: brokers, banks, but also systematic asset managers who often develop their own execution algorithms. For assets traded in centralized markets, trading automation is nowadays massive. For instance, in the cash equity world, a vast majority of the execution is now carried out using algorithms. For assets traded in dealer markets, the automation of the market making process has been at the agenda for a few years and more and more banks are developing market making algorithms for various asset classes (currencies, bonds, etc.).

In the academic literature, many market making models have been proposed since the eighties. In the early literature on market making, the two main references are the paper of Ho and Stoll [19] and the paper of Grossman and Miller [10]. Ho and Stoll introduced indeed a very relevant framework to tackle the main problem faced by market makers: inventory management. Grossman and Miller, who were more interested in capturing the essence of liquidity, proposed a very simple model with 3 periods that encompassed both market makers and final customers, enabled to understand what happens at equilibrium, and contributed to the important literature on the price formation process. If the latter paper belongs to a strand of literature that is extremely important to go beyond the simple Walrasian view of markets, it is of little help to build market making algorithms. The former paper however, after more than 25 years, has paved the way to a recent mathematical literature on algorithmic market making.

The seminal reference of the new literature on market making is the paper of Avellaneda and Stoikov [1] who revived the dynamic approach proposed by Ho and Stoll. They indeed showed how the quoting and inventory management problems of market makers could be addressed using the tools of stochastic optimal control. Since then, many models have been proposed, most of them to tackle the same problem of single-asset market making as that of Avellaneda and Stoikov. For instance, Guéant, Lehalle, and Fernandez-Tapia provided in [14] a rigourous analysis of the Avellaneda-Stoikov stochastic optimal control problem and proved that the problem could be simplified into a system of linear ordinary differential equations (ODE) in the case of exponential intensity functions. Cartea, Jaimungal, and coauthors contributed a lot to the literature and added many features to the initial models: alpha signals, ambiguity aversion, etc. (see [5, 6, 7]). They also considered a different objective function: a risk-adjusted expectation instead of a Von Neumann-Morgenstern expected utility.

The models proposed in the above papers all share the same characteristics: (i) they are agnostic with respect to the market structure but are in fact more adapted to OTC markets,555The models are also suited for assets traded in limit order books if the bid-ask spread to tick size ratio is large. (ii) they only deal with single-asset market making, and (iii) they do not deal with the market making of options.

In fact, models have been specifically developed by Guilbaud and Pham (see [16, 17]) for assets traded through limit order books (e.g. most stocks) and for assets traded on platforms with a pro-rata microstructure (e.g. some currency pairs). Interestingly, these models enable the use of aggressive orders by market makers, which is – surprisingly – a standard behavior on equity markets (see [22]).

As far as multi-asset market making is concerned, models have been developed recently to account for the correlation structure between asset price changes. Guéant extended to a multi-asset framework both models à la Avellaneda-Stoikov and models à la Cartea-Jaimungal (see [11], [12], and [13]) and showed that the problem boils down, for general intensity functions, to solving a system of (a priori nonlinear) ODEs. The associated question of the numerical methods to approximate the solution of the equations characterizing the optimal quotes of a multi-asset market maker is addressed in [2] using a factorial approach and in [15] using reinforcement learning, both with applications to corporate bond markets.

Finally, as far as asset classes are concerned, there have been few attempts to address market making problems outside of the cash world. Market making models for derivative contracts are indeed intrinsically more complicated because they must account for the strategies on both the market for the underlying asset and the market for the derivatives, and usually for numerous contracts (e.g. options for lots of strikes and maturities). Option market making is only addressed in a paper by El Aoud and Abergel (see [8]) and in a paper by Stoikov and Sağlam [23]. In the former, the authors consider a single-option market driven by a stochastic volatility model and assume that the position is always $\Delta$-hedged. They provide optimal bid and ask quotes for the option and focus on the risk of model misspecification. In the latter, the authors consider three different settings, but all with only one option: (i) a market maker in a complete market where continuous trading in the perfectly liquid underlying stock is allowed, (ii) a market maker who may not trade continuously in the underlying stock, but rather sets bid and ask quotes in the option and the stock, and (iii) a market maker in an incomplete market with residual risks due to stochastic volatility and overnight jumps.

In this paper, we consider the case of a market maker in charge of a book of options whose prices are driven by a stochastic volatility model. We assume that trading in continuous time can be carried out in the underlying asset so that the residual risk is only that of the vega associated with the inventory. Using a constant-vega approximation, we show that the problem of an option market maker boils down to solving a low-dimensional functional equation of the Hamilton-Jacobi-Bellman type that can be tackled numerically using a simple Euler scheme along with interpolation techniques. In particular, in spite of the large number of assets, the market making problem is tractable.

In Section 1, we describe the model and present the optimization problem of the option market maker. In Section 2, we show how that problem can be simplified under the constant-vega approximation. In particular, we show that solving the high-dimensional stochastic optimal control problem of the market maker boils down to solving a low-dimensional functional equation. In Section 3, we consider the example of a book of options with several strikes and maturities and provide numerical results obtained through interpolation techniques and an explicit Euler scheme.

2 Description of the problem

We consider a probability space (\Omega,\mathcal{F},\mathbb{P}\big{)} with a filtration $(\mathcal{F}_{t})_{t\in\mathbb{R}_{+}}$ satisfying the usual conditions. Throughout the paper, we assume that all stochastic processes are defined on \big{(}\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\in\mathbb{R}_{+}},\mathbb{P}\big{)}.

2.1 The market

We consider an asset whose price dynamics is described by a one-factor stochastic volatility model of the form

[TABLE]

where $\mu\in\mathbb{R},\xi\in\mathbb{R}_{+}^{*}$, $(W_{t}^{S},W_{t}^{\nu})_{t\in\mathbb{R}_{+}}$ is a couple of Brownian motions with quadratic covariation given by $\rho=\frac{d\langle W^{S},W^{\nu}\rangle}{dt}\in(-1,1)$, and $a_{\mathbb{P}}$ is such that the processes are well defined (in particular, we assume that the process $(\nu_{t})_{t\in\mathbb{R}_{+}}$ stays positive almost surely).

Remark 1**.**

*A classical example for the function $a_{\mathbb{P}}$ is that of the Heston model (see [18]), i.e. $a_{\mathbb{P}}:(t,\nu)\mapsto\kappa_{\mathbb{P}}(\theta_{\mathbb{P}}-\nu)$ where $\kappa_{\mathbb{P}},\theta_{\mathbb{P}}\in\mathbb{R}^{*}_{+}$ satisfy the Feller condition $2\kappa_{\mathbb{P}}\theta_{\mathbb{P}}>\xi^{2}$.

Remark 2**.**

*For the sake of simplicity, we consider throughout this paper a one-factor model where the instantaneous variance is the main variable of interest. Similar results could be obtained with a one-factor model focused on forward variances, such as the classical one-factor Bergomi model (see [3, 4]). Moreover, it is noteworthy that our approach can easily be extended to two-factor stochastic volatility models such as the celebrated two-factor Bergomi model (see [3, 4]), up to an increase – by $1$ – of the dimension of the equation to solve.

Assuming interest rates are equal to [math], we introduce an equivalent risk-neutral/pricing probability measure666For references, see for instance [9]. $\mathbb{Q}$ under which the price and volatility processes become

[TABLE]

where $(\widehat{W}_{t}^{S},\widehat{W}_{t}^{\nu})_{t\in\mathbb{R}_{+}}$ is another couple of Brownian motions, this time under $\mathbb{Q}$, with quadratic covariation given by $\rho=\frac{d\langle\widehat{W}^{S},\widehat{W}^{\nu}\rangle}{dt}\in(-1,1)$, and where $a_{\mathbb{Q}}$ is such that the processes are well defined.

We consider $N\geq 1$ European options written on the above asset (hereafter, the underlying asset). For each $i\in\left\{1,\dots,N\right\}$, the maturity date of the $i$-th option is denoted by $T^{i}$ and we denote by $(\mathcal{O}_{t}^{i})_{t\in[0,T^{i}]}$ the price process associated with the $i$-th option.

Remark 3**.**

*In applications, the options under consideration will always be call and/or put options. However, our setting enables to consider any European payoff.

In the above one-factor model, we know that for all $i\in\left\{1,\dots,N\right\}$, and all $t\in[0,T^{i}]$, $\mathcal{O}_{t}^{i}=O^{i}(t,S_{t},\nu_{t})$ where $O^{i}$ is solution on $[0,T^{i})\times\mathbb{R}^{2}_{+}$ of the following partial differential equation (PDE):

[TABLE]

Remark 4**.**

*Options prices are also characterized by a terminal condition corresponding to the payoff. However, we will only consider short-term optimization problems for which the time horizon is before the maturity of all the options under consideration. Therefore, we shall never use the final condition associated with Eq. (1).

2.2 The optimization problem of the market maker

We consider an option market maker in charge of providing bid and ask quotes for the $N$ above options over the period $[0,T]$ where $T<\text{min}_{i\in\{1,\ldots,N\}}T^{i}$ (see Remark 4). For all $i\in\left\{1,\dots,N\right\}$, we denote by $\mathcal{O}^{i}_{t}-\delta_{t}^{i,b}(z)$ and $\mathcal{O}^{i}_{t}+\delta_{t}^{i,a}(z)$ the bid and ask prices (per contract) proposed by the market maker for a transaction corresponding to $z$ contracts of the $i$-th option, where (\delta^{i}_{t}(.))_{t\in[0,T]}:=\big{(}\delta_{t}^{i,b}(.),\delta_{t}^{i,a}(.)\big{)}_{t\in[0,T]} is $\mathbb{F}$-predictable and bounded from below by a given constant $\delta_{\infty}$.777In applications, we always choose $\delta_{\infty}$ negative enough so that this lower bound is never binding. Hereafter, we denote by $\mathcal{A}$ the set of $\mathbb{F}$-predictable $\mathbb{R}^{2N}$-valued maps that are bounded from below by $\delta_{\infty}$, and, by abuse of language, we call these maps admissible control processes. The dynamics of the inventory process $(q_{t})_{t\in[0,T]}:=(q_{t}^{1},\dots,q_{t}^{N})^{\prime}_{t\in[0,T]}$ of the market maker is given by

[TABLE]

where, $\forall i\in\{1,\ldots,N\}$, $N^{i,b}(dt,dz)$ and $N^{i,a}(dt,dz)$ are two right-continuous $\mathbb{R}_{+}^{*}$-marked point processes, with almost surely no simultaneous jumps,888See Appendix A.3 for more details on the construction of those processes. modelling the transactions of the $i$-th option on the bid and ask side, whose respective intensity processes $(\lambda_{t}^{i,b}(dz))_{t\in\mathbb{R}_{+}}$ and $(\lambda_{t}^{i,a}(dz))_{t\in\mathbb{R}_{+}}$ are given by

[TABLE]

with $(e^{i})_{i\in\{1,\ldots,N\}}$ the canonical basis of $\mathbb{R}^{N}$, $\mathcal{Q}$ the set of authorized inventories999The frontier of this set defines the risk limits of the market maker. for the market maker, and $(\mu^{i,b},\mu^{i,a})$ a couple of probability measures on $\mathbb{R}_{+}^{*}$ modelling the distributions of transaction sizes. For $i\in\left\{1,\dots,N\right\}$, $\Lambda^{i,b}$ and $\Lambda^{i,a}$ are positive functions satisfying the following classical hypotheses (see [11, 12] for similar assumptions):

• •

$\Lambda^{i,b}$ and $\Lambda^{i,a}$ are twice continuously differentiable.

• •

$\Lambda^{i,b}$ and $\Lambda^{i,a}$ are strictly decreasing, with $\Lambda^{i,b^{\prime}}<0$ and $\Lambda^{i,a^{\prime}}<0$.

• •

$\underset{\delta\rightarrow+\infty}{\text{lim}}\Lambda^{i,b}(\delta)=\underset{\delta\rightarrow+\infty}{\text{lim}}\Lambda^{i,a}(\delta)=0$.

• •

\underset{\delta\in\mathbb{R}}{\text{sup}}\frac{\Lambda^{i,b}(\delta)\Lambda^{i,b^{\prime\prime}}(\delta)}{\big{(}\Lambda^{i,b^{\prime}}(\delta)\big{)}^{2}}<2 and \underset{\delta\in\mathbb{R}}{\text{sup}}\frac{\Lambda^{i,a}(\delta)\Lambda^{i,a^{\prime\prime}}(\delta)}{\big{(}\Lambda^{i,a^{\prime}}(\delta)\big{)}^{2}}<2.

The above conditions are sufficiently general to allow for the vast majority of relevant forms of intensities: the exponential intensities initially introduced in [1] and used in most of the literature, logistic intensities as in [2], or many SU Johnson intensities as in [15].

In addition to quoting prices for the $N$ options, the market maker can buy and sell the underlying asset. We assume that the market for that asset is liquid enough to ensure a perfect $\Delta$-hedging.

Remark 5**.**

*In practice, for a portfolio that is not vega-hedged, it is usually suboptimal to perfectly $\Delta$-hedge the portfolio because of the correlation between the spot process and the instantaneous variance process. Nevertheless, we assume here for the sake of simplicity that $\Delta$-hedging is carried out in continuous time. A study of the optimal position in the underlying asset and its consequence on our problem is carried out in Appendix A.1.

In what follows, we denote by $(\Delta_{t})_{t\in[0,T]}$ the $\Delta$ of the portfolio:

[TABLE]

The resulting dynamics for the cash process $(X_{t})_{t\in[0,T]}$ of the market maker is:

[TABLE]

We denote by $(V_{t})_{t\in[0,T]}$ the process for the Mark-to-Market (MtM) value of the market maker’s portfolio (cash, shares, and options), i.e.,

[TABLE]

The dynamics of that process is given by

[TABLE]

For all $i\in\left\{1,\dots,N\right\}$, the vega of the $i$-th option is defined as

[TABLE]

Therefore, we can rewrite the dynamics of the portfolio as

[TABLE]

Following the academic literature on market making, we can consider two objective functions. As in the initial Avellaneda and Stoikov setting [1] (see also [11, 12, 14]), we can consider the following expected utility objective function:

[TABLE]

where $\gamma>0$ is the risk–aversion parameter of the market maker. Instead, as in [5, 6, 7], but also in [12], we can consider a risk–adjusted expectation for the objective function, i.e.

[TABLE]

The second objective function in our case writes

[TABLE]

where

[TABLE]

These two objective functions are close to one other in practice. Guéant showed in [12] that they give similar optimal quotes in practical examples. Furthermore, in many cases, the expected utility framework with exponential utility function can be reduced to the maximization of the expected PnL minus a quadratic penalty of the above form, up to a change in the intensity functions (see [20]).

In what follows, we consider the second framework. Therefore, we define the value function

[TABLE]

associated with (2) as

[TABLE]

where $\mathcal{A}_{t}$ is the set of admissible controls ($\mathbb{F}$-predictable $\mathbb{R}^{2N}$-valued maps that are bounded from below by $\delta_{\infty}$) defined on $[t,T]$.

2.3 Assumptions and approximations

The above stochastic optimal control problem can be addressed from a theoretical point of view using an approach similar to that of [12]. However, when it comes to approximating the optimal quotes a market maker should set for the $N$ options, classic numerical methods are of no help because the value function $u$ has $N+2$ variables (in addition to the time variable). In order to beat the curse of dimensionality we propose a method based on the following assumptions/approximations:

Assumption 1**.**

We approximate the vega of each option over $[0,T]$ by its value at time $t=0$, namely

[TABLE]

Assumption 2**.**

We assume that the set of authorized inventories is associated with vega risk limits, i.e.

[TABLE]

*where $\overline{\mathcal{V}}\in\mathbb{R}_{+}^{*}$ is the vega risk limit of the market maker.

The first assumption is acceptable if $T$ is not too large. This raises in fact the deep question of the reasonable value of $T$, as there is no natural choice for the horizon of the optimization problem. In practice, $T$ has to be sufficient large to allow for several transactions in many options and small enough for the constant-vega approximation to be relevant (and smaller than the maturities of the options). It is also noteworthy, although it is time-inconsistent, that one can use the output of the model (with the constant-vega approximation) over a short period of time and then run the model again with updated vegas. This is a classical practice in applied optimal control when the parameters are estimated online.101010In Appendix A.2 we propose a method to relax the constant-vega assumption. This method is based on a Taylor expansion around the constant-vega case. The curse of dimensionality is tamed by the reduction of the problem to a Monte-Carlo simulation.

The second assumption states that risk limits are related to the only source of risk (as the portfolio is $\Delta$-hedged). This is a natural assumption. The only drawback is that no risk limit can be set to individual options.

3 An approximate solution to the problem

3.1 Change of variables: beating the curse of dimensionality

Under the above assumptions, the $N+2$ state variables can be replaced by only two: the instantaneous variance and the vega of the portfolio. This portfolio vega, defined by $\mathcal{V}_{t}^{\pi}:=\sum_{i=1}^{N}q_{t}^{i}\mathcal{V}^{i}$, has the following dynamics:

[TABLE]

It is clear then that the value function $u$ verifies

[TABLE]

where

[TABLE]

Because $(\nu_{t},\mathcal{V}^{\pi}_{t})_{t\in[0,T]}$ is a Markov process, the problem boils down, under the two above assumptions, to a low-dimensional optimal control problem where the two state variables are driven by $2N$ controlled point processes and a standard Brownian motion.

3.2 Hamilton-Jacobi-Bellman equation and optimal controls

Following [21], the Hamilton-Jacobi-Bellman equation associated with (3) is given by

[TABLE]

with final condition $v(T,\nu,\mathcal{V}^{\pi})=0$, where

[TABLE]

We end up therefore with a low-dimensional functional equation of the Hamilton-Jacobi-Bellman type.

Once the value function is known, the optimal controls, which are the optimal mid-to-bid and ask-to-mid associated with the $N$ options, are given by the following formula (see [2, 12]):

[TABLE]

Remark 6**.**

In the case where $a_{\mathbb{P}}=a_{\mathbb{Q}}$, it is evident that $v$ does not depend on $\nu$. In that case indeed, we have $v(t,\nu,\mathcal{V}^{\pi})=w(t,\mathcal{V}^{\pi})$ where $w$ is solution of the simpler Hamilton-Jacobi-Bellman equation

[TABLE]

with final condition $w(T,\mathcal{V}^{\pi})=0$.

4 Numerical results

4.1 Model parameters

In this section we consider a book of options and derive the optimal quotes using the above approach.

For this purpose, we consider an underlying stock with the following characteristics:

Stock price at time $t=0$: $S_{0}=10\ \text{€}.$

• -

Instantaneous variance at time $t=0$: $\nu_{0}=0.0225\text{ year}^{-1}$.

• -

Heston model with $a_{\mathbb{P}}(t,\nu)=\kappa_{\mathbb{P}}(\theta_{\mathbb{P}}-\nu)$ where $\kappa_{\mathbb{P}}=2\text{ year}^{-1}$ and $\theta_{\mathbb{P}}=0.04\text{ year}^{-1}$, and $a_{\mathbb{Q}}(t,\nu)=\kappa_{\mathbb{Q}}(\theta_{\mathbb{Q}}-\nu)$ where $\kappa_{\mathbb{Q}}=3\text{ year}^{-1}$ and $\theta_{\mathbb{Q}}=0.0225\text{ year}^{-1}$.

• -

Volatility of volatility parameter: $\xi=0.2\text{ year}^{-1}$.

• -

Spot-variance correlation: $\rho=-0.5$.

We consider the case of a market maker dealing with 20 European call options written on that stock where the strike$\times$maturity couples are the elements $(K^{i},T^{i})_{i=1,\ldots,20}$ of the set $\mathcal{K}\times\mathcal{T}$, where

[TABLE]

The associated implied volatility surface is plotted in Figure 1.111111This plot has been computed using $10^{5}$ Monte-Carlo simulations for each option.

As far as the liquidity parameters of these options are concerned, we consider the following intensity functions:

[TABLE]

where $\lambda^{i}=\frac{252\times 30}{1+0.7\times|S_{0}-K^{i}|}\text{ year}^{-1}$, $\alpha=0.7$, and $\beta=150\text{ year}^{\frac{1}{2}}$. The choice of $\lambda^{i}$ corresponds to $30$ requests per day for at-the-money options, and decreases to 12.5 for the most in- and out-the-money options. The choice of $\alpha$ corresponds to a probability of $\frac{1}{1+e^{0.7}}\approx 33\%$ to trade when the answered quote is the mid-price. The choice of $\beta$ corresponds to a probability of $\frac{1}{1+e^{-0.8}}\approx 69\%$ to trade when the answered quote corresponds to an implied volatility $1\%$ better for the client and a probability of $\frac{1}{1+e^{2.2}}\approx 10\%$ to trade when the answered quote corresponds to an implied volatility $1\%$ worse for the client.

Moreover, we assume transactions of constant size with $z^{i}=\frac{5\cdot 10^{5}}{\mathcal{O}_{0}^{i}}\text{ contracts}$ for option $i$. This corresponds approximately to $500000\text{ \text{€}}$ per transaction.121212This is only an approximation as trade sizes are in number of options and option prices move. In other words, the measures $\mu^{i,b}$ and $\mu^{i,a}$ are Dirac masses at $z^{i}$.

Regarding the vega risk limits, we consider $\overline{\mathcal{V}}=10^{7}\text{ \text{€}}\cdot\text{year}^{\frac{1}{2}}$.

As far as the time horizon is concerned, we consider $T=0.0012$ year (i.e. $0.3$ day). This short time horizon surprisingly ensures convergence towards stationary quotes at time $t=0$ (see Figure 3 below).

Eventually, we consider a risk aversion parameter $\gamma=1\cdot 10^{-3}\text{ \text{€}}^{-1}$.

4.2 Optimal quotes

Using a monotone explicit Euler scheme with linear interpolation on a grid of size $180\times 30\times 40$, we approximate the value function solution to (4) (with Neumann conditions at the boundaries in $\nu$) on the domain $[0,T]\times[0.0144,0.0324]\times\left[-\overline{\mathcal{V}},\overline{\mathcal{V}}\right]$. This value function is plotted in Figure 2. We see that the dependence in $\nu$ is low. This is not surprising given that, at the time scale of the market making problem, the discrepancy between $a_{\mathbb{P}}$ and $a_{\mathbb{Q}}$ could be ignored.

From that value function, we deduce the optimal bid and ask quotes of the market maker (in fact mid-to-bid and ask-to-mid) for each option as a function of the portfolio vega. As mentioned above, we chose $T=0.0012$ year (i.e. $0.3$ day) – a choice that ensures convergence of the optimal quotes to their stationary values (see Figure 3).

Focusing on the asymptotic values, we now present in Figures 4, 5, 6, 7, and 8, the optimal mid-to-bid quotes as a function of the portfolio vega for each strike and maturity. More precisely, as the options we consider can have very different prices, we consider instead of the optimal bid quotes themselves the ratio between each optimal mid-to-bid quote and the price (at time $t=0$) of the corresponding option.

We clearly see that the optimal mid-to-bid quotes are increasing with the vega of the portfolio. This is expected since the options all have a positive vega: the higher the vega of the portfolio, the lower the willingness to buy another option, and therefore the lower the bid price (and the higher the mid-to-bid) proposed by the market maker.

Since option market makers usually reason in terms of implied volatility, we plot the same optimal bid quotes for the 20 options in terms of implied volatility (divided by the implied volatility at time $t=0$) in Figures 9, 10, 11, 12, and 13. We obtain of course decreasing curves.

Conclusion

In this article, we tackled the problem of an option market maker dealing with options on a single underlying asset.131313As noted when presenting this paper to a wide audience, our method could easily be extended to the case of multiple underlying assets using the same method as in [2] if the instantaneous variance processes of the different assets are driven by a few factors. Using a constant-vega approximation, we showed how to reduce the problem to a low-dimensional functional equation whose solution can easily be approximated using an explicit Euler scheme and linear interpolation techniques. Furthermore, our method scales linearly in the number of options and can therefore be used with large books of options. Our method is illustrated by an example involving 20 European calls, but our model can be used with any European options.

Appendix A Appendices

A.1 An alternative to the $\Delta$-hedging assumption

Throughout the body of this paper, we assumed that the market maker ensured $\Delta$-hedging. In this appendix, we show that this assumption can be relaxed without much change in the reasoning.

Let us introduce the process $(q_{t}^{S})_{t\in[0,T]}$ representing the inventory of the market maker in the underlying asset. The dynamics of the cash process of the market maker $(X_{t})_{t\in[0,T]}$ rewrites as

[TABLE]

The Mark-to-Market value of the portfolio writes

[TABLE]

and its dynamics is

[TABLE]

Denoting by $\Delta^{\pi}_{t}:=\sum_{i=1}^{N}q_{t}^{i}\partial_{S}O^{i}(t,S_{t},\nu_{t})$ the $\Delta$ of the market maker’s portfolio at time $t$, our mean-variance optimization problem becomes

[TABLE]

where

[TABLE]

Noticing that

[TABLE]

we easily see that the variance term is minimized for $q^{S}={q^{S}}^{*}$ where

[TABLE]

and that its minimum value is

[TABLE]

Therefore, the optimization problem boils down to

[TABLE]

[TABLE]

and we recover the same optimization problem as in the body of the paper, except that the volatility of volatility parameter $\xi$ is multiplied by $\sqrt{1-\rho^{2}}$ to account for the reduction of risk made possible by the optimal trading strategy in the underlying asset in presence of vol-spot correlation.

A.2 Beyond the constant-vega assumption

In this appendix we propose a method to relax our main assumption: the constant-vega approximation.

If, for all $i\in\{1,\ldots,N\}$, the process $\left(\mathcal{V}^{i}_{t}\right)_{t\in[0,T]}$ stays close to its initial value $\mathcal{V}^{i}_{0}=:\mathcal{V}^{i}$, then it is reasonable to consider a perturbative approach around the constant-vega approximation. In particular, instead of assuming that $\sum_{i=1}^{N}q^{i}_{t}\partial_{\sqrt{\nu}}O^{i}(t,S_{t},\nu_{t})\simeq\sum_{i=1}^{N}q^{i}_{t}\partial_{\sqrt{\nu}}O^{i}(0,S_{0},\nu_{0})=\sum_{i=1}^{N}q^{i}_{t}\mathcal{V}^{i}=\mathcal{V}_{t}^{\pi}$, we consider the expansion

[TABLE]

and we consider an expansion of the value function $u$ of the following form:

[TABLE]

Assuming that $\mathcal{Q}=\mathbb{R}^{N}$ and noting that the Hamilton-Jacobi-Bellman equation associated with $u$ is

[TABLE]

with terminal condition equal to [math], the first-order term in $\varepsilon$ in the Taylor expansion gives

[TABLE]

with terminal condition equal to [math].

This equation is linear and therefore $\varphi(t,S,\nu,q)$ admits a Feynman-Kac representation that tames the curse of dimensionality for practical applications:

[TABLE]

where, for each $i\in\{1,\ldots,N\},$ the processes $N^{i,b}$ and $N^{i,a}$ have respective intensities

[TABLE]

with as before $dq^{i}_{t}=\int_{\mathbb{R}_{+}^{*}}z\left(N^{i,b}(dt,dz)-N^{i,a}(dt,dz)\right)$.

Subsequently, the function $\varphi$ can be computed using a Monte-Carlo algorithm and quotes accounting for the variation of the vegas can therefore be computed (to the first order in $\varepsilon$).

A.3 On the construction of the processes $N^{i,b}$ and $N^{i,a}$

Let us consider a new filtered probability space \big{(}\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\in\mathbb{R}_{+}},\tilde{\mathbb{P}}\big{)}. For the sake of simplicity, assume that there is only one option (the generalization is straightforward). Let us introduce $\bar{N}^{b}$ and $\bar{N}^{a}$ two independent compound Poisson processes of intensity 1 whose increments follow respectively the distributions $\mu^{b}(dz)$ and $\mu^{a}(dz)$ with support on $\mathbb{R}_{+}^{*}$. We denote by $\bar{N}^{b}(dt,dz)$ and $\bar{N}^{a}(dt,dz)$ the associated random measures. Let $N^{b}$ and $N^{a}$ be two processes, starting at 0, solutions of the coupled stochastic differential equation:

[TABLE]

Then, under $\tilde{\mathbb{P}},$ $N^{b}$ and $N^{a}$ are two $\mathbb{R}^{*}_{+}$-marked point processes with respective intensity kernels

[TABLE]

where $q_{t}:=N^{b}_{t}-N^{a}_{t}.$ We denote by $N^{b}(dt,dz)$ and $N^{a}(dt,dz)$ the associated random measures. For each $\delta\in\mathcal{A}$, we introduce the probability measure $\tilde{\mathbb{P}}^{\delta}$ given by the Radon-Nikodym derivative

[TABLE]

where $\left(L_{t}^{\delta}\right)_{t\geq 0}$ is the unique solution of the stochastic differential equation

[TABLE]

with $L_{0}^{\delta}=1$, where $\tilde{N}^{b}$ and $\tilde{N}^{a}$ are the compensated processes associated with $N^{b}$ and $N^{a}$ respectively.

We then know from Girsanov theorem that under $\tilde{\mathbb{P}}^{\delta}$, the jump processes $N^{b}$ and $N^{a}$ have respective intensity kernels

[TABLE]

as in the body of the paper.