Size matters for OTC market makers: general results and dimensionality reduction techniques

TL;DR

This paper introduces a dimensionality reduction method using factor models for multi-asset OTC market making, incorporating variable transaction sizes and price responses, to overcome computational challenges.

Contribution

It proposes a novel factor-based dimensionality reduction technique for multi-asset market making models with variable transaction sizes.

Findings

Effective reduction of model complexity in multi-asset settings

Inclusion of variable transaction sizes enhances model realism

Improved computational efficiency for market making optimization

Abstract

In most OTC markets, a small number of market makers provide liquidity to other market participants. More precisely, for a list of assets, they set prices at which they agree to buy and sell. Market makers face therefore an interesting optimization problem: they need to choose bid and ask prices for making money while mitigating the risk associated with holding inventory in a volatile market. Many market making models have been proposed in the academic literature, most of them dealing with single-asset market making whereas market makers are usually in charge of a long list of assets. The rare models tackling multi-asset market making suffer however from the curse of dimensionality when it comes to the numerical approximation of the optimal quotes. The goal of this paper is to propose a dimensionality reduction technique to address multi-asset market making by using a factor model. Moreover, we generalize existing market making models by the addition of an important feature: the existence of different transaction sizes and the possibility for the market makers in OTC markets to answer different prices to requests with different sizes.

Tables

Table 1. Table 1: Statistics associated with our 2000 simulations starting from zero inventory (with optimal quotes).

Mean PnL	Stdev PnL	Stdev coming from RFQs	Objective function
72081	80432	5959	69293

Table 2. Table 2: Statistics associated with our 2000 simulations starting from zero inventory (with myopic quotes).

Mean PnL	Stdev PnL	Stdev coming from RFQs	Objective function
73410	265906	6211	43953

Table 3. Table 3: Statistics associated with our 2000 simulations starting from zero inventory (with the optimal quotes of the one-factor case).

Mean PnL	Stdev PnL	Stdev coming from RFQs	Objective function
72523	96746	6033	68567

Table 4. Table 4: Statistics associated with our 2000 simulations starting from zero inventory (with the two-factor optimal quotes).

Mean PnL	Stdev PnL	Stdev coming from RFQs	Objective function
61471	64911	5338	59765

Equations

$$d{S}_t^i={\sigma }^{i}d{W}_t^i,$$

$$\begin{array}{r}{\nu }_t^{i,b}(dz)={\Lambda }^{i,b}({\delta }^{i,b}(t,z)){\mu }^{i,b}(dz),\\ {\nu }_t^{i,a}(dz)={\Lambda }^{i,a}({\delta }^{i,a}(t,z)){\mu }^{i,a}(dz),\end{array}$$

$$\forall i\in \{1,\dots ,d\},d{q}_t^i={\int }_{{\mathbb{R}}_{+}^{\ast }}z{J}^{i,b}(dt,dz)-{\int }_{{\mathbb{R}}_{+}^{\ast }}z{J}^{i,a}(dt,dz),$$

$$d{X}_t$$

$$\mathcal{A}$$

$$\mathbb{E}\left[X_T+\sum _{i=1}^d{q}_T^i{S}_T^i-{\ell }_d(q_T)-{\int }_0^T\psi (q_t)dt\right],$$

$$\mathbb{E}\Bigg{[}\int\limits_{0}^{T}\Bigg{\{}\sum_{i=1}^{d}\int_{\mathbb{R}_{+}^{*}}\left(\delta^{i,b}(t,z)z\Lambda^{i,b}(\delta^{i,b}(t,z))\mu^{i,b}(dz)+\delta^{i,a}(t,z)z\Lambda^{i,a}(\delta^{i,a}(t,z))\mu^{i,a}(dz)\right)-\psi(q_{t})\Bigg{\}}dt-\ell_{d}(q_{T})\Bigg{]},$$

$$\begin{split}\mathcal{J}(t,q,(\delta^{i})_{i\in\{1,\ldots,d\}})=\mathbb{E}\Bigg{[}\int\limits_{t}^{T}\Bigg{\{}\sum_{i=1}^{d}\int_{\mathbb{R}_{+}^{*}}\bigg{(}&\delta^{i,b}(s,z)z\Lambda^{i,b}(\delta^{i,b}(s,z))\mu^{i,b}(dz)+\delta^{i,a}(s,z)z\Lambda^{i,a}(\delta^{i,a}(s,z))\mu^{i,a}(dz)\bigg{)}\\ &-\psi\big{(}q^{t,q,(\delta^{i})_{i\in\{1,\ldots,d\}}}_{s}\big{)}\Bigg{\}}ds-\ell_{d}\big{(}q^{t,q,(\delta^{i})_{i\in\{1,\ldots,d\}}}_{T}\big{)}\Bigg{]},\end{split}$$

$$\theta (t,q)=\underset{({\delta }^i{)}_{i\in \{1,\dots ,d\}}\in \mathcal{A}}{\sup }\mathcal{J}(t,q,({\delta }^i{)}_{i\in \{1,\dots ,d\}}),\forall (t,q)\in [0,T]\times {\mathbb{R}}^d.$$

$$0$$

$$H^{i,b}:p\in \mathbb{R}\mapsto \underset{\delta \ge -{\delta }_{\infty }}{\sup }{\Lambda }^{i,b}(\delta )(\delta -p)\text{ and }{H}^{i,a}:p\in \mathbb{R}\mapsto \underset{\delta \ge -{\delta }_{\infty }}{\sup }{\Lambda }^{i,a}(\delta )(\delta -p),$$

$$h_p^i:\delta \in \mathbb{R}⟼{\Lambda }^{i,a}(\delta )(\delta -p).$$

$$p={\tilde{\delta }}^{i,a\ast }(p)+\frac{{\Lambda }^{i,a}({\tilde{\delta }}^{i,a\ast }(p))}{{\Lambda }^{i,a^{\prime }}({\tilde{\delta }}^{i,a\ast }(p))}.$$

$${\tilde{\delta }}^{i,a{\ast }^{\prime }}(p)=\frac{1}{2-\frac{{\Lambda }^{i,a}({\tilde{\delta }}^{i,a\ast }(p)){\Lambda }^{i,a^{\prime \prime }}({\tilde{\delta }}^{i,a\ast }(p))}{({\Lambda }^{i,a^{\prime }}({\tilde{\delta }}^{i,a\ast }(p)){)}^2}}>0,\forall p\in \mathbb{R}.$$

$${\tilde{H}}^{i,a^{\prime }}(p)=-{\Lambda }^{i,a}({\tilde{\delta }}^{i,a\ast }(p))<0.$$

$${\tilde{\delta }}^{i,a\ast }(p)={\left({\Lambda }^{i,a}\right)}^{-1}\left(-{\tilde{H}}^{i,a^{\prime }}(p)\right).$$

$$H^{i,a}(p)=h_p^i({\tilde{\delta }}^{i,a\ast }(p)).$$

$$H^{i,a}(p)=h_p^i(-{\delta }_{\infty }).$$

$${\delta }^{i,a\ast }(p)=\max ({\tilde{\delta }}^{i,a\ast }(p),-{\delta }_{\infty }).$$

$${H^{i,a}}^{\prime }(p)=-{\Lambda }^{i,a}({\tilde{\delta }}^{i,a\ast }(p))=-{\Lambda }^{i,a}({\delta }^{i,a\ast }(p)).$$

$${H^{i,a}}^{\prime }(p)=-{\Lambda }^{i,a}(-{\delta }_{\infty })=-{\Lambda }^{i,a}({\delta }^{i,a\ast }(p)).$$

$${\mathcal{C}}_{\pi }=\left\{u\in C^0({\mathbb{R}}^d,\mathbb{R})|\underset{q\in {\mathbb{R}}^d}{\sup }\left|\frac{u(q)}{1+\pi (q)}\right|<+\infty \right\}.$$

$$F(u):q\in {\mathbb{R}}^d\mapsto \psi (q)-\underset{i=1}{\sum ^d}{\int }_{{\mathbb{R}}_{+}^{\ast }}z{H}^{i,b}\left(\frac{u(q)-u(q+z{e}^i)}{z}\right){\mu }^{i,b}(dz)-\underset{i=1}{\sum ^d}{\int }_{{\mathbb{R}}_{+}^{\ast }}z{H}^{i,a}\left(\frac{u(q)-u(q-z{e}^i)}{z}\right){\mu }^{i,a}(dz)$$

$$\underset{n\to +\infty }{\lim }z{H}^{i,b}\left(\frac{u(q_n)-u(q_n+z{e}^i)}{z}\right)=z{H}^{i,b}\left(\frac{u(q)-u(q+z{e}^i)}{z}\right).$$

$$z{H}^{i,b}\left(\frac{u(q_n)-u(q_n+z{e}^i)}{z}\right)$$

$$\left|\frac{F(u)(q)}{1+\pi (q)}\right|$$

$$\displaystyle\qquad\qquad\qquad\qquad\quad-\underset{i=1}{\overset{d}{\sum}}\int_{\mathbb{R}_{+}^{*}}\frac{z}{1+\pi(q)}H^{i,a}\left(\frac{u(q)-u(q-ze^{i})}{z}\right)\mu^{i,a}(dz)\Bigg{|}$$

$$\le \Vert \psi {\Vert }_{\pi }+\underset{i=1}{\sum ^d}{\int }_{{\mathbb{R}}_{+}^{\ast }}\frac{z}{1+\pi (q)}{H}^{i,b}\left(\frac{u(q)-u(q+z{e}^i)}{z}\right){\mu }^{i,b}(dz)$$

$$+\underset{i=1}{\sum ^d}{\int }_{{\mathbb{R}}_{+}^{\ast }}\frac{z}{1+\pi (q)}{H}^{i,a}\left(\frac{u(q)-u(q-z{e}^i)}{z}\right){\mu }^{i,a}(dz)$$

$$\le \Vert \psi {\Vert }_{\pi }+\underset{i=1}{\sum ^d}{\int }_{{\mathbb{R}}_{+}^{\ast }}\frac{1}{1+\pi (q)}\left(z{H}^{i,b}(0)+L^{i,b}\left|u(q)-u(q+z{e}^i)\right|\right){\mu }^{i,b}(dz)$$

$$+\underset{i=1}{\sum ^d}{\int }_{{\mathbb{R}}_{+}^{\ast }}\frac{1}{1+\pi (q)}\left(z{H}^{i,a}(0)+L^{i,a}\left|u(q)-u(q-z{e}^i)\right|\right){\mu }^{i,a}(dz)$$

$$\le \Vert \psi {\Vert }_{\pi }+\underset{i=1}{\sum ^d}{\int }_{{\mathbb{R}}_{+}^{\ast }}\left(z{H}^{i,b}(0)+L^{i,b}\Vert u{\Vert }_{\pi }+C{L}^{i,b}\Vert u{\Vert }_{\pi }\left(1+z^p\right)\right){\mu }^{i,b}(dz)$$

$$+\underset{i=1}{\sum ^d}{\int }_{{\mathbb{R}}_{+}^{\ast }}\left(z{H}^{i,a}(0)+L^{i,a}\Vert u{\Vert }_{\pi }+C{L}^{i,a}\Vert u{\Vert }_{\pi }\left(1+z^p\right)\right){\mu }^{i,a}(dz).$$

$$F(u)(q)=\psi (q)-\underset{i=1}{\sum ^d}{\int }_{{\mathbb{R}}_{+}^{\ast }}z{H}^{i,b}\left(\frac{u(q)-u(q+z{e}^i)}{z}\right){\mu }^{i,b}(dz)-\underset{i=1}{\sum ^d}{\int }_{{\mathbb{R}}_{+}^{\ast }}z{H}^{i,a}\left(\frac{u(q)-u(q-z{e}^i)}{z}\right){\mu }^{i,a}(dz).$$

Full text

Key words: Market making, Stochastic optimal control, Curse of dimensionality, Integro-differential equations, Risk factor models.

1 Introduction

The electronification of financial markets has changed the traditional role played by market makers. This is evident in the case of most order-driven markets, such as many stock markets, where the traditional market makers in charge of maintaining fair and orderly markets now often compete with high-frequency market making companies. Surprisingly maybe, many OTC markets organized around dealers have also undergone upheaval linked to electronification over the last ten years. This is the case of the corporate bond markets on both sides of the Atlantic ocean where the electronification process is dominated by Multi-dealer-to-client (MD2C) platforms enabling clients to send the same request for quote (RFQ) to several dealers simultaneously and therefore instantly put them into competition with one another. Electronification is also in progress inside investment banks as most of them replace their traders by algorithms to be able to provide clients with quotes for a large set of assets and automate their market making business, at least for small tickets.

Building market making algorithms is a difficult task as the optimization problem faced by a market maker involves both static and dynamic components. A market maker faces indeed a first (static) trade-off: high margin and low volume versus low margin and high volume. A market maker quoting a large bid-ask spread (with no skew) trades indeed rarely, but each trade is associated with large Mark-to-Market (MtM) gain. Conversely, a market maker who quotes a narrow bid-ask spread (with no skew) trades often, but each trade is associated with a small MtM gain. In addition to this simple static trade-off, market makers face a dynamic problem: in a volatile market, they must quote in a dynamic way to mitigate their market risk exposure and, in particular, skew their quotes as a function of their inventory. For example, a single-asset market maker with a long inventory should price in a conservative manner on the bid side and rather aggressively on the ask side, if she wants – a reasonable behaviour – to decrease her probability to buy and increase her probability to sell.

The optimization problem faced by market makers has been addressed in a long list of academic papers. The first two references commonly cited in the market making literature are two economic papers: Grossman and Miller [7] and Ho and Stoll [15]. If the former is a classic from a theoretical point of view, the latter was revived in 2008 by Avellaneda and Stoikov [1] to build the first practical model of single-asset market making. Since then, many models have been proposed, most of them to tackle the same problem of single-asset market making. For instance, [12] provides a rigourous analysis of the stochastic optimal control problem introduced by Avellaneda and Stoikov and proves that the problem boils down to a system of linear ordinary differential equations (ODE) in the case of exponential intensity functions. Cartea et al. ([3, 4, 5]) contributed a lot to the literature and added many features to the initial models: alpha signals, ambiguity aversion, etc. They also considered a different objective function: a risk-adjusted expectation instead of the Von Neumann-Morgenstern expected utility of [1] and [12]. Multi-asset market making has been considered in [9, 10] for both kinds of objective functions and the author shows that the problem boils down, for general intensity functions, to solving a system of (a priori nonlinear) ODEs. Most of the above models are well suited to tackle market making in OTC markets or in order-driven markets when the tick/spread ratio is large. For major stock markets or for some foreign-exchange platforms, other models are better suited such as those of Guilbaud and Pham who really took the microstructure into account (see [13, 14]).333Option market making has also been addressed, see for instance [2, 6, 16].

In spite of a large and growing literature on market making, several problems are rarely addressed. A first example is that of trade sizes: in markets organized around requests for quotes, quotes can and should depend on the size of the requests. A second and more general problem is that of the numerical approximation of the optimal quotes. If optimal quotes can theoretically be computed through the solution of a system of ODEs, the size of that system (which grows exponentially with the number of assets) prevents any concrete computation with grid methods when it comes to portfolios with more than 4 or 5 assets. To our knowledge, the only attempt to approximate the solution of the Hamilton-Jacobi equations associated with market making models in high dimension is [8] in which the authors propose a method – inspired by reinforcement learning techniques – that uses neural networks instead of grids.

In this paper, our goal is twofold. Our first goal is to generalize existing models to introduce a distribution of trade size. This extension is not straightforward as the optimal controls cannot be modeled anymore with real-valued stochastic processes, but must instead be modeled with predictable maps. A consequence, in terms of mathematics, is that the problem does not anymore boil down to a finite system of ODEs but instead to an integro-differential equation of the Hamilton-Jacobi type that can be regarded as an ordinary differential equation in an infinite-dimensional space. Our second goal is to propose a numerical method for approximating the optimal bid and ask quotes of a market maker over a large universe of assets. For that purpose, we show that the real dimension of the problem is not that of the number of assets, but rather that of the rank of the correlation matrix of asset prices. Then, by using a factor model, we show how to approximate the optimal quotes of a market maker. Indeed, if market risk is projected on a low-dimensional space of factors, solving the market making problem boils down to solving a low-dimensional Hamilton-Jacobi equation. In particular, if the number of factors is lower than 3, grid methods can be applied independently of the number of assets. Furthermore, we suggest a Monte-Carlo method to approximate the influence of the residual risk not taken into account when the risk is projected on the space of factors.

In Section 2, we present our market making model with distributed request sizes. We characterize the value function associated with the stochastic optimal control problem as the solution of an integro-differential equation of the Hamilton-Jacobi type by using ODE techniques in a well-chosen infinite-dimensional space and a verification argument. We subsequently provide expressions for the optimal quotes as a function of time, inventory, and request size. In Section 3, we show how the equations can be simplified when the dependence structure between the prices of the different assets can be modeled by risk factors. We then show how this simplification leads to an approximation that helps to tackle the curse of dimensionality by solving a low-dimensional equation on a grid. We also explain how Monte-Carlo simulations could be used to account for the part of the risk not accounted by the factors. We apply these techniques in Section 4 to portfolios of 2 and 30 bonds and discuss the results.

