Low-rank tensor approximation for Chebyshev interpolation in parametric option pricing

TL;DR

This paper introduces a low-rank tensor approximation method using tensor train format to efficiently perform high-dimensional Chebyshev interpolation in parametric option pricing, overcoming the curse of dimensionality.

Contribution

It extends Chebyshev interpolation for high-dimensional problems by exploiting low-rank tensor structures, enabling efficient computation in up to 25-dimensional parameter spaces.

Findings

The method effectively handles high-dimensional parameter spaces in option pricing.

Numerical results show the low-rank structure and efficiency of the approach.

Compared to existing techniques, the proposed method demonstrates superior performance.

Abstract

Treating high dimensionality is one of the main challenges in the development of computational methods for solving problems arising in finance, where tasks such as pricing, calibration, and risk assessment need to be performed accurately and in real-time. Among the growing literature addressing this problem, Gass et al. [14] propose a complexity reduction technique for parametric option pricing based on Chebyshev interpolation. As the number of parameters increases, however, this method is affected by the curse of dimensionality. In this article, we extend this approach to treat high-dimensional problems: Additionally exploiting low-rank structures allows us to consider parameter spaces of high dimensions. The core of our method is to express the tensorized interpolation in tensor train (TT) format and to develop an efficient way, based on tensor completion, to approximate the interpolation coefficients. We apply the new method to two model problems: American option pricing in the Heston model and European basket option pricing in the multi-dimensional Black-Scholes model. In these examples we treat parameter spaces of dimensions up to 25. The numerical results confirm the low-rank structure of these problems and the effectiveness of our method compared to advanced techniques.

Figures

Tables

Table 1. Table 1: Completion results on 𝒫 𝒫 {\mathcal{P}} for the parametric American put option pricing problem in the Heston model.

final $\left|\Omega \right|$	rel err on last ${\Omega }_c^{new}$	rel err on full $𝒫$	completion time (s)
8050 (5 %)	$2.56\cdot {10}^{-5}$	$2.75\cdot {10}^{-5}$	366.12

Table 2. Table 2: Results on American put option pricing via combined methodology and reference method.

time reference method (s)	time interpolation (s)	max abs error
$3.65\cdot {10}^{-2}$	$4.89\cdot {10}^{-4}$	$1.95\cdot {10}^{-4}$

Table 3. Table 3: Completion parameters for constructing 𝒫 𝒫 {\mathcal{P}} . Case of uncorrelated assets.

	$d$	$\rho$	$tol$	$to{l}^{\prime }$	$r_{\max }$	initial $\left|\Omega \right|$	$\left|{\Omega }_C\right|$	$p$
$n=4$	5	0	${10}^{-2}$	${10}^{-8}$	$5$	31	31	${10}^{-1}$
	10	0	${10}^{-2}$	${10}^{-8}$	$5$	78	78	${10}^{-2}$
	15	0	${10}^{-2}$	${10}^{-8}$	$5$	214	214	${10}^{-5}$
	20	0	${10}^{-2}$	${10}^{-8}$	$5$	763	763	${10}^{-8}$
	25	0	${10}^{-2}$	${10}^{-8}$	$5$	2086	2086	${10}^{-11}$
$n=6$	5	0	${10}^{-3}$	${10}^{-8}$	$7$	17	17	${10}^{-1}$
	10	0	${10}^{-3}$	${10}^{-8}$	$7$	282	141	${10}^{-3}$
	15	0	${10}^{-3}$	${10}^{-8}$	$7$	475	475	${10}^{-6}$
	20	0	${10}^{-3}$	${10}^{-8}$	$7$	798	798	${10}^{-10}$
	25	0	${10}^{-3}$	${10}^{-8}$	$7$	1341	1341	${10}^{-15}$

Table 4. Table 4: Completion results on 𝒫 𝒫 {\mathcal{P}} for the basket option pricing problem in the Black and Scholes model. Case of uncorrelated assets.

	$d$	final $\left|\Omega \right|$	rel err on last ${\Omega }_C^{new}$	completion time (s)
$n=4$	5	124	$3.42\cdot {10}^{-3}$	9.90
	10	546	$2.54\cdot {10}^{-6}$	67.44
	15	1712	$3.55\cdot {10}^{-8}$	171.14
	20	2289	$5.03\cdot {10}^{-8}$	193.90
	25	4172	$3.96\cdot {10}^{-9}$	226.38
$n=6$	5	204	$2.40\cdot {10}^{-4}$	52.55
	10	987	$1.20\cdot {10}^{-6}$	198.27
	15	1900	$2.28\cdot {10}^{-7}$	429.39
	20	3192	$2.97\cdot {10}^{-7}$	732.49
	25	4023	$1.35\cdot {10}^{-7}$	999.25

Table 5. Table 5: Basket option prices computed via Chebyshev interpolation (combined methodology) versus MC reference method with 10 4 superscript 10 4 10^{4} simulations for n = 4 𝑛 4 n=4 and 10 5 superscript 10 5 10^{5} simulations for n = 6 𝑛 6 n=6 . Case of uncorrelated assets.

	$d$	time reference method (s)	time interpolation (s)	max abs error
$n=4$	5	$0.18$	$0.45\cdot {10}^{-3}$	$3.75\cdot {10}^{-3}$
	10	$0.19$	$0.64\cdot {10}^{-3}$	$5.21\cdot {10}^{-4}$
	15	$0.20$	$0.73\cdot {10}^{-3}$	$4.38\cdot {10}^{-4}$
	20	$0.20$	$1.09\cdot {10}^{-3}$	$3.16\cdot {10}^{-4}$
	25	$0.21$	$0.97\cdot {10}^{-3}$	$2.08\cdot {10}^{-4}$
$n=6$	5	1.84	$0.40\cdot {10}^{-3}$	$5.20\cdot {10}^{-4}$
	10	1.91	$0.61\cdot {10}^{-3}$	$1.42\cdot {10}^{-4}$
	15	1.99	$0.78\cdot {10}^{-3}$	$1.02\cdot {10}^{-4}$
	20	2.04	$0.93\cdot {10}^{-3}$	$1.01\cdot {10}^{-4}$
	25	2.10	$1.04\cdot {10}^{-3}$	$9.36\cdot {10}^{-5}$

Table 6. Table 6: Completion results on 𝒫 𝒫 {\mathcal{P}} for the basket option pricing problem in the Black and Scholes model. Case of correlated assets.

	$d$	final $\left|\Omega \right|$	rel err on last ${\Omega }_C^{new}$	completion time (s)
$n=4$	5	124	$1.86\cdot {10}^{-3}$	8.95
	10	390	$2.19\cdot {10}^{-4}$	65.73
	15	1284	$1.72\cdot {10}^{-7}$	118.73
	20	1526	$2.49\cdot {10}^{-8}$	168.20
	25	4172	$7.52\cdot {10}^{-9}$	215.44
$n=6$	5	255	$4.40\cdot {10}^{-4}$	66.54
	10	987	$2.06\cdot {10}^{-4}$	200.15
	15	1900	$1.79\cdot {10}^{-7}$	432.58
	20	3990	$1.82\cdot {10}^{-8}$	852.13
	25	5364	$2.88\cdot {10}^{-7}$	1335.76

Table 7. Table 7: Basket option prices computed via Chebyshev interpolation (combined methodology) versus MC reference method with 10 4 superscript 10 4 10^{4} simulations for n = 4 𝑛 4 n=4 and 10 5 superscript 10 5 10^{5} simulations for n = 6 𝑛 6 n=6 . Case of correlated assets.

	$d$	time reference method (s)	time interpolation (s)	max abs error
$n=4$	5	0.18	$0.50\cdot {10}^{-3}$	$1.39\cdot {10}^{-3}$
	10	0.20	$0.56\cdot {10}^{-3}$	$4.82\cdot {10}^{-4}$
	15	0.20	$0.70\cdot {10}^{-3}$	$2.82\cdot {10}^{-4}$
	20	0.23	$0.91\cdot {10}^{-3}$	$2.93\cdot {10}^{-4}$
	25	0.23	$1\cdot {10}^{-3}$	$4.30\cdot {10}^{-4}$
$n=6$	5	1.85	$0.38\cdot {10}^{-3}$	$3.55\cdot {10}^{-4}$
	10	1.90	$0.57\cdot {10}^{-3}$	$5.58\cdot {10}^{-4}$
	15	1.99	$0.74\cdot {10}^{-3}$	$1.39\cdot {10}^{-4}$
	20	2.06	$0.90\cdot {10}^{-3}$	$1.41\cdot {10}^{-4}$
	25	2.15	$0.96\cdot {10}^{-3}$	$9.28\cdot {10}^{-5}$

Equations

$${\mathsf{Price}}^{\mathbf{p}}\approx \sum _{j_1=0}^{n_1}\cdots \sum _{j_d=0}^{n_d}{c}_{j_1,\dots ,j_d}{T}_{j_1,\dots ,j_d}(\mathbf{p}).$$

$$I_{\overline{\mathbf{n}}}({\mathsf{Price}}^{(\cdot )})(\mathbf{p})=\sum _{j_1=0}^{n_1}\cdots \sum _{j_d=0}^{n_d}{c}_{j_1,\dots ,j_d}{T}_{j_1,\dots ,j_d}(\mathbf{p}).$$

$$T_{j_1,\dots ,j_d}(\mathbf{p})=\prod _{i=1}^d{T}_{j_i}(p_i),T_{j_i}(p_i)=\cos (j_i\arccos (p_i)),$$

$$c_{j_{1},\dots,j_{d}}=\Big{(}\prod_{i=1}^{d}\frac{2^{\mathbbm{1}_{n_{i}>j_{i}>0}}}{n_{i}}\Big{)}\sideset{}{{}^{\prime\prime}}{\sum}_{k_{1}=0}^{n_{1}}\dots\sideset{}{{}^{\prime\prime}}{\sum}_{k_{d}=0}^{n_{d}}{\mathcal{P}}(k_{1},\dots,k_{d})\prod_{i=1}^{d}\cos\big{(}j_{i}\pi\frac{k_{i}}{n_{i}}\big{)},$$

$$\mathcal{P}(k_1,\dots ,k_d)={\mathsf{Price}}^{q_{k_1,\dots ,k_d}},$$

$$X^{<\mu >}\in {\mathbb{R}}^{(n_1{n}_2\cdots n_{\mu })\times (n_{\mu +1}\cdots n_d)},$$

$${\text{rank}}_{\text{TT}}(\mathcal{X})=(r_0,r_1,\cdots ,r_d):=(1,\text{rank}(X^{<1>}),\cdots ,\text{rank}(X^{<d-1>}),1).$$

$$\mathcal{X}(i_1,i_2,\cdots ,i_d)=U_1(i_1)U_2(i_2)\cdots U_d(i_d),$$

$$\mathcal{X}(i_1,i_2,\cdots ,i_d)=\sum _{k_1=1}^{r_1}\cdots \sum _{k_{d-1}=1}^{r_{d-1}}{\mathbf{U}}_1(1,i_1,k_1){\mathbf{U}}_2(k_1,i_2,k_2)\cdots {\mathbf{U}}_d(k_{d-1},i_d,1).$$

$$⟨\mathcal{X},\mathcal{Y}⟩=⟨\text{vec}(\mathcal{X}),\text{vec}(\mathcal{Y})⟩=\sum _{i_1=1}^{n_1}\cdots \sum _{i_d=1}^{n_d}\mathcal{X}(i_1,\dots ,i_d)\mathcal{Y}(i_1,\dots ,i_d),$$

$$\mathcal{Z}(i_1,\cdots ,i_{\mu -1},j,i_{\mu +1}\cdots ,i_d)=\sum _{i_{\mu }=1}^{n_{\mu }}\mathcal{X}(i_1,\cdots ,i_d)M(j,i_k),j=1,\dots ,m.$$

$$\begin{array}{rl}\underset{\mathcal{X}}{\min }& ||P_{\Omega }\mathcal{X}-P_{\Omega }\mathcal{A}|{|}^2\\ \text{subject to}& \mathcal{X}\in {\mathcal{M}}_{\mathbf{r}}:=\{\mathcal{X}\in {\mathbb{R}}^{n_1\times \cdots \times n_d}|{\text{rank}}_{\text{TT}}=\mathbf{r}\},\end{array}$$

$${ϵ}_{\Omega }({\mathcal{X}}_k):=\frac{\Vert P_{\Omega }\mathcal{A}-P_{\Omega }{\mathcal{X}}_k\Vert }{\Vert P_{\Omega }\mathcal{A}\Vert },{ϵ}_{{\Omega }_C}({\mathcal{X}}_k):=\frac{\Vert P_{{\Omega }_C}\mathcal{A}-P_{{\Omega }_C}{\mathcal{X}}_k\Vert }{\Vert P_{{\Omega }_C}\mathcal{A}\Vert }.$$

$$\frac{|{ϵ}_{\Omega }({\mathcal{X}}_k)-{ϵ}_{\Omega }({\mathcal{X}}_{k+1})|}{|{ϵ}_{\Omega }({\mathcal{X}}_k)|}<\delta \text{and}\frac{|{ϵ}_{{\Omega }_C}({\mathcal{X}}_k)-{ϵ}_{{\Omega }_C}({\mathcal{X}}_{k+1})|}{|{ϵ}_{{\Omega }_C}({\mathcal{X}}_k)|}<\delta ,$$

$${\mathrm{err}}_{\mathrm{new}}\leftarrow {ϵ}_{\Gamma }({\mathcal{X}}_c).$$

$$f:[0,1{]}^4\to \mathbb{R},f(\mathbf{x})=\exp (-\Vert \mathbf{x}\Vert )$$

$$\mathcal{C}(i_1,i_2,\cdots ,i_d)=c_{i_1-1,i_2-1,\cdots ,i_d-1},$$

$${\mathcal{C}}(j+1)=\frac{2^{\mathbbm{1}_{n_{1}>j>0}}}{n_{1}}\sideset{}{{}^{\prime\prime}}{\sum}_{k=0}^{n_{1}}{\mathcal{P}}(k+1)\cos\big{(}j\pi\frac{k}{n_{1}}\big{)},\quad j=0,\cdots,n_{1},$$

$$\mathcal{C}=\frac{2}{n_1}\left(\begin{array}{ccccc}\frac{1}{4}& \frac{1}{2}& \dots & \frac{1}{2}& \frac{1}{4}\\ \frac{1}{2}& \cos (\frac{\pi }{n_1})& \dots & \cos (\frac{\pi (n_1-1)}{n_1})& \frac{1}{2}\cos (\pi )\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ \frac{1}{2}& \cos (\frac{\pi (n_1-1)}{n_1})& \dots & \cos (\frac{\pi (n_1-1{)}^2}{n_1})& \frac{1}{2}\cos (\pi (n_1-1))\\ \frac{1}{4}& \frac{1}{2}\cos (\pi )& \dots & \frac{1}{2}\cos (\pi (n_1-1))& \frac{1}{4}\cos (\pi n_1)\end{array}\right)\mathcal{P},$$

$$I_{\overline{\mathbf{n}}}({\mathsf{Price}}^{(\cdot )})(\mathbf{p})=⟨\mathcal{C},{\mathcal{T}}_{\mathbf{p}}⟩.$$

$$d{S}_t=r{S}_{t}dt+\sqrt{v_t}{S}_{t}d{W}_t^1,$$

$$d{v}_t=\kappa (\theta -v_t)dt+\sigma \sqrt{v_t}d{W}_t^2.$$

$$\mathsf{Price}=\underset{t<\tau <T}{\sup }\mathbb{E}[e^{-r\tau }f(S_{\tau })|S_t=s,v_t=v],$$

$$f(x)=(K-x{)}^{+},$$

$$\begin{array}{r}\{\begin{array}{l}{\partial }_t\mathsf{Price}\ge \mathcal{G}\mathsf{Price}\\ \mathsf{Price}\ge f\\ (\mathsf{Price}-f)({\partial }_t\mathsf{Price}-\mathcal{G}\mathsf{Price})=0,\end{array}\end{array}$$

$$\mathcal{G}g(s,v)=\frac{1}{2}{s}^{2}v{\partial }_{ss}^{2}g+\rho \sigma sv{\partial }_{sv}^{2}g+\frac{1}{2}{\sigma }^{2}v{\partial }_{vv}^{2}g+rs{\partial }_{s}g+\kappa (\theta -v){\partial }_{v}g-rg.$$

$$S_0=2,v_0=0.0175,r=0.1,T=0.25,$$

$$(K,\rho ,\sigma ,\kappa ,\theta )\in [2;4]\times [-1;1]\times [0.2;0.5]\times [1;2]\times [0.05;0.2]$$

$$\rho =0,tol=10^{-3},to{l}^{\prime }=10^{-8},r_{\max }=10,|\Omega |=805,|{\Omega }_C|=805,p=\mathrm{0.2.}$$

$$\max (|P_{\text{Int}}-P_{\text{Ref}}|),$$

Full text

Key words

Chebyshev interpolation, parametric option pricing, high-dimensional problem, tensor train format, low-rank tensor approximation, tensor completion

1 Introduction

Financial problems are, by their nature, multi- and high-dimensional, because a large number of risk factors contribute to the prices of each financial asset. Moreover, the banking, insurance and hedge fund industry draws on investments in large portfolios. The interdependencies of both the risk factors and the assets make basic computational tasks such as model calibration, pricing, and hedging as well as more global tasks such as uncertainty quantification, risk assessment and capital reserve calculation computationally extremely challenging, see for instance [5].

Automatic and high-speed trading challenge the computational methods in that the results need to be available fast and with minimal storage requirement. Moreover, we observe rising regulatory requirements. On the one hand, more realistic modeling demands more prudent considerations, which leads to rising computational complexity. On the other hand, the availability of requested performance characteristics is expected to be delivered within shorter periods of time. This poses a high challenge for traditional approaches, which typically suffer from low convergence rates in higher dimensions, see for instance [9, 12].

For the reasons explained above, the development of efficient computational methods for high-dimensional problems in finance is an utmost active field of research in both academia and industry. For example, further developments of the Monte Carlo method have been very successfully applied to financial problems; we refer to [35, 16] for the quasi Monte Carlo method and to [15] for the multilevel Monte Carlo method. Besides stochastic integration, deterministic numerical integration has been exploited using sparse grid techniques, see [19, 27, 6]. Also PDE methods have been extended to multivariate problems in finance. For instance using operator splitting methods as in [28], principal component analysis and expansions as in [42], and wavelet compression techniques proposed in [36, 25, 26].

Exploiting the particular structure of a problem, complexity reduction techniques exhibit great potential to save run-time and storage capacity while maintaining the required accuracy. In numerical analysis and a large variety of applications, for example in engineering and medicine, complexity reduction techniques have been developed and implemented with great success. For instance the field of reduced basis methods to efficiently solve parametric partial differential equations (PDEs) has experienced a tremendous development over the last decade, see, e.g., [23, 40, 41] and the references therein. Pioneered by [43, 10] the potential of reduced basis methods is also increasingly exploited for problems in finance; see [8, 37, 7] for examples. These methods can be viewed as high-dimensional interpolation methods that are trained in an offline step to solve a specific class of parametric PDEs. In this article we explore direct interpolation of multivariate functions as a unified approach to complexity reduction for finance.

Our starting point is the tensorized Chebyshev interpolation of conditional expectations in the parameter and state space, as introduced in [14]. Having observed for a large set of applications that these functions are highly regular, admitting sensitivities of high order or even being analytic, and that the domain of interest can be restricted to a hyperrectangular, Chebyshev interpolation is a promising choice: Its convergence is subexponential for multivariate analytic functions, its implementation is numerically stable, and the coefficients are simply given by a linear transformation of the function values at the nodal points. In this article we exploit this favorable structure further for high dimensionality. In passing, we point out that, while we choose Chebyshev interpolation for the reasons listed in this paragraph, the technique presented in this paper extends to other tensorized interpolation techniques. Also, our approach is applicable beyond option pricing and finance.

The basis of our approach is the following. In an offline phase, the price as function of parameters ${\mathbf{p}}\in[-1,1]^{d}$, ${\mathbf{p}}\mapsto\mathsf{Price}^{\mathbf{p}}$ is evaluated at selected parameter samples ${\mathbf{p}}$ to prepare an approximation by tensorized Chebyshev polynomials $T_{j_{1},\dots,j_{d}}$ with pre-computed Fourier coefficients $c_{j_{1},\dots,j_{d}}$, as follows,

[TABLE]

To evaluate the function in the online phase, only the multivariate polynomials on the right-hand side need to be evaluated. However, implementing (1) in a straightforward manner exposes the method to the curse of dimensionality in both the offline and the online phase: In the offline phase, the prices need to be evaluated on a tensorized grid of Chebyshev nodes, amounting to $O(n^{d})$ parameter samples when $n$ nodes are required for each parameter. This is computationally costly, especially if the underlying pricing method is already computationally demanding. In the online phase alike, $O(n^{d})$ operations are needed for evaluating the approximating multivariate polynomial. Even for a number as low as $n=3$, corresponding to quadratic polynomials, a problem with $d=20$ parameters becomes infeasible.

One approach to breaking the curse of dimensionality that has already proven effective in a number of areas is to exploit low-rank structures of high-dimensional tensors; see [18, 20, 29] and the references therein. These techniques reduce, sometimes dramatically, memory requirements and the cost of operating with tensors. In the context of parametric PDEs, low-rank tensor structures have been successfully exploited in, e.g., [1, 4, 30, 34, 45]. As option prices are characterized as solutions of parabolic PDEs, this gives hope that low-rank structures can be exploited in finance as well. The following questions arise:

Can we detect low-rank structures for the problem of form (1)? Existing theoretical studies only provide partial answers to this question, either not reflecting the observed effectiveness of low-rank techniques or being limited to rather specific function classes; see [11, 20, 44] for examples. We therefore approach the question from an experimental perspective and analyze examples of different nature and different dimensionality in Section 3. The results clearly indicate an approximate low-rank structure of the tensor $\mathcal{P}$ containing the prices evaluated at the nodes of the tensorized Chebyshev grid. In the specific case of the interpolation of American option prices in the Heston model in five parameters we can explicitly compare the full tensor $\mathcal{P}$ with the one resulting from low-rank approximation. We perform this comparison in Section 3.1, which confirms the low-rank structure of $\mathcal{P}$. In Section 3.2 we consider prices of basket options in the Black-Scholes model with up to $25$ underlyings and interpolate in the initial values of the underlyings. Although the resulting full tensor $\mathcal{P}$ is too large to be explicitly computed and compared with, we provide a structural analysis in Section 3.2.3 that explains why $\mathcal{P}$ is expected to exhibit low-rank structure.

How can we exploit low-rank structures for the problem of form (1)? Expressing the problem in a tensor format reveals that exploiting the tensor structure itself (even without low-rank structure) leads to a considerable efficiency gain in both the offline and the online phase. Next, we explore existing low-rank tensor techniques. In order to efficiently exploit these techniques for problem (1), we need to introduce several new components resulting in the new method. We detail these steps below.

In order to construct the interpolation coefficients $c_{j_{1},\dots,j_{d}}$ in the offline phase, it is first required to compute or approximate all values of the tensor $\mathcal{P}$, containing the prices in the tensorized Chebyshev grid. Evaluating $\mathcal{P}$ explicitly is too costly for larger $d$, especially when the underlying pricing procedure is computationally expensive. Instead we only compute part of the entries of $\mathcal{P}$ and then need to deal with an incomplete tensor. This leads us to the following first step:

We start by computing the prices for a small portion of the Chebyshev grid points only. Then, we adapt a completion algorithm (in Section 2.3) which allows us to approximate the tensor of prices for the complete Chebyshev grid by fitting tensors of pre-specified low rank to the provided data points. As it is not reasonable to assume a priori knowledge of low-rank structure, the completion procedure needs to be combined with an adaptive rank and sampling strategy. Specifically, we repeat the process of adding new samples and increasing the pre-specified rank until an adequate stopping criterion is fulfilled. This completion algorithm is designed to work with tensors built and stored in tensor train (TT) format.

With the low-rank approximation of the tensor ${\mathcal{P}}$ in TT format at hand, we can then approximate efficiently the Fourier coefficients $c_{j_{1},\dots,j_{d}}$. This is the last step of the offline phase:

The computation of the tensor $\mathcal{C}$, containing the Fourier coefficients $c_{j_{1},\dots,j_{d}}$, is computed by a sequence of $d$ tensor-matrix multiplications. The particular structure of the involved matrices facilitates the use of the fast Fourier transform, leading to a complexity of $O(dnr^{2}\log(n))$, where $r$ is determined by the ranks of $\mathcal{P}$. This step is explained in Section 2.4.2.

Suppose now that, in the online phase, we want to compute the interpolated price (1) for a new set of parameter samples. Given the tensor $\mathcal{C}$ in TT format, the evaluation of (1) for a price ${\mathbf{p}}$ is performed efficiently as follows:

First, each of the Chebyshev polynomials involved in the tensorized Chebyshev basis is evaluated in ${\mathbf{p}}$. It turns out that (1) can be viewed as inner product between $\mathcal{C}$ and a rank-one tensor. Thanks to the TT format, the complexity of computing this inner product is $O(dnr^{2})$; see Section 2.2. As long as $r$ is reasonably small, this compares favorably with the $O(n^{d})$ operations needed by the standard approach.

In Section 3, we test the performance of the new method for two different option pricing problems, the interpolation of

• –

American option prices in the Heston model in $d=5$ parameters, and of

• –

prices of basket options in the Black-Scholes model in up to $d=25$ underlyings.

At comparable accuracy, the interpolation in American option prices reveals a promising gain in efficiency when compared to an ADI-based PDE solver. The efficiency gain for the basket option prices is shown in comparison to a Monte Carlo simulation with variance reduction.

2 TT format and tensor completion for Chebyshev interpolation

This section describes the methodology proposed in this work. We start with recalling the tensorized Chebyshev interpolation method from [14]. After introducing the TT format [39], we present and extend the tensor completion approach from [45]. Finally, we explain how to combine these algorithms in order to efficiently price parametric options for a large number of parameters.

2.1 Chebyshev interpolation for parametric option pricing

We consider an option price that depends on a vector of $d$ parameters $\mathbf{p}$ contained in $[-1,1]^{d}$; general hyperrectangular parameter domains can be addressed by a suitable affine transformation. The basic idea developed in [14] consists of using tensorized Chebyshev interpolation in the parameters (model and payoff parameters) to increase the efficiency of computing option prices, while maintaining satisfactory accuracy. Writing $\mathsf{Price}^{\mathbf{p}}$ for the price evaluated in ${\mathbf{p}}$, the Chebyshev interpolation of order $\overline{\mathbf{n}}:=(n_{1},\dots,n_{d})$ with $n_{i}\in\mathbb{N}_{0}$ is given by

[TABLE]

The basis functions $T_{j_{1},\dots,j_{d}}$ are constructed from Chebyshev polynomials by

[TABLE]

and the coefficients $c_{j_{1},\dots,j_{d}}$ are defined as

[TABLE]

where $\sideset{}{{}^{\prime\prime}}{\sum}$ indicates that the first and the last summand are halved. The tensor ${\mathcal{P}}$ contains the prices on the tensorized Chebyshev grid:

[TABLE]

where $q_{k_{1},\dots,k_{d}}:=(q_{k_{1}},\dots,q_{k_{d}})$ is defined via Chebyshev nodes $q_{k_{i}}:=\cos(\pi\frac{k_{i}}{n_{i}})$ for $k_{i}=0,\dots,n_{i}$ and $i=1,\dots,d$. A convergence analysis of the tensorized Chebyshev interpolation in the setting of option pricing is given in [14].

The tensor $\mathcal{P}$ in equation (4) is of order $d$ and size $(n_{1}+1)\times\dots\times(n_{d}+1)$. The interpolation procedure first requires to compute each entry of this tensor with the reference method. This becomes expensive when the interpolation order and the dimension $d$ increase. We will use tensor completion to lower this cost.

Remark 2.1** (Choice of interpolation order).**

In our numerical experiments, the interpolation order $\overline{\mathbf{n}}$ is chosen a priori for simplicity. However, this choice can be made adaptively as explained in [22] for the case $d=3$ (the extension to general $d$ is straightforward).

2.2 TT format

For recalling the TT format introduced in [39], we consider a general tensor $\mathcal{X}\in\mathbb{R}^{n_{1}\times n_{2}\times\cdots\times n_{d}}$ of order $d$. For each $\mu=1,\ldots,d-1$, the entries of $\mathcal{X}$ can be rearranged into a matrix

[TABLE]

which is called the $\mu$th unfolding of $\mathcal{X}$. For this purpose, the first $\mu$ indices of $\mathcal{X}$ are merged into the row index and the last $n-\mu$ indices into a column index; see [39] for a formal definition. The TT ranks of $\mathcal{X}$ form an integer tuple

[TABLE]

Every entry $\mathcal{X}(i_{1},i_{2},\cdots,i_{d})$ can be expressed as a product of $d$ matrices

[TABLE]

with $U_{\mu}(i_{\mu})$ a matrix of size $r_{\mu-1}\times r_{\mu}$. For each $\mu=1,\cdots,d$, one can then collect the $n_{\mu}$ matrices $U_{\mu}(i_{\mu})$, $i_{\mu}=1,2,\cdots,n_{\mu}$ into a third order tensor $\mathbf{U}_{\mu}$ of size $r_{\mu-1}\times n_{\mu}\times r_{\mu}$. These tensors are called TT cores and, by construction, we have

[TABLE]

Figure 1 illustrates this so-called TT decomposition by a tensor network diagram [38].

Provided that the TT ranks remain moderate, a significant memory reduction is obtained by storing instead of $\mathcal{X}$ the TT cores: from $O(n^{d})$ to $O(dnr^{2})$, where $r=\max\{r_{0},\ldots,r_{d}\}$ and $n=\max\{n_{1},\ldots,n_{d}\}$.

Some operations can be effected quite cheaply in the TT format for tensors of low TT ranks. Let us first consider the inner product of two tensors $\mathcal{X},\mathcal{Y}\in\mathbb{R}^{n_{1}\times\cdots\times n_{d}}$ defined as

[TABLE]

where ${\text{vec}}(\cdot)$ stacks the entries of a tensor into a long vector. The corresponding tensor network diagram when $\mathcal{X}$ and $\mathcal{Y}$ are both in TT decomposition is shown in Figure 2. It can be seen that the summations in (7) become contractions between the TT cores of $\mathcal{X}$ and $\mathcal{Y}$. By carrying out these contractions of cores from the left to right, the cost of evaluating the inner product reduces from $O(n^{d})$ to $O(dnr^{3})$, where $r$ denotes the maximum of all involved TT ranks.

The mode-$\mu$ matrix multiplication between a tensor $\mathcal{X}\in\mathbb{R}^{n_{1}\times\cdots\times n_{d}}$ and a matrix $M\in\mathbb{R}^{m\times n_{\mu}}$ results in a tensor $\mathcal{Z}\in\mathbb{R}^{n_{1}\times\cdots n_{\mu-1}\times m\times n_{\mu+1}\cdots\times n_{d}}$ defined by

[TABLE]

We will denote this operation by $\mathcal{Z}=\mathcal{X}\times_{k}M$. If $\mathcal{X}$ is in TT decomposition (6) then it is straightforward to obtain a TT decomposition for $\mathcal{Z}$, by performing a mode-$2$ matrix multiplication of $\mathbf{U}_{\mu}$ with $M$. Letting $c_{M}$ denote the cost of multiplying $M$ with a vector, this requires $O(c_{M}nr)$ operations instead of the $O(c_{M}n^{d-1})$ operations needed when $\mathcal{X}$ is a general tensor.

2.3 Completion algorithm

The goal of completion algorithms is to reconstruct a given data set from a small fraction of its entries. As this is clearly an ill-posed task, one needs to additionally impose some regularization, such as smoothness conditions. In this work, we impose low TT ranks on the tensor ${\mathcal{P}}$ containing the prices and reconstruct ${\mathcal{P}}$ using the completion algorithm proposed in [45].

In the following, we briefly summarize the approach from [45]. Let ${\mathcal{A}}\in\mathbb{R}^{n_{1}\times\cdots\times n_{d}}$ denote the original data tensor for which only the entries in a (small) training set $\Omega\subset\{1,n_{1}\}\times\cdots\times\{1,n_{d}\}$ are known. When aiming at fitting a tensor of fixed (low) TT ranks $\mathbf{r}=(r_{0},\ldots,r_{d})$ to this data, completion takes the form of the constrained optimization problem

[TABLE]

where $P_{\Omega}{\mathcal{X}}$ denotes the orthogonal projection onto $\Omega$ and $\|\cdot\|$ is the norm induced by the inner product (7). It is known that ${\mathcal{M}}_{\mathbf{r}}$ is a smooth embedded submanifold, which enables one to apply Riemannian optimization techniques to (8). Specifically, in [45] it is proposed to employ a Riemannian conjugate gradient (CG) method (see Algorithm 1 in [45]). This method produces iterates that stay on the manifold and, in turn, can be stored and manipulated efficiently in the TT format. One iteration requires $O(dnr^{3}+d|\Omega|r^{2})$ operations, where $|\Omega|$ denotes the cardinality of $\Omega$.

Our stopping criterion of Riemannian CG is designed to attain a level of accuracy warranted by the data and the chosen TT ranks. Following [45], we choose a test set ${\Omega_{C}}$ of, say, $100$ additional parameter samples not in the training set $\Omega$. Letting ${\mathcal{X}}_{k}$ denote the $k$th iterate of Riemannian CG algorithm, we measure the errors on the training and the test set:

[TABLE]

The algorithm is stopped once these errors stagnate, that is,

[TABLE]

holds for some small $\delta>0$.

2.3.1 Adaptive rank and adaptive sampling strategy

To set up the optimization problem (8), two issues remain to be discussed: The choice of the TT ranks $\mathbf{r}$ and a suitable training set $\Omega$. For our application, these are not known a priori and thus need to be chosen adaptively.

Concerning the choice of TT ranks, we follow the adaptive strategy proposed in [45]. We start by solving (8) for the smallest sensible choice of TT ranks, $\mathbf{r}=(1,\dots,1)$. Most likely, this choice will not suffice to obtain satisfactory accuracy and the error on the test set will be relatively large. To decrease it, the obtained solution is used as starting value for Riemannian CG applied again to (8), but this time with the increased TT ranks $\mathbf{r}=(1,2,1,\dots,1)$ as discussed in [45]. See also [47] for a greedy rank update procedure in the context of matrix completion. The described procedure is repeated by increasing cyclically every TT rank $r_{\mu}$. The overall algorithm stops as soon as increasing any of the TT ranks does not improve the test set error anymore or the maximal possible rank $r_{\max}$ is reached; see Algorithm 1.

For the adaptive choice of the sampling set $\Omega$, which has not been addressed in [45], we present two different strategies. The core idea is to gradually increase the size of $\Omega$ in order to improve the approximation of the tensor. Both strategies are also combined with Algorithm 1 and they differ only in the measurement of the error.

The steps of the first adaptive sampling strategy are as follows.

Start with a sample set $\Omega$ of small size and a test set ${\Omega_{C}}$ of a certain prescribed size $|\Omega_{C}|$. Run Algorithm 1. 2. 2.

Measure the relative error on the test set ${\Omega_{C}}$ and stop if the stopping criterion is satisfied. If not satisfied, add the test set ${\Omega_{C}}$ to the sample set $\Omega$ and create a new test set of size $|\Omega_{C}|$. In our applications, this corresponds to computing new option prices on the Chebyshev grid using the reference method. 3. 3.

Run again Algorithm 1 from line 2 to the end, by using a rank $\mathbf{r}=(1,\dots,1)$ approximation of the result from the previous step as initial guess for the CG algorithm. 4. 4.

Repeat 1-3 until a maximal sampling percentage is reached or an a priori chosen stopping criterion is satisfied.

The pseudo-code in Algorithm 2 summarizes this first strategy.

The second adaptive sampling strategy that we propose is designed in a similar way. The only difference is that the error is measured on an a priori defined fixed set $\Gamma$ and not on $\Omega_{C}$, which changes at each step. Therefore, this strategy follows the same steps as the first one, with the only difference that in Step 2 we measure the error on the set $\Gamma$, which has been previously defined. The algorithm summarizing this second strategy can be obtained by replacing line 3 and line 12 in Algorithm 2 with

[TABLE]

The stopping criterion of line 13 can be also defined in different ways. We choose to stop the algorithm if one of the following criteria is satisfied:

if $\mathrm{err}_{\mathrm{new}}<\mathrm{tol}$, where $\mathrm{tol}$ is a prescribed tolerance; 2. 2.

if $|\mathrm{err}_{\mathrm{new}}-\mathrm{err}_{\mathrm{old}}|<\mathrm{tol}^{\prime}$ where $\mathrm{tol}^{\prime}$ is a prescribed tolerance; 3. 3.

if $\exists\mu$ such that $r_{\mu}({\mathcal{X}}_{c})==r_{\max}$.

The first criterion allows us to stop as soon as the error goes below a certain level, the second stops the algorithm whenever the error stagnates and the last one when the TT rank has reached the maximal allowed rank at least in one mode $\mu$.

We test the new adaptive sampling strategies on a problem with known solution in the next section.

2.3.2 Numerical test for adaptive sampling strategies

We consider the problem of Chapter 5.4.2 in [45] and we apply our adaptive sampling strategies to it in order to compare them and to investigate their advantages and disadvantages. We expect a similar performance of both strategies in terms of accuracy and compression. In this numerical example, as well as in the rest of the paper, we choose $\|\cdot\|$ to be the 2-norm and $\delta=10^{-4}$ in (9). The problem consists of discretizing the function

[TABLE]

using $n=20$ equally spaced discretization points on $[0,1]$ in each mode. We aim at reconstructing the tensor containing the function values in the grid. In Algorithm 2 we set the maximum rank to $\mathbf{r}_{\max}=(1,7,7,7,1)$ and we start with an initial sampling set $\Omega$ satisfying $|\Omega|/n^{4}=0.01$. Moreover, we set the acceptance parameter $\rho$ of Algorithm 1 to $\rho=10^{-4}$. In order to analyze the behavior of the error, we do not impose any stopping criterion, but we let our adaptive sampling strategies run until $|\Omega|/{\rm size}({\mathcal{A}})>0.25$. The size of each $\Omega_{C}$ is set to $2000$ and $|\Gamma|=3000$ for the second strategy. Figures 4 and 4 show the results for the two different strategies.

First, we observe that both strategies eventually reach the same accuracy and the same final TT ranks, which makes both of them valid. We observe an oscillatory behavior in Figure 4. This non-smooth decay can be expected since in each step the error is measured on a different test set $\Omega_{C}$. We observe that the amplitude of the oscillations becomes smaller as $|\Omega|$ increases. This indicates an error stagnation over the whole tensor which cannot be improved by enlarging $\Omega_{C}$ further. On the other hand, the error in the second strategy behaves almost monotonically and stagnates much earlier than in the previous case. This is due to the fact that we measure it on the fixed set $\Gamma$. In practice, the earlier error stagnation of the second strategy is preferable as it triggers the stopping criterion 2. However, the second strategy has the disadvantage of the initial additional cost of evaluating the tensor in the set $\Gamma$. In our numerical experiments in Section 3 we choose the first strategy, which turned out to be more favorable since the stopping criterion 1 was triggered.

2.4 Combined methodology

We are now in the position to combine the concepts and the algorithms in order to develop an efficient procedure for high-dimensional tensorized Chebyshev interpolation.

We would like to price options that depend on a vector ${\mathbf{p}}=(p_{1},\cdots,p_{d})$ of $d$ varying parameters. It is reasonable to assume that every combination of parameters ${\mathbf{p}}$ belongs to a compact hyper-rectangular $[\underline{p}_{1},\overline{p}_{1}]\times[\underline{p}_{2},\overline{p}_{2}]\times\cdots\times[\underline{p}_{d},\overline{p}_{d}]$. For example, if time-to-maturity $T$ belongs to the set of varying parameters, we can assume that $T\in[0.05,2]$; similarly for the other payoff or model parameters. The combined methodology consists of two phases: offline phase and online phase, as already introduced in [14].

2.4.1 Offline phase - Computation of ${\mathcal{P}}$

The offline phase starts by performing following operations:

Fix an interpolation order $\overline{\mathbf{n}}=(n_{1},\dots,n_{d})$ and compute the entries of the tensor ${\mathcal{P}}$ (as defined in (4)) from an a priori chosen subset $\Omega$ of Chebyshev nodes, using the reference pricing technique. 2. 2.

Apply tensor completion with adaptive sampling strategy (Algorithm 2) in order to get a low-rank approximation of the tensor ${\mathcal{P}}$ in TT format.

For simplicity, we denote the obtained low-rank approximation of ${\mathcal{P}}$ again by ${\mathcal{P}}$. In the last step of the offline phase we construct the interpolation coefficients, defined in (4). We denote the tensor of coefficients by $\mathcal{C}\in\mathbb{R}^{(n_{1}+1)\times(n_{2}+1)\times\cdots\times(n_{d}+1)}$. Its entries are therefore given by (adjusting the ordering according the Sections 2.1 and 2.2)

[TABLE]

for $i_{j}=1,\cdots,n_{j}+1$ and $j=1,\cdots,d$. The tensor ${\mathcal{C}}$ can be efficiently computed in TT format, as explained in the following subsection.

2.4.2 Offline phase - Efficient computation of $\mathcal{C}$

In order to explain the algorithm we first consider the simple case $d=1$. In this case ${\mathcal{P}}$ and ${\mathcal{C}}$ are in $\mathbb{R}^{(n_{1}+1)\times 1}$, where $n_{1}$ is the chosen interpolation order. The entries of ${\mathcal{C}}$ are given by

[TABLE]

so that the whole vector ${\mathcal{C}}$ can be computed via the matrix-vector multiplication

[TABLE]

and we denote by $F_{n_{1}}\in\mathbb{R}^{(n_{1}+1)\times(n_{1}+1)}$ the matrix multiplying ${\mathcal{P}}$ in (11).

For a general dimension $d>1$, the same reasoning can be applied and the tensor ${\mathcal{C}}$ of interpolation coefficients can be computed by sub-sequentially multiplying ${\mathcal{P}}$ with $F_{n_{i}}$ ($i=1,\cdots,d$) via the mode-$\mu$ multiplication, defined in Section 2.2. The final procedure for an efficient computation of ${\mathcal{C}}$ is given in Algorithm 3.

Note that if $n_{1}=\cdots=n_{d}=\mathrel{\mathop{:}}n$ (as for example in our numerical experiments in Section 3), Algorithm 3 can be further simplified by computing the matrix $F_{n}$ only once. The particular structure of the matrices $F_{n_{i}}$ allows us to apply a Fast-Fourier-Transform based algorithm which computes each mode multiplication in $O(r^{2}n\log(n))$ (instead of $O(r^{2}n^{2}$) as mentioned in Section 3). Therefore, the total complexity for computing ${\mathcal{C}}$ is $O(dnr^{2}\log(n))$.

The offline phase can be finally completed by performing the step

Construct the tensor $\mathcal{C}$ as explained in Algorithm 3.

2.4.3 Online phase

Once we have stored ${\mathcal{C}}$ in TT format, we can use it to compute every option price via interpolation during the online phase. For any particular choice of parameters ${\mathbf{p}}$, we first perform the step

Evaluate the Chebyshev tensor basis (3) in ${\mathbf{p}}$.

This step returns a tensor $\mathcal{T}_{\mathbf{p}}\in\mathbb{R}^{(n_{1}+1)\times(n_{2}+1)\times\cdots\times(n_{d}+1)}$ of TT rank $(1,\cdots,1)$, that we store in TT format. The interpolated price, defined in (2), can now be rewritten as the inner product

[TABLE]

The final step of our combined methodology is then defined as

Compute the interpolated price (12) in TT format as in (7).

If we consider a fixed interpolation order $n$ in each dimension and if the TT ranks of ${\mathcal{P}}$ and ${\mathcal{C}}$ are approximately $r$, then the total cost for performing both Step 3 and Step 5 is given by $O(dnr^{2}+dnr^{2}\log(n))$. These two steps are represented via a tensor network diagram in Figure 5 (for $d=5$), where we denoted by $\mathbf{P}_{i}$ the core tensors of ${\mathcal{P}}$ and by $\mathbf{T}_{i}$ the ones of $\mathcal{T}_{\mathbf{p}}$.

Finally, we summarize our complete methodology in Algorithm 4.

In the next section we see how this combined methodology performs on concrete examples.

3 Financial applications and numerical experiments

Putting the new approach to test, we implement the method described in Section 2 for two different types of applications. In the first one, we tackle computational intense option pricing methods in a parametric model. We treat option prices as functions in the parameter space which consists of model and option parameters. We then approximate the price function by Chebyshev interpolation in the parameter space. This approach has been successfully tested in cases where the parameter space is low-dimensional. In various applications, several varying parameters are of interest. If the interpolation is even efficient in the full parameter space, it is indeed a new pricing methodology. Here, we combine Chebyshev interpolation and low-rank approximation to cope with higher dimensionality in the parameter space. Already for pricing single asset options, it is promising to tackle medium and high-dimensional parameters spaces in this approach. As a generic example, we choose to approximate American put option prices in the Heston model with the varying parameters $K,\rho,\sigma,\kappa$ and $\theta$. It turns out that the computational complexity reduces significantly in this case.

As second type of application we examine the interpolation of basket option prices in the $d$-variate Black-Scholes model as function of the initial stock prices. This is a prototypical example for the computation of generalized conditional moments of high-dimensional Markov processes.

All algorithms have been implemented in Matlab and run on a standard laptop (Intel Core i7, 2 cores, 256kB/4MB L2/L3 cache). In order to deal with tensors, we used the toolboxes [39] by Oseledets and [2, 3], while for the completion algorithm we used the TT completion toolbox described in [32, 33, 45, 46]. Note that in this toolbox the most expensive steps have been implemented in C using the Mex-function capabilities of Matlab.

3.1 Pricing American options in Heston’s model

We consider pricing single asset American put options in the Heston model. As introduced by Heston in [24], the price dynamics of the financial asset under the risk neutral measure are given by

[TABLE]

where the square of the volatility $v_{t}$ is modeled by the square root process

[TABLE]

Here, the two Brownian motions $W^{1}$ and $W^{2}$ are correlated with correlation parameter $\rho$, mean-reversion rate $\kappa>0$, long-term mean $\theta>0$, volatility of the variance $\sigma>0$ and, finally, fixed and deterministic continuously compounding interest rate $r$.

The price of an American option at time $t<T$, maturing at $T$, with initial underlying price $s\geq 0$ and initial volatility $v\geq 0$ is given by

[TABLE]

where the $\sup$ is taken over all stopping times $\tau$ in $[t,T]$. Here, $f$ denotes the payoff function of the European put option, i.e.

[TABLE]

where $K$ denotes the strike price.

It is well-known (see e.g. [13]) that the price (13) of the American option satisfies the following partial differential complementarity problem (PDCP):

[TABLE]

where $\mathcal{G}$ is the infinitesimal generator of $(s,v)$ in the Heston model, defined as

[TABLE]

The problem (14) has been well studied in the literature and different pricing algorithms have been developed so far. In our example we consider, as reference method for our combined methodology, the pricing algorithm explained in [21]. More precisely, the authors propose different schemes for the time discretization and we consider the Hundsdorfer Verwer - Ikonen Toivanen (HV-IT) scheme, explained at page 219 of [21].

Solving the discretized PDCP yields an approximate price for all values of $S_{0},v_{0}$ and $T$ in each grid point of the pre-specified domain. For many applications we would like to have the solution at hand for other parameters (as well). In calibration, for instance, we observe $S_{0}$ and $r$, and one could estimate $v_{0}$ from historical stock price data. Then the calibration problem reduces to fitting the parameters $(K,\rho,\sigma,\kappa,\theta)$ to the observed option price data. To do so one needs to solve an optimization problem where prices need to be computed for large sets of parameters $(K,\rho,\sigma,\kappa,\theta,T)$. Since the price for different maturities can be obtained by rescaling $\kappa$ and $\sigma$, effectively we need the prices for combinations of the parameters $K,\rho,\sigma,\kappa$ and $\theta$. This motivates the following set up, where we fix the model and payoff parameters

[TABLE]

and we let vary the five parameters

[TABLE]

in their corresponding domain.

In order to compute the reference prices we consider $50$ equidistant spatial grid points in both directions $s$ and $v$ with $s_{\min}=0,s_{\max}=5,v_{\min}=0,v_{\max}=1$, $40$ time steps and the Crank-Nicholson time stepping scheme.

We start by performing the offline phase of Algorithm 4. We consider an interpolation order $n_{1}=\cdots=n_{5}=\mathrel{\mathop{:}}n=10$ in each direction and we construct the tensor ${\mathcal{P}}$ by tensor completion as explained in Section 2.3. We apply the first adaptive sampling strategy as in Algorithm 2. We choose the completion parameters as

[TABLE]

For this particular example, we were also able to explicitly construct the full tensor (in more than 1 hour and 40 minutes!). In Table 1 we show the size of the final set $\Omega$ (first column), the relative error of the completed tensor on the last $\Omega_{c}^{new}$ (second column), the relative error between the obtained completed tensor and the full one (third column), the runtime of the completion, Algorithm 2, in seconds (fourth column), the TT-rank of ${\mathcal{P}}$ (fifth column), the storage needed to save ${\mathcal{P}}$ in TT format, denoted by store(TT) and measured in bytes (sixth column) and finally, the storage needed to save the full tensor, denoted by store(full) and again measured in bytes. Matlab requires $8$ bytes to store a floating-point number of type double, which gives us the formula store(full)$=8\cdot(n+1)^{d}$ for the storage of the full tensor and store(TT)$=8\cdot(n+1)(r_{1}r_{2}+\cdots+r_{d-2}r_{d-1})+8\cdot(n+1)(r_{1}+r_{d-1})$ for the storage of the tensor in TT format, see [39].

Table 1 shows that a sample set of $5\%$ is sufficient for the algorithm to reach the prescribed accuracy. Furthermore, the relative error of the completed tensor ${\mathcal{P}}$ in the 2-norm over the last test sample parameter space $\Omega^{new}_{c}$ and the relative error over the full tensor, i.e. over all Chebyshev nodes, is only in the $6$th digit. This is one order of magnitude smaller than the relative error on the full ${\mathcal{P}}$. This is a good indication that the approach can be extended to more complex cases, where the computation of the full tensor $\mathcal{P}$ is not feasible any more (see Section 3.2). The completion time was about $6$ minutes. Finally, the rank properties together with its storage reduction of a factor of $115$ confirm the low-rank structure of the problem.

For constructing the tensor $\mathcal{C}$ (last step of the offline phase) we applied Algorithm 3 and the computation time was $0.0037$ seconds, which is negligible compared to the completion time. Hence, almost all the computation time in the offline phase is spent in the construction of the tensor ${\mathcal{P}}$.

Next, we compute American put option prices for the online phase in both ways using our methodology and the reference algorithm. We compute $243$ prices with random model parameters uniformly drawn from the reference set $[2;4]\times[-1;1]\times[0.2;0.5]\times[1;2]\times[0.05;0.2]$. We measure the maximal absolute error over the computed options prices, i.e. we report the quantity

[TABLE]

where $P_{\text{Int}}$ is a vector containing all interpolated prices for the different choices of model parameters; analogously is $P_{\text{Ref}}$ for the reference method. In Table 2 we also report the computation time for computing one single option price for both methods. One can notice that the online phase of the interpolation compared to the reference method accelerates the procedure by a factor of $75$. The accuracy of the reference method is reported in part C of Figure 1 in [21] for one specific parameter set to be of the order $10^{-3}$ in the maximum norm. The interpolation error is one order smaller, making the new procedure at least as accurate as the reference method. Therefore, we can conclude that the methodology strongly outperforms the reference method in the online phase while keeping the same accuracy.

We would like to emphasize that this approach can be further extended to an interpolation in the full set of parameters $(S_{0},v_{0},r,T,K,\rho,\sigma,\kappa,\theta)$. Since then the offline phase needs to be performed only once, this would result in a new pricing method. Here, in the offline phase one could explore the fact that the PDCP solver returns the price for all $(S_{0},v_{0},T)$ in the grid to make the sampling steps more efficient. This opens up an interesting topic for future research.

3.2 Basket options in multivariate Black-Scholes model

In the $d$-variate Black-Scholes model with $d$ assets $S^{1},\cdots,S^{d}$, the risk neutral dynamics are given by

[TABLE]

where $r$ is a fixed deterministic interest rate, $(\sigma_{1},\cdots,\sigma_{d})$ is the vector of volatilities and $(W^{1},\cdots,W^{d})$ is a vector of correlated Brownian motions with correlation matrix $\Sigma$. The solution to (15) is given by

[TABLE]

In this section we apply the new methodology in order to price basket options with payoff function $f:\mathbb{R}^{d}\to\mathbb{R}$ defined as

[TABLE]

where $K$ is the strike and $(w_{1},\cdots,w_{d})$ is a vector of weights satisfying $\sum_{n=1}^{d}w_{n}=1$. The risk neutral price at time $t=0$ of the basket option with maturity $T$ is, as usual, given by

[TABLE]

From now on, we consider the parameters $r,\sigma_{i}$ $(i=1,\cdots,d)$ and the correlation matrix $\Sigma$ to be fixed, and we let the vector $\mathbf{S_{0}}\in\mathbb{R}^{d}$ of initial asset prices be the varying parameter. The reference pricing algorithm will be of Monte Carlo (MC) type combined with a variance reduction technique. In particular, we use the control variates method presented in [17], where the control variate is given by

[TABLE]

Since the only varying parameter is the vector of initial asset prices, it is very convenient to split the Monte Carlo simulation in two parts in order to make the completion more efficient. More precisely, in a pre-computation phase (Algorithm 5) we simulate a certain number of realizations (e.g. $10^{4}$) of

[TABLE]

and in a second moment we multiply the vector $\mathbf{S_{0}}$ (for all required parameter combinations) with all the realizations and we compute the Monte Carlo price by applying the chosen variance reduction technique (Algorithm 6). In order to generate the correlated random variables $W_{T}^{i}$, we use the Cholesky factorization of the correlation matrix, which is then multiplied by a vector of independently generated standard normal distributed random variates. Note that $\circ$ in Algorithm 6 represents the Hadamard (component-wise) product between vectors.

Algorithm 5 is executed at the beginning of the whole procedure and Algorithm 6 whenever needed in later stages. The advantage of splitting the MC algorithm is twofold. Firstly, it supports a considerable gain in efficiency in the performance of the completion algorithm: When we adaptively increment the sampling set $\Omega$ (which consists of sampling Chebyshev nodes in $\mathbf{S}_{0}$) in Algorithm 2, we need to compute new prices in the Chebyshev grid, which can be done by using Algorithm 6 only. The second advantage regards the analysis of the methodology and the completion accuracy: Since we use the same set of simulations for every Chebyshev price, the MC simulation does not introduce any further error to the completion. Moreover, we will see in Section 3.2.3 that this splitting procedure allows for a qualitative analysis of the rank structure of ${\mathcal{P}}$.

Next, we perform numerical experiments for different settings of model parameters, first for uncorrelated then for correlated assets.

3.2.1 Basket options of uncorrelated assets

In this example we consider the special case of uncorrelated assets. We investigate the performance of the proposed method for two different interpolation orders $n_{1}=\cdots=n_{d}=\mathrel{\mathop{:}}n=4$ and $n_{1}=\cdots=n_{d}=\mathrel{\mathop{:}}n=6$. We apply the combined methodology (Algorithm 4) to portfolios consisting of $d\in\{5,10,15,20,25\}$ assets. The set of fixed parameters is given by

[TABLE]

where $I_{d}$ denotes the $d\times d$ identity matrix. We let $\mathbf{S}_{0}$ vary in the hyper-rectangular

[TABLE]

so that we consider ITM options and ATM options as well.

For each value of $d$, we start by performing Algorithm 5 with $NumberSim=10^{3}$ for $n=4$ and with $NumberSim=10^{4}$ for $n=6$. In a second moment we construct the tensor ${\mathcal{P}}$ by applying the tensor completion with the adaptive sampling strategy of Algorithm 2 (first strategy). Table 3 shows the completion parameters for each value of $d$ and each interpolation order. The results of the tensor completions are displayed in Table 4. As in the previous subsection, we report the final size of the set $\Omega$, the relative error measured on the last set $\Omega_{C}^{new}$, the completion time and the memory needed to store both the obtained tensor in TT format and the full tensor. For the TT ranks of the completed tensor, we do not report the full tuple $(r_{0},\cdots,r_{d})$ (see Definition (5)) but only the quantity $\max_{\mu\in\{0,\cdots,d\}}r_{\mu}$.

It is interesting to analyze the size of the finally obtained set $\Omega$ in Algorithm 2 for different values of $d$ and $n$ (different sizes of ${\mathcal{P}}$). Figure 6 shows a plot of $|\Omega|$ (final) against $d$ for the two chosen interpolation orders. The graphical representation clearly suggests that the number of sampled entries, i.e. $|\Omega|$, required for the chosen tolerance $tol=10^{-2}$ for a fixed interpolation order $n=4$ and $tol=10^{-3}$ for a fixed $n=6$ is roughly of $O(d^{2})$, whereas the size of the full tensor is $n^{d}$. On the practical side, this means that by the completion algorithm we can reduce the complexity of the first step of the offline phase from an exponential growth down to a quadratic growth in the dimensionality. The exponential growth typically is referred to as curse of dimensionality. The reduction in absolute numbers is already tremendous for $d=5$ and $n=4$, where we observe $|\Omega|=124$ and the full tensor size equals $(n+1)^{d}=3125$. The compression is dramatic for $n=6$ and $d=25$, namely the numbers of required entries shrinks by a factor of more than $3\times 10^{17}$.

As in the previous numerical example, the computation time to build the tensor $\mathcal{C}$ of interpolation coefficients is negligible in the offline phase. Indeed, for all choices of $d$ and $n$ it is less than $0.01$ seconds, for instance $0.0045$ seconds for $n=4$ and $d=5$, and $0.0095$ seconds for $n=6$ and $d=25$.

We now perform the online phase of Algorithm 4 in order to see how efficient becomes pricing basket options in the new setting. We start by computing 100 basket option prices via Chebyshev interpolation (combined methodology), choosing random initial asset prices $\mathbf{S_{0}}$ in the reference hypercube $[1;1.5]^{d}$. We then compare the obtained prices with reference prices computed by applying the reference method (Monte Carlo with control variates) with $10^{4}$ new simulations for $n=4$ and $10^{5}$ new simulations for $n=6$. In particular, we measure again the maximal absolute error over all computed prices

[TABLE]

where $P_{\text{Int}}$ is a vector containing all 100 interpolated prices for the different choices of $\mathbf{S_{0}}$; analogously is $P_{\text{Ref}}$ for the reference method. The errors together with the computational times are shown in Table 5. Note that we report again the computational time to compute one single option price.

One can see that the online phase of the new procedure compared to the MC reference method accelerates the computation of a factor between 200 and 400 for $n=4$ and of a factor between 2000 and 4000 for $n=6$. Note that the difference in the acceleration between the two chosen interpolation orders is given by the different numbers of simulations chosen in the MC reference method ($10^{4}$ for $n=4$ and $10^{5}$ for $n=6$). Therefore, for both interpolation orders and for all choices of $d$, the acceleration is dramatic. In order to judge the accuracy of our method we have computed the $95\%$ confidence interval of the reference method, which results to be of a size between $10^{-4}$ and $5\cdot 10^{-4}$ for all choices of $\mathbf{S}_{0}$ and $d$ or $n$. This, together with the last column of Table 5, leads us to the conclusion that the new method is as accurate as the reference MC algorithm.

Finally, in Figure 7 we show the gain in efficiency of the new method when computing basket option prices for $d=25$ and both choices of interpolation orders. In particular, on the x-axis we consider a possible number of computed prices and on the y-axis we present

the computational time of the reference MC method, 2. 2.

the computational time of the new combined methodology (offline phase + online phase ),

required to compute the corresponding amount of prices.

The plots in Figure 7 show that after an initial investment the computational time grows very slowly in the number of computed prices for the new method. This is due to the fact that the online phase in Algorithm 4 is very cheap, as shown in the numerical experiments. This proves that the method is useful whenever one can split the task in a pre-computational phase during idle times and a run-time phase where execution is required to be fast. Moreover, it will outperform the reference methods if a large number of prices needs to be computed. The first plot in Figure 7 indicates that for the case $n=4$ it is convenient to use the reference MC method if we want to compute up to $1000$ option prices. For the case $n=6$ the break-even point is already reached with $500$ prices.

3.2.2 Basket options of correlated assets

In this second numerical experiment we repeat the test of the previous subsection but, this time, we consider correlated assets. In particular, we choose again the interpolation orders $n=4$, $n=6$ and the other parameters are given by

[TABLE]

where $R_{d}$ denotes a random correlation matrix. The free parameters $S_{0}^{i}$, $i=1,\cdots,d$ are again contained in [1;1.5]. We perform the offline phase by considering again the set of completion parameters listed in Table 3. The obtained results of the completion are now in Table 6 and Figure 8 shows the required size of $\Omega$ to go below the tolerance $tol=10^{-2}$ for $n=4$ and $tol=10^{-3}$ for $n=6$. We notice that the completion results are similar to the case of uncorrelated assets and that $|\Omega|$ scales again like $O(d^{2})$. The computational time to construct $\mathcal{C}$ was again measured to be less than $0.01$ seconds for all choices of $d$ and $n$.

The online phase is performed similarly to the previous chapter, in particular we compute again 100 prices using the new method and the reference one. The MC parameters are set as before and the results are shown in Table 7.

The performance of the new method in terms of accuracy and computational efficiency is similar to the one observed in the case of uncorrelated assets. To summarize, the new methodology achieves a very good performance for uncorrelated as well as for correlated assets.

3.2.3 Rank structure of ${\mathcal{P}}$

In this section we qualitatively analyze the rank structure of the tensor ${\mathcal{P}}$. For simplicity, we perform this analysis for the standard Monte Carlo approach (without any variance reduction technique). Assume that we have already simulated the realizations of the correlated geometric Brownian motions stored in the matrix $M$ (Algorithm 5). Then, the price in the point $\mathbf{S}_{0}$ is given by the function

[TABLE]

where $D$ is the hyper-rectangular domain for the interpolation, $M(n,:)$ is the $n$-th row of $M$ and $N_{S}$ is the number of Monte Carlo simulations. This expression can be rewritten in the form

[TABLE]

where the $\alpha_{i}(n)$’s are coefficients multiplying $S_{0}^{i}$ depending on the $n$-th simulation and on the $i$-th weight $\omega_{i}$. The function $p$ is piecewise affine in the variables $S_{0}^{i}$.

To explore the rank structure of ${\mathcal{P}}$ let us consider the case of a single Monte Carlo simulation $N_{S}=1$. Then $p$ is of the form

[TABLE]

Now we analyze three different cases. First, consider the case where the price is positive for any $\mathbf{S}_{0}$ in the hyper-rectangular $D$. Here, $p$ is affine. This implies that the TT ranks are bounded by $d$. This follows from the fact that the CP rank (rank of the Canonical Polyadic Decomposition, see [31]) of ${\mathcal{P}}$, which is an upper bound for each $r_{\mu}$ in the TT ranks (see [20]), is equal to $d$. Second, if we observe a vanishing price for all $\mathbf{S}_{0}$ in the hyper-rectangular, then ${\mathcal{P}}$ is the zero-tensor, which has rank [math]. These two cases obviously yield a low-rank structure of ${\mathcal{P}}$, a favorable case for the new combined methodology.

In the third case where $p$ is only piecewise affine the situation is more complex and to gain an intuition we consider the case $d=2$, where $p$ is of the form

[TABLE]

on a squared domain $D$. Now, define the set

[TABLE]

When $L$ intersects the domain $D$ it cuts it in two regions. Only if $\alpha_{1},\alpha_{2}$ and $K$ are of a specific form that leads L to be the diagonal of $D$, the rank of ${\mathcal{P}}$ is almost full. In the Monte Carlo simulation context, this special case is very unlikely. In all other cases, ${\mathcal{P}}$ exhibits a lower rank structure. In particular, we expect the rank to be the lower the more the sizes of the two regions differ.

In order to visualize these findings we consider three different pairs $(\alpha_{1},\alpha_{2})$ together with $r=0$, $K=1$ and evaluate the corresponding $p$ on the discretized $D=[1;1.5]^{2}$ using 50 equidistant points in each direction. Figure 9 shows the sparsity pattern and the rank of the obtained matrices ${\mathcal{P}}$.

This qualitative explanation indicates that the rank structure of ${\mathcal{P}}$ depends on $D$. We expect the rank to be lower for domains $D$ with an asymmetry with respect to the strike $K$. Next we construct ${\mathcal{P}}$ as in the experiments of Section 3.2 for $K=1$, $d=2$ and different interpolation orders $n$ for both $D=[0.5;1.5]^{2}$ and $D=[1;1.5]^{2}$. In particular, we first construct the matrix $M$ via Algorithm 5 with $10^{5}$ simulations and subsequently compute ${\mathcal{P}}$ using Algorithm 6. In Figure 10 we display the decay of the singular values for all treated cases. As expected, the decay is faster for $D=[1;1.5]^{2}$. However, also for $D=[0.5;1.5]^{2}$ the decay of the singular values is reasonably fast. This implies that the new methodology would be still beneficial in this case.

4 Summary and future work

We have presented a unified approach to efficiently compute parametric option prices. The starting point of our methodology was the Chebyshev interpolation technique developed in [14], which we briefly summarized in Section 2.1. We refined both the offline and the online phase to treat high-dimensional problems with parameter spaces up to dimension $25$. We have exploited the low-rank structure of the tensors involved in the interpolation procedure, which have been stored in TT format (summarized in Section 2.2). In particular, we have developed a completion technique (explained in Section 2.3) which allows us the construct the tensor ${\mathcal{P}}$, containing the option prices in the Chebshev tensor grid. All ingredients have been efficiently assembled to finally build a combined methodology, explained in Section 2.4.

In the second part of the paper, Section 3, we have tested our approach in two different concrete option pricing settings: We have treated the American option pricing problem in the Heston model (Section 3.1) and the European basket option pricing problem in the $d$-dimensional Black and Scholes model (Section 3.2). Both examples show that our approach allows for a substantial gain in efficiency, while maintaining very accurate results, whose precision is comparable to the one of the considered reference methods. For instance, the interpolation of American option prices in $5$ parameters accelerates the procedure by a factor of $75$, when compared to the FD reference method [21]. For basket option pricing with $25$ underlyings the efficiency gain reaches factors up to $4000$. See Tables 2, 5, 7 and Figure 7 for further results. Finally, for both examples we qualitatively investigated the rank structure of ${\mathcal{P}}$, which confirmed that our initial low-rank assumption was indeed reasonable. For instance, for the American put, we obtain a compression factor of $115$ of the completed tensor ${\mathcal{P}}$ with respect to the full one, with a relative error in the $5$th digit only, see Table 1. For the basket option the full tensor containing prices in the Chebyshev grid is too large to be computed, however in Section 3.2.3 it is qualitatively explained why ${\mathcal{P}}$ is expected to have a low-rank structure. This is also confirmed by the compression rates observed in Tables 4 and 6 that go up to $3\times 10^{17}$.

Seen the promising performance of this new approach and considering the fact that this methodology can be easily tailored to different problem settings, we expect it to be applicable in several domains in finance. For instance, pricing, calibration and sensitivity analysis in equity markets, fixed income and credit, and parameter uncertainty quantification are some of the possible domains of application.