An approximate solution for options market-making in high dimension
AAn approximate solution for options market-making in highdimension ∗ Bastien
Baldacci † Joffrey
Derchu ‡ Iuliia
Manziuk § September 3, 2020
Abstract
Managing a book of options on several underlying involves controlling positions of severalthousands of financial assets. It is one of the most challenging financial problems involving bothpricing and microstructural modeling. An options market maker has to manage both long- andshort-dated options having very different dynamics. In particular, short-dated options inventoriescannot be managed as a part of an aggregated inventory, which prevents the use of dimension-ality reduction techniques such as a factorial approach or first-order Greeks approximation. Inthis paper, we show that a simple analytical approximation of the solution of the market maker’sproblem provides significantly higher flexibility than the existing algorithms designing optionsmarket making strategies.
Keywords:
Option market making, stochastic control, partial differential equations.
After the electronification of delta-one trading, where high-frequency trading companies provide thevast majority of the liquidity on several thousands of assets, systematic options trading seems to bethe next main challenge in quantitative trading. For assets listed in a central limit order book, asin the equity world, execution and market making are carried out using algorithms. However, forless mature markets such as a great proportion of fixed income securities, systematic market makingactivities are driven by request-for-quote (RFQ for short) systems: the client sends a request to obtaina buy or sell price, for a given quantity of a security, to one or several market makers, who proposeprices based on their current positions. Given the prices, the client accepts or refuses one or several ∗ This work benefits from the financial support of the Chaires Analytics and Models for Regulation, Financial Risk andFinance and Sustainable Development. Bastien Baldacci and Joffrey Derchu gratefully acknowledge the financial sup-port of the ERC Grant 679836 Staqamof. The authors would like to thank Mathieu Rosenbaum (Ecole Polytechnique),Dylan Possamai (ETH Zurich), Greg Sidier (G-Research) and Olivier Guéant (Université Paris-1 Panthéon-Sorbonne)for numerous fruitful discussions. In particular, Mathieu Rosenbaum deserves warm thanks for his careful reading ofthe paper and his many suggestions to improve its quality. † École Polytechnique, CMAP, 91128, Palaiseau, France, [email protected]. ‡ École Polytechnique, CMAP, 91128, Palaiseau, France, joff[email protected]. § École Polytechnique, CMAP, 91128, Palaiseau, France, [email protected]. a r X i v : . [ q -f i n . T R ] S e p ransactions. On OTC markets, such as the corporate bonds market, the proportion of the volumetraded with electronic market makers is increasing.For more than three decades, the optimal market making problem on cash markets has been theobject of many academic studies. The two primary references are [14, 18]. In [14], the authors pro-posed a simple three-period economic model representing the interaction between market makers andmarket-takers and analyzed the equilibrium state. In [18], the authors studied the behavior of a mar-ket maker facing a stochastic demand and an inventory risk and obtained his optimal strategy usingthe stochastic optimal control theory. In the well-known paper of Avellaneda and Stoikov [1] inspiredby this framework, they proposed a model applicable for market making on the order-driven marketat the high-frequency. However, due to the continuous nature of the market maker’s spreads, and theassumption that the underlying asset is a diffusion process, this model is more suited to quote-drivenmarkets such as corporate bonds market.By providing a rigorous analysis of the stochastic control problem of [1], the authors of [15] show, inthe case of a CARA utility function, that the market maker’s problem boils down to a system of linearordinary differential equations. A large part of the contribution to the market making literature comesfrom works of Cartea and Jaimungal, who enriched the initial model by introducing alpha signals,ambiguity aversion, competition with other agents, see, for example, [8, 9, 10, 11]. In these works,they consider a risk-adjusted expectation maximization. As shown in [19], the solution of such for-mulation can also be obtained through CARA utility maximization after a suitable intensity functiontransformation. More recently, multi-asset market making, still on linear markets, has been addressedthrough reinforcement learning techniques, see [16], and dimensionality reduction techniques, as in [5].Regardless of how rich is the part of academic literature considering linear markets, the part study-ing optimal market making on options is far less extensive. A reasonable market making modelfor options has to take into account a lot of stylized facts. First, option market makers trade si-multaneously derivatives and the corresponding underlying, which implies the construction of morecomplex trading strategies taking into account, for example, the Delta-Vega hedging. Consequently,one needs to impose a factorial stochastic volatility model, possibly with jumps, on the underlyingasset. Second, option market makers need to manage several thousands of positions, which lead tovery high-dimensional problems that cannot be solved using classical numerical schemes. Even if ma-chine learning techniques are used, involving, for example, deep reinforcement learning methods (see[16, 22]), the computation time can still be an obstacle. The market maker has to answer a requestfrom a client in a given time, which can be insufficient to recalibrate the model if some parameterchanges need to be applied (for example, the correlation structure). Finally, when dealing with shortmaturity options, the market maker has to manage the positions individually to avoid sudden highexposure due to the Gamma of a specific position. This specificity prevents the use of some dimen-sionality reduction techniques.In the existing academic literature, options market making is addressed in [2, 12, 21]. In [21], theauthors consider three different settings for a market maker managing a single option and its underly-ing. The first setting is a complete market with continuous trading in the perfectly liquid underlying.The second is a complete market with an illiquid underlying where the market maker sets bid andask quotes in the option and the stock. The third is an incomplete market with residual risks due to2tochastic volatility and overnight jumps in the stock price. In [12], the authors consider a marketmaker in charge of a single option in a framework à la Avellaneda-Stoikov, where an underlying followsa one-factor stochastic volatility model, and the market maker is always Delta-hedged. They provideoptimal bid and ask quotes for the option taking into account the risk of model misspecification.Finally, in [2], the authors consider a perfectly Delta-hedged market maker in charge of a book ofoptions with long maturities, whose prices are driven by a stochastic volatility model. The only riskfactor comes from the Brownian motion driving the volatility of the underlying. Using a first-orderapproximation of the Vega of the portfolio, they show that the problem of an options market makerboils down to a three-dimensional Hamilton-Jacobi-Bellman (HJB) equation, which can be solvedusing classical finite difference schemes. By linearizing the value function of the market maker aroundthe Vegas at the initial time, they provide a way to relax the constant Vega assumption. However,the disadvantage of this approach is its time-consumption due to the necessity to simulate inventorytrajectories. Moreover, the constant Vega assumption, making the control problem time-inconsistent,is only valid for a market maker in charge of long-dated options where possible jumps in the under-lying do not influence the global risk position drastically. Finally, if one adds other Greeks such asVanna and Vomma, the model becomes hardly tractable as the HJB equation is in dimension 5.In this article, our goal is to propose a market making algorithm that considers the three specificitiesmentioned above, more flexible and applicable in practice. To this end, we consider a market makerin charge of a book of options on different underlyings. The assets follow a one-factor stochasticvolatility model with jumps, and the Brownian motions driving the underlying and the volatility ofeach asset are correlated. We first consider the case of a perfectly Delta-hedged market maker whomanages his volatility Greeks, namely the Vega, the Vanna, and the Vomma, for all his positions.Inspired by [13], we approximate the jump-diffusion HJB equation corresponding to the optimizationproblem of the market maker with an elliptic Partial Differential Equation (PDE for short). Usingan ansatz quadratic in the inventories, we approximate the value function by a system of non-linearPDEs, which can be easily solved via classical numerical methods for a small number of assets. For anumber of underlyings above two, we recast the ansatz by adding a non-local term, enabling the use ofthe Deep Galerkin method as in [17] to solve the system of PDEs rapidly due to its simple non-linearity.The method presented in this paper has several advantages. First, contrary to [12] and similarlyto [2], the market maker can design trading strategies on a high number of options. Contrary to [2],the market maker controls each position individually, which is particularly important for short-datedoptions that must be managed one by one. Moreover, it enables us to reproduce classic option marketmaking behavior where one option is hedged with another. Second, we allow continuous updates ofthe Greeks (Delta, Vega, Vanna, Vomma) of each option, and the dependence of the intensities oforders arrival on the dynamics of the underlying and its stochastic volatility. This is a major im-provement compared to [2], as the quotes of the market maker are adjusted dynamically with respectto the evolution of both an underlying and stochastic volatility, allowing the problem to be solvedin a time-consistent way. Third, we can use a model for the underlying dynamics with an arbitrarynumber of factors without increasing the computation time. We show numerically how this algorithmoutperforms the one in [2] in terms of average PnL for a portfolio of options, where Vegas vary sig-nificantly.The paper has the following structure: in Section 2, we present the framework of options market3aking and the corresponding optimization problem faced by the market maker. In Section 3, weshow how to simplify the problem by approximating the value function. Finally, Section 4 is devotedto numerical experiments. We consider a filtered probability space (Ω , F , P ) where all stochastic processes are defined, and atime horizon T >
0. We consider d > ( dS it = b i P ( t, S it ) dt + σ i ( t, S it , ν it ) dW i,St + R R Z i ( dt, dz ) dν it = a i P ( t, ν it ) dt + v i P ( t, ν it ) dW i,νt , (2.1)where ( W i,St , W i,νt ) t ∈ R + is a couple of Brownian motions with quadratic covariation given by thecoefficients ρ i = d h W i,S ,W i,ν i dt ∈ ( − , a i P , b i P , v i P , σ i are such that the SDEs (2.1) admit a uniquestrong solution . The processes Z i ( dt, dz ) are marked point processes independent from the Brownianmotions, with intensity kernels κ it ( dz ). We also assume that there exists covariance matrices Σ S , Σ ν ∈M d ( R ) which correspond to the correlation structure of the stocks and the stochastic volatility in theoption book. There also exists a risk-neutral probability measure Q such that ( dS it = σ i ( t, S it , ν it ) d ˆ W i,St + R R Z i ( dt, dz ) dν it = a i Q ( t, ν it ) dt + v i Q ( t, ν it ) d ˆ W i,νt , where ( ˆ W i,St , ˆ W i,νt ) , i ∈ { , . . . , d } are Q − Brownian motions.
Remark 2.1.
As the reader will see in the following, by applying the ansatz detailed in Section 3, onecan use a multi-factor stochastic volatility model for the underlying without increasing the complexityof the algorithm. For example, one can work with the well-known two-factor Bergomi model easily,see [6, 7].
On every underlying i ∈ { , . . . , d } we consider a set of N i European options O i,j of maturity T i,j , for j ∈ { , . . . , N i } . In the above one-factor model, we know that for all ( i, j ) ∈ { , . . . , d } × { , . . . , N i } ,and all t ∈ [0 , T i,j ] such that T < min i,j T i,j , O i,jt = O i,j ( t, S it , ν it ) where O i,j is a solution on [0 , T i,j ) × R of the following partial differential equation under the probability Q :0 = ∂ t O i,j ( t, S i , ν i ) + a i Q ( t, ν i ) ∂ ν i O i,j ( t, S i , ν i ) + 12 (cid:16) σ i ( t, S i , ν i ) (cid:17) ∂ S i S i O i,j ( t, S i , ν i )+ ρ i,i ν i Q ( t, ν i ) σ i ( t, S i , ν i ) ∂ ν i S i O i,j ( t, S i , ν i ) + 12 (cid:16) v i Q ( t, ν i ) (cid:17) ∂ ν i ν i O i,j ( t, S i , ν i )+ Z R (cid:18) O i,j ( t, S i + γ i ( t, z ) , ν i ) − O i,j ( t, S i , ν i ) (cid:19) κ i ( dz ) . As the time horizon T is small compared to the maturity of the options (which can be from one dayup to several years), the terminal condition of the PDEs does not have to be specified. In Section 4,numerical experiments are addressed using European call options but any other option with a path-independent payoff can be considered. We now define the market maker’s problem. In particular, for the sake of readability, we assume that there is no correlation between the volatility process of anasset and the variations of another asset. This assumption can be directly relaxed. .2 The market maker’s problem on OTC markets We consider a market maker in charge of providing bid and ask quotes for the P i ∈{ ,...,d } N i optionsover the period [0 , T ] where T < min i,j T i,j . The bid and ask prices on the option j ∈ { , . . . , N i } ofstock i ∈ { , . . . , d } are defined, for transaction size z , by P i,j,bt = O i,jt − δ i,j,bt ( z ) , P i,j,at = O i,jt + δ i,j,at ( z ) , where (cid:16) δ i,j,bt ( · ) , δ i,j,at ( · ) (cid:17) ∈ A , where A is the set of uniformly bounded F -predictable processes. Theyrepresent the spread on the bid or ask side of the option O i,j . The number or transactions onthese options are defined by marked point processes N i,j,b ( dt, dz ) , N i,j,a ( dt, dz ), with almost surely nosimultaneous jumps, whose respective intensity processes are given byΛ i,j,bt ( S, ν, dz ) = λ i,j,b (cid:16) S, ν, δ i,j,bt ( z ) (cid:17) µ i,j,b ( dz ) , Λ i,j,at ( S, ν, dz ) = λ i,j,a (cid:16) S, ν, δ i,j,at ( z ) (cid:17) µ i,j,a ( dz ) . The couples ( µ i,j,b , µ i,j,a ) are probability measures on R ? + modeling the distribution of transactionsizes for the options. Note that, in our framework, the intensities are allowed to depend on both theunderlying and the stochastic volatility of the assets.The market maker manages his inventory process on each option, that is dq i,jt = Z R ? + z (cid:16) N i,j,b ( dt, dz ) − N i,j,a ( dt, dz ) (cid:17) . For the sake of simplicity, we represent the vector of inventories as follows: q T = (cid:16) q , , . . . , q ,N , . . . , q d,N d (cid:17) ∈ M P dl =1 N l , ( R ) . Assuming perfect Delta-hedging , the ∆ of the portfolio on the i -th asset, i ∈ { , . . . , d } , is given by∆ it = X j ∈{ ,...,N i } ∂ S i O i,j ( t, S it , ν it ) q i,jt , ∆ t = X i ∈{ ,...,d } ∆ it . The cash process of the market maker at time t is defined by dX t = X ( i,j ) ∈{ ,...,d }×{ ,...,N i } (cid:18) Z R ? + z (cid:16) δ i,j,bt ( z ) N i,j,bt ( dt, dz ) + δ i,j,at ( z ) N i,j,at ( dt, dz ) (cid:17) − O i,jt dq i,jt (cid:19) + X i ∈{ ,...,d } (cid:16) S it d ∆ it + d h ∆ i , S i i t (cid:17) . We finally define the Mark-to-Market value of the portfolio of the market maker as V t = X t − X i ∈{ ,...,d } ∆ it S it + X ( i,j ) ∈{ ,...,d }×{ ,...,N i } q i,jt O i,jt . For all ( i, j ) ∈ { , . . . , d } × { , . . . , N i } , the Vega, the Vomma and the Vanna of the option O i,jt aredefined as V i,jt = ∂ √ ν i O i,j ( t, S it , ν it ) = 2 √ ν i ∂ ν i O i,j ( t, S it , ν it ) , This assumption can be relaxed by assuming that the market maker acts on the stock market. This way, themean-variance objective function will take into account the Delta of the portfolio. VO ) i,jt = ∂ √ ν i √ ν i O i,j ( t, S it , ν it ) = 4 ν i ∂ ν i ν i O i,j ( t, S it , ν it ) , ( VA ) i,jt = ∂ S √ ν i O i,j ( t, S it , ν it ) = 2 √ ν i ∂ Sν i O i,j ( t, S it , ν it ) . We also define the vectors e i,j ∈ R P dl =1 N l where e i,jk = { k = P i − l =1 N l + j } and ( e , . . . , e d ) as the canonicalbasis of R d . If we denote by Γ it = v i P ( t,ν it )2 √ ν it P j ∈{ ,...,N i } q i,jt V i,jt , we can write the market maker’s problemas sup δ ∈A E (cid:20) V T − γ X ( i,k ) ∈{ ,...,d } Z T Γ it Γ kt Σ ν,i,k dt (cid:21) . (2.2)Here we penalize the portfolio’s total Vega. Any other penalization could be used, as long as it isquadratic in q . For example, this includes more complicated penalties linked to another position tohedge, or some target for the Greeks. We define the Hamiltonians H i,j,a ( S, ν, p ) = sup δ λ i,j,a ( S, ν, δ ) (cid:16) δ − p (cid:17) , H i,j,b ( S, ν, p ) = sup δ λ i,j,b ( S, ν, δ ) (cid:16) δ − p (cid:17) , and the following processes G ( t, S, ν ) ∈ R P dl =1 N l such that G j ( t, S, ν ) = V k j ,j − (cid:16)P kj − l =1 N l (cid:17) t a k j P ( t, ν k j ) − a k j Q ( t, ν k j )2 √ ν k j + ρ k j ( VA ) k j ,j − (cid:16)P kj − l =1 N l (cid:17) t v k j P ( t, ν k j ) − v k j Q ( t, ν k j )2 √ ν k j σ k j ( t, S k j , ν k j )+ ( VO ) k j ,j − (cid:16)P kj − l =1 N l (cid:17) t (cid:16) v k j P ( t, ν k j ) (cid:17) − (cid:16) v k j Q ( t, ν k j ) (cid:17) ν k j , where k j = i if j ∈ { P i − l =1 N l , . . . , P il =1 N l } , for i ∈ { , . . . , d } . We also define R ( t, S, ν ) ∈ M P dl =1 N l ,d ( R )such that R j,i ( t, S, ν ) = v i P ( t, ν i )2 q ν it V i,j − (cid:16)P kj − l =1 N l (cid:17) t , for j ∈ { i − X l =1 N l , . . . , i X l =1 N l } , i ∈ { , . . . , d } , and 0 otherwise. Finally, denote the diffusion part of the HJB equation as L ( t, S, ν, q, u ) = X i ∈{ ,...,d } b i P ( t, S i ) ∂ S i u ( t, S, ν, q ) + X i ∈{ ,...,d } a i P ( t, ν i ) ∂ ν i u ( t, S, ν, q )+ 12 X ( i,k ) ∈{ ,...,d } ∂ S i S k u ( t, S, ν, q ) σ i ( t, S i , ν i ) σ j ( t, S k , ν j )Σ S,i,k + X i ∈{ ,...,d } Z R κ i ( dz ) (cid:18) u (cid:16) t, S + e i γ i ( t, z ) , ν, q (cid:17) − u ( t, S, ν, q ) (cid:19) + 12 X ( i,k ) ∈{ ,...,d } ∂ ν i ν j u ( t, S, ν, q ) v i P ( t, ν i ) v k P ( t, ν k )Σ ν,i,k + X i ∈{ ,...,d } ∂ ν i S i u ( t, S, ν, q ) ρ i v i P ( t, ν i ) σ i ( t, S i , ν i ) . ∂ t u ( t, S, ν, q ) + L ( t, S, ν, q, u ) + q T G ( t, S, ν ) − γ q T R ( t, S, ν )Σ ν R T ( t, S, ν ) q + X ( i,j ) ∈{ ,...,d }×{ ,...,N i } Z R + zH i,j,b (cid:18) S, ν, u ( t, S, ν, q ) − u ( t, S, ν, q + ze i,j ) z (cid:19) µ i,j,b ( dz )+ X ( i,j ) ∈{ ,...,d }×{ ,...,N i } Z R + zH i,j,a (cid:18) S, ν, u ( t, S, ν, q ) − u ( t, S, ν, q − ze i,j ) z (cid:19) µ i,j,a ( dz ) . (2.3)with terminal condition u ( T, S, ν, q ) = 0. The proof of existence and uniqueness of a viscosity solutionto (2.3) associated to the control problem (2.2) relies on classic arguments of second order viscositysolutions with jumps, see for example [3, 4, 5].
Equation (2.3) is intractable with classical numerical methods when dealing with several optionson several underlyings. Notably, the method proposed in [2] to overcome the constant Vega as-sumption requires Monte-Carlo simulations of high-dimensional inventory trajectories, which is verytime-consuming. In this section, inspired by [13], we propose an approximation of the value functionof the market maker, quadratic with respect to the vector of inventories to reduce the dimensionalityof the problem.A Taylor expansion at 0 on the third variable with respect to (cid:15) gives H i,j,b (cid:18) S, ν, u ( t, S, ν, q ) − u ( t, S, ν, q + (cid:15)ze i,j ) z (cid:19) + H i,j,a (cid:18) S, ν, u ( t, S, ν, q ) − u ( t, S, ν, q − (cid:15)ze i,j ) z (cid:19) = H i,j,b ( S, ν,
0) + H i,j,a ( S, ν,
0) + (cid:15) (cid:16) H i,j,a ( S, ν, − H i,j,b ( S, ν, (cid:17) ∂ q u ( t, S, ν, q )+ (cid:15) (cid:18) H i,j,a ( S, ν, (cid:16) ∂ q u ( t, S, ν, q ) (cid:17) − zH i,j,a ( S, ν, ∂ qq u ( t, S, ν, q ) (cid:19) + (cid:15) (cid:18) H i,j,b ( S, ν, (cid:16) ∂ q u ( t, S, ν, q ) (cid:17) − zH i,j,b ( S, ν, ∂ qq u ( t, S, ν, q ) (cid:19) + o ( (cid:15) ) , and by taking (cid:15) = 1, Equation (2.3) becomes0 = ∂ t u ( t, S, ν, q ) + L ( t, S, ν, q, u ) + q T G ( t, S, ν ) − γ q T R ( t, S, ν )Σ ν R T ( t, S, ν ) q + X ( i,j ) ∈{ ,...,d }×{ ,...,N i } Z R + H i,j,b ( S, ν, − H i,j,b ( S, ν, ∂ q u ( t, S, ν, q )7 12 (cid:18) H i,j,b ( S, ν, (cid:16) ∂ q u ( t, S, ν, q ) (cid:17) − zH i,j,b ( S, ν, ∂ qq u ( t, S, ν, q ) (cid:19) µ i,j,b ( dz ) (3.1)+ X ( i,j ) ∈{ ,...,d }×{ ,...,N i } Z R + H i,j,a ( S, ν,
0) + H i,j,a ( S, ν, ∂ q u ( t, S, ν, q )+ 12 (cid:18) H i,j,a ( S, ν, (cid:16) ∂ q u ( t, S, ν, q ) (cid:17) − zH i,j,a ( S, ν, ∂ qq u ( t, S, ν, q ) (cid:19) µ i,j,a ( dz ) . In the following we will show how a simple ansatz, quadratic with respect to the vector of inventories,leads to significant simplifications. For the sake of the simplicity of the notation, assume that H i,j,a = H i,j,b = H i,j (extension to asymmetric intensities is straightforward). By setting u ( t, S, ν, q ) = θ ( t, S, ν ) + q T θ ( t, S, ν ) − q T θ ( t, S, ν ) q, where θ ∈ R , θ ∈ R P dl =1 N l , θ ∈ M P dl =1 N l ( R ) are solutions of the following system of non-linearPDEs: ∂ t θ ( t, S, ν ) + L ( t, θ , ν, S ) + 2 P ( i,j ) ∈{ ,...,d }×{ ,...,N i } H i,j ( S, ν, R R + (cid:18) zH i,j ( S, ν, θ j,j ( t, S, ν ) + H i,j ( S, ν, (cid:16) θ j ( t, S, ν ) (cid:17) (cid:19) µ i,j ( dz )0 = ∂ t θ ( t, S, ν ) + L ( t, θ , ν, S ) + G ( t, S, ν ) + 4 θ ( t, S, ν )diag (cid:16) H ( S, ν, (cid:17) θ ( t, S, ν )0 = ∂ t θ ( t, S, ν ) + L ( t, θ , ν, S ) − γ R ( t, S, ν )Σ ν R T ( t, S, ν )+ 4 θ ( t, S, ν )diag (cid:16) H ( S, ν, (cid:17) θ ( t, S, ν ) , (3.2)where L ( t, θ, ν, S ) = X i ∈{ ,...,N } b i P ( t, S i ) ∂ S i θ ( t, S, ν ) + X i ∈{ ,...,N } a i P ( t, ν i ) ∂ ν i θ ( t, S, ν )+ 12 X ( i,k ) ∈{ ,...,d } (cid:16) ∂ S i S j θ ( t, S, ν ) σ i ( t, S i , ν i ) σ k ( t, S k , ν k )Σ S,i,k + ∂ ν i ν k θ ( t, S, ν ) v i P ( t, ν i ) v k P ( t, ν k )Σ ν,i,k (cid:17) + X i ∈{ ,...,d } ∂ ν i S i θ ( t, S, ν ) ρ i v i P ( t, ν i ) σ i ( t, S i , ν i )+ Z R κ i ( dz ) (cid:18) θ (cid:16) t, S + e i γ i ( t, z ) , ν (cid:17) − θ ( t, S, ν ) (cid:19)! . and θ ( T, S, ν ) = 0 , θ ( T, S, ν ) = P dl =1 N l , , θ ( T, S, ν ) = P dl =1 N l . In system (3.1), one can note thatthe PDE with respect to θ is independent from the two others, which reduces the overall complexity.It can easily be solved for a small number of underlyings and a large number of options using finitedifference schemes. Note that a higher order expansion does not yield a polynomial solution. How-ever, it is possible to truncate the high degree terms to obtain a polynomial solution. This does notlead to a significant change of the value function or the controls if the penalty term is at most quadratic.We now show some numerical applications of the methodology.8 Numerical results
To perform a comparison with respect to the existing methods, we first recall the methodology of [2].In this article, the authors consider a market maker managing a book of options on a single underlying,and they suppose he is perfectly delta-hedged. We have the following set of market parameters: • d = 1 , N = N = 20: there are 20 call options on a single underlying. • Stock price at time t = 0: S = 100 e . • Instantaneous variance at time t = 0: ν = 0 .
04 year − . • Heston model parameters: b P ( t, S ) = µS , σ ( t, S, ν ) = S √ ν , v P ( t, ν ) = v Q ( t, ν ) = ξ √ ν , with ξ = 0 . − . • a P ( t, ν ) = κ P ( θ P − ν ) , a Q ( t, ν ) = κ Q ( θ Q − ν ), with κ P = κ Q = 2 year − , θ P = θ Q = 0 .
04 year − . • Z ( dt, dz ) = 0: there is no jump in the dynamics of the underlying. • Spot-variance correlation: ρ = − . × maturity couples are the elements ( K j , T j ) , j ∈ { , ..., } of the set K × T , where K = { , , , } , T = { . , . , . , . , . } . These market parameters provide the implied volatility surface as in Figure 1. T i m e t o M a t u r i t y S t r i k e I m p li ed V o l a t ili t y Figure 1: Implied volatility surface associated with the market parameters.We consider mainly in-the-money options with maturity ranging from 3 to 6 months so that, due tothe influence of both Vanna and Vomma, the Vega of the portfolio changes noticeably and the pricesof options are non negligible.We define the following intensity functions:Λ j,a ( S, ν, δ ) = Λ j,b ( S, ν, δ ) = λ j (cid:16) α + β V jt δ (cid:17) , j ∈ { , . . . , N } , where λ j = × . ×| S − K j | year − , α = − .
7, and β = 10 year . The choice of λ j corresponds to 50 requests per day for at-the-money options, and decreases to 13.2 for the mostin-the-money options. The choice of α corresponds to a probability of e − . ≈
66% to trade whenthe answered quote is the mid-price (i.e δ = 0). The choice of β corresponds to a probability of e − . ≈
68% to trade when the answered quote corresponds to an implied volatility 1% better forthe client and a probability of e − . ≈
64% to trade when the answered quote corresponds to animplied volatility 1% worse for the client.We assume transactions of constant size with z j = × O j contracts for option j , in other words, themeasures µ j,b , µ j,a are Dirac masses at z j . This corresponds approximately to 500000 e per transaction.We finally set T = 0 .
004 year (i.e 1 day), and a risk aversion parameter γ = 2 ˙10 − e − .The HJB equation using the constant Vega assumption of [2] is0 = ∂ t u ( t, ν, V π ) + a P ( t, ν ) ∂ ν u ( t, ν, V π ) + 12 νξ ∂ νν u ( t, ν, V π ) + V π a P ( t, ν ) − a Q ( t, ν )2 √ ν − γξ V π ) + X j ∈{ ,...,N } z j H j,b (cid:18) u ( t, ν, V π ) − u ( t, ν, V π + z j V j ) z j (cid:19) + X j ∈{ ,...,N } z j H j,a (cid:18) u ( t, ν, V π ) − u ( t, ν, V π − z j V j ) z j (cid:19) , with terminal condition u ( T, ν, V π ) = 0, and V πt = X j ∈{ ,...,N } z j V j q jt ,H j,a/b ( p ) = sup δ j,a/b Λ j,a/b ( δ j,a/b )( δ j,a/b − p ) . In the case where Vega are not constant, we use the following ansatz: u ( t, S, ν, q ) = θ ( t, S, ν ) + q T θ ( t, S, ν ) + q T θ ( t, S, ν ) q, where θ ∈ R , θ ∈ R N , θ ∈ M N ( R ). Define˜ L ( t, S, ν, θ ) = a P ( t, ν ) ∂ ν θ ( t, S, ν ) + 12 νξ ∂ νν θ ( t, S, ν ) + 12 νS ∂ SS θ ( t, S, ν ) + ρνSξ∂ νS θ ( t, S, ν ) , and assume symmetry of intensity functions, that is H j,b = H j,a = H j , we obtain the following systemof coupled PDEs: ∂ t θ ( t, S, ν ) + ˜ L ( t, S, ν, θ ) + 2 P j ∈{ ,...,N } H j ( S, ν,
0) + 2 P j ∈{ ,...,N } z j H j ( S, ν, θ j,j ( t, S, ν )+ P j ∈{ ,...,N } H j ( S, ν, (cid:16) θ j ( t, S, ν ) (cid:17) ∂ t θ ( t, S, ν ) + ˜ L ( t, S, ν, θ ) + V t a P ( t,ν ) − a Q ( t,ν )2 √ ν + 4 θ ( t, S, ν )diag (cid:18) H ( S, ν, (cid:19) θ ( t, S, ν )0 = ∂ t θ ( t, S, ν ) + ˜ L ( t, S, ν, θ ) − γξ diag (cid:16) V t (cid:17) T N diag (cid:16) V t (cid:17) + 4 θ ( t, S, ν )diag (cid:16) H ( S, ν, (cid:17) θ ( t, S, ν ) , (4.1)10here V t = (cid:16) ∂ √ ν O ( t, S, ν ) , . . . , ∂ √ ν O N ( t, S, ν ) (cid:17) T , = (1 , . . . , T ∈ R N . We first show in Figures 2 and 3 some plots of the value function obtained by solving (4.1). −100000−50000050000100000−150000−100000−500000 50000100000150000−150000−100000−50000050000
Optimal value function with ν = 0.02 −100000−50000050000 −100000−50000050000100000−150000−100000−500000 50000100000150000−400000−300000−200000−1000000
Optimal value function with ν = 0.04 −400000−300000−200000−1000000−100000−50000050000100000−150000−100000−500000 50000100000150000−800000−600000−400000−2000000
Optimal value function with ν = 0.06 −600000−400000−2000000
Figure 2: Value function for different inventories in (97 , .
3) and (98 , .
3) options, inventories in otheroptions assumed to be equal 0, for different values of ν . −100000−50000050000100000−150000−100000−500000 50000100000150000−600000−500000−400000−300000−200000−1000000 Optimal value function with ν = 0.02 −600000−400000−2000000 −100000−50000050000100000−150000−100000−500000 50000100000150000−1750000−1500000−1250000−1000000−750000−500000−2500000
Optimal value function with ν = 0.04 −1500000−1000000−5000000−100000−50000050000100000−150000−100000−500000 50000100000150000−3000000−2500000−2000000−1500000−1000000−5000000
Optimal value function with ν = 0.06 −2000000−10000000
Figure 3: Value function for different inventories in (97 , .
3) and (100 , .
7) options, inventories inother options assumed to be equal 0, for different values of ν .11he value function often has higher values on the diagonals. The market maker can compensate along position in an option with a short position in another one. The values are noticeably lower forhigher values of the volatility.We present in Figure 4 the evolution of the optimal ask quotes with respect to the stochastic volatilityfor the spot S = 100. ν O p t i m a l a sk quo t e T = 0.3
K= 97K= 98K= 99K= 100 0.020 0.025 0.030 0.035 0.040 0.045 0.050 0.055 0.060 ν O p t i m a l a sk quo t e T = 0.4
K= 97K= 98K= 99K= 1000.020 0.025 0.030 0.035 0.040 0.045 0.050 0.055 0.060 ν O p t i m a l a sk quo t e T = 0.5
K= 97K= 98K= 99K= 100 0.020 0.025 0.030 0.035 0.040 0.045 0.050 0.055 0.060 ν O p t i m a l a sk quo t e T = 0.6
K= 97K= 98K= 99K= 1000.020 0.025 0.030 0.035 0.040 0.045 0.050 0.055 0.060 ν O p t i m a l a sk quo t e T = 0.7
K= 97K= 98K= 99K= 100
Figure 4: Optimal ask quotes with respect to ν for different options maturities.We observe the usual increasing behavior of the optimal quotes with respect to both maturity andvolatility of the underlying.In Figure 5, we plot the evolution of the optimal ask quotes with respect to the underlying asset forthe volatility ν = 0 .
04. 12 S O p t i m a l a sk quo t e T = 0.3
K= 97K= 98K= 99K= 100 90.0 92.5 95.0 97.5 100.0 102.5 105.0 107.5 110.0 S O p t i m a l a sk quo t e T = 0.4
K= 97K= 98K= 99K= 10090.0 92.5 95.0 97.5 100.0 102.5 105.0 107.5 110.0 S O p t i m a l a sk quo t e T = 0.5
K= 97K= 98K= 99K= 100 90.0 92.5 95.0 97.5 100.0 102.5 105.0 107.5 110.0 S O p t i m a l a sk quo t e T = 0.6
K= 97K= 98K= 99K= 10090.0 92.5 95.0 97.5 100.0 102.5 105.0 107.5 110.0 S O p t i m a l a sk quo t e T = 0.7
K= 97K= 98K= 99K= 100
Figure 5: Optimal ask quotes with respect to S for different options maturities.The behavior of the optimal quotes with respect to the strike depends on the expiry. We can see thatthe quotes are of the U-shaped nature, the quotes are decreasing in the spot price until some pointdepending on the strike and the expiry, and then become increasing. The inflection point decreaseswith the strike decreasing, and conversely for the expiry date. This way we can see that, for exam-ple, the quote for the option ( K, T ) = (97 , .
7) is monotonously increasing in the spot price for theconsidered grid, which is fairly representative of the possible prices during one day. Conversely, forthe option (
K, T ) = (100 , .
3) the quote is decreasing for almost all values of the grid.In Figure 6, we show the average PnL per request of the trader during the day over 1000 simulations,using the constant Greek approximation of [2] and our algorithm.At the beginning of the trading day, both methods yield a similar PnL per request. Notice that thePnL per request for the method with constant Greek approximation is slightly higher. Indeed theparameters at the beginning of the day correspond to the calibration parameters, and our algorithm ismore conservative as it takes into account the risk that the underlying price could change. However,13fter roughly a tenth of the trading day the method with constant Greek approximations starts tounderperform our algorithm. This underperformance increases along the day as the constant Vegaapproximation becomes less accurate. On the contrary, with our method the PnL per request remainsconstant: there is no need for recalibration.Figure 6: Average PnL per request over the trading day using constant and non-constant Greekapproximations.In Figure 7, we show one of 1000 simulation examples of the trajectories for the Vega of each option.We see that Vegas for this set of options are changing considerably during the day. i ,T i ) = (97.0,0.3)(K i ,T i ) = (98.0,0.3)(K i ,T i ) = (99.0,0.3)(K i ,T i ) = (100.0,0.3)(K i ,T i ) = (97.0,0.4)(K i ,T i ) = (98.0,0.4)(K i ,T i ) = (99.0,0.4)(K i ,T i ) = (100.0,0.4)(K i ,T i ) = (97.0,0.5)(K i ,T i ) = (98.0,0.5) 0 5 10 15 20 25 30 3560708090100110120 (K i ,T i ) = (99.0,0.5)(K i ,T i ) = (100.0,0.5)(K i ,T i ) = (97.0,0.6)(K i ,T i ) = (98.0,0.6)(K i ,T i ) = (99.0,0.6)(K i ,T i ) = (100.0,0.6)(K i ,T i ) = (97.0,0.7)(K i ,T i ) = (98.0,0.7)(K i ,T i ) = (99.0,0.7)(K i ,T i ) = (100.0,0.7) Figure 7: Example of Vega trajectories.Figure 8: Cumulative distribution functions of the PnL over the trading day using both methods.14inally, we present in Figure 8 the cumulative distribution function of the PnL of the trader usingthe constant Greek approximation of [2] and our algorithm. We observe that the tail distribution ofthe PnL using our non-constant Greek approximation is higher compared to the method in [2].
Appendix A The market maker’s problem for large numberof underlyings
In this appendix, we present the system of low-dimensional PDEs analogous to (3.2) for more com-plex cases such as the market making problem on several underlyings or the case where a number ofdifferent options’ parameters is large (over one hundred).We can rewrite the system of ( P i ∈{ ,...,d } N i ) equations (3.2) on θ as a set of d equations by addingthe strike and the maturity to the state variables. The same can be applied for the θ equation. Thisway we obtain a smaller set of equations, though having more dimensions and some non-local terms.Let O i = n ( T i,j , K i,j ) , j ∈ { , ...N i } o be the set of parameters of options on the underlying i ∈{ , ..., d } and let us define ˆ θ i : [0 , T ] × R × R + × O i → R such that, for all j ∈ { , ...N i } ,ˆ θ i ( t, S, ν, ( T i,j , K i,j )) = θ P i − l =1 N l + j ( t, S, ν ) . Similarly for i , i ∈ { , ..., d } , define ˆ θ i ,i : [0 , T ] × R × R + × O i × O i → R such that, for any j ∈ { , ...N i } and l ∈ { , ...N i } ,ˆ θ i ,i ( t, S, ν, ( T i ,j , K i ,j ) , ( T i ,l , K i ,l )) = θ P i − l =1 N l + j, P i − l =1 N l + l ( t, S, ν ) . Then the system of non-linear PDEs (3.2) can be rewritten as ∂ t θ ( t, S, ν ) + L ( t, S, ν, θ ) + 2 P i ∈{ ,...,d } P ( T,K ) ∈ O i H i ( S, ν,
T, K )+ 2 P i ∈{ ,...,d } P ( T,K ) ∈ O i R R + zH i ( S, ν,
T, K )ˆ θ i,i (cid:16) t, S, ν, ( T, K ) , ( T, K ) (cid:17) µ i, ( T,K ) ( dz )+ P i ∈{ ,...,d } P ( T,K ) ∈ O i H i ( S, ν,
T, K ) (cid:18) ˆ θ i (cid:16) t, S, ν, ( T, K ) (cid:17)(cid:19) ∂ t ˆ θ i ( t, S, ν, ( T , K )) + L (cid:16) t, S, ν, ˆ θ i , ( T , K ) (cid:17) + G i ( t, S, ν, ( T , K ))+ 4 P i ∈{ ,...,d } P ( T,K ) ∈ O i ˆ θ i,i ( t, S, ν, ( T , K ) , ( T, K )) H i ( S, ν,
T, K )ˆ θ i (cid:16) t, S, ν, ( T, K ) (cid:17) ∂ t ˆ θ i ,i (cid:16) t, S, ν, ( T , K ) , ( T , K ) (cid:17) + L (cid:16) t, S, ν, ˆ θ i ,i , ( T , K ) , ( T , K ) (cid:17) − γ R i (cid:16) t, S, ν, ( T , K ) (cid:17) Σ ν,i ,i R i (cid:16) t, S, ν, ( T , K ) (cid:17) + 4 P i ∈{ ,...,d } P ( T,K ) ∈ O i ˆ θ i ,i (cid:16) t, S, ν, ( T , K ) , ( T, K ) (cid:17) ˆ H i ( S, ν,
T, K )ˆ θ i ,i (cid:16) t, S, ν, ( T, K ) , ( T , K ) (cid:17) , where (cid:16) ( T , K ) , ( T , K ) (cid:17) ∈ (cid:16) Q i ∈{ ,...,d } O i (cid:17) and, for j ∈ { , ...N i } , l ∈ { , ...N i } , H i ( S, ν, T i,j , K i,j ) = H i,j ( S, ν, , i ( t, S, ν, ( T i,j , K i,j )) = G ( t, S, ν ) P i − l =1 N l + j , R i ( t, S, ν, ( T i,j , K i,j )) = R ( t, S, ν ) P i − l =1 N l + j,i ,µ i, ( T i,j ,K i,j ) = µ i,j , L (cid:16) t, S, ν, ˆ θ i , ( T i,j , K i,j ) (cid:17) = L ( t, S, ν, θ ) P i − l =1 N l + j , L (cid:16) t, S, ν, ˆ θ i ,i , ( T i ,j , K i ,j ) , ( T i ,l , K i ,l ) (cid:17) = L ( t, S, ν, θ ) P i − l =1 N l + j, P i − l =1 N l + l . In particular, if d = 1, ˆ θ and ˆ θ are solutions of non-local PDEs in dimensions 5 and 7 respectively.The observed regularity of the solution with respect to the strike and expiry implies that the high-dimensional PDEs can be solved, for example, by a non-local variant of the Deep Galerkin Method,see [17, 20]. References [1] M. Avellaneda and S. Stoikov. High-frequency trading in a limit order book.
QuantitativeFinance , 8(3):217–224, 2008.[2] B. Baldacci, P. Bergault, and O. Guéant. Algorithmic market making for options. arXiv , pagesarXiv–1907, 2019.[3] G. Barles and C. Imbert. Second-order elliptic integro-differential equations: viscosity solutions’theory revisited. In
Annales de l’IHP Analyse non linéaire , volume 25, pages 567–585, 2008.[4] B. Bastien, B. Philippe, D. Joffrey, and R. Mathieu. On bid and ask side-specific tick sizes. arXivpreprint arXiv:2005.14126 , 2020.[5] P. Bergault and O. Guéant. Size matters for otc market makers: viscosity approach and dimen-sionality reduction technique. arXiv preprint arXiv:1907.01225 , 2019.[6] L. Bergomi. Smile dynamics iii.
Available at SSRN 1493308 , 2008.[7] L. Bergomi.
Stochastic volatility modeling . CRC press, 2015.[8] A. Cartea, R. Donnelly, and S. Jaimungal. Algorithmic trading with model uncertainty.
SIAMJournal on Financial Mathematics , 8(1):635–671, 2017.[9] A. Cartea, R. Donnelly, and S. Jaimungal. Enhancing trading strategies with order book signals.
Applied Mathematical Finance , 25(1):1–35, 2018.[10] A. Cartea and S. Jaimungal. Incorporating order-flow into optimal execution.
Mathematics andFinancial Economics , 10(3):339–364, 2016.[11] Á. Cartea, S. Jaimungal, and J. Penalva.
Algorithmic and high-frequency trading . CambridgeUniversity Press, 2015.[12] S. El Aoud and F. Abergel. A stochastic control approach to option market making.
Marketmicrostructure and liquidity , 1(01):1550006, 2015.1613] D. Evangelista and D. Vieira. New closed-form approximations in multi-asset market making. arXiv preprint arXiv:1810.04383 , 2018.[14] S. J. Grossman and M. H. Miller. Liquidity and market structure. the Journal of Finance ,43(3):617–633, 1988.[15] O. Guéant, C.-A. Lehalle, and J. Fernandez-Tapia. Dealing with the inventory risk: a solutionto the market making problem.
Mathematics and financial economics , 7(4):477–507, 2013.[16] O. Guéant and I. Manziuk. Deep reinforcement learning for market making in corporate bonds:beating the curse of dimensionality. arXiv preprint arXiv:1910.13205 , 2019.[17] A. Hirsa and W. Fu. An unsupervised deep learning approach in solving partial-integro differentialequations, 2020.[18] T. Ho and H. R. Stoll. Optimal dealer pricing under transactions and return uncertainty.
Journalof Financial economics , 9(1):47–73, 1981.[19] I. Manziuk.
Optimal control and machine learning in finance: contributions to the literature onoptimal execution, market making, and exotic options . PhD thesis, PhD dissertation, UniversitéParis 1 Panthéon-Sorbonne, 2019.[20] J. Sirignano and K. Spiliopoulos. DGM: A deep learning algorithm for solving partial differentialequations.
Journal of Computational Physics , 375:1339–1364, Dec 2018.[21] S. Stoikov and M. Sağlam. Option market making under inventory risk.
Review of DerivativesResearch , 12(1):55–79, 2009.[22] E. Weinan, J. Han, and A. Jentzen. Deep learning-based numerical methods for high-dimensionalparabolic partial differential equations and backward stochastic differential equations.