2 Market making with marked point processes

In all this paper, we consider a filtered probability space $\left(\Omega,\mathcal{F},\mathbb{P};\mathbb{F}=(\mathcal{F}_{t})_{t\geq 0}\right)$ satisfying the usual conditions. We assume this probability space to be large enough to support all the processes we introduce.

2.1 Modeling framework and notations

We consider a market maker in charge of $d$ assets. For $i\in\{1,\ldots,d\}$, the reference price of asset $i$ is modeled by a process $(S^{i}_{t})_{t\geq 0}$ with the following dynamics:

[TABLE]

with $S^{i}_{0}$ given, $\sigma^{i}>0$, and $\left((W^{1}_{t},\ldots,W^{d}_{t})^{\prime}\right)_{t\geq 0}$ a $d$-dimensional444The sign ′ designates the transpose operator. It transforms here a line vector into a column vector. Brownian motion with correlation matrix $(\rho^{i,j})_{1\leq i,j\leq d}$, adapted to the filtration $(\mathcal{F}_{t})_{t\geq 0}$. We denote by $\Sigma=\left(\rho^{i,j}\sigma^{i}\sigma^{j}\right)_{1\leq i,j\leq d}$ the variance-covariance matrix associated with the process $(S_{t})_{t\geq 0}=\left((S^{1}_{t},\ldots,S^{d}_{t})^{\prime}\right)_{t\geq 0}$.

Those assets are traded with requests for quote (RFQ): the market maker first receives a request for quote from a client wishing to buy or sell a given asset and she then proposes a price to the client who finally decides whether she accepts to trade at that price or not.

At any time, the market maker must be ready to propose bid and ask quotes to buy and sell any of the $d$ assets. These bid and ask quotes depend on the size $z\in\mathbb{R}_{+}^{*}$ of the RFQ (in all this paper, we use the notation $\mathbb{R}_{+}^{*}:=(0,+\infty)$). For a given asset $i$, they are modeled by maps $S^{i,b},S^{i,a}:\Omega\times[0,T]\times\mathbb{R}_{+}^{*}\rightarrow\mathbb{R}$ which are $\mathcal{P}\otimes\mathcal{B}(\mathbb{R}_{+}^{*})$-measurable, where $\mathcal{P}$ denotes the $\sigma$-algebra of $\mathbb{F}$-predictable subsets of $\Omega\times[0,T]$ and $\mathcal{B}(\mathbb{R}_{+}^{*})$ denotes the Borelian sets of $\mathbb{R}_{+}^{*}$.

For each $i\ \in\ \{1,\ldots,d\}$, we introduce $J^{i,b}(dt,dz)$ and $J^{i,a}(dt,dz)$ two càdlàg $\mathbb{R}_{+}$-marked point processes.555These processes are explicitly constructed in the Appendix. Note that in our model, as in most real OTC markets, there are no simultaneous RFQs in multiple assets.

For $i\in\{1,\ldots,d\}$, we denote by $\left(\nu^{i,b}_{t}(dz)\right)_{t\geq 0}$ and $\left(\nu^{i,a}_{t}(dz)\right)_{t\geq 0}$ the intensity kernels of $J^{i,b}(dt,dz)$ and $J^{i,a}(dt,dz)$, respectively. In addition, we assume that $\left(\nu^{i,b}_{t}(dz)\right)_{t\geq 0}$ and $\left(\nu^{i,a}_{t}(dz)\right)_{t\geq 0}$ verify:

[TABLE]

where for all $i\in\{1,\ldots,d\}$, $\left(\mu^{i,b},\mu^{i,a}\right)$ is a couple of probability measures on $\mathbb{R}_{+}^{*}$, $\delta^{i,b}(t,z)=S^{i}_{t}-S^{i,b}(t,z)$, $\delta^{i,a}(t,z)=S^{i,a}(t,z)-S^{i}_{t}$, and $\left(\Lambda^{i,b},\Lambda^{i,a}\right)$ is a couple of functions satisfying the following hypotheses (H):

• •

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

• •

$\Lambda^{i,b}$ and $\Lambda^{i,a}$ are decreasing, with $\forall\delta\in\mathbb{R}$, ${\Lambda^{i,b}}^{\prime}(\delta)<0$ and ${\Lambda^{i,a}}^{\prime}(\delta)<0$,

• •

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

• •

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

For all $i\in\{1,\ldots,d\}$, $J^{i,b}(dt,dz)$ and $J^{i,a}(dt,dz)$ model respectively the volumes of transactions at the bid and at the ask for asset $i$. The inventory of the market maker, modeled by the $d$-dimensional process $(q_{t})_{t\geq 0}=\left((q^{1}_{t},\ldots,q^{d}_{t})^{\prime}\right)_{t\geq 0}$, has therefore the following dynamics:

[TABLE]

with $q_{0}$ given.

Remark 1**.**

*For a given asset $i$, $\Lambda^{i,.}$ typically has the form $\Lambda^{i,.}(\delta)=\lambda^{i,.}_{RFQ}f^{i,.}(\delta)$, where $\lambda^{i,.}_{RFQ}$ is the (constant) intensity of arrival of requests for quote and $f^{i,.}(\delta)$ gives the probability that a request will result in a transaction given the quote $\delta$ proposed by the market maker. Furthermore, $\mu^{i,.}$ should be seen as the distribution of sizes.

Finally, the process $(X_{t})_{t\geq 0}$ modeling the market maker’s cash account has the dynamics

[TABLE]

We fix $\delta_{\infty}\geq 0$ and define the set $\mathcal{A}$ of admissible controls666We introduce here a unique lower bound for the quotes, independently of the asset, the side, the size, and the time. Generalizations are straightforward. by

[TABLE]

As proved in [10], under assumptions (H), the functions $\delta\in\mathbb{R}\mapsto\delta\Lambda^{i,b}(\delta)$ and $\delta\in\mathbb{R}\mapsto\delta\Lambda^{i,a}(\delta)$ have a unique maximum on $\mathbb{R}$. It is also easy to see that on $[-\delta_{\infty},+\infty)$, they are bounded from below by $-\delta_{\infty}\Lambda^{i,b}(-\delta_{\infty})$ and $-\delta_{\infty}\Lambda^{i,a}(-\delta_{\infty})$, respectively.

For two given continuous penalty functions $\psi:\mathbb{R}^{d}\rightarrow\mathbb{R}_{+}$ and $\ell_{d}:\mathbb{R}^{d}\rightarrow\mathbb{R}_{+}$, modeling the risk aversion of the market maker, we aim at maximizing the objective function

[TABLE]

over the set $\mathcal{A}$ of admissible controls.

Remark 2**.**

*For instance, we can choose $\psi(q)=\frac{\gamma}{2}q^{\prime}\Sigma q$ or $\psi(q)=\gamma\sqrt{q^{\prime}\Sigma q}$ (for $\gamma>0$) and $\ell_{d}(q)=0$, $\ell_{d}(q)=\frac{\zeta}{2}q^{\prime}\Sigma q$ or $\ell_{d}(q)=\zeta\sqrt{q^{\prime}\Sigma q}$ (for $\zeta>0$), as done in [4], [5], [8], and [10].

After applying Itô’s formula to $\left(X_{t}+\sum_{i=1}^{d}q^{i}_{t}S^{i}_{t}\right)_{t\geq 0}$ between [math] and $T$, it is easy to see that the problem is equivalent to maximizing

[TABLE]

over the set of admissible controls $\mathcal{A}$.

We introduce the function $\mathcal{J}:[0,T]\times\mathbb{R}^{d}\times\mathcal{A}\rightarrow\mathbb{R}$ such that, $\forall t\in[0,T]$, $\forall q=(q^{1},\ldots,q^{d})^{\prime}\in\mathbb{R}^{d}$, $\forall(\delta^{i})_{i\in\{1,\ldots,d\}}\in\mathcal{A}$,

[TABLE]

where \big{(}q^{t,q,(\delta^{i})_{i\in\{1,\ldots,d\}}}_{s}\big{)}_{s\geq t} is the inventory process starting in state $q$ at time $t$ and controlled by $(\delta^{i})_{i\in\{1,\ldots,d\}}$.

The value function $\theta:[0,T]\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ of the problem is then defined as follows:

[TABLE]

We will show that $\theta$ is the unique (in a large class of functions) classical solution to the following integro-differential Hamilton-Jacobi (HJ) equation:

[TABLE]

with terminal condition $w(T,q)=-\ell_{d}(q),\forall q\in\mathbb{R}^{d}$, where

[TABLE]

and where $\left(e^{1},\ldots,e^{d}\right)$ denotes the canonical basis of $\mathbb{R}^{d}$.

2.2 Existence and uniqueness of a solution to (2)

Lemma 1**.**

$\forall i\in\{1,\ldots,d\}$*, $H^{i,b}$ and $H^{i,a}$ are two globally Lipschitz continuously differentiable decreasing functions. Moreover, the supremum in the definition of $H^{i,b}(p)$ (respectively $H^{i,a}(p)$) is reached at a unique $\delta^{i,b*}(p)$ (respectively $\delta^{i,a*}(p)$). Furthermore, $\delta^{i,b*}$ and $\delta^{i,a*}$ are continuous and nondecreasing functions.

Proof.

We prove the result only for the ask side. The proof is similar for the bid side.

Let $i\in\{1,\ldots,d\}$. For $p\in\mathbb{R}$, we define

[TABLE]

$h^{i}_{p}$ is a continuously differentiable function, positive for $\delta\in(p,+\infty)$ and nonpositive otherwise. It is easy to prove (see [10]) that there is a unique maximizer $\tilde{\delta}^{i,a*}(p)$ of $h^{i}_{p}$ on $\mathbb{R}$ characterized by

[TABLE]

By the implicit function theorem, $p\in\mathbb{R}\mapsto\tilde{\delta}^{i,a*}(p)$ is continuously differentiable and

[TABLE]

In particular, $\tilde{\delta}^{i,a*}$ is increasing.

We introduce $\tilde{H}^{i,a}:p\in\mathbb{R}\rightarrow\underset{\delta\in\mathbb{R}}{\sup}\ h^{i}_{p}(\delta)$. Then $\forall p\in\mathbb{R}$, we have $\tilde{H}^{i,a}(p)=h^{i}_{p}(\tilde{\delta}^{i,a*}(p))$ and

[TABLE]

So $\tilde{H}^{i,a}$ is decreasing and

[TABLE]

Let us now recall that $\forall p\in\mathbb{R}$, $H^{i,a}(p)=\underset{\delta\geq-\delta_{\infty}}{\sup}\ h^{i}_{p}(\delta)$.

For all $p\in\mathbb{R}$ such that $-\delta_{\infty}\leq\tilde{\delta}^{i,a*}(p)$, we clearly have

[TABLE]

Otherwise, if $\tilde{\delta}^{i,a*}(p)<-\delta_{\infty}$, we can easily see that $h_{p}^{i}(.)$ is increasing on $]-\infty,\tilde{\delta}^{i,a*}(p)]$ and decreasing on $[\tilde{\delta}^{i,a*}(p),+\infty[$, which implies

[TABLE]

This means that the supremum in $H^{i,a}(p)$ is reached at a unique $\delta^{i,a*}(p)$ given by

[TABLE]

In particular, $\delta^{i,a*}$ is continuous and nondecreasing, so $H^{i,a}$ is continuous. Moreover, for all $p\in\mathbb{R}$ such that $\tilde{\delta}^{i,a*}(p)>-\delta_{\infty}$, we have $H^{i,a}(p)=\tilde{H}^{i,a}(p)$ so $H^{i,a}$ is decreasing on $]\tilde{\delta}^{i,a*-1}(-\delta_{\infty}),+\infty[$ and its derivative on this interval is

[TABLE]

On $]-\infty,\tilde{\delta}^{i,a*-1}(-\delta_{\infty})[$, $H^{i,a}$ is affine and its derivative is

[TABLE]

Thus, by continuity of $\delta^{i,a*}$, $H^{i,a}$ is continuously differentiable and decreasing on $\mathbb{R}$. In particular, $|{H^{i,a}}^{\prime}(p)|\leq\Lambda^{i,a}(-\delta_{\infty})$ for all $p\in\mathbb{R},$ so $H^{i,a}$ is Lipschitz.

∎

In what follows, we denote by $L^{i,a}$ the Lipschitz constant of $H^{i,a}$ for all $i\in\{1,\ldots,d\},$ and we define similarly $L^{i,b}$ the Lipschitz constant of $H^{i,b}$ for all $i\in\{1,\ldots,d\}.$

For $\pi\in C^{0}(\mathbb{R}^{d},\mathbb{R}_{+})$, let us consider $\mathcal{C}_{\pi}$ the following vector space:

[TABLE]

Equipped with the norm $u\in\mathcal{C}_{\pi}\mapsto\|u\|_{\pi}=\underset{q\in\mathbb{R}^{d}}{\sup}\left|\frac{u(q)}{1+\pi(q)}\right|$, $\mathcal{C}_{\pi}$ is a Banach space.

We now consider for the rest of the paper that there exists $p\in\mathbb{N}^{*}$ and $C>0$ such that:

• •

$\forall q\in\mathbb{R}^{d},\pi(q)\leq C\left(1+\|q\|^{p}\right)$,

• •

$\forall q,y\in\mathbb{R}^{d},\frac{1+\pi(q+y)}{1+\pi(q)}\leq C\left(1+\|y\|^{p}\right)$,

• •

$\forall i\in\{1,\ldots,d\}$, $\int_{\mathbb{R}_{+}^{*}}\left(z^{p}\mu^{i,b}(dz)+z^{p}\mu^{i,a}(dz)\right)<+\infty$,

where $\|.\|$ denotes the Euclidean norm on $\mathbb{R}^{d}$.

Moreover, we assume that $\psi,\ell_{d}\in\mathcal{C}_{\pi}$.777This assumption implies in particular that $\psi$ and $\ell_{d}$ have, at most, polynomial growth at infinity.

Remark 3**.**

*For the examples of Remark 2, it is natural to choose a quadratic function $\pi$ such that $\psi,\ell_{d}\leq\pi$. Then, the above assumptions are satisfied for $p=2$ whenever $\mu^{i,b}$ and $\mu^{i,a}$ have a finite second moment.

Proposition 1**.**

For all $u\in\mathcal{C}_{\pi}$, the function

[TABLE]

*is in $\mathcal{C}_{\pi}.$

Proof.

Let $u\in\mathcal{C}_{\pi}$.

Let us consider $q\in\mathbb{R}^{d}$ and a sequence $(q_{n})_{n}$ converging towards $q$.

From the continuity of $\psi$, we have $\lim_{n\to+\infty}\psi(q_{n})=\psi(q)$.

Also, $\forall i\in\{1,\ldots,d\},\forall z\in\mathbb{R}_{+}^{*}$, from the continuity of $H^{i,b}$ and $u$, we have

[TABLE]

Now, we write $H^{i,b}(p)\leq H^{i,b}(0)+L^{i,b}|p|$ so that we get

[TABLE]

which is integrable by assumption. Using the same technique for the terms associated with the ask side and Lebesgue’s dominated convergence theorem, we conclude that $\lim_{n\to+\infty}F(u)(q_{n})=F(u)(q)$, hence the continuity of $F(u)$.

Moreover, for all $q\in\mathbb{R}^{d},$ we have

[TABLE]

We conclude that $\underset{q\in\mathbb{R}^{d}}{\sup}\left|\frac{F(u)(q)}{1+\pi(q)}\right|<+\infty$ and therefore that $F(u)\in\mathcal{C}_{\pi}.$

∎

We can therefore define a functional $F:\mathcal{C}_{\pi}\rightarrow\mathcal{C}_{\pi}$ such that, for all $u\in\mathcal{C}_{\pi}$ and for all $q\in\mathbb{R}^{d},$

[TABLE]

We now come to the main property of the function $F$.

Proposition 2**.**

$F$* is Lipschitz on $\mathcal{C}_{\pi}.$

Proof.

Let $u,v\in\mathcal{C}_{\pi}.$ For all $q\in\mathbb{R}^{d}$, we have

[TABLE]

Therefore

[TABLE]

We obtain therefore,

[TABLE]

By taking the supremum over $q$, we get that there exists a constant $K>0$ such

[TABLE]

We conclude that $F$ is Lipschitz continuous.

∎

The Lipschitz property of $F$ allows to obtain the following existence and uniqueness theorem:

Theorem 1**.**

*There exists a unique function $\mathcal{W}\in C^{1}([0,T],\mathcal{C}_{\pi})$ such that $w:(t,q)\in[0,T]\times\mathbb{R}^{d}\mapsto\mathcal{W}(t)(q)$ is solution to (2) with terminal condition $w(T,q)=-\ell_{d}(q),\forall q\in\mathbb{R}^{d}$.

Proof.

Let us observe that $\mathcal{W}\in C^{1}([0,T],\mathcal{C}_{\pi})$ is solution of the Cauchy problem

[TABLE]

if and only if $w:(t,q)\in[0,T]\times\mathbb{R}^{d}\mapsto\mathcal{W}(t)(q)$ is solution to (2) with terminal condition $w(T,q)=-\ell_{d}(q),\forall q\in\mathbb{R}^{d}$.

As $\left(\mathcal{C}_{\pi},\|.\|_{\pi}\right)$ is a Banach space and $F:\mathcal{C}_{\pi}\rightarrow\mathcal{C}_{\pi}$ is Lipschitz continuous, we know by Cauchy-Lipschitz theorem that there exists a unique maximal solution $\mathcal{W}$ to the above equation, and that this solution is in fact global, meaning in particular that $\mathcal{W}$ is defined on $[0,T]$.

∎

2.3 Verification theorem

We now want to prove that $\theta$ is in fact the function $w$ defined in Theorem 1 and deduce the optimal controls associated with the problem (1) using a verification argument.

Theorem 2**.**

*Let $w$ be the function defined in Theorem 1.

*Let $(t,q)\in[0,T)\times\mathbb{R}^{d}$.

Let us define $(\bar{\delta}^{i})_{i\in\{1,\ldots,d\}}=(\bar{\delta}^{i,b},\bar{\delta}^{i,a})_{i\in\{1,\ldots,d\}}\in\mathcal{A}$ such that $\forall i\in\{1,\ldots,d\},\forall s\in[t,T],\forall z>0$:

[TABLE]

*where $\delta^{i,b*}$ and $\delta^{i,a*}$ are the functions defined in Lemma 1 and $(q_{s})_{t\leq s\leq T}=(q^{t,q,(\bar{\delta}^{1},\ldots,\bar{\delta}^{d})}_{s})_{t\leq s\leq T}$.

*Then, $\theta(t,q)=w(t,q)$ and $(\bar{\delta}^{1},\ldots,\bar{\delta}^{d})$ is an optimal control for our stochastic control problem starting at time $t$ with $q_{t}=q$.

Proof.

Let $(\delta^{i})_{i\in\{1,\ldots,d\}}=(\delta^{i,b},\delta^{i,a})_{i\in\{1,\ldots,d\}}\in\mathcal{A}$ be an arbitrary control and let us denote by $(q_{s})_{s\in[t,T]}$ the process $\left(q_{s}^{t,q,(\delta^{1},\ldots,\delta^{d})}\right)_{s\in[t,T]}$.

Let us first prove that for all $i\in\left\{1,\ldots,d\right\}$,

[TABLE]

Denoting by $M^{w}$ the quantity $\underset{t\in[0,T]}{\sup}\|w(t,\cdot)\|_{\pi}$, we have

[TABLE]

Therefore

[TABLE]

Subsequently, we just have to prove that

[TABLE]

Since $\|q_{s}\|\leq\|q\|+\|q_{s}-q\|$, $\|q_{s}\|^{p}\leq 2^{p-1}\left(\|q\|^{p}+\|q_{s}-q\|^{p}\right)$, and we need to prove that

[TABLE]

As we are working in $\mathbb{R}^{d}$, it is equivalent to prove that

[TABLE]

where $\|(x_{1},\ldots,x_{d})^{\prime}\|_{p}=\left(\sum_{i=1}^{d}|x_{i}|^{p}\right)^{1/p}$.

For that purpose, we introduce for each $j\in\{1,\ldots,d\},$ two independent Poisson processes $N^{j,b}$ and $N^{j,a}$ with respective intensities $\Lambda^{j,b}(-\delta_{\infty})$ and $\Lambda^{j,a}(-\delta_{\infty})$, and $(\xi^{j,b}_{k})_{k\geq 1}$ and $(\xi^{j,a}_{k})_{k\geq 1}$ two sequences of i.i.d. random variables with respective distributions $\mu^{j,b}$ and $\mu^{j,a}.$ Then, we have

[TABLE]

Using the above, we have, for all $i\in\left\{1,\ldots,d\right\}$,

[TABLE]

Of course, we can similarly prove that, for all $i\in\left\{1,\ldots,d\right\}$,

[TABLE]

Now, by applying Itô’s formula, we get

[TABLE]

By taking expectation, we get

[TABLE]

which, by definition of $w$, gives us the following inequality:

[TABLE]

with equality when $(\delta^{i})_{i\in\{1,\ldots,d\}}=\left(\bar{\delta}^{i}\right)_{i\in\{1,\ldots,d\}}$.

In other words,

[TABLE]

with equality when $(\delta^{i})_{i\in\{1,\ldots,d\}}=\left(\bar{\delta}^{i}\right)_{i\in\{1,\ldots,d\}}$.

By taking the supremum over $(\delta^{i})_{i\in\{1,\ldots,d\}}\in\mathcal{A}$, we get $\theta(t,q)=w(t,q)$ and the fact that $\left(\bar{\delta}^{i}\right)_{i\in\{1,\ldots,d\}}$ is optimal.

∎

3 Solving the multi-asset market making problem with factors

Let us now consider the particular case of problem (1) where $\forall q\in\mathbb{R}^{d}$, $\psi(q)=\bar{\psi}\left(q^{\prime}\Sigma q\right)$ and $\ell_{d}(q)=\bar{\ell}_{d}(q^{\prime}\Sigma q)$ for some continuous functions $\bar{\psi}$ and $\bar{\ell}_{d}$ with, at most, polynomial growth at infinity. This particular case covers the examples of the literature (see Remark 2).

If the prices of the $d$ assets are modeled using a small number $k$ of factors, as it is the case in most econometric models of financial asset prices, then the variance-covariance matrix $\Sigma$ takes the form

[TABLE]

where $\beta$ is a $d$-by-$k$ matrix of real coefficients, $V$ the $k$-by-$k$ variance-covariance matrix of the factors, and $R$ the $d$-by-$d$ variance-covariance matrix of the residuals.

If the explanatory power of the factors is high, $R$ should be small compared to $\Sigma$ (in Frobenius norm for instance). Our approach consists in ignoring the residuals, i.e. setting $R$ to [math]. In other words, we project the market risk on a space of factors of dimension $k$. As we shall see in Section 4, this approach provides very good results as measured by the objective function (1).

In what follows, we also discuss an approximation method based on Monte-Carlo simulations to account for the influence of $R$ once one has computed the optimal quotes in the case with no residual risk. The advantages and drawbacks of this additional approximation method will be discussed in Section 4.

3.1 A low-dimensional approximation

Let us now assume that $\Sigma=\beta V\beta^{\prime}$, i.e. $R=0$. Under this assumption, we can write problem (1) as the maximization of

[TABLE]

Using the same ideas as in Section 2, this expression can be written as

[TABLE]

Let us introduce $(f_{t})_{t\in[0,T]}=(\beta^{\prime}q_{t})_{t\in[0,T]}$. Then, the problem of maximizing (3) is equivalent to that of maximizing

[TABLE]

The state process of our problem is now the Markov process $(f_{t})_{t\in[0,T]}$ instead of $(q_{t})_{t\in[0,T]}$: we have reduced the dimension of the problem from $d$ to $k$.

Let us introduce $\tilde{\mathcal{J}}:[0,T]\times\mathbb{R}^{k}\times\mathcal{A}\rightarrow\mathbb{R}$ such that, $\forall t\in[0,T]$, $\forall f=(f^{1},\ldots,f^{k})^{\prime}\in\mathbb{R}^{k}$, $\forall(\delta^{i})_{i\in\{1,\ldots,d\}}\in\mathcal{A}$

[TABLE]

where $(f_{s})_{s\in[t,T]}=(f^{t,f,(\delta^{i})_{i\in\{1,\ldots,d\}}}_{s})_{s\in[t,T]}$ is here the state process starting in state $f$ at time $t$ and controlled by $(\delta^{i})_{i\in\{1,\ldots,d\}}$.

The value function $\tilde{\theta}:[0,T]\times\mathbb{R}^{k}\rightarrow\mathbb{R}$ of the problem is then defined as follows:

[TABLE]

By using the same arguments as in Section 2, we can show that $\tilde{\theta}$ is the unique (in a large class of functions) smooth solution to the following integro-differential Hamilton-Jacobi equation:

[TABLE]

with terminal condition $\tilde{\theta}(T,f)=-\bar{\ell}_{d}(f^{\prime}Vf),\forall f\in\mathbb{R}^{k}$, where $\forall i\in\{1,\ldots,d\},\tilde{e}^{i}=\beta^{\prime}e^{i}$.

Furthermore, the optimal controls are now given by:

[TABLE]

When $R=0$, the problem boils down therefore to finding the solution $\tilde{\theta}$ of (4) with the appropriate terminal condition. In particular, from a numerical point of view, we need to approximate the solution of an equation involving time plus $k$ space dimensions, and this is doable with grid methods if $k$ is small.

3.2 A Monte-Carlo method to take the residual risk into account

As we shall see in Section 4, the above approximation method provides very good results as measured by the value of the objective function (1). Nevertheless, when market risk is projected on a low-dimensional space of factors, there are linear combinations of assets that falsely appear to be risk-free. To prevent trajectories of the inventory visiting too often regions that are falsely associated with low risk, it makes sense to look for methods that account for the residual risk measured by the matrix $R$.

In what follows, we propose an approximation method to take the residual risk into account. The idea consists in considering the first-order expansion in $\varepsilon$ where

[TABLE]

The rationale behind this idea is that, for a factor model with high explanatory power, $R$ should be small and it makes sense therefore to use a perturbative approach.

When $\varepsilon=0,$ we know how to solve the problem, and the value function is given by $\tilde{\theta}.$ To approximate the value function $\theta$ of the problem for $\varepsilon>0$, we consider a first-order expansion of the form

[TABLE]

By plugging this expression into equation (2), we formally get

[TABLE]

and

[TABLE]

Assuming that $\bar{\psi}$ and $\bar{\ell_{d}}$ are $C^{1}$ and performing a Taylor expansion, we obtain

[TABLE]

and

[TABLE]

As $\tilde{\theta}$ verifies (4), we get

[TABLE]

and $\eta(T,q)=-\bar{\ell}^{\prime}_{d}\left((\beta^{\prime}q)^{\prime}V(\beta^{\prime}q)\right)q^{\prime}Rq.$

This equation, although in space-dimension $d$, is linear. Therefore, by the Feynman-Kac representation theorem, we have the following formula:

[TABLE]

where under $\tilde{\mathbb{P}}$, for all $i\in\{1,\ldots,d\}$, $J^{i,b}$ and $J^{i,a}$ have their respective intensity kernels given by

[TABLE]

Remark 4**.**

*It is noteworthy that the dynamics of $(q_{s})_{s\in[t,T]}$ under $\tilde{\mathbb{P}}$ is that associated with the use of the optimal quotes when $R=0$.

Thanks to this probabilistic representation, we can easily compute $\eta(t,q)$ for a given time $t$ and inventory $q$ using a Monte-Carlo method, and therefore easily compute both an approximation of the value function and an approximation of the optimal quotes that account, to the first order, for the residual risk. Of course, in practice, it would be prohibitively expensive in terms of computation time to carry out a Monte-Carlo simulation for all possible values of the inventory, but this method can alternatively be used (online) upon receiving a request for quote for a specific asset and given the current time and inventory (this will be discussed in Section 4.1.4).

Remark 5**.**

In the computation of the optimal quotes associated with asset $i$, one relies on the approximation

[TABLE]

*To compute $\eta(t,q_{t-})-\eta(t,q_{t-}\pm ze^{i})$, the same sample paths should be used in the estimations of $\eta(t,q_{t-})$ and $\eta(t,q_{t-}\pm ze^{i})$. This is the same remark as for the computation of the Greeks of derivatives contracts with Monte-Carlo techniques.

4 Numerical results

4.1 The case of two assets: one factor vs. two factors

4.1.1 Model parameters

In this section, we apply our multi-asset market making model to the case of two highly-correlated assets (here bonds). Our goal is to show that, in this case, the reduced one-factor model gives very similar results to the complete two-factor model. For this purpose, we consider two assets with the following characteristics:

• •

Asset prices: $S^{1}_{0}=S^{2}_{0}=100\ \textrm{€}$.

• •

Volatility of asset 1: $\sigma^{1}=1.2\ \textrm{€}\cdot\textrm{day}^{-\frac{1}{2}}$.

• •

Volatility of asset 2: $\sigma^{2}=0.6\ \textrm{€}\cdot\textrm{day}^{-\frac{1}{2}}$.

• •

Correlation: $\rho=0.9$.

• •

Intensity functions:

[TABLE]

with $\lambda_{RFQ}=30\ \textrm{day}^{-1}$, $\alpha_{\Lambda}=0.7$, and $\beta_{\Lambda}=30\ \textrm{€}^{-1}$. This corresponds to 30 RFQs per day for each asset, a probability of $\frac{1}{1+e^{0.7}}\simeq 33\%$ to trade when the answered quote is the reference price and a probability of $\frac{1}{1+e^{-0.2}}\simeq 55\%$ to trade when the answered quote is the reference price improved by 3 cents.

• •

Request sizes are distributed according to a Gamma distribution $\Gamma(\alpha_{\mu},\beta_{\mu})$ with $\alpha_{\mu}=4$ and $\beta_{\mu}=4\cdot 10^{-4}$. This corresponds to an average request size of $10000$ assets (i.e. approximately $1000000\ \textrm{€}$) and a standard deviation equal to half the average.

The variance-covariance matrix is therefore given by

[TABLE]

We can see that the second eigenvalue is very small in comparison to the first. This justifies that it is reasonable to approximate the two-dimensional problem with a one-dimensional problem using the result of Section 3, i.e. by considering $\beta\simeq\begin{bmatrix}0.906\\ 0.424\end{bmatrix}$ and $V\simeq 1.744$.

Regarding the objective function, we consider the following:

• •

Time horizon given by $T=12\ \textrm{days}$. This horizon ensures convergence towards stationary quotes at time $t=0$ – see Figures 5 and 20 below.

• •

$\psi:q\in\mathbb{R}^{2}\mapsto\frac{\gamma}{2}q^{\prime}\Sigma q$ with $\gamma=8\cdot 10^{-7}\ \textrm{€}^{-1}$.

• •

$\ell_{d}=0$.

4.1.2 Results with 2 factors

Since $\theta$ and $\tilde{\theta}$ are defined on $[0,T]\times\mathbb{R}^{2}$, a first step for approximating the value functions consists in restricting the state space to a compact set. A traditional way to proceed consists in setting boundary conditions. In what follows, we equivalently impose risk limits in the sense that no trade that would result in an inventory $q\in\mathbb{R}^{2}$ such that $q^{\prime}\Sigma q>B$ is admitted, where $B=2.4\cdot 10^{10}$.888These risk limits can be expressed in the space of factors instead of being expressed with the inventory vector. We then approximate the solution $\tilde{\theta}$ to (4) with two factors using a monotone explicit Euler scheme with linear interpolation on a grid of size $141\times 141$ for the factors and a discretization of the RFQ size distribution with 4 sizes: $z^{1}=6250$, $z^{2}=12500$, $z^{3}=18750$, and $z^{4}=25000$ assets – thereafter respectively designated by very small, small, large, and very large size – with respective probability $p^{1}=0.53$, $p^{2}=0.35$, $p^{3}=0.10$, and $p^{4}=0.02$.999When there are as many factors as assets, one could directly consider the problem with inventory variables.

The value function (at time $t=0$) as a function of the factors is plotted in Figure 1. The corresponding value function as a function of the inventory, obtained through linear interpolation is plotted in Figure 2.

From the value function, we obtain the optimal bid and ask quotes of the market maker as a function of inventories and request size. The optimal bid quotes (at time $t=0$) for asset 1 and asset 2 (in the case of the smallest RFQ size) are plotted in Figures 3 and 4. The ask quotes are of course symmetric and are not plotted.

We see in Figures 3 and 4 that the optimal bid quotes for both assets are increasing functions of both the inventory in asset 1 and asset 2, as expected given the positive value of the correlation parameter $\rho$.

As discussed above, we chose $T=12$ days to ensure convergence of the optimal quotes to their stationary values. This is illustrated in Figure 5.

To see the impact of the RFQ size on the optimal quotes, we plot in Figure 6 the four functions

[TABLE]

in Figure 7 the four functions

[TABLE]

and in Figure 8 the four functions

[TABLE]

where $(q^{1},q_{\text{major axis}}^{2}(q^{1}))$ spans the major axis of the ellipse of authorized inventories.

Likewise for asset 2: we plot in Figure 9 the four functions

[TABLE]

in Figure 10 the four functions

[TABLE]

and in Figure 11 the four functions

[TABLE]

where $(q^{1},q_{\text{major axis}}^{2}(q^{1}))$ spans the major axis of the ellipse of authorized inventories.

Remark 6**.**

*The reader should be aware that the $x$-axis does not span the same range of inventory on all graphs even when the variable is the same. This is due to the elliptic geometry of the set of authorized values for the inventory.

We see that accounting for the size of RFQs significantly impacts the optimal quotes of asset 1. This is less the case for asset 2 (this difference is due to the fact that the volatility of asset 1 is twice that of asset 2). In all cases, the monotonicity is unsurprising.101010Boundary effects related to impossible interpolation explain the surprising position of some extreme points.

In order to check that the value of $B$ defining the risk limits does not have a significant impact on our numerical approximation, we carried out a Monte-Carlo simulation with 2000 trajectories starting from zero inventory, using the optimal quotes. The distribution of inventory is plotted in Figure 12.111111The shades of gray are in logarithmic scale. We clearly see that the ellipse of authorized inventory is wide enough to have little influence on the outcome.

The statistics associated with our simulations are documented in Table 1: the average PnL at time $T$, the standard deviation of that PnL at date $T$, the part of that standard deviation not related to market risk – i.e. only related to the randomness of RFQs –,121212Using the law of total variance, it is easy to see that

$\mathbb{V}(\text{PnL}_{T})=\mathbb{E}\left[\int_{0}^{T}q_{t}^{\prime}\Sigma q_{t}dt\right]+\mathbb{V}\left(\int_{0}^{T}\sum_{i=1}^{2}\int_{\mathbb{R}_{+}^{*}}\left(\delta^{i,b}(t,z)zJ^{i,b}(dt,dz)+\delta^{i,a}(t,z)zJ^{i,a}(dt,dz)\right)\right).$

This formula enables to distinguish the part of the variance of the PnL at time $T$ coming from market risk (the first term) from that coming from the randomness of RFQs (the second term). and the estimated value of the objective function, i.e. the empirical mean of $\text{PnL}_{T}-\frac{\gamma}{2}\int_{0}^{T}q_{t}^{\prime}\Sigma q_{t}dt$.

These figures have to be compared with those associated with a basic naive strategy. A basic strategy consists, for each asset and side, in always quoting the same “myopic” quote that maximizes the expected instantaneous PnL. In other words, these myopic quotes are defined, for all $i\in\{1,\ldots,d\}$, by

[TABLE]

In our case, the myopic quotes are all equal to $0.03854\ \textrm{€}$.

We carried out 2000 simulations using these quotes with the same source of randomness as above. The distribution of inventory is plotted in Figure 13, and the statistics associated with our simulations are documented in Table 2.

With these figures, we clearly see that the main source of risk is market risk and not the risk associated with the randomness of RFQs. This justifies our choice of objective function that only penalizes the part of the variance coming from market risk.

We also clearly see that the use of the optimal quotes drastically reduces the variance of the PnL and results in a high value of the objective function. More precisely, although the use of optimal quotes reduces the average PnL from around $73410$ to around $72081$, it enables a reduction by a factor $3$ of the standard deviation of the PnL from around $265906$ to around $80432$, hence a major increase of the objective function from around $43953$ to around $69293$ (a figure consistent with the maximum of the value function plotted in Figure 1 which is approximately 69174).

4.1.3 Results with the one-factor model and comparison

Let us now compare the results with two factors to the results with one factor, i.e. when the smallest eigenvalue of $\Sigma$ is replaced by [math].

As above, we start with an approximation of the solution $\tilde{\theta}$ to (4) with one factor. We used a monotone explicit Euler scheme on a grid of size $141$ for the unique factor and the same discretization (with 4 sizes) as in the previous paragraphs for the RFQ size distribution. The set of authorized inventory $\{q\in\mathbb{R}^{2}|q^{\prime}\Sigma q\leq B\}$ is of course replaced by the set $\{q\in\mathbb{R}^{2}|(\beta^{\prime}q)^{\prime}V(\beta^{\prime}q)\leq B\}$ which corresponds, in terms of the unique factor, to the interval $\left[-\sqrt{\frac{B}{V}},\sqrt{\frac{B}{V}}\right]$.

The value function (at time $t=0$) as a function of the inventory, obtained through linear interpolation is plotted in Figure 14. The difference between the one-factor value function and the two-factor one is plotted in Figure 15.

We see that the value function in the one-factor case is above that of the two-factor case. This comes from the fact that not all the risk is taken into account in the one-factor case. We also see that the difference between the two value functions is very large at the two extremes of the major axis of the ellipse. This comes from the fact that the market maker using the one-factor model believes that positions close to the major axis of the ellipse are associated with low risk whereas this is less and less the case as the inventory in each asset increases in absolute value.

As in the two-factor case, we deduce from the value function the optimal bid and ask quotes of the market maker (at time $t=0$) as a function of inventory and request size. The optimal bid quotes for asset 1 and asset 2 (in the case of the smallest RFQ size) are plotted in Figures 16 and 17. The differences between the optimal bid quotes in the one-factor case and two-factor case (in the case of the smallest RFQ size) are plotted in Figure 18 for asset 1 and Figure 19 for asset 2.

We clearly see that the larger (in absolute value) the inventory in each asset, the larger the difference in optimal quotes between the exact model and the one-factor approximation. This is in line with expectation.

As in the two-factor case, we checked the convergence of the quotes towards their asymptotic values. We see in Figure 20 that the optimal quotes reach their stationary values in less than $3$ days, far before the $12$ days.

To better compare the quotes and see the impact of RFQ size, we plot in Figures 21, 22, and 23 the optimal bid quotes of asset 1 when $q_{2}=0$ for different values of $q_{1}$, when $q_{1}=0$ for different values of $q_{2}$, and when $(q_{1},q_{2})$ spans the major axis of the ellipse of authorized inventory.

Likewise, we plot in Figures 24, 25, and 26 the optimal bid quotes of asset 2 when $q_{1}=0$ for different values of $q_{2}$, when $q_{2}=0$ for different values of $q_{1}$, and when $(q_{1},q_{2})$ spans the major axis of the ellipse of authorized inventory.

These plots confirm that the optimal quotes in the one-asset model are good approximations of the true optimal ones whenever the inventory in each asset is not too large, all the more for inventories that are not close to the major axis of the ellipse of authorized inventory. Moreover in our example, the one-factor model seems to return quotes closer to the true optimal ones for asset 1 than for asset 2: this is due to the fact that here, the factor obtained through PCA explains better the risk of asset 1 than that of asset 2 (the residual variance of the latter is four times that of the former).

Comparing quotes is important but what really matters is to compare the distribution of the PnL at time $T$ when using the quotes obtained within the one-factor model with the distribution of the PnL at time $T$ when using the optimal quotes (of the two-factor model). For that purpose, we carried out a Monte-Carlo with 2000 simulations using the same source of randomness as in Section 4.1.2. The distribution of inventory when using the optimal quotes in the one-factor model is plotted in Figure 27. The statistics associated with our simulations are documented in Table 3.

We clearly see that the performance of the one-factor approximation is very good. The value of the objective function is indeed around $68567$ when using the quotes obtained in the one-factor model, a bit smaller than the value of approximately $69293$ obtained with the same source of randomness when using the optimal quotes. In fact, the average PnL is higher with the one-factor optimal quotes but, since the distribution of inventory is denser in areas that are falsely believed to be risk-free or wrongfully associated with low risk, the standard deviation of the PnL is in fact higher (around $96746$ versus around $80432$), resulting in a lower value of the objective function.

4.1.4 Taking the residual risk into account with our Monte-Carlo method

We have seen above that the use of the optimal quotes of the one-factor model provides very good results in terms of the value of the objective function. Nevertheless, the distribution of the inventory plotted in Figure 27 differs from the distribution associated with the true optimal quotes plotted in Figure 12 because the major axis of the ellipse of authorized inventory is associated with zero risk in the one-factor model. In this section, we illustrate the Monte-Carlo method proposed in Section 3.2 in order to account for the residual risk in the approximation of the value function and the optimal quotes.

It is noteworthy that the Monte-Carlo method of Section 3.2, unlike the grid method of Section 3.1, does not allow to compute the optimal quotes for all assets, sides, sizes, and values of the inventory (we ignore time by focusing on $t=0$) at once. Instead, it requires a different Monte-Carlo simulation for each desired quote. In particular, should it be used by practitioners, its use should be online. In other words, the computations should only be carried out upon receiving an RFQ or slightly beforehand if one wants to prepare the quotes (given the current inventory) for the most probable RFQs.

A related point is that, even for illustration and even with two assets, it is too time-consuming to compute the Monte-Carlo adjustment for all assets, sides, sizes, and possible inventories. As a consequence, it is too time-consuming to carry out simulations of the PnL with the quotes amended by the Monte-Carlo method of Section 3.2. Instead of a complete analysis, we focus on a sectional analysis by looking at the cases $q^{2}=0$, $q^{1}=0$, and $(q^{1},q^{2})$ on the major axis of the ellipse of authorized inventory.

In Figures 28, 29, and 30, we compare the value function obtained in the two-factor case, i.e. the true value function, to the value function of the one-factor case and to its adjustment through the Monte-Carlo technique of Section 3.2 – we use 50 simulations for each point (with the same source of randomness for all points). Figure 28 deals with the comparison of the values on the section $\{q^{2}=0\}$, Figure 29 deals with the comparison of the values on the section $\{q^{1}=0\}$, and Figure 30 deals with the comparison of the values on the major axis of the ellipse of authorized inventory.

We clearly see that, unsurprisingly, the Monte-Carlo adjustment goes in the right direction. However, the Monte-Carlo method leads to (i) an overestimation of the gap between the value function of the one-factor case and the true value function and (ii) an overestimation of the degree of concavity of the value function (this is particularly the case for the third section).

It is noteworthy that the quality of the approximation is the best around [math]. This point is interesting if one wants to estimate the degree of suboptimality of a quoting strategy in a scenario starting with zero inventory. The poor approximation of the concavity is however a limitation since quotes are based on finite differences of the value function. This is well illustrated by Figures 31, 32, and 33. Nevertheless, even though the Monte-Carlo adjustments of quotes are too large, especially on the major axis of the ellipse (see Figure 33), this drawback of the Monte-Carlo method of Section 3.2 should be qualified as the quotes obtained with the Monte-Carlo technique naturally lead – because of the overestimated slope of the quotes – to trajectories of the inventory more concentrated around [math] and therefore to a very rare use of the quotes that are too different from the optimal ones.

Before we go on with an example including 30 assets, let us conclude on the two-asset case. The method we propose to tackle the curse of dimensionality is based on the projection of market risk on a low-dimensional space of factors. It works very well in the two-asset case as the quotes permit to reach a value of the objective function close to the optimal one. The Monte-Carlo adjustment we suggested in Section 3.2 allows to approximate the true value function at the point of zero inventory, which is quite useful when one does not have access to the true value function, as is the case in high dimension. However, it overestimates the changes one must make to the quotes computed with the low-dimensional approximation.

4.2 Dealing with 30 assets

We now consider the more challenging case of a market maker in charge of 30 assets (here bonds) with the following characteristics:

• •

Asset prices: $S^{i}_{0}=100\ \textrm{€},\quad\forall i\in\{1,\ldots,30\}$.

• •

Volatility of assets: $\sigma^{i}=1.2\ \textrm{€}\cdot\textrm{day}^{-\frac{1}{2}},\forall i\in\{1,\ldots,15\}$, and $\sigma^{i}=0.6\ \textrm{€}\cdot\textrm{day}^{-\frac{1}{2}},\forall i\in\{16,\ldots,30\}$.

• •

Correlation matrix: $\begin{bmatrix}R_{11}&R_{12}\\ R_{21}&R_{22}\\ \end{bmatrix},$ where

[TABLE]

• •

Intensity functions:

[TABLE]

with $\lambda_{RFQ}=10\ \textrm{day}^{-1}$, $\alpha_{\Lambda}=0.7$, and $\beta_{\Lambda}=30\ \textrm{€}^{-1}$. This corresponds to 10 RFQs per day for each asset, a probability of $\frac{1}{1+e^{0.7}}\simeq 33\%$ to trade when the answered quote is the reference price and a probability of $\frac{1}{1+e^{-0.2}}\simeq 55\%$ to trade when the answered quote is the reference price improved by 3 cents.

• •

Request sizes are distributed according to a Gamma distribution $\Gamma(\alpha,\beta)$ with $\alpha=4$ and $\beta=4\cdot 10^{-4}$. This corresponds to an average request size of $10000$ assets (i.e. approximately $1000000\ \textrm{€}$) and a standard deviation equal to half the average.

The variance-covariance matrix $\Sigma$ has two eigenvalues equal to $19.895060$ and $4.584941$, and 28 eigenvalues below $0.15$. The first eigenspace is spanned by a vector with all coordinates of the same sign. The associated factor – the first factor – represents an index of the 30 assets. The second eigenspace is spanned by a vector with the first 15 coordinates of the same sign and the next 15 of the same, but opposite, sign. The associated factor – the second factor – allows to separate the two groups of assets.

We can therefore legitimately approximate our 30-asset problem by a two-factor one and solve the corresponding Hamilton-Jacobi equation (4) to approximate the optimal quotes.

Regarding the objective function, we consider the following:

• •

Time horizon given by $T=2\ \textrm{days}$. This horizon ensures convergence towards stationary quotes at time $t=0$ – see Figure 37.

• •

$\psi:q\in\mathbb{R}^{2}\mapsto\frac{\gamma}{2}q^{\prime}\Sigma q$ with $\gamma=8\cdot 10^{-7}\ \textrm{€}^{-1}$.

• •

$\ell_{d}=0$.

We approximate the solution $\tilde{\theta}$ to (4) with two factors using a monotone explicit Euler scheme with linear interpolation on a grid of size $71\times 71$ for the factors and a discretization of the RFQ size distribution with the same 4 sizes as in the above two-asset example.131313We considered risk limits similar to those of the above two-asset example. Here, no trade that would result in an inventory $q\in\mathbb{R}^{30}$ such that $q^{\prime}\Sigma q>B$ was admitted, where $B=5\cdot 10^{10}$.

The value function (at time $t=0$) is plotted in Figure 34. The associated optimal bid quotes for asset 1 and asset 16 (for the smallest RFQ size) are plotted in Figures 35 and 36. We see on these graphs that the optimal quotes depend monotonously on the two factors.141414Exceptions to this monotonicity property are related to boundary effects.

As discussed above, we chose $T=2$ days to ensure convergence of the quotes to their stationary values. This is illustrated in Figure 37.

To see the role of RFQ size on the optimal quotes, we plot in Figure 38 the four functions $f^{1}\mapsto\bar{\delta}^{1,b}(0,f^{1},0,z^{k}),k\in\{1,\ldots,4\}$ and in Figure 39 the four functions $f^{2}\mapsto\bar{\delta}^{1,b}(0,0,f^{2},z^{k}),k\in\{1,\ldots,4\}.$

Likewise, we plot in Figure 40 the four functions $f^{1}\mapsto\bar{\delta}^{16,b}(0,f^{1},0,z^{k}),k\in\{1,\ldots,4\}$ and in Figure 41 the four functions $f^{2}\mapsto\bar{\delta}^{16,b}(0,0,f^{2},z^{k}),k\in\{1,\ldots,4\}.$

We see, especially in Figures 39 and 40, that the size of the RFQ significantly impacts the quotes that should be answered (as computed with our two-factor approximation).

Unlike what we did in the two-asset case, it is impossible in our 30-asset case to know how far from real optimality are the optimal quotes computed with the two-factor approximation. Nevertheless, we can use Monte-Carlo simulations to estimate the value of the objective function associated with the optimal quotes computed with the two-factor approximation in a scenario starting from zero inventory, and compare such an estimation to an approximation of the value function at $(t,q)=(0,0)$ computed through the Monte-Carlo approximation of Section 3.2.

We carried out 2000 trajectories starting from zero inventory, using the optimal quotes computed with the two-factor approximation. These 2000 simulations enable to illustrate the distribution of the PnL at time $T$. The statistics associated with our simulations are documented in Table 4.

The value of the objective function documented in Table 4 has to be compared to an approximation of the true value function at $(t,q)=(0,0)$. Given the value $\tilde{\theta}(0,0)=60156$ obtained with our numerical scheme and given an estimation of $\eta(0,0)$ equal to $-643$ obtained using the Feynman-Kac representation of Section 3.2 – with 500 trajectories –, we obtain an approximation of the true value function at $(t,q)=(0,0)$ equal to $59513$. From the very small value of $\frac{\eta(0,0)}{\tilde{\theta}(0,0)}$, we deduce that our two-factor approximation is quite satisfactory. The near-optimality of the quotes obtained with our two-factor approximation is confirmed by the value $59765$ obtained with our 2000 trajectories (see Table 4) which is even slightly above $59513$.

Conclusion

In this paper, we generalized existing market making models to introduce trade size variability. This extension led to an integro-differential equation of the Hamilton-Jacobi type that can be solved using ODE techniques in an infinite-dimensional space. Then, we introduced a numerical method for approximating the optimal bid and ask quotes of a market maker over a large set of assets using a dimensionality reduction technique based on a factor decomposition of the risk. To exemplify our findings, and show that they contribute to beating the curse of dimensionality, we considered two cases of market making with respectively 2 and 30 assets. Our method scales linearly in the number of assets and exponentially in the number of factors, and can therefore be used on large markets driven by a few number of factors.

Appendix: On the construction of the processes $J^{i,b}$ and $J^{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 $d=1$ and let us omit the superscript $i$ (the generalization is straightforward). Let us introduce $N^{b}$ and $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 $J^{b}(dt,dz)$ and $J^{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{J}^{b}(dt,dz)$ and $\tilde{J}^{a}(dt,dz)$ stand for the compensated measures.

We then know from Girsanov theorem (in its Brémaud-Jacod version) that under $\tilde{\mathbb{P}}^{\delta}$, $J^{b}(dt,dz)$ and $J^{a}(dt,dz)$ have respective intensity kernels

[TABLE]

as in the body of the paper.

Data Availability Statement

Data sharing is not applicable to this article as no new data were created or analyzed in this study.