Performance Enhancing Market Risk Calculation Through Gaussian Process Regression and Multi-Fidelity Model

1. Introduction

Market risk refers to the risk of losses resulting from adverse movements in market prices, encompassing both on- and off-balance sheet positions. From a regulatory standpoint, market risk emanates from positions held in a bank trading book, as well as from commodity and foreign exchange risk positions within the entire balance sheet. Counterparty credit risk, on the other hand, relates to the risk that a counterparty in a transaction may get default before all the cash flows have been settled. Credit valuation adjustment (CVA) risk pertains to the risk of incurring losses due to fluctuations in CVA values triggered by changes in counterparty credit spreads and market risk factors. Financial risk management has always been a central activity in banking industry, but since the 2008 global financial crisis, it has become even more important: there is an increasing focus on how financial risks are being assessed, monitored and managed. In the last three decades, extensive research, spanning both academic and financial industry domains [13,15,16,48], has been dedicated to advancing banking practices and refining risk calculations to meet existing and forthcoming challenges.

The European Banking Authority (EBA) plays a crucial role in ensuring consistent implementation of the new regulations across the European Union, particularly with respect to the Fundamental Review of the Trading Book (FRTB) as part of the Basel III reforms [3], by developing technical standards, guidelines, and reports. Specifically, the Basel Committee on Banking Supervision introduced the FRTB framework to enhance accuracy and consistency in trading book capital requirements. This initiative has compelled banks with trading operations to swiftly implement significant changes. Banks had to reinforce their models, data management, information technology infrastructure, and processes. Additionally, capital markets businesses had to optimize product structures, adjust hedging strategies, and reprice products. For some banks, unprofitable business lines due to FRTB led to downsizing or exiting their operations. The impact of FRTB on banks is multifaceted, affecting methodologies, data management, processes, and systems. These changes pose significant challenges, particularly in areas like non-modellable risk factors and the development of default risk charge (DRC) methodologies. Simultaneously, the implementation of FRTB necessitates substantial upgrades in banking systems, including computational capacities and alignment of risk factors. The FRTB also emphasizes the importance of market and reference data quality, demanding a coherent framework for sourcing and mapping positions accurately. These interconnected challenges highlight the intricate demands faced by banks in adapting to FRTB regulatory requirements. Banks must increase computational capacities, revalue positions daily, and develop new platforms to compute standardized capital charges: significantly higher calculation capacity will be required compared to that for current value-at-risk (VaR) calculations and banks estimate an increase in number of full revaluation runs by a factor of 5 to 10 [4]. Simultaneously, the banking industry response to these complex demands involves exploring the potential of machine learning, which has gained prominence for big data management, function interpolation and exploration in various sectors. In finance, machine learning applications in, such as pricing, fraud detection, credit scoring, and portfolio management, have become increasingly significant [17,19,36]. This technological advancement offers a potential solution to the intricate challenges posed by the FRTB, offering innovative ways to adapt and streamline banking processes in the face of regulatory requirements.

A recent study by the [19] highlights the potential of machine learning in enhancing risk models by detecting complex patterns in large datasets. The neural network family, famous for its universal approximation and its customization flexibility, is applied to various financial tasks ranging from derivative pricing and hedging [10,29] to calibration [26] and partial differential equations solving [28]. [50] show improvements in the volatility forecast using neural network and GELM, where the latter is a combination of the generalized autoregressive conditional heteroskedasticity (GARCH) model and extreme machine learning (ELM) model [27], this then enhances the accuracy and efficiency in value-at-risk calculations. [46] make a survey on the application of neural network models in quantitative finance. Bayesian models such as the Gaussian process regression can be explored for the yield curve forecast and derivatives pricing [40], reducing computation time. [14] apply GP models to the fast derivatives pricing and hedging under the assumption that the market is modelled by a limited set of parameters. They provide the numerical test of European call, American put and barrier option in the variance gamma and Heston models. [23,38] replace the linear regression model in the [37] model by a Gaussian process regression to price Bermudan options. In particular, at each backward evaluation time step, they train a GP regression to estimate the time value of the option. [12] focus on the approximation of European call in the Heston model, they furthermore study the multi-output Gaussian processes (GPs) structure to simultaneously estimate various option prices underwritten on the same asset. In [34] the market value of the whole trading portfolio of different options on a common underlying is efficiently approximated by a Gaussian process regression. In a nutshell, the choice between neural networks and Gaussian process regression depends on the specific task requirements. The former is The latter is preferred when uncertainty analysis in predictions is crucial, or when there are limitations in the number of training points available. Despite these advancements, machine learning methods have not fully integrated into the financial risk measurement framework yet.

This work aligns with the research flow of [34], which produced encouraging results for the $\mathrm{VaR}$ and $\mathrm{ES}$ calculation of the portfolio consisting of financial instruments written on a single asset. However, one drawback of the approach is the curse of dimensionality. As the input space dimension of the pricing function increases, the number of sampling points required to train the machine learning algorithm, used as a proxy pricing function, increases substantially. This issue becomes practically challenging when evaluating multi-asset options such as best-of, worst-of and basket options. The GP techniques applied to fixed income derivatives like Bermudan swaptions also suffer from the curse of dimensionality. At this stage, the main difficulty lies in the limitation of training data rather than the scalability of Gaussian processes model for learning large datasets (cf. Section 3.2). While the direct use of machine learning model may be efficient enough, some enhancement techniques can be applied to achieve optimal results. For example, when performing the risk calculation of an interest rate derivatives portfolio, [47] apply dimensional reduction techniques, e.g. principal component analysis (PCA) or learning from diffusive risk factors [1], to reduce the number of market risk factors before interpolating the portfolio value by Chebyshev tensors. However, when products are highly nonlinear and depend on many diffusive risk factors, e.g. multi-asset derivatives, the dimensional reduction techniques proposed in [47] may be no longer appropriate. This leads us to explore other enhancement techniques such as multi-fidelity modeling. Fidelity in modeling refers the degree of detail and precision offered by a predictive model or simulation. High-fidelity models provide highly accurate outcomes but require substantial computational resources, for instance, the pricing engine in our financial case. Due to its high cost, the use of high-fidelity models is often limited in practice. By contrast, low-fidelity models, while resulting to less accurate outcomes, are cheaply accessible and often offer some useful insights of the learning problem. By combining models of different fidelity levels, multi-fidelity modeling exploits all the information along these models and their correlation structure to provide quick yet accurate predictions. [18] do a survey of the multi-fidelity modeling. This technique can be used along with any machine learning models, such as neural networks [35,39] or Gaussian process regression [6,30,32]. In the research, we investigate multi-fidelity Gaussian process regression (mGPR). While the latter is well-known and widely used in geostatistics and physics under the name cokriging (see lecture notes in https://geostatisticslessons.com and references therein), to the best of our knowledge, we have not seen any of its application in quantitative finance. This motivates us to study this modeling for the financial application, where high-fidelity model is the derivative pricing engine and lower fidelity models can be either its simpler approximations or fast pricing engines of correlated products.

The rest of the chapter is organized as follows. Section 2 introduces the principle of market risk according to the regulation, highlighting the computational challenge in the practical implementation of VaR and ES model. Section 3 reviews the Gaussian process regression and its application in the risk calculation of the trading portfolios. Section 4 describes multi-fidelity modeling and how it works with GP model, coming along with some ideas of their application in quantitative finance. Section 5 sketches the setup and the training specification for our numerical experiments, which are reported in Section 6. Section 7 concludes our findings and provides perspectives for future research.

2. Market Risk Assessment and Computational Challenge

According to [4], market risk is defined as "the risk of losses (in on- and off-balance sheet positions) arising from movements in market prices". Market risks include default risk, interest rate risk, credit spread risk, equity risk, foreign exchange (FX) risk commodities risk, and so on. Regulatory separates market risk with other kinds of financial risk, for example, credit risks referring to the risk of loss resulting from counterparty (the party with whom the bank makes engagement in the contract) default or operational risk referring the risk of loss resulting from the failures of internal banking system. Under the FRTB framework, banks can choose either the standardized measurement method (SMM) or the internal model-based approach (IMA) to quantify their market risks. The first approach, which can be easily implemented, is not appealing to banks due to its overestimation of the capital required to cover market risks. The second approach more accurately reflects the risk sensitivities and better postulates the economic risk carried by the banking balance sheet. However, to use their own proposed internal model, banks must gain the regulatory approval through a rigorous backtesting procedure.

Under the internal model approach the calculation of capital charge is based on tail risk measures, namely value-at-risk ($\mathrm{VaR}$) and expected shortfall ($\mathrm{ES}$), of the possible loss of the holding portfolio during a certain trading period. $\mathrm{VaR}$ and $\mathrm{ES}$ are measured at typically high confidence levels $\alpha =97.5\%$ (for $\mathrm{ES}$) or $99\%$ (for $\mathrm{VaR}$) in a short liquidity horizon, i.e. $h=1$ day or 10 days. Assume that the value of a considered trading portfolio $V_t$ at date t is a function of risk factors $R{F}_t$, denoted by $p_t\left(R{F}_t\right)$. At date $t+h$, the possible loss of the portfolio (or profit and loss), denoted by $L_t\left(h\right)$, is defined by

$$\begin{array}{cc}\hfill L_t\left(h\right)=-(V_{t+h}-V_t)& =p_t\left(R{F}_t\right)-p_{t+h}\left(R{F}_{t+h}\right)\approx p_t\left(R{F}_t\right)-p_t(R{F}_t+\Delta R{F}_h),\hfill \end{array}$$

(1)

where $V_{t+h}=p_{t+h}\left(R{F}_{t+h}\right)$ is the value of portfolio at date $t+h$ and the approximation is probably the convention used by the bank, considering time horizon h is small, e.g. one day or ten days. By the approximation in (1), $L_t\left(h\right)$ can be interpreted as the potential loss of the portfolio resulted from the shock $\Delta R{F}_h$ of the risk factors $R{F}_t$. The value-at-risk at a confidence level $\alpha$ represents the minimum capital required to cover the market risk, while the expected shortfall means the average loss given the latter exceeds the $\mathrm{VaR}$. These two risk measures are defined by

$$\begin{array}{cc}\hfill \mathrm{VaR}(L_t,h,\alpha )& =min\{q\in \mathbb{R}:\mathrm{P}[L_t\left(h\right)\le q]=\alpha \},\hfill \\ \hfill \mathrm{ES}(L_t,h,\alpha )& =\mathrm{E}\left[L_t\left(h\right)\right|L_t\left(h\right)\ge \mathrm{VaR}(L_t,h,\alpha )]={\displaystyle \frac{1}{1-\alpha }}{\int }_{\alpha }^1\mathrm{VaR}(L_t,h,\gamma )d\gamma .\hfill \end{array}$$

(2)

Statistical methods for calculating (2) include analytical, historical and Monte Carlo approaches ([45] Section 2.2 page 61). According to the the Monte Carlo method, $\mathrm{VaR}(L_t,h,\alpha )$ is the $\alpha$-percentile of the set of possible losses simulated for time horizon h, while $\mathrm{ES}(L_t,h,\alpha )$ is the empirical average of losses above the $\mathrm{VaR}(L_t,h,\alpha )$. Shocks of risk factors at date $t+h$ are generated by synthetic dynamic models calibrated from the historical market data. The market risk calculation procedure by Monte Carlo is summarized in Algorithm 1, see also [25] for a review of risk calculation by the Monte Carlo methods.

Algorithm 1:Calculation of $\mathrm{VaR}$ and $\mathrm{ES}$ by full repricing approach

input:  Calibrated diffusion model, pricing engine $p_t$, an actual vector of risk factors $R{F}_t$, a confidence level $\alpha$, time step h and a large number N output:  estimated $\mathrm{VaR}(L_t,h,\alpha )$ and $\mathrm{ES}(L_t,h,\alpha )$ Compute $p_t\left(R{F}_t\right)$ Simulate N scenarios of shock $\Delta R{F}_h$ using calibrated diffusion model Compute N corresponding prices $p_t(R{F}_t+\Delta R{F}_h)$ Compute N scenarios of losses $L_t\left(h\right)=-\left(p_t(R{F}_t+\Delta R{F}_h)-p_t\left(R{F}_t\right)\right)$ Compute $\mathrm{VaR}$ and $\mathrm{ES}$ by Monte Carlo

Remark 1. 

In the FRTB, the calculation of market risk for a trading book by the internal model-based approach needs to be implemented for various liquidity horizons associated with different risk factors. To be more specific, for a considered portfolio, banks need to classify involved risk factors into five categories of liquidity horizons: 10, 20, 40, 60 and 120 days. The FRTB requires banks to assess market risk for these liquidity horizons by taking into account shocks to risk factors at the same and longer horizons. For example, the risk calculation for 10-day horizon must consider shocks to risk factors from all five categories, while the risk calculation for 20-day horizon needs to consider shocks at 20 days and longer. This framework provides a more comprehensive and accurate assessment of market risk, but requires a complex and time-consuming calculation process.

2.1. Computational Challenge and Applications of Machine Learning

The market risk calculation appears to be doable by Algorithm 1. Nevertheless, we have not discussed yet about the computational complexity of the pricing engine. In banks, the pricing engine of a trading derivative portfolio involves various numerical algorithms depending on the nature of the portfolio. Some products or portfolios can be rapidly priced, such as using analytical formulas or accurately fast approximations. However, when financial models and/or products of the portfolio are complex as usually the case in practice, some heavily computational algorithms need to be referenced including Monte Carlo schemes, tree methods, see for example ( [11] Chapter 6 and 7). The risk calculation using any of the algorithms mentioned above is referred to as the full repricing approach (i.e. full revaluation approach). The latter is the most accurate but very costly (for a portfolio of complex products) as the measurement of the tail risk, i.e. VaR and ES, of the portfolio needs to recall intensively the pricing engine.

In the context of implementing the $\mathrm{VaR}$ and $\mathrm{ES}$ calculation using Monte Carlo simulation, banking institutions tend to avoid direct use of pricers, as it is time consuming and makes them inefficient. Certain banks prefer to use proxies based on sensitivity calculation and second-order Taylor developments to compute the possible losses of derivative portfolios [8], see Section B. Although these approaches offer performance gains, their precision is often questioned. Recently, alternative approaches have emerged that aim to strike a balance between performance and precision in risk calculation, such as value-at-risk and expected shortfall. Among these methods, the works of [12], and [14] explore the application of Gaussian process regression in finance, and [46] discuss the use of neural networks in the same context. The common idea behind these studies is to use machine learning models as interpolation tools, also known as surrogate models, to estimate derivative prices. More recently, [47] use Chebyshev polynomials as interpolators, combined with principal component analysis techniques (PCA), for financial risk calculation, aligning with the same line of though. Figure 1 illustrates the risk calculation process that employs these innovative machine learning approaches, involving only a few complete revaluations of a derivative portfolio.

2.2. Equity Options

In the remainder of this chapter, the objective is to evaluate and manage the market risk associated with a trading portfolio consisting of equity derivatives. These products are written on a stock whose the price at date t is denoted by $S_t$. Among these derivative products, vanilla options, barrier options, and American options on a single underlying asset can be mentioned. In the Black-Scholes model, the values of these options are simply determined by analytical formulas depending on the underlying asset price $S_t$, except the American option values for which a numerical method, e.g. binomial tree, is mandatory. The portfolio may also include multi-underlying options such as basket options, best-of options and worst-off options. For this case, we denote by $S_t=[S_t^1,\cdots S_t^d]$ the prices of d underlying assets that are governed by lognormal processes:

$$\begin{array}{cc}& d{S}_t^j=r{S}_t^{j}dt+{\sigma }^j{S}_t^{j}d{W}_t^j,forj=1\cdots d,\hfill \\ & d{\langle W^j,W^{j^{\prime }}\rangle }_t={\rho }_{j{j}^{\prime }}dt,\hfill \end{array}$$

(3)

where r is a risk-free rate, ${\sigma }^j$ denote the asset volatilities, and $W^j$ are Brownian motions with a correlation matrix $\mathbf{R}={\left({\rho }_{j{j}^{\prime }}\right)}_{j,j^{\prime }\in \{1,\cdots ,d\}}$. In what follows, we assume the constant risk-free rate and volatilities.

Let $\eta \left(S_T\right)$ denote the value of the European multi-asset derivative at the maturity T, also called the payoff. Subject to the above assumptions, the value of the European multi-asset option $V_t$ at date t is the conditional expectation under the risk neutral probability of the discounted payoff, denoted by ${\mathrm{E}}^{*}$,

$$\begin{array}{c}\hfill p_t\left(S_t\right)={\mathrm{E}}^{*}\left[e^{-r(T-t)}\eta \left(S_T\right)\right|S_t].\end{array}$$

(4)

We are particularly interested in four options whose the payoff can be written as a vanilla option payoff against a strike K given in Table 1.

Given (3), the price of geometric average options are fast referenced by the Black-Scholes formula, and the price of basket options can be accurately approximated by the Black-Scholes formula thanks to the convention that the sum of lognormal variables is approximately a lognormal variable [7]. Otherwise, the price of other options in Table 1 are evaluated by Monte Carlo simulation.

3. Gaussian Process Regression for Option Pricing

Gaussian process regression is the canonical method for Bayesian modeling of spatial functions. The method is highly recommended in cases where data is sparse because of its generalization of the Gaussian distribution from finite-dimensional vector spaces to infinite-dimensional functional spaces. Due to this principle, the model is considered as non-parametric. In this section we briefly introduce the Gaussian process regression and its application to model financial derivative products. The more background of GP model can be referenced in [44] and Chapter 15 [41].

In what follows, bold lowercase letters, e.g. $\mathbf{x}$, stand for vector notation, whereas bold uppercase letters, e.g. $\mathbf{X}$, are used for matrix notation.

3.1. Gaussian Processes Regression and Prediction

Let $(\Omega ,\mathcal{A},\mathrm{P})$ be a probability space, consider a random vector $\mathbf{x}:\Omega \to \mathcal{X}\subset {\mathbb{R}}^d$ for a positive integer d corresponding to the dimension of inputs, and a response variable $y\in \mathbb{R}$. In supervised learning, the goal is to learn a target function $f:\mathcal{X}\to \mathbb{R}$ portraying the relation between inputs and the response variable. To this end, we assume the following regression

$$\begin{array}{c}\hfill y_i=f\left({\mathbf{x}}_i\right)+{ϵ}_i,\end{array}$$

(5)

where ${ϵ}_i$ are i.i.d Gaussian white noise with zero mean and variance ${\sigma }^2$. Different from other regression models which usually assume a parameterization form of f, the Gaussian processes place directly a multivariate normal prior on the space of functions. Implicitly, for any set of input points $[{\mathbf{x}}_1,\cdots ,{\mathbf{x}}_N]$ in $\mathcal{X}$, the random vector ${[f\left({\mathbf{x}}_1\right),\cdots ,f\left({\mathbf{x}}_N\right)]}^{\top }$ is multivariate normally distributed. More precisely, given a data set of N observations $\mathcal{D}=(\mathbf{X},\mathbf{y})=\{({\mathbf{x}}_i,y_i)\mid i=1,\cdots ,N\}$, where ${\mathbf{x}}_i$ is taken from set $\mathcal{X}$, $\mathbf{X}$ is the concatenated matrix in which each row is the row-vector of one observation of inputs. Let ${[f\left({\mathbf{x}}_1\right)\cdots f\left({\mathbf{x}}_N\right)]}^{\top }$. In the Gaussian process regression,

$$\begin{array}{c}\hfill \mathcal{N}({\mu }_{\mathbf{X}},K_{\mathbf{X},\mathbf{X}})implies\mathbf{y}\mid \mathbf{X}\sim \mathcal{N}({\mu }_{\mathbf{X}},K_{\mathbf{X},\mathbf{X}}+{\sigma }^2{\mathbf{I}}_N),\end{array}$$

(6)

where ${\mathbf{I}}_N$ is the identity matrix of size N, ${\mu }_{\mathbf{X}}={[\mu \left({\mathbf{x}}_1\right),\cdots ,\mu \left({\mathbf{x}}_N\right)]}^{\top }$ is a prior mean vector parameteried by a mean function $\mu \mathcal{X}\mapsto \mathbb{R}$, and $K_{\mathbf{X},\mathbf{X}}$ denotes the covariance matrix characterized by a kernel function $k:\mathcal{X}\times \mathcal{X}\mapsto \mathbb{R}$ such that $K_{\mathbf{X},\mathbf{X}}={\left[k({\mathbf{x}}_i,{\mathbf{x}}_{i^{\prime }})\right]}_{i,i^{\prime }=1\cdots N}$.

Unless some extra knowledge about the prior mean function, for simplicity one can choose $\mu (·)=0$. Note that the convention is related only to the prior distribution and does not imply that the posterior distribution (the prediction) has zero mean. Similarly, we can impose some knowledge of the target function on the choice of the kernel function. For example, the exponential sine squared kernel takes into account the periodic characteristic of function ([44] pp. 92). The Matern kernel and the squared exponential kernel are the most frequently used stationary kernels ([44] pp. 83-84). The Matern kernel reads

$$\begin{array}{c}\hfill k\left({\mathbf{x}}_i,{\mathbf{x}}_{i^{\prime }}\right)={\displaystyle \frac{1}{\Gamma \left(\nu \right)2^{\nu -1}}}{\left({\displaystyle \frac{\sqrt{2\nu }}{l}}∥{\mathbf{x}}_i-{\mathbf{x}}_{i^{\prime }}∥\right)}^{\nu }{K}_{\nu }\left({\displaystyle \frac{\sqrt{2\nu }}{l}}∥{\mathbf{x}}_i-{\mathbf{x}}_{i^{\prime }}∥\right),\end{array}$$

(7)

where $\Gamma$ is the gamma function and $K_{\nu }$ is the modified Bessel function of the second kind. The smoothness of the approximate function is controllable through a positive parameter $\nu$, i.e. (7) is once differentiable when $\nu =3/2$ and twice differentiable when $\nu =5/2$. As $\nu$ tends to infinity, the Matern kernel (7) becomes infinitely differentiable and converges to the squared exponential kernel. The length-scale parameter $l>0$ can be a d positive vector and each component will be used to standardize a corresponding component of inputs. If the variances by column of the input $\mathbf{X}$ are equal, the vector length-scale parameter will measure the impact of input variables into the predictive conditional variance, yielding an insight about the variable importance ([44] Section 5.1 and 6.6). Other kernel functions can be found in ( [44] Section 4). When $y_i$ in $\mathcal{D}$ is observed directly from the groundtruth function f, i.e. $y_i=f\left({\mathbf{x}}_i\right)$, the white noise ${ϵ}_i$ in (5) and its variance $\sigma$ in above equations will be removed and this refers to the interpolation with noise-free observations. In brief, Gaussian processes provide a lot of flexibility to integrate extra knowledge of the learning problem, apart from data set $\mathcal{D}$, into the prior model hypothesis $p\left(f\right)$, which is particularly useful for financial applications.

For a matrix of $N^{\prime }$ test points, denoted by ${\mathbf{X}}^{*}={[{\mathbf{x}}_1^{*},\cdots ,{\mathbf{x}}_{N^{\prime }}^{*}]}^T$, the joint prior distribution of the response reads

$$\left[\begin{array}{c}\mathbf{y}\end{array}\right]\sim \mathcal{N}\left(\begin{array}{c}{\mu }_{\mathbf{X}}\\ {\mu }_{{\mathbf{X}}^{*}}\end{array},\left[\begin{array}{cc}{K}_{\mathbf{X},\mathbf{X}}+{\sigma }^2{\mathbf{I}}_n& K_{\mathbf{X},{\mathbf{X}}^{*}}\\ K_{{\mathbf{X}}^{*},\mathbf{X}}& K_{{\mathbf{X}}^{*},{\mathbf{X}}^{*}}\end{array}\right]\right),$$

where $*={[f\left({\mathbf{x}}_1^{*}\right),\cdots ,f\left({\mathbf{x}}_{N^{\prime }}^{*}\right)]}^{\top }$, ${\mu }_{{\mathbf{X}}^{*}}={[\mu \left({\mathbf{x}}_1^{*}\right),\cdots ,\mu \left({\mathbf{x}}_{N^{\prime }}^{*}\right)]}^{\top }$; and $K_{{\mathbf{X}}^{*},{\mathbf{X}}^{*}}={\left[k({\mathbf{x}}_i^{*},{\mathbf{x}}_{i^{\prime }}^{*})\right]}_{i,i^{\prime }=1\cdots N^{\prime }}$ and $K_{\mathbf{X},{\mathbf{X}}^{*}}=K_{{\mathbf{X}}^{*},\mathbf{X}}^{\top }$ are respectively covariance matrix of ${\mathbf{X}}^{*}$, and of $\mathbf{X}$ and ${\mathbf{X}}^{*}$. Once the model is learned (see Section 3.2), the predictive (posterior) distribution of * is

$$\begin{array}{c}\hfill *\mid \mathbf{X},\mathbf{y},{\mathbf{X}}^{*}\sim \mathcal{N}(\mathrm{E}[*\mid \mathbf{X},\mathbf{y},{\mathbf{X}}^{*}],\mathrm{Var}[*\mid \mathbf{X},\mathbf{y},{\mathbf{X}}^{*}]),\end{array}$$

(8)

where the posterior mean and variance are given as follows

$$\begin{gathered} \begin{array}{c}\hfill {\mathbf{f}}^{*}:=\mathrm{E}[*\mid \mathbf{X},\mathbf{y},{\mathbf{X}}^{*}]={\mu }_{{\mathbf{X}}^{*}}+K_{{\mathbf{X}}^{*},\mathbf{X}}{\left[K_{\mathbf{X},\mathbf{X}}+{\sigma }^2{\mathbf{I}}_n\right]}^{-1}\mathbf{y}, \\ \mathrm{Var}[*\mid \mathbf{X},\mathbf{y},{\mathbf{X}}^{*}]=K_{{\mathbf{X}}^{*},{\mathbf{X}}^{*}}-K_{{\mathbf{X}}^{*},\mathbf{X}}{\left[K_{\mathbf{X},\mathbf{X}}+{\sigma }^2{\mathbf{I}}_n\right]}^{-1}{K}_{\mathbf{X},{\mathbf{X}}^{*}}.¯\end{array} \end{gathered}$$

(9)

Remark 2. 

(i)As we will discuss in Section 3.2, the repeated computation of the inverse matrix ${\left[K_{\mathbf{X},\mathbf{X}}+{\sigma }^2{\mathbf{I}}_n\right]}^{-1}$ in the learning procedure is a bottleneck of the GP model. However, in the prediction phase (cf. (9)) the matrix inversion is required only one time and it is even already computed during the learning phase. In particular, the posterior mean is numerically computed using the following linear expression, for $i^{\prime }=1\cdots N^{\prime }$,

$$\begin{array}{c}\hfill {\mathbf{f}}_{i^{\prime }}^{*}=\mu \left({\mathbf{x}}_{i^{\prime }}^{*}\right)+{\sum }_{i=1}^N{\omega }_{i}k({\mathbf{x}}_i,{\mathbf{x}}_{i^{\prime }}^{*}),¯\end{array}$$

(10)

$where\omega ={[{\omega }_1,\cdots ,{\omega }_N]}^{\top }={\left[K_{\mathbf{X},\mathbf{X}}+{\sigma }^2{\mathbf{I}}_n\right]}^{-1}\mathbf{y}.$

(ii)The posterior mean can be also written by the following linear combination, without the loss of generality, suppose the zero mean prior,

$$\begin{array}{c}\hfill {\mathbf{f}}_{i^{\prime }}^{*}={\sum }_{i=1}^N{\gamma }_i{y}_i,¯\end{array}$$

(11)

where$\gamma ={[{\gamma }_1,\cdots ,{\gamma }_N]}^{\top }=K_{{\mathbf{X}}^{*},\mathbf{X}}{\left[K_{\mathbf{X},\mathbf{X}}+{\sigma }^2{\mathbf{I}}_n\right]}^{-1}.$ (11) interprets the GP prediction of unction f at a new point, say${\mathbf{x}}_{i^{\prime }}^{*}$ as a linear combination from a set of N observable points${\left(y_i\right)}_{i=1\cdots N}.$ However, the GPs model is not a linear interpolation since the non-linear term is incorporated in y_i. Moreover, in the general case where exact observations are not accessible, the model will interpolate from the noisy observations taking into account the variance of noise. Hence GP model is considered as a powerful non-linear interpolation tool

3.2. Estimation of Model Parameters

The learning procedure of Gaussian process model involves finding the best suited parameters of the mean function ${\mu }_{\mathbf{X}}$, the covariance kernel function k, jointly denoted by θ and the white noise variance σ. The set of GP parameters can be estimated either by maximum likelihood or maximum a posteriori, or cross-validation method (see ( [44] Section 5.4) and ([41] Section 15.2.4)). By mean of computational facility, maximum likelihood is the the most used among these and this method is chosen to be a default method for the standard Gaussian processes. The criterion that we aim to maximize is the marginal likelihood, also called by the (model) evidence reads

$$\begin{array}{c}\hfill (\mu ,\sigma ,\theta )=log\mathrm{P}\left(\mathbf{y}\right|\mathbf{X})=-\left[{(\mathbf{y}-{\mu }_{\mathbf{X}})}^{\top }{\left(K_{\mathbf{X},\mathbf{X}}+{\sigma }^2{\mathbf{I}}_N\right)}^{-1}(\mathbf{y}-{\mu }_{\mathbf{X}})\right.\\ \hfill \left.+logdet\left(K_{\mathbf{X},\mathbf{X}}+{\sigma }^2{\mathbf{I}}_N\right)\right]+constant.\end{array}$$

Notice that the first term of (12) correspond to the data-fit can be seen as the the negative squared loss with a kernel-based weight in linear regression whereas the second term relates to the penalty of the model complexity (see [44] Section 2.4 and 5.4.1).

Maximizing the evidence (or minimizing the negative evidence) is usually solved by a global optimization method, namely BFGS ([42] Chapter 6), as in Scikit-learn Python package or by any gradient-based method such as Adam [31] in Gpytorch Python package. In any case, the numerical implementation must be done by full batch and involves reversing the covariance matrix of size $N\times N$ whose the complexity is $\mathcal{O}\left(N^3\right)$ (using Cholesky decomposition). Contrarily, the prediction step (cf. (9)) costs only $\mathcal{O}\left(N^2\right)$. Hence, even Gaussian process model is a great interpolation tool, the classical model is not appropriate for high-dimensional problem that requires an important amount of training points N to explore the target surface. Many innovative approaches are proposed to tackle this issue and in general they are developed through the use of low-rank approximation (or possibly decomposition) of the covariance matrix [22,49] and/or the use of other learning methods such as (stochastic) variational inference [5,24,49]. We refer the reader to the documentation of Gpytorch package on Python [21] for the detail as well as numerical implementation of these approaches.

3.3. Application to Derivative Portfolio Valuation

The application of Gaussian process regression is not restricted to evaluating mono-asset derivatives as illustrated below, but also applies to multi-asset derivatives in the same manner. However, when the number of underlying assets is high, one needs to use Monte Carlo sampling instead of the grid point sampling technique used in low-dimensional problems. The regression model can directly approximate the value of the entire portfolio without any additional computational cost for learning and prediction, regardless of the portfolio size and complexity. This is hence favorable for risk calculation applications ([9] Remark 3).

When the considered portfolio $V_t$ consists of a possibly large number of mono-asset derivatives, one can still use efficiently GP regression by the decomposition technique provided in [34]. More precisely, the portfolio may depend on a vector of d assets (risk factors) $S=[S^1,\cdots ,S^d]$ but each of its derivatives is underwritten on single stock. The portfolio price $V_t$ is the conditional expectation of the total derivative payoffs given the stock asset vector $S_t$, hence, a function of $S_t$. By the nature of the considered portfolio, we are able to categorize $V_t$ by risk factor

$$\begin{array}{c}\hfill V_t=\sum _{j=1}^d{V}_t^j\left(S_t^j\right),\end{array}$$

(12)

where $V_t^j$ is the sub-portfolio corresponding to the j-th stock. We build the approximate of each $V_t^j$ by one-dimensional Gaussian process model. To this end, for $j=1,\cdots ,d$, we sample N equidistant training points of $S_t^j$, say ${\mathbf{s}}^j=[s_1^j,\cdots ,s_N^j]$, in a determined range $[{\underline{s}}^j,{\overline{s}}^j]$, and the corresponding price ${\mathbf{v}}^j=[V_t^j\left(s_1^j\right),\cdots ,V_t^j\left(s_N^j\right)]$. Since this is the one-dimensional interpolation case and the price functions in finance are usually well-behaved, we only need to recall the pricing engine few times for sampling training data. For example, in the mono-asset portfolio experiments [34], the GPs can reach the good results using only $N=10$. Once the learning is done, the GP price of a sub-portfolio at a new point $s^{*}$ is the predictive mean of the posterior distribution

$$\begin{array}{c}\hfill \mathrm{E}[V_t^j\left(s^{*}\right)\mid {\mathbf{s}}^j,{\mathbf{v}}^j,s^{*}]=\sum _{i=1}^N{\gamma }_j{V}_t^j\left(s_i^j\right),\end{array}$$

(13)

where ${\gamma }_j$ is determined in Remark 2 (ii) (with the adaptation of notations). The GP price of the whole portfolio is then deduced by (13). The market risk calculation is straightforward by using the GP pricing model instead of the pricing engine in the Monte Carlo method (cf. Algorithm 1).

Remark 3.(i)For a derivative which depends on several market risk factors such as asset volatility, dividend and interest rate and so on, the bank sometimes uses the Black model for the pricing engine. In this case, we can train bring it back to the two-dimensional interpolation problem, i.e. the discounted future price and the volatility, and efficiently apply GP.

(ii)In the case of basket options, one can use the current price of ${\sum }_{j=1}^d{\alpha }_j{S}_t^j$ (cf. Table 1) as the unique learning feature, leading to one-dimension Gaussian process regression.

(ii)When the derivative depends on many risk factors, i.e. multi-asset products. The algorithm may require a large number of training data for a good approximation. In our financial application, the time to get pricing training points (sampled by costly pricing engine) is more relevant than the model learning difficulty of GPR mentioned in Section 3.2.

4. Multi-fidelity Gaussian Processes Regression

In the previous section, we describe how to apply GPR to interpolate the surface of price function using a few number of price points. However, the number of required points increases dramatically with the number of risk factors. Suppose that we have access to other kinds of information, e.g. a weak pricing engine based on a lower number of Monte Carlo simulations or using few discretization time steps in the tree model, the question arises how to take into benefit of this knowledge in modeling. For this purpose, we investigate multi-fidelity Gaussian process regression (mGPR) that allows us to learn the ground truth target from several information sources, each with a different level of fidelity. Any detail of this model can be found in [6,30,32].

Assuming that we have s data sets, denoted by ${\mathcal{D}}_1=({\mathbf{X}}_1,{\mathbf{y}}_1),\cdots ,{\mathcal{D}}_s=({\mathbf{X}}_s,{\mathbf{y}}_s)$, sorted by increasing level of fidelity (i.e. ${\mathcal{D}}_s$ is the most accurate and expensive data whereas ${\mathcal{D}}_1$ is the most simplified and abundant data). In [30], the autoregressive model is used to model the cross-correlation structure between different levels:

$$\begin{array}{cc}\hfill f_j\left(\mathbf{x}\right)& ={\rho }_{j-1}{f}_{j-1}\left(\mathbf{x}\right)+{\xi }_j\left(\mathbf{x}\right),j=2\cdots s,\hfill \\ \hfill {\mathbf{y}}_1{]}_i& =f_1\left({\mathbf{x}}_i\right)+ϵ,{\mathbf{x}}_i\in {\mathbf{X}}_1,\hfill \end{array}$$

(14)

where ${\rho }_{j-1}$ are scalar correlation parameters between ${\mathbf{y}}_{j-1}$ and ${\mathbf{y}}_j$, ${\xi }_j\left(\mathbf{x}\right)$ measures the mismatch between j and ${\rho }_{j-1}{f}_{j-1}\left(\mathbf{x}\right)$ and ϵ is the while noise of the lowest fidelity model. In mGPR, ${\xi }_j$ is supposed independent of $f_{j-1},\cdots ,f_1$. The goal of the model is to lean $f_s$ which is the relationship between the input and the highest fidelity response.

4.1. Two-fidelity Gaussian Process Regression Model

In this section, we focus on the mGPR with two fidelity levels. Let assume a data set of $N_l$ low-fidelity points, denoted by ${\mathcal{D}}_l=({\mathbf{X}}_l,{\mathbf{y}}_l)={({\mathbf{x}}_i,y_{l_i})}_{i=1\cdots N_l}$ and a second set of $N_h(\ll N_l)$ high-fidelity points, denoted by ${\mathcal{D}}_h=({\mathbf{X}}_h,{\mathbf{y}}_h)={({\mathbf{x}}_i,y_{h_i})}_{i=N_l-N_h+1\cdots N_l}$, such that $N_h$ last points of low-fidelity data coincides with $N_h$ high-fidelity data, i.e.

$$\begin{array}{c}\hfill {\mathbf{X}}_l=\left(\begin{array}{c}{\mathbf{x}}_1^{\top }\\ ⋮\\ {\mathbf{x}}_{N_l-N_h+1}^{\top }\\ ⋮\\ {\mathbf{x}}_{N_l}^{\top }\end{array}\right)and{\mathbf{X}}_h=\left(\begin{array}{c}{\mathbf{x}}_{N_l-N_h+1}^{\top }\\ ⋮\\ {\mathbf{x}}_{N_l}^{\top }\end{array}\right).\end{array}$$

(15)

The multi-fidelity modeling involves the following regression functions

$$\left\{\begin{array}{c}{f}_h\left(\mathbf{x}\right)=\rho f_l\left(\mathbf{x}\right)+\xi \left(\mathbf{x}\right),\mathbf{x}\in \mathcal{X},\hfill \\ f_l\left(\mathbf{x}\right)\perp \xi \left(\mathbf{x}\right)\mid \mathbf{x},\hfill \end{array}\right.$$

(16)

where the parameter ρ describes the correlation of the high- and low-fidelity data, $f_l$ is a Gaussian process model of ${\mathbf{y}}_l$ against ${\mathbf{X}}_l$ and δ is the second GP model conditionally independent with the former. In what follows, we consider GP models with constant prior mean [20,30] parameterized by coefficient β, Marten kernel covariance function $k_{\theta }$ with $\nu =2.5$ and parameterized by θ as per (7), and a multiplicative form of the marginal variance ${\sigma }^2$. In particular, the prior model is

$$\begin{array}{c}\hfill l\sim \mathcal{N}({\beta }_l,{\sigma }_l^2{K}_{{\theta }_l}({\mathbf{X}}_l,{\mathbf{X}}_l)),\end{array}$$

(17)

and

$$\begin{array}{c}\hfill \xi \sim \mathcal{N}({\beta }_h,{\sigma }_h^2{K}_{{\theta }_h}({\mathbf{X}}_h,{\mathbf{X}}_h)),\end{array}$$

(18)

where l and ξ are mean vectors ${[f_l\left({\mathbf{x}}_1\right)\cdots f_l\left({\mathbf{x}}_{N_l}\right)]}^{\top }$ and ${[\xi \left({\mathbf{x}}_{N_l-N_h+1}\right)\cdots \xi \left({\mathbf{x}}_{N_l}\right)]}^{\top }$, and the covariance matrices $K_{{\theta }_{l\left(h\right)}}({\mathbf{X}}_{l\left(h\right)},{\mathbf{X}}_{l\left(h\right)})={\left[k_{{\theta }_{l\left(h\right)}}({\mathbf{x}}_i,{\mathbf{x}}_{i^{\prime }})\right]}_{{\mathbf{x}}_i,{\mathbf{x}}_{i^{\prime }}\in {\mathbf{X}}_{l\left(h\right)}}$. Hence the multi-fidelity model is characterized by the parameters of the first and second Gaussian processes: $({\beta }_l,{\sigma }_l,{\theta }_l)$ and $({\beta }_h,{\sigma }_h,{\theta }_h)$, and the correlation parameter ρ.

4.2. Illustrative Application of Multi-Fidelity Model in Option Pricing

For the pricing and risk calculation problem, our proposed approach is to learn the pricing engine, cf. the high fidelity model, with the less precise price values, cf. lower fidelity data. For instance, if the pricing engine involves Monte Carlo simulation of $100,000$ paths, we can sample weak price values calculated by only, e.g. 100 paths as low-fidelity data (see Section 6.2 for numerical experiments). In the American option pricing case, where the accurate price value is calculated using a binomial tree with 100 time steps, the weaker price can be generated by a sparser tree with only, e.g. 10 time steps. Conversely, if the target is the American option, which is expensively computed by tree model, its European counterpart, which is quickly evaluated using the Black-Scholes formula, can serve as lower fidelity data.

Remark 4.[2] propose an analytic approximation of an American option which is equal to the value of its European counterpart plus an add-on or their intrinsic values with a switch condition detailed in Section A.

Figure 2 and Table 2 demonstrate the numerical test of an American put of strike 62. The configuration of this experiment is detailed in Section 5 and Section 6.1. In particular, the GPR is trained with only 5 American put prices, while the mGPR integrates additionally 10 European prices as low-fidelity data. All training points are drawn from a regular grid of the underlying asset within an interval$[1,140]$, and test points are uniformly sampled within the same interval. We observe a significant mismatch of the GPR estimate (cf. blue curve), notably in the in-the-money region, whereas the mGPR estimate (cf. yellow curve) satisfactorily fits the American put price computed using a binomial tree with 100 time steps. Additionally, the mGPR provides smaller prediction uncertainty represented by the light yellow interval, which closely aligns with the yellow curve and is visually indistinguishable in Figure 2. For further comparison, we also implement a second standard GPR, learning the discrepancy between American and European prices, and the approximation introduced in [2] (see Section A). Their results are then reported in Table 2.

4.3. Conditional Distribution of the Estimate

We describe now the posterior mean and variance of the prediction once the model is completely learned. Following the demonstration in Appendix A [20], for a new point ${\mathbf{x}}^{*}\in \mathcal{X}$, the prediction of the high-fidelity model follows normal distribution:

$$\begin{array}{c}\hfill f_h\left({\mathbf{x}}^{*}\right)\mid l={\mathbf{y}}_l,h={\mathbf{y}}_h,({\beta }_l,{\sigma }_l,{\theta }_l),(\rho ,{\beta }_h,{\sigma }_h,{\theta }_h)\sim \mathcal{N}(m_h\left({\mathbf{x}}^{*}\right),s_h^2\left({\mathbf{x}}^{*}\right)),\end{array}$$

(19)

with mean and variance functions:

$$\begin{array}{cc}\hfill m_h\left({\mathbf{x}}^{*}\right)& =\beta +K{\left({\mathbf{x}}^{*}\right)}^{\top }{\mathbf{K}}^{-1}(\mathbf{y}-\beta ),\hfill \\ \hfill s_h^2\left({\mathbf{x}}^{*}\right)& ={\rho }^2{\sigma }_l^2+{\sigma }_h^2-K{\left({\mathbf{x}}^{*}\right)}^{\top }{\mathbf{K}}^{-1}K\left({\mathbf{x}}^{*}\right),\hfill \end{array}$$

(20)

where

$$\begin{array}{c}\hfill \beta ={\displaystyle \frac{{\mathbf{1}}^{\top }{\mathbf{K}}^{-1}\mathbf{y}}{{\mathbf{1}}^{\top }{\mathbf{K}}^{-1}\mathbf{1}}},with\mathbf{y}=\left[\begin{array}{c}{\mathbf{y}}_l\\ ⋮\\ {\mathbf{y}}_h\end{array}\right].\end{array}$$

(21)

In addition, $K{\left({\mathbf{x}}^{*}\right)}^{\top }$ is the covariance function between $f_h\left({\mathbf{x}}^{*}\right)$ and $[f_h\left({\mathbf{x}}_1\right)\cdots f_h\left({\mathbf{x}}_{N_l}\right),f_h\left({\mathbf{X}}_h\right)]$

$$\begin{array}{c}\hfill K{\left({\mathbf{x}}^{*}\right)}^{\top }=\left[\rho {\sigma }_l^2{k}_{{\theta }_l}({\mathbf{x}}^{*},{\mathbf{X}}_l){\rho }^2{\sigma }_l^2{k}_{{\theta }_l}({\mathbf{x}}^{*},{\mathbf{X}}_h)+{\sigma }_h^2{k}_{{\theta }_h}({\mathbf{x}}^{*},{\mathbf{X}}_h)\right],\end{array}$$

(22)

and the covariance matrix $\mathbf{K}$ of vector $(l,h)$ reads

$$\begin{array}{c}\hfill \mathbf{K}=\left(\begin{array}{cc}{\sigma }_l^2{K}_{{\theta }_l}({\mathbf{X}}_l,{\mathbf{X}}_l)& \rho {\sigma }_l^2{K}_{{\theta }_l}({\mathbf{X}}_l,{\mathbf{X}}_h)\\ \rho {\sigma }_l^2{K}_{{\theta }_l}({\mathbf{X}}_h,{\mathbf{X}}_l)& {\rho }^2{\sigma }_l^2{K}_{{\theta }_l}(X_h,X_h)+{\sigma }_h^2{K}_{{\theta }_h}({\mathbf{X}}_h,{\mathbf{X}}_h)\end{array}\right).\end{array}$$

(23)

We can notice that the computation of (21) involves reversing the massive covariance matrix $\mathbf{K}$ of size $(N_l+N_h)\times (N_l+N_h)$. To reduce the computational complexity, one can apply the Schur complement of block matrix $\mathbf{K}$ which is thus presented by Proposition 3.1 in [32].

prop 4.1.By sorting the input data such that ${\mathbf{X}}_l=({\mathbf{X}}_l\setminus {\mathbf{X}}_h,{\mathbf{X}}_h)$ (cf. (16)) and with all above notations, the inverse matrix of $\mathbf{K}$ in (24) has the form:

$$\begin{array}{c}\hfill {\mathbf{K}}^{-1}=\left(\begin{array}{cc}{\sigma }_l^{-2}{K}_{{\theta }_l}{({\mathbf{X}}_l,{\mathbf{X}}_l)}^{-1}+\left(\begin{array}{cc}0& 0\\ 0& {\rho }^2{\sigma }_h^{-2}{K}_{{\theta }_h}{({\mathbf{X}}_h,{\mathbf{X}}_h)}^{-1}\end{array}\right)& -\left(\begin{array}{c}0\\ \rho {\sigma }_h^{-2}{K}_{{\theta }_h}{({\mathbf{X}}_h,{\mathbf{X}}_h)}^{-1}\end{array}\right)\\ -\left(\begin{array}{cc}0& \rho {\sigma }_h^{-2}{K}_{{\theta }_h}{({\mathbf{X}}_h,{\mathbf{X}}_h)}^{-1}\end{array}\right)& {\sigma }_h^{-2}{K}_{{\theta }_h}{({\mathbf{X}}_h,{\mathbf{X}}_h)}^{-1}\end{array}\right).\end{array}$$

(24)

Instead of performing the inversion matrix in (24) which has a complexity of $\mathcal{O}\left({(N_l+N_h)}^3\right)$, the computation of ${\mathbf{K}}^{-1}$ in (25), involving the inverse of $K_{{\theta }_l}{({\mathbf{X}}_l,{\mathbf{X}}_l)}^{-1}$ and $K_{{\theta }_h}{({\mathbf{X}}_h,{\mathbf{X}}_h)}^{-1}$, has a much lower complexity of $\mathcal{O}(N_l^3+N_h^3)$. By substituting (25) into (21) and doing some calculus, one can simplify the above formulas as follows

$$\begin{array}{c}\hfill \beta ={\displaystyle \frac{{\sigma }_l^{-2}{\mathbf{1}}^{\top }{K}_{{\theta }_l}{({\mathbf{X}}_l,{\mathbf{X}}_l)}^{-1}{\mathbf{y}}_l+(1-\rho ){\sigma }_h^{-2}{K}_{{\theta }_h}{({\mathbf{X}}_h,{\mathbf{X}}_h)}^{-1}({\mathbf{y}}_h-\rho {[y_{l_{N_l-N_h+1}}\cdots y_{l_{N_l}}]}^T)}{{\sigma }_l^{-2}{\mathbf{1}}^{\top }{K}_{{\theta }_l}{({\mathbf{X}}_l,{\mathbf{X}}_l)}^{-1}\mathbf{1}+{(1-\rho )}^2{\sigma }_h^{-2}{\mathbf{1}}^{\top }{K}_{{\theta }_h}{({\mathbf{X}}_h,{\mathbf{X}}_h)}^{-1}\mathbf{1}}},\end{array}$$

(25)

yielding the posterior mean and variance

$$\begin{array}{cc}\hfill m_h\left({\mathbf{x}}^{*}\right)& =\beta +\rho K_{{\theta }_l}({\mathbf{x}}^{*},{\mathbf{X}}_l)K_{{\theta }_l}{({\mathbf{X}}_l,{\mathbf{X}}_l)}^{-1}({\mathbf{y}}_l-\beta )+K_{{\theta }_h}({\mathbf{x}}^{*},{\mathbf{X}}_h)K_{{\theta }_h}{({\mathbf{X}}_h,{\mathbf{X}}_h)}^{-1}({\mathbf{y}}_h-\beta ),\hfill \\ \hfill s_h\left({\mathbf{x}}^{*}\right)& ={\sigma }_h^2+\rho {\sigma }_l^2-{\sigma }_h^2{K}_{{\theta }_h}({\mathbf{x}}^{*},{\mathbf{X}}_h)K_{{\theta }_h}{({\mathbf{X}}_h,{\mathbf{X}}_h)}^{-1}{K}_{{\theta }_h}({\mathbf{X}}_h,{\mathbf{x}}^{*})\hfill \\ & -{\rho }^2{\sigma }_l^2{K}_{{\theta }_l}({\mathbf{x}}^{*},{\mathbf{X}}_l)K_{{\theta }_l}{({\mathbf{X}}_l,{\mathbf{X}}_l)}^{-1}{K}_{{\theta }_l}({\mathbf{X}}_l,{\mathbf{x}}^{*}).\hfill \end{array}$$

(26)

4.4. Bayesian Estimation of Model Parameters

Due to the conditional independence assumption in the second line of (17), the first set $({\beta }_l,{\sigma }_l,{\theta }_l)$ can be separately learnt by, for example, the maximum likelihood method using data set ${\mathcal{D}}_l$. While the parameters of the second Gaussian processes $({\beta }_h,{\sigma }_h,{\theta }_h)$ and the correlation parameter ρ need to be jointly estimated to carry out the idea of formulation (17) ([32] Section 3.3.2). Similarly to (12), the parameters of the low-fidelity model are estimated by maximizing the following log-likelihood

$$\begin{array}{c}\hfill -{\displaystyle \frac{1}{2}}\left[N_{l}log\left({\sigma }_l^2\right)+log(det\left(K_{{\theta }_l}({\mathbf{X}}_l,{\mathbf{X}}_l)\right)+{\displaystyle \frac{{({\mathbf{y}}_l-{\beta }_l)}^{\top }{K}_{{\theta }_l}{({\mathbf{X}}_l,{\mathbf{X}}_l)}^{-1}({\mathbf{y}}_l-{\beta }_l)}{{\sigma }_l^2}}\right]\end{array}$$

(27)

Differentiating the above equation w.r.t. ${\beta }_l$ and ${\sigma }_l^2$ and solving it at 0 yield

$$\begin{array}{c}\hfill {\widehat{\beta }}_l={\displaystyle \frac{{\mathbf{1}}^{\top }{K}_{{\theta }_l}{({\mathbf{X}}_l,{\mathbf{X}}_l)}^{-1}{\mathbf{y}}_l}{{\mathbf{1}}^{\top }{K}_{{\theta }_l}{({\mathbf{X}}_l,{\mathbf{X}}_l)}^{-1}\mathbf{1}}},\end{array}$$

(28)

and

$$\begin{array}{c}\hfill {\widehat{\sigma }}_l^2={\displaystyle \frac{1}{N_l}}{({\mathbf{y}}_l-{\widehat{\beta }}_l)}^{\top }{K}_{{\theta }_l}{({\mathbf{X}}_l,{\mathbf{X}}_l)}^{-1}({\mathbf{y}}_l-{\widehat{\beta }}_l).\end{array}$$

(29)

By substituting (29) and (30) into (28), we have the following optimization problem

$$\begin{array}{c}\hfill {\widehat{\theta }}_l=\underset{{\theta }_l}{\arg \min }\left[N_{l}log\left({\widehat{\sigma }}_l^2\right)+log(det\left(K_{{\theta }_l}({\mathbf{X}}_l,{\mathbf{X}}_l)\right))\right].\end{array}$$

(30)

To learn the second Gaussian process model, we define

$$\begin{array}{c}\hfill \mathbf{d}={\mathbf{y}}_h-\rho \left(\begin{array}{c}{y}_{l_{N_l-N_h+1}}\\ ⋮\\ y_{l_{N_l}}\end{array}\right).\end{array}$$

(31)

The log-likelihood of the model reads

$$\begin{array}{c}\hfill -{\displaystyle \frac{1}{2}}\left[N_{h}log\left({\sigma }_h^2\right)+log(det\left(K_{{\theta }_h}({\mathbf{X}}_h,{\mathbf{X}}_h)\right)+{\displaystyle \frac{{(\mathbf{d}-{\beta }_h)}^{\top }{K}_{{\theta }_h}{({\mathbf{X}}_h,{\mathbf{X}}_h)}^{-1}(\mathbf{d}-{\beta }_h)}{{\sigma }_h^2}}\right],\end{array}$$

(32)

yieling the MLEs of

$$\begin{array}{c}\hfill {\widehat{\beta }}_h={\displaystyle \frac{{\mathbf{1}}^{\top }{K}_{{\theta }_h}{({\mathbf{X}}_h,{\mathbf{X}}_h)}^{-1}\mathbf{d}}{{\mathbf{1}}^{\top }{K}_{{\theta }_h}{({\mathbf{X}}_h,{\mathbf{X}}_h)}^{-1}\mathbf{1}}},and\end{array}$$

(33)

$$\begin{array}{c}\hfill {\widehat{\sigma }}_h^2={\displaystyle \frac{1}{N_l}}{(\mathbf{d}-{\widehat{\beta }}_h)}^{\top }{K}_{{\theta }_h}{({\mathbf{X}}_h,{\mathbf{X}}_h)}^{-1}(\mathbf{d}-{\widehat{\beta }}_h).\end{array}$$

(34)

To compute the above parameters, one needs to estimate ${\widehat{\theta }}_h,\rho$ by solving the following optimization

$$\begin{array}{c}\hfill {\widehat{\theta }}_h,\rho =\underset{{\theta }_h,\rho }{\arg \min }\left[N_{h}log\left({\widehat{\sigma }}_h^2\right)+log(det\left(K_{{\theta }_h}({\mathbf{X}}_h,{\mathbf{X}}_h)\right))\right].\end{array}$$

(35)

Note that (36) can be directly solved by numerical methods mentioned in Section 3.2 since the inversion of the small matrix of size $N_h\times N_h$ is not problematic. While the same procedure may not be applied for system (31) when the number of low-fidelity $N_l$ is large. We then present in the next subsection the sparse Gaussian process model for low-fidelity model.

5. Experiment Design and Model Specification

In this section, we present numerical tests dedicated to evaluating the accuracy and performance of the algorithms studied earlier in this chapter. More specifically, these experiments are conducted on two trading portfolios (cf. Section 6.1 and Section 6.3), interspersed by an examination of individual multi-asset calls (cf. Section 6.2). The first one concerns a portfolio of single underlying derivatives, including plain vanilla (call and put), barrier and American options. The second one consists of multi-asset options, such as geometric, basket, best-of- and worst-of- options. The comparison is based on benchmark methods and scores explained below. We conclude this section by specifying the model of the implemented algorithms.

First, we analyze the case of the single-asset options portfolio described in [34]. This portfolio, with a time maturity of ten years, consists of 100 distinct calls and puts issued on each of the four considered underlying assets. In particular, each asset is used to write 10 vanilla, 10 barrier and 5 American options with their strikes randomly generated from 60 to 140. More details about market and product parameters in the portfolio are presented in Tables 1 and A.1 [34]. We then evaluate the price of each multi-asset option, as indicated in Table 1. Each option expires after two years and is based on the same five assets. The actual values of geometric call options can be calculated using the Black-Scholes formula and serve as a reference. In the case of high dimensionality, the considered portfolio comprises 500 derivatives, depending on the performance of 20 stocks, generating a total of 20 market risk factors. Each derivative is characterized by various parameters, including a random selection of underlying assets from 2 to 20, a strike price ranging from 80 to 160, and a number of options from 100 to 500. It is noteworthy that each derivative can be either a call or a put option among the four types introduced in Section 2.2. The portfolio maturity is set at 5 years, meaning that all its derivatives have a maturity ranging from 6 months to 5 years. In the Black-Scholes model framework, we choose a risk-free rate of $1\%$. Initial asset prices are arbitrarily chosen within a range of 100 to 120, and the volatility of these assets is set to vary between 0.2 and 0.5 with a random correlation matrix, refer to Table 3 and Figure 3. Finally, we examine the pricing and risk calculation problem for different time horizons, namely one day, ten days, and one month. At longer horizons, the magnitude of shocks and their volatilities increase, leading to a significant challenge in adapting the number of interpolation points. Financially, the price function becomes less linear, and the tail of the potential loss distribution becomes heavier.

5.1. Benchmarking

For the numerical implementation of market risk metrics calculation, the full pricing approach uses $N=100,000$ paths in any case. This level is chosen to ensure a satisfactory level of convergence of high quantiles estimation. However, a limited number ($N_{train}$) of these prices is chosen to train the machine learning models.

Regarding the exact pricing engine, the value of vanilla and barrier options in the mono-asset derivatives portfolio are calculated by analytical formulas, while the binomial tree with 100 time steps is used to evaluate American options. For multi-asset options, the exact pricing engine is defined by a Monte Carlo simulation of $100,000$ paths, except for geometric options whose the price can be computed by the Black-Scholes formula.

We implement the sensitivity-based pricing approximation (SxS) in Section B. To get sensitivities data, we shock relatively $0.1\%$ of each risk factor and use the central finite difference, yielding totally 40 portfolio evaluations for the first order sensitivity-based approach. The second order sensitivities are numerical unstable and lead to a worse approximation of the price. Hence they are not considered in our numerics.

Regarding machine learning methods, we implement a two-hidden-layer neural network (see Section C) (NN) and two Gaussian process models introduced in Section 3 (GPR) and Section 4 (mGPR). In the latter, the high-fidelity data consists of the $100,000$ Monte Carlo paths prices, while the low-fidelity comprises weak prices based on 100 MC paths in Section 6.2 and 500 MC paths in Section 6.3. For measuring the pricing mismatch, we use the mean absolute percentage error (MAPE) over $N=100,000$ scenarios in the $\mathrm{VaR}$ and $\mathrm{ES}$ computation:

$$\begin{array}{c}\hfill \mathrm{MAPE}={\displaystyle \frac{1}{N}}\sum _{i=1}^N{\displaystyle \frac{|{\widehat{p}}_0\left(S_h\right)-p_0\left(S_h\right)|}{p_0\left(S_0\right)}},\end{array}$$

(36)

where $p_0\left(S_h\right)$ and ${\widehat{p}}_0\left(S_h\right)$ are respectively actual price (i.e. analytical or Monte Carlo price) and predicted price (e.g. GPR approximative price), and $p_0\left(S_0\right)$ is the actual price without shock. The MAPE evaluating the average relative error of the estimate and the true observation is a common choice for benchmarking in finance. For risk assessment purpose, we measure the error of value-at-risk and expected shortfall calculated based on surrogate models (i.e. $\widehat{\mathrm{VaR}}(·,h,\alpha )$ and $\widehat{\mathrm{ES}}(·,h,\alpha )$) against the full revaluation ones calculated based on actual prices (i.e. $\mathrm{VaR}(·,h,\alpha )$ and $\mathrm{ES}(·,h,\alpha )$). For the relative measure, the error is then standardised by today price $p_0\left(S_0\right)$

$$\begin{array}{c}\hfill err.\mathrm{VaR}\alpha ={\displaystyle \frac{|\widehat{\mathrm{VaR}}(·,h,\alpha )-\mathrm{VaR}(·,h,\alpha )|}{p_0\left(S_0\right)}},err.\mathrm{ES}\alpha ={\displaystyle \frac{|\widehat{\mathrm{ES}}(·,h,\alpha )-\mathrm{ES}(·,h,\alpha )|}{p_0\left(S_0\right)}}.\end{array}$$

(37)

Finally, the computation time is also reported to compare the resource efficiency.

5.2. Model Specification

For neural network model, we train a two-hidden-layer Relu network of 50 hidden units (cf. Section C). We observe numerically that more complex architectures do not significantly improve the approximation performance. Network parameters are calibrated by the stochastic gradient descent algorithm, i.e. Adam optimizer [31], using 2000 epochs. The learning rate is set at $0.01$ and each batch takes one fifth of the training data size.

The standard Gaussian process regression is trained using the scikit-learn Python package. We set a zero mean prior distribution and use Matern kernels with $\nu =2.5$ (see (7)) to model the covariance matrix. A single length-scale is selected; however, the input is normalized by subtracting its empirical mean and then dividing by its empirical standard error. The optimization is restarted 5 times for the numerical stability. To implement multi-fidelity Gaussian process regression, we use emukit package [43] which is built upon the Gpy framework developed by machine learning researchers from the University of Sheffield. The model configuration is either set as the standard Gaussian processes or left at its default. The optimization is rerun 5 times.

6. Numerical Results

The main aim of this section is to provide a numerical analysis of algorithms based on Gaussian process regression (GPR) and the multi-fidelity model (mGPR) when applied to the repeated valuation and market risk calculation of an equity trading portfolio. It is worth emphasizing that a key advantage of the proposed new methods is that valuation and risk assessment are performed at the portfolio level, rather than on a deal-by-deal basis, which is the case with traditional approaches.

6.1. Mono-Asset Options Portfolio Case

Figure 4 and Figure 5 and Table 4 outline the performance of Gaussian Process Regression (GPR) in the full revaluation and risk calculation of a mono-asset portfolio with the sub-portfolio techniques outlined in Section 3.3. The accuracy of the GPR is evident in both pricing and risk calculation. For this mono-asset experiment, we train the regression models by a regular grid and backtest with more (uniformly) random data [14,34]. In particular, the training and test stock prices are shocked around their initial values using their diffusion models, e.g.

$$\left[S_0^{i}exp\left((r-{\displaystyle \frac{{\left({\sigma }^i\right)}^2}{2}})h-3\sqrt{h}{\sigma }^i\right),S_0^{i}exp\left((r-{\displaystyle \frac{{\left({\sigma }^i\right)}^2}{2}})h+3\sqrt{h}{\sigma }^i\right)\right],$$

where $S_0^i,{\sigma }^i,r$ and h denote respectively the initial stock price and volatility of the i-th asset, and the common interest rate and time step.

In Figure 4 and Figure 5, a precise fit of the price function is achieved when employing 5 points or more for training each sub-portfolio. The risk measures estimated by the GPR approach consistently converge to the benchmark measures obtained by the full revaluation approach, i.e. the intensive Monte Carlo method with 100,000 paths. The latter takes approximately 2 hours and 25 minutes on our machine. When learning from 10 points (resp. 5 points), the GPR achives a speed up of 1,000 times (resp. 2,000 points) compared to the full revaluation approach. The high runtime of the full revaluation approach in this experiment is primarily due to the American pricing calculations using a binomial tree with 100 time steps. This section concludes the performance of the Gaussian process regression in the mono-asset portfolio pricing and risk calculation.

6.2. Multi-Asset Options Case

Figure 6 and Table 5 display the results of the pricing and risk calculation for a five-asset geometric average call. Both GPR and mGPR models are trained using Monte Carlo prices with 100,000 paths, while the analytical prices are used as a reference for backtesting. Additionally, mGPR uses 200 prices obtained through the Monte Carlo method with 100 paths as low-fidelity data, which is still less expensive than one price point used 100,000 paths. The top left plot in Figure 6 shows that two GPR models trained with 50 points can outperform Monte Carlo prices. The error decreases, reaching 10 basis points at 200 training points. Notably, mGPR outperforms GPR at lower numbers of data points before their error curves converge.

The quantile-quantile (q-q) plots in Figure 6 trace the left tail of the simulated price distribution below the 10th percentile, corresponding to the right tail of the loss in the risk calculation in Table 5. Despite the Monte Carlo price having a higher MAPE error compared to GPR prices, its tail perfectly matches the true price. The tails of GPR and mGPR prices also converge as the number of training points increases. We stop the convergence curve of mGPR at 200 training points (corresponding to high-fidelity data) because we only use 200 low-fidelity points. In Table 5, the one-day $\mathrm{VaR}99\%$ and $\mathrm{ES}97.5\%$ estimated by GPRs using 200 data points have an error of around 7 basis points. Notably, GPR models significantly save time in risk calculation, including sampling and training time with 200 points (13 seconds for GPR compared to 7488 seconds for the brute force Monte Carlo approach). The mGPR model remains an effective approximation for risk measurement, especially in cases of high pricing engine costs, such as with complex portfolios, and when using minimal price points is preferred.

Figure 7, Figure 8 and Figure 9 and Table 6, Table 7 and Table 8 show the results for other options, with similar computational times to those in Table 5. Most of the remarks made in the geometric average call case are still valid for these other cases. In the arithmetic average basket call case, the risk measures computed by two GPR models have an error of about 23 basis points, which is less than 10 basis points in the two other cases. These above errors can be computed by applying results in Table 6, Table 7 and Table 8 into (37).

6.3. Multi-Asset Options Portfolio Case

Figure 10 and Table 9 illustrate the results of a multi-asset options portfolio . The time horizon of the value-at-risk and expected shortfall is set at one day. The second-order sensitivity-based approach (SxS) and the neural network learning from 40 points are not reported here due to their low precision. The learning of mGPR incorporates 500 prices computed by 500 Monte Carlo paths as low-fidelity data. Despite the portfolio complexity, both GPR models (GPR and mGPR) achieve impressive results with only 40 points. The GPR model has a MAPE of $0.61\%$, slightly outperformed by the mGPR model with a MAPE of $0.57\%$. While SxS performs about 10 times worse with an error of $6.13\%$. Figure 10 highlights this difference, showing that the predicted left tail of SxS prices diverges significantly from the true price represented by the red diagonal line. While the predicted left tails of GPR and mGPR prices are not perfect, but they still concentrate around the tail of the true price, i.e. the red line.

In Table 9, mGPR predictions offer more accurate estimates of the $\mathrm{VaR}99\%$ and $\mathrm{ES}97.5\%$ compared to GPR predictions, with errors of 34bps and 32bps, respectively, versus 54bps for both errors in GPR. As more data becomes available, the performance of all models improves. Generally, two GPR models outperform NN, with a MAPE of $0.42\%$ ($0.38\%$) compared to the one of NN’s $0.5\%$ ($0.43\%$) when learning with $N_{train}=100$ ($N_{train}=500$). The difference is more pronounced in the prediction of the left tail, as seen in the two bottom plots in Figure 10. The left tails predicted by GPR and mGPR closely align with the true one, while NN underestimates the left tail. When trained with $N_{train}=100$, mGPR’s errors in risk measures, around 30bps, are only half of GPR’s errors. However, this difference disappears when training with $N_{train}=500$.

Table 10 displays the runtime for various risk calculation methods on our server with an Intel(R) Xeon(R) Gold 5217 central processing unit (CPU) and a Nvidia Tesla V100 graphics processing unit (GPU). The full repricing approach involves simulating $100,000$ scenarios of shock, with each scenario being a Monte Carlo simulation of $100,000$ paths computed over 500 derivatives. Without advanced programming on the GPU, the sampling procedure in banks may take up to one day. Alternative approaches show significant computational performance gains. For instance, standard Gaussian processes regression achieves a speedup of 230 times compared to the full repricing approach when learning with $N_{train}=100$ points, i.e. a total of 2.6 comparing to 600 seconds. Despite mGPR’s longer learning times (i.e. 56 when learning with $N_{train}=100$), its computational efficiency compared to full repricing remains remarkable. It is important to note that, under realistic constraints of the more time-consuming pricing engine in banks, the mGPR approach may take a bit longer than the GPR approach for the learning procedure, but it provides better precision in risk calculation, especially when limiting the use of the pricing engine. Table 11 and Table 12 and Figure 11 and Figure 12 correspond to the risk assessment results for a ten-day and one-month horizon, and observations from the one-day horizon case hold true for other time horizons.

7. Conclusions

We have introduced statistical and machine learning tools to handle the repeated valuation of extensive portfolios comprising linear and non-linear derivatives. To be more precise, we applied a Bayesian non-parametric technique, known as Gaussian process regression (GPR). To evaluate the accuracy and performance of these algorithms, we have considered a multi-asset options portfolio and a portfolio of non-linear derivatives, including vanilla and barrier/American options on an equity asset. The numerical tests demonstrated that the GPR algorithm outperforms pricing models in efficiently conducting repeated valuations of large portfolios. It is noteworthy that the GPR algorithm eliminates the need to separately revalue each derivative security within the portfolio, leading to a significant speed-up in calculating value-at-risk (VaR) and expected shortfall (ES). This independence from the size and composition of the trading portfolio is remarkable. Consequently, we argue that it is more advantageous for banks to construct their risk models using the power of GPR techniques in terms of calculation accuracy and speed-up.

Moreover, we have investigated multi-fidelity modeling technique and have explored its applications in pricing and risk calculation, a relatively novel concept in quantitative finance. The multi-fidelity modeling approach leverages limited access to the price engine while incorporating information from other related, more affordable resources to enhance the approximation. Trained regression models offer price estimates at the portfolio level, making them more favorable for risk calculation compared to classical pricing approaches that evaluate products individually. Our numerical findings indicate that machine learning models such as neural networks and Gaussian process regression provide more accurate results and significant time savings compared to traditional sensitivity-based approaches, maintaining precision in risk calculation akin to the full repricing approach. With the limitation of training price data, GPR models yield more impressive results than neural network. Despite the longer learning time required for mGPR, especially when low-fidelity data is available, it yields better or at least equal precision compared to the standard GPR model. Notably, as multi-fidelity modeling is not yet well-known in quantitative finance, it presents an intriguing direction for research, allowing scientists to explore traditional relationships in finance to define low-fidelity models and enhance the learning process.

Future research delves into other topics and applications of algorithms such that Gaussian process regression and multi-fidelity modeling, with a particular focus on the fixed income portfolio. We also envisage the use of the replicating portfolio value, assumed to be accessible, as low-fidelity data to leverage the learning of the target portfolio value. Additionally, the price of European basket options can serve as low-fidelity data to approximate the price of their American/Bermudian counterparts. The Black-Scholes price of derivatives can be used to improve the approximation of their Heston price. We leave these interesting applications there for future study.