A gradient-based calibration method for the Heston model

Abstract

The Heston model is a well-known two-dimensional financial model. Because the Heston model contains implicit parameters that cannot be determined directly from real market data, calibrating the parameters to real market data is challenging. In addition, some of the parameters in the model are non-linear, which makes it difficult to find the global minimum of the optimization problem within the calibration. In this paper, we present a first step towards a novel space mapping approach for parameter calibration of the Heston model. Since the space mapping approach requires an optimization algorithm, we focus on deriving a gradient descent algorithm. To this end, we determine the formal adjoint of the Heston PDE, which is then used to update the Heston parameters. Since the methods are similar, we consider a variation of constant and time-dependent parameter sets. Numerical results show that our calibration of the Heston PDE works well for the various challenges in the calibration process and meets the requirements for later incorporation into the space mapping approach. Since the model and the algorithm are well known, this work is formulated as a proof of concept.

Mathematics Subject Classifications:

• 91G50
• 65M06

1. Introduction

In finance, calibrating model parameters to fit real market data is challenging because most model parameters are implicit in the real market data [4,8,12,14,15,17]. We consider the well-known two-dimensional Heston model, which contains at least four parameters implicit in the market data. The Heston proposed a two-dimensional stochastic differential equations (SDE) model to simulate the behaviour of the stock price [6] and presented a closed-form valuation formula for his model. Some calibration techniques are based on this formula [4,15]. The closed-form equation has some restrictions w.r.t. improving the model by considering non-constant parameters [15]. Another strategy to improve the model is to introduce additional processes, which leads to an increase in the difficulty of the calibration process [17]. Since the stock price and variance are stochastic processes in the Heston model, Monte Carlo optimization methods can be applied [17]. New approaches also consider deep learning strategies and neural networks to calibrate the Heston model [8,11–13]. Another recent approach is to improve the calibration by using faster methods, such as a parallel GPU implementation [5]. We consider the well-known two-dimensional Heston model, which contains at least four parameters implicit in the market data. The parameter calibration is formulated as a constrained optimization problem to minimize a cost functional. The cost functional describes the difference between the reference data, the subsequent market data, and the data obtained by numerically solving our model. In the space mapping approach [1], a coarse and a fine optimization solver are used. This paper presents the first step ”a gradient descent algorithm” for the Heston model towards a space mapping approach in finance. Later, within space mapping, the gradient descent algorithm will compute a coarse and cheap approximation for the calibration problem of the Heston model.

Our goal is to introduce the space mapping approach to financial applications, and therefore we use the log-transformed normalized Heston model formulated as a partial differential equation (PDE) and the gradient descent algorithm as a pre-step, since both are well known and thus we can focus on the novel aspects. During our research, we didn't find any other PDE-based calibration approach for the Heston model. Of course, we are aware that there are already algorithms to compute the exact solution of the Heston calibration problem [4]. However, these methods are limited to the assumption that the parameters are constant, whereas our approach considers time-dependent parameters and constant parameter calibration. Furthermore, if the parameters are assumed to be time dependent, an analytical solution can't be derived and thus an exact solution is not available. Therefore, as mentioned before, the paper is considered as a proof of concept to introduce a new calibration method based on financial research.

We focus on deriving a gradient descent algorithm for the Heston model that satisfies the requirements for its use in the space mapping approach. In addition to our financial model, the gradient descent algorithm requires the adjoint of the model. Therefore, we formally derive the adjoint of the Heston model and construct the gradient using well-known techniques from optimization of partial differential equations (PDEs) [7,18]. The gradient descent algorithm has previously been applied to the Heston model using neural networks, so it has not been explicitly computed [13]. In this work, we focus on the calibration rather than the numerical solution of the model, i.e. our obtained results can be further improved by using more accurate numerical methods. Nevertheless, the calibration results obtained are remarkable. Our attempt can be easily extended to include time-dependent parameters for the variance, thus providing an improvement over current calibration strategies based on the semi-analytical solution of the Heston model.

The rest of the article is organized as follows. In Section 2, we introduce the Heston model and our approach to solving the resulting log-transformed Heston PDE. We then focus on parameter calibration. In Section 3, we derive the adjoint of the log-transformed Heston model and the corresponding gradient. We then discretize the method in Section 4 Finally, in Section 5 we focus on numerical results for the application of the gradient descent algorithm. We analyse the behaviour of the algorithm in various problems for constant and time-dependent parameter calibration. The paper ends with a conclusion and an outlook in Section 6.

2. The Heston model

The Heston model was proposed by Heston in 1993 [6] and describes the dynamics of the underlying asset $S\in {\mathbb{R}}^{>0}$ and the variance $\nu \in {\mathbb{R}}^{>0}$ by a two-dimensional SDE. This model is an extension of the well-known Black-Scholes model, which considers only a stochastic process for the asset and leads to (1) $\mathrm{d}{S}_t=r{S}_t\mathrm{d}t+{\sigma }_S{S}_t\mathrm{d}{W}_t^S,S_0>0,$(1) where $r\in {\mathbb{R}}^{>0}$ is the risk-free rate, ${\sigma }_S\in {\mathbb{R}}^{>0}$ is the volatility, and $W_t^S$ is the Brownian motion. Heston introduced a second stochastic process for volatility by introducing ${\sigma }_S=\sqrt{\nu }$. He assumed that volatility follows an Ornstein–Uhlenbeck process, which can be rewritten as a Cox–Ingersoll process for the variance ν. This leads to the SDE system of Heston's model under the risk-neutral measure given by (2) $\{\begin{array}{ll}\mathrm{d}{S}_t=r{S}_t\mathrm{d}t+\sqrt{{\nu }_t}{S}_t\mathrm{d}{W}_t^S,& S_0>0,\\ \mathrm{d}{\nu }_t={\kappa }_{\nu }({\mu }_{\nu }-{\nu }_t)\mathrm{d}t+{\sigma }_{\nu }\sqrt{{\nu }_t}\mathrm{d}{W}_t^{\nu },& {\nu }_0>0,\end{array},$(2) where ${\kappa }_{\nu }\in {\mathbb{R}}^{>0}$ is the mean reversion rate, ${\mu }_{\nu }\in {\mathbb{R}}^{>0}$ is the long term mean, and ${\sigma }_{\nu }\in {\mathbb{R}}^{>0}$ is the volatility of variance. The Brownian motions $W_t^S$ and $W_t^{\nu }$ are correlated by the constant $\rho \in [-1,1]$ [6]. For the variance process to be positive, the Feller condition $2{\kappa }_{\nu }{\mu }_{\nu }\ge {\sigma }_{\nu }^2$ must be satisfied. If the Feller condition is violated, the (2(2) $\{\begin{array}{ll}\mathrm{d}{S}_t=r{S}_t\mathrm{d}t+\sqrt{{\nu }_t}{S}_t\mathrm{d}{W}_t^S,& S_0>0,\\ \mathrm{d}{\nu }_t={\kappa }_{\nu }({\mu }_{\nu }-{\nu }_t)\mathrm{d}t+{\sigma }_{\nu }\sqrt{{\nu }_t}\mathrm{d}{W}_t^{\nu },& {\nu }_0>0,\end{array},$(2) ) becomes complex which leads to computational problems. Using Kolmogorov's backward equation, we derive the Heston PDE under the risk-neutral measure (3) $0=V_t+\frac{\nu }{2}{S}^2{V}_{\mathrm{SS}}+\frac{1}{2}{\sigma }_{\nu }^2\nu V_{\mathrm{\nu \nu }}+\mathrm{rS}{V}_S+{\kappa }_{\nu }({\mu }_{\nu }-\nu )V_{\nu }+{\sigma }_{\nu }\mathrm{\nu \rho S}{V}_{\mathrm{S\nu }}-\mathrm{rV},$(3) where $V(S,\nu ,t)$ denotes the fair price of the option. If we look at the European plain vanilla put option, the terminal condition (‘pay-off’) is as follows (4) $V(S,\nu ,T)=(K-S,0{)}^{+},$(4) where $K\in {\mathbb{R}}^{>0}$ denotes the predefined strike price. Heston proposed the following boundary conditions (5) $\begin{array}{rl}V(0,\nu ,t)& =\mathrm{Kexp}(-r(T-t)),\end{array}$(5) (6) $\begin{array}{rl}V(S,\nu ,t)& =0,\text{for }S\to \mathrm{\infty },\end{array}$(6) (7) $\begin{array}{rl}V(S,0,t)& =\mathrm{rs}{V}_S(S,0,t)+{\kappa }_{\nu }{\mu }_{\nu }{V}_{\nu }(S,0,t)-\mathrm{rV}(S,0,t),\end{array}$(7) (8) $\begin{array}{rl}V(S,v,t)& \sim \mathrm{Kexp}(-r(T-t)),\text{ for }\nu \to \mathrm{\infty }.\end{array}$(8) The boundary conditions for the underlying follow directly from the assumptions on the financial market. We reverse the time direction by introducing $\tau =T-t$ and thus the payoff condition (4(4) $V(S,\nu ,T)=(K-S,0{)}^{+},$(4) ) becomes an initial condition. Next, it is advantageous to use the variable transformation $x=\log (S/K)$ for the asset, as it simplifies the numerical implementation, i.e. we obtain a log-transformed variable with normalization to the strike price K. We rewrite the option price equation as $\tilde{V}(x,\nu ,\tau )=V(S,\nu ,t)/K$ and obtain the so-called log-transformed normalized Heston PDE (9) ${\tilde{V}}_{\tau }=\frac{\nu }{2}{\tilde{V}}_{\mathrm{xx}}+\frac{1}{2}{\sigma }_{\nu }^2\nu {\tilde{V}}_{\mathrm{\nu \nu }}+(r-q-\frac{\nu }{2}){\tilde{V}}_x+{\kappa }_{\nu }({\mu }_{\nu }-\nu ){\tilde{V}}_{\nu }+{\sigma }_{\nu }\mathrm{\nu \rho }{\tilde{V}}_{\mathrm{x\nu }}-r\tilde{V},$(9) defined on the semi-unbounded domain $x\in \mathbb{R}$, $\nu \ge 0$, $0\le \tau \le T$ and supplied with the initial condition (10) $\tilde{V}(x,\nu ,\tau )=(1-\exp (x),0{)}^{+}$(10) and the following boundary conditions (11) $\begin{array}{rl}\tilde{V}(x,\nu ,\tau )& \sim \exp (-\mathrm{r\tau }),\text{ for }x\to -\mathrm{\infty },\end{array}$(11) (12) $\begin{array}{rl}\tilde{V}(x,\nu ,\tau )& =0,\text{ for }x\to \mathrm{\infty },\end{array}$(12) (13) $\begin{array}{rl}\tilde{V}(x,\nu ,\tau )& \sim \exp (-\mathrm{r\tau }),\text{ for }\nu \to \mathrm{\infty }.\end{array}$(13) The boundary conditions for the underlying asset are the normalized log-transformation of Heston's proposed boundary conditions (5(5) $\begin{array}{rl}V(0,\nu ,t)& =\mathrm{Kexp}(-r(T-t)),\end{array}$(5) ), (6(6) $\begin{array}{rl}V(S,\nu ,t)& =0,\text{for }S\to \mathrm{\infty },\end{array}$(6) ) and (8(8) $\begin{array}{rl}V(S,v,t)& \sim \mathrm{Kexp}(-r(T-t)),\text{ for }\nu \to \mathrm{\infty }.\end{array}$(8) ). At $\nu =0$ the parabolic PDE degenerates to a first-order hyperbolic PDE, and therefore we need to consider the Fichera theory [2,10] to assess whether it is necessary to provide these analytic boundary conditions or not.

From the Fichera condition at $\nu =0$ given by (14) $b(\nu )={\kappa }_{\nu }({\mu }_{\nu }-\nu )-\frac{1}{2}{\sigma }_{\nu }^2$(14)

it follows

• if $\underset{\nu \to 0^{+}}{lim}b(\nu )\ge 0$ (outflow) we must not supply any boundary conditions at $\nu =0$.

• if $\underset{\nu \to 0^{+}}{lim}b(\nu )<0$ (inflow) we have to supply boundary conditions at $\nu =0$.

In the sequel, we assume that the Feller condition $2{\kappa }_{\nu }{\mu }_{\nu }\ge {\sigma }_{\nu }^2$ is satisfied. Therefore we obtain an outflow boundary at $\nu =0$ and must not impose any analytic boundary conditions on this boundary. But we obviously need a numerical boundary condition to complete the scheme, which will be later discussed in Section 4.

In addition to the Heston PDE with constant parameters for the variance process, we also consider the time-dependent parameters ${\tilde{\kappa }}_{\nu }$, ${\tilde{\mu }}_{\nu }$, ${\tilde{\sigma }}_{\nu }$, as well as $\tilde{\rho }$. Then the corresponding PDE formulation is given by (15) ${\tilde{V}}_{\tau }=\frac{\nu }{2}{\tilde{V}}_{\mathrm{xx}}+\frac{1}{2}{\tilde{\sigma }}_{\nu }^2\nu {\tilde{V}}_{\mathrm{\nu \nu }}+(r-q-\frac{\nu }{2}){\tilde{V}}_x+{\tilde{\kappa }}_{\nu }({\tilde{\mu }}_{\nu }-\nu ){\tilde{V}}_{\nu }+{\tilde{\sigma }}_{\nu }\nu \tilde{\rho }{\tilde{V}}_{\mathrm{x\nu }}-r\tilde{V}.$(15) In the numerical simulations 5, we will also discuss variations of the constant and time-dependent parameter.

3. Derivation of a gradient-based optimization strategy

For a given data set $V_{\mathrm{data}}$ and reference parameters ${\xi }_{\mathrm{ref}}$, we formulate the calibration as optimization task with the cost functional given by $J(\tilde{V},\xi )=\frac{1}{2}{\int }_0^T\Vert K\tilde{V}-V_{\mathrm{data}}{\Vert }^2\mathrm{d}\tau +\frac{\lambda }{2}\Vert \xi -{\xi }_{\mathrm{ref}}{\Vert }^2.$As our aim is to fit the parameter to real market data and therefore no ${\xi }_{\mathrm{ref}}$ exists, we set $\lambda =0$ and the cost functional reduces to (16) $J(\tilde{V},u)=\frac{1}{2}{\int }_0^T\Vert K\tilde{V}-V_{\mathrm{data}}{\Vert }^2\mathrm{d}\tau .$(16) In the following, we derive a gradient-based algorithm that allows us to calibrate the parameters numerically. To this end we use a Lagrangian approach.

3.1. First-order optimality conditions for the Heston model

Let us denote the Lagrange multipliers by $\psi =(\phi ,{\phi }^a,{\phi }^b,{\phi }^c,{\phi }^d)$, set $\mathrm{\Omega }=(0,{\nu }_{max})\times (x_{min},x_{max})$ and split the boundary $\mathrm{\partial \Omega }$ into $\begin{array}{rl}& {\mathrm{\Gamma }}_a=\mathrm{\partial \Omega }\cap \{x=-\mathrm{\infty }\},{\mathrm{\Gamma }}_b=\mathrm{\partial \Omega }\cap \{x=\mathrm{\infty }\},\\ & {\mathrm{\Gamma }}_c=\mathrm{\partial \Omega }\cap \{\nu =0\},{\mathrm{\Gamma }}_d=\mathrm{\partial \Omega }\cap \{\nu =\mathrm{\infty }\}.\end{array}$First, we focus on the log-transformed normalized Heston equation with constant parameters (9(9) ${\tilde{V}}_{\tau }=\frac{\nu }{2}{\tilde{V}}_{\mathrm{xx}}+\frac{1}{2}{\sigma }_{\nu }^2\nu {\tilde{V}}_{\mathrm{\nu \nu }}+(r-q-\frac{\nu }{2}){\tilde{V}}_x+{\kappa }_{\nu }({\mu }_{\nu }-\nu ){\tilde{V}}_{\nu }+{\sigma }_{\nu }\mathrm{\nu \rho }{\tilde{V}}_{\mathrm{x\nu }}-r\tilde{V},$(9) ) and define $\begin{array}{r}A=\frac{1}{2}\nu \left(\begin{array}{cc}{\sigma }_{\nu }^2& {\sigma }_{\nu }\rho \\ {\sigma }_{\nu }\rho & 1\end{array}\right),b=\left(\begin{array}{c}{\kappa }_{\nu }({\mu }_{\nu }-\nu )-\frac{1}{2}{\sigma }_{\nu }^2\\ r-q-\frac{\nu }{2}-\frac{1}{2}{\sigma }_{\nu }\rho \end{array}\right)\tilde{b}=\left(\begin{array}{c}{\kappa }_{\nu }{\mu }_{\nu }\\ r-q\end{array}\right).\end{array}$ Then it can be written as ${\tilde{V}}_{\tau }-\mathrm{\nabla }\cdot A\mathrm{\nabla }\tilde{V}-b\cdot \mathrm{\nabla }\tilde{V}+r\tilde{V}=0.$

Next, we define the operator e which will represent the constraint in the Lagrangian. Since at ${\mathrm{\Gamma }}_c$ no boundary condition has to be given, we introduce $\tilde{\mathrm{\Omega }}=\mathrm{\Omega }\cap \mathrm{\Gamma }$ and the operator e is implicitly defined by (17) $\begin{array}{rl}⟨e(\tilde{V},\xi ),\psi ⟩& ={\int }_0^T{\int }_{\tilde{\mathrm{\Omega }}}[{\tilde{V}}_{\tau }-\mathrm{\nabla }\cdot A\mathrm{\nabla }\tilde{V}-b\cdot \mathrm{\nabla }\tilde{V}+r\tilde{V}]\phi \mathrm{d}z\mathrm{d}\tau \\ & +{\int }_0^T{\int }_{{\mathrm{\Gamma }}_a}[\tilde{V}-\exp (-\mathrm{r\tau })]{\phi }^a\mathrm{d}s\mathrm{d}\tau +{\int }_0^T{\int }_{{\mathrm{\Gamma }}_b}\tilde{V}{\phi }^b\mathrm{d}s\mathrm{d}\tau \\ & +{\int }_0^T{\int }_{{\mathrm{\Gamma }}_d}[\tilde{V}-\exp (-\mathrm{r\tau })]{\phi }^d\mathrm{d}s\mathrm{d}\tau \\ & =:T_1+T_2+T_3+T_4.\end{array}$(17) The Lagrangian for the constrained parameter calibration problem is then given by $L(\tilde{V},\xi ,\psi )=J(\tilde{V},\xi )-⟨e(\tilde{V},\xi ),\psi ⟩.$We formally compute the first-order optimality conditions by setting dL = 0. For details on the method we refer to [7,18]. Before computing the Gâteaux derivatives of L in arbitrary directions [7], we note that by Green's first identity, it holds (18) ${\int }_{\mathrm{\Omega }}(b\cdot \mathrm{\nabla }\tilde{V})\phi \mathrm{d}z={\int }_{\mathrm{\partial \Omega }}(b\cdot \overrightarrow{n})\tilde{V}\phi \mathrm{d}s-{\int }_{\mathrm{\Omega }}\tilde{V}\mathrm{\nabla }\cdot (\mathrm{b\phi })\mathrm{d}z.$(18) Therefore, we can rewrite $\begin{array}{rl}{T}_1& ={\int }_0^T{\int }_{\tilde{\mathrm{\Omega }}}\phi {\tilde{V}}_{\tau }+A\mathrm{\nabla }\tilde{V}\cdot \mathrm{\nabla \phi }-\frac{1}{2}b\cdot \mathrm{\nabla }\tilde{V}\phi +\frac{1}{2}\tilde{V}b\cdot \mathrm{\nabla \phi }+(r+\frac{1}{2}\mathrm{\nabla }\cdot b)\tilde{V}\phi \mathrm{d}z\\ & -{\int }_{\mathrm{\partial \Omega }}(A\mathrm{\nabla }\tilde{V})\cdot \overrightarrow{n}\phi \mathrm{d}s-\frac{1}{2}{\int }_{\mathrm{\partial \Omega }}(b\cdot \overrightarrow{n})\tilde{V}\phi \mathrm{d}s\mathrm{d}\tau \\ & ={\left[{\int }_{\mathrm{\Omega }}\phi \tilde{V}\mathrm{d}z\right]}_{\tau =0}^{\tau =T}+{\int }_0^T{\int }_{\tilde{\mathrm{\Omega }}}\tilde{V}[-{\phi }_{\tau }-\mathrm{\nabla }\cdot A^{\mathrm{\top }}\mathrm{\nabla \phi }+b\cdot \mathrm{\nabla \phi }+(r+\mathrm{\nabla }\cdot b)\phi ]\mathrm{d}z\\ & +{\int }_{\mathrm{\partial \Omega }}[(A^{\mathrm{\top }}\mathrm{\nabla \phi })\cdot \overrightarrow{n}-(b\cdot \overrightarrow{n})\phi ]\tilde{V}\mathrm{d}s-{\int }_{\mathrm{\partial \Omega }}(A\mathrm{\nabla }\tilde{V})\cdot \overrightarrow{n}\phi \mathrm{d}s\mathrm{d}\tau .\end{array}$As e is linear in V, the Gâteaux derivative in some arbitrary direction h reads $\begin{array}{rl}{\mathrm{d}}_V⟨e(V,\xi ),\psi ⟩[h]& ={\left[{\int }_{\mathrm{\Omega }}\mathrm{\phi h}\mathrm{d}z\right]}_{\tau =0}^{\tau =T}+{\int }_0^T{\int }_{\mathrm{\Omega }}h[-{\phi }_{\tau }-\mathrm{\nabla }\cdot A^{\mathrm{\top }}\mathrm{\nabla \phi }+b\cdot \mathrm{\nabla \phi }+(r+\mathrm{\nabla }\cdot b)\phi ]\mathrm{d}z\\ & +{\int }_{\mathrm{\partial \Omega }}[(A^{\mathrm{\top }}\mathrm{\nabla \phi })\cdot \overrightarrow{n}-(b\cdot \overrightarrow{n})\phi ]h\mathrm{d}s-{\int }_{\mathrm{\partial \Omega }}(A\mathrm{\nabla h})\cdot \overrightarrow{n}\phi \mathrm{d}s\mathrm{d}\tau \\ & +{\int }_0^T{\int }_{{\mathrm{\Gamma }}_a}h{\phi }^a\mathrm{d}s\mathrm{d}\tau +{\int }_0^T{\int }_{{\mathrm{\Gamma }}_b}h{\phi }^b\mathrm{d}s\mathrm{d}\tau +{\int }_0^T{\int }_{{\mathrm{\Gamma }}_d}h{\phi }^d\mathrm{d}s\mathrm{d}\tau .\end{array}$For the cost functional we have ${\mathrm{d}}_{\tilde{V}}J(\tilde{V},\xi )[h]={\int }_0^T{\int }_{\mathrm{\Omega }}\mathrm{hK}(K\tilde{V}-V_{\mathrm{data}})\mathrm{d}z\mathrm{d}\tau .$

To identify the adjoint equation, we consider $\begin{array}{rl}0& =d_{\tilde{V}}L(\tilde{V},u,\psi )[h]\\ & ={\int }_0^{T}h[{\int }_{\mathrm{\Omega }}K(K\tilde{V}-V_{\mathrm{data}})+{\phi }_{\tau }+\mathrm{\nabla }\cdot A^{\mathrm{\top }}\mathrm{\nabla \phi }-b\cdot \mathrm{\nabla \phi }-(r+\mathrm{\nabla }\cdot b)\phi ]\mathrm{d}z\mathrm{d}\tau .\end{array}$for arbitrary h. Note that we are not allowed to vary $\tilde{V}$ at $\tilde{V}(x,\nu ,0)$ as the initial condition is fixed. Therefore we have $h(0,\nu ,x)\equiv 0$.

For choosing $h\equiv 0$ on $\mathrm{\partial \Omega }$ and $h(T,\nu ,x)=0$, we find with the Variational Lemma ${\phi }_{\tau }+\mathrm{\nabla }\cdot A^{\mathrm{\top }}\mathrm{\nabla }\phi -b\cdot \mathrm{\nabla \phi }-(r+\mathrm{\nabla }\cdot b)\phi =-K(K\tilde{V}-V_{\mathrm{data}})\mathrm{on}\mathrm{\Omega }.$Now, choosing $h(T,\nu ,x)\ne 0$, we then obtain the terminal condition $\phi (T,\nu ,x)=0$.

We consider the four boundary conditions separately. At ${\mathrm{\Gamma }}_c$, also the parabolic adjoint PDE degenerates to a first-order hyperbolic PDE, and thus we have to consider the Fichera theory [2,10] for the variance again.

The Fichera condition w.r.t. the variance at $\nu =0$ of the adjoint is the same as before. Therefore no analytic boundary condition is supplied for this boundary, as we assume that the Feller condition holds. On ${\mathrm{\Gamma }}_a$ we have (19) $0={\int }_{{\mathrm{\Gamma }}_a}[(A^{\mathrm{\top }}\mathrm{\nabla }\phi )\cdot \overrightarrow{n}-(b\cdot \overrightarrow{n})\phi ]h-(A\mathrm{\nabla h})\cdot \overrightarrow{n}\phi +h{\phi }^a\mathrm{d}s.$(19) Choosing $h\equiv \mathrm{const}\ne 0$ yields $0={\int }_{{\mathrm{\Gamma }}_a}h[(A^{\mathrm{\top }}\mathrm{\nabla \phi })\cdot \overrightarrow{n}-(b\cdot \overrightarrow{n})\phi +{\phi }^a]\mathrm{d}s,$hence $(A^{\mathrm{\top }}\mathrm{\nabla }\phi )\cdot \overrightarrow{n}-(b\cdot \overrightarrow{n})\phi +{\phi }^a=0$. On the other hand, choosing $\mathrm{\nabla }h\ne 0$ (19(19) $0={\int }_{{\mathrm{\Gamma }}_a}[(A^{\mathrm{\top }}\mathrm{\nabla }\phi )\cdot \overrightarrow{n}-(b\cdot \overrightarrow{n})\phi ]h-(A\mathrm{\nabla h})\cdot \overrightarrow{n}\phi +h{\phi }^a\mathrm{d}s.$(19) ) must still hold. This yields $\phi =0$ on ${\mathrm{\Gamma }}_a$ and $\begin{array}{r}{\phi }^a=-(A^{\mathrm{\top }}\mathrm{\nabla }\phi )\cdot \overrightarrow{n}=-(A^{\mathrm{\top }}\mathrm{\nabla \phi })\cdot \left(\begin{array}{c}0\\ -1\end{array}\right)=\frac{1}{2}\nu {\sigma }_{\nu }\rho {\phi }_{\nu }+\frac{1}{2}\nu {\phi }_x=\frac{1}{2}\nu {\phi }_x\end{array}$As $\tilde{V}(x_{min},\nu ,0)=\exp (-\mathrm{r\tau })$ is given and independent of x, we obtain ${\phi }_x=0$ there and thus ${\phi }^a=\phi =0$ at this boundary. Similarly, we find on ${\mathrm{\Gamma }}_b$ that (20) $0={\int }_{{\mathrm{\Gamma }}_b}[(A^{\mathrm{\top }}\mathrm{\nabla }\phi )\cdot \overrightarrow{n}-(b\cdot \overrightarrow{n})\phi ]h-(A\mathrm{\nabla h})\cdot \overrightarrow{n}\phi +h{\phi }_b\mathrm{d}s.$(20) With the same arguments, we obtain $\phi =0$ and ${\phi }^b=0$ on ${\mathrm{\Gamma }}_b$ as well. Following the same arguments for ${\mathrm{\Gamma }}_d$ we obtain $(A^{\mathrm{\top }}\mathrm{\nabla }\phi )\cdot \overrightarrow{n}-(b\cdot \overrightarrow{n})\phi +{\phi }^d=0$ with $\phi =0$ and thus it reduces to $\begin{array}{r}{\phi }^d=-(A^{\mathrm{\top }}\mathrm{\nabla }\phi )\cdot \overrightarrow{n}=-(A^{\mathrm{\top }}\mathrm{\nabla \phi })\cdot \left(\begin{array}{c}1\\ 0\end{array}\right)=-\frac{1}{2}\nu {\sigma }_{\nu }^2{\phi }_{\nu }+\frac{1}{2}\nu {\sigma }_{\nu }\rho {\phi }_x=-\frac{1}{2}\nu {\sigma }_{\nu }^2{\phi }_{\nu }.\end{array}$As $\tilde{V}(x,\mathrm{\infty },\tau )=\exp (-\mathrm{r\tau })$ is given and independent of ν, we obtain ${\phi }_{\nu }=0$ there and thus ${\phi }^d=\phi =0$ at this boundary.

Altogether, the adjoint equation reads (21) $\frac{\mathrm{\partial }\phi }{\mathrm{\partial }\tau }+\mathrm{\nabla }\cdot A^{\mathrm{\top }}\mathrm{\nabla }\phi -b\cdot \mathrm{\nabla \phi }-(r+\mathrm{\nabla }\cdot b)\phi =-K(K\tilde{V}-V_{\mathrm{data}})\mathrm{on}\mathrm{\Omega },$(21)

which is equivalent to (22) $\begin{array}{rl}& {\phi }_{\tau }+\frac{1}{2}\nu {\sigma }_{\nu }^2{\phi }_{\mathrm{\nu \nu }}+\nu {\sigma }_{\nu }\rho {\phi }_{\mathrm{x\nu }}+\frac{1}{2}\nu {\phi }_{\mathrm{xx}}+({\sigma }_{\nu }^2-{\kappa }_{\nu }({\mu }_{\nu }-\nu )){\phi }_{\nu }\\ & +(q-r+\frac{\nu }{2}+{\sigma }_{\nu }\rho ){\phi }_x+({\kappa }_{\nu }-r)\phi =-K(K\tilde{V}-V_{\mathrm{data}})\mathrm{on}\mathrm{\Omega }\end{array}$(22) with terminal condition $\phi (T)=0$ and $\phi =0$ on the boundaries ${\mathrm{\Gamma }}_a$, ${\mathrm{\Gamma }}_b$ and ${\mathrm{\Gamma }}_d$ and the outflow boundary at $\nu =0$.

3.2. Derivation of the gradient

Let $\xi =({\sigma }_{\nu },\rho ,{\kappa }_{\nu },{\mu }_{\nu })$ be the parameters to be identified, as r and q are given by the data. We compute the optimality condition by setting $d_{\xi }L(V,\xi ,\psi )=0$. Since the boundaries ${\mathrm{\Gamma }}_a$, ${\mathrm{\Gamma }}_b$ and ${\mathrm{\Gamma }}_d$ are zero, we focus on $\tilde{\mathrm{\Omega }}$. In the following, we state the derivatives w.r.t. the different parameters separately. For ${\sigma }_{\nu }$ we get ${\mathrm{d}}_{{\sigma }_{\nu }}⟨e(\tilde{V},\xi ),\psi ⟩={\int }_0^T{\int }_{\tilde{\mathrm{\Omega }}}\tilde{V}[-{\sigma }_{\nu }\nu {\phi }_{\mathrm{\nu \nu }}-2{\sigma }_{\nu }{\phi }_{\nu }-\rho {\phi }_x-\mathrm{\rho \nu }{\phi }_{\mathrm{x\nu }}]\mathrm{d}z\mathrm{d}\tau .$Similarly, we obtain for the other derivatives $\begin{array}{rl}{\mathrm{d}}_{\rho }⟨e(\tilde{V},\xi ),\psi ⟩& ={\int }_0^T{\int }_{\tilde{\mathrm{\Omega }}}\tilde{V}[-{\sigma }_{\nu }{\phi }_x-{\sigma }_{\nu }\nu {\phi }_{\mathrm{x\nu }}]\mathrm{d}z\mathrm{d}\tau ,\\ {\mathrm{d}}_{{\kappa }_{\nu }}⟨e(\tilde{V},\xi ),\psi ⟩& ={\int }_0^T{\int }_{\tilde{\mathrm{\Omega }}}\tilde{V}[({\mu }_{\nu }-\nu ){\phi }_{\nu }-\phi ]\mathrm{d}z\mathrm{d}\tau ,\\ {\mathrm{d}}_{{\mu }_{\nu }}⟨e(\tilde{V},\xi ),\psi ⟩& ={\int }_0^T{\int }_{\tilde{\mathrm{\Omega }}}{\kappa }_{\nu }\tilde{V}{\phi }_{\nu }\mathrm{d}z\mathrm{d}\tau .\end{array}$Note that $d_{\xi }L(\tilde{V},\xi ,\psi )[h_{\xi }]=0$ needs to hold for arbitrary directions $h_{\xi }$. Therefore, we can read off the gradient from the above expressions.

We extend this gradient formulation for time-dependent parameter $\tilde{\xi }=({\tilde{\sigma }}_{\nu },\tilde{\rho },{\tilde{\kappa }}_{\nu },{\tilde{\mu }}_{\nu })$. The gradient is then time-dependent as well and given by $\begin{array}{rl}{\mathrm{d}}_{{\sigma }_{\nu }}⟨e(\tilde{V},\tilde{\xi }),\psi ⟩& ={\int }_{\tilde{\mathrm{\Omega }}}\tilde{V}[-{\tilde{\sigma }}_{\nu }\nu {\phi }_{\mathrm{\nu \nu }}-2{\tilde{\sigma }}_{\nu }{\phi }_{\nu }-\tilde{\rho }{\phi }_x-\tilde{\rho }\nu {\phi }_{\mathrm{x\nu }}]\mathrm{d}z\mathrm{d}\tau \\ {\mathrm{d}}_{\rho }⟨e(\tilde{V},\tilde{\xi }),\psi ⟩& ={\int }_{\tilde{\mathrm{\Omega }}}\tilde{V}[-{\tilde{\sigma }}_{\nu }{\phi }_x-{\tilde{\sigma }}_{\nu }\nu {\phi }_{\mathrm{x\nu }}]\mathrm{d}z\mathrm{d}\tau ,\\ {\mathrm{d}}_{{\kappa }_{\nu }}⟨e(\tilde{V},\tilde{\xi }),\psi ⟩& ={\int }_{\tilde{\mathrm{\Omega }}}\tilde{V}[({\tilde{\mu }}_{\nu }-\nu ){\phi }_{\nu }-\phi ]\mathrm{d}z\mathrm{d}\tau ,\\ {\mathrm{d}}_{{\mu }_{\nu }}⟨e(\tilde{V},\tilde{\xi }),\psi ⟩& ={\int }_{\tilde{\mathrm{\Omega }}}{\tilde{\kappa }}_{\nu }\tilde{V}{\phi }_{\nu }\mathrm{d}z.\end{array}$

3.3. Gradient descent algorithm for the parameter calibration

Solving the first-order optimality condition all at once is difficult due to the forward-backward structure. Therefore, we propose a gradient descent algorithm in the following.

For a given initial parameter set ${\xi }_0$, we can solve the state equation for the Heston model with constant control variable ξ (9(9) ${\tilde{V}}_{\tau }=\frac{\nu }{2}{\tilde{V}}_{\mathrm{xx}}+\frac{1}{2}{\sigma }_{\nu }^2\nu {\tilde{V}}_{\mathrm{\nu \nu }}+(r-q-\frac{\nu }{2}){\tilde{V}}_x+{\kappa }_{\nu }({\mu }_{\nu }-\nu ){\tilde{V}}_{\nu }+{\sigma }_{\nu }\mathrm{\nu \rho }{\tilde{V}}_{\mathrm{x\nu }}-r\tilde{V},$(9) ) or time-dependent parameter (15(15) ${\tilde{V}}_{\tau }=\frac{\nu }{2}{\tilde{V}}_{\mathrm{xx}}+\frac{1}{2}{\tilde{\sigma }}_{\nu }^2\nu {\tilde{V}}_{\mathrm{\nu \nu }}+(r-q-\frac{\nu }{2}){\tilde{V}}_x+{\tilde{\kappa }}_{\nu }({\tilde{\mu }}_{\nu }-\nu ){\tilde{V}}_{\nu }+{\tilde{\sigma }}_{\nu }\nu \tilde{\rho }{\tilde{V}}_{\mathrm{x\nu }}-r\tilde{V}.$(15) ). With the state solution at hand, we can compute the corresponding adjoint equation (21(21) $\frac{\mathrm{\partial }\phi }{\mathrm{\partial }\tau }+\mathrm{\nabla }\cdot A^{\mathrm{\top }}\mathrm{\nabla }\phi -b\cdot \mathrm{\nabla \phi }-(r+\mathrm{\nabla }\cdot b)\phi =-K(K\tilde{V}-V_{\mathrm{data}})\mathrm{on}\mathrm{\Omega },$(21) ) or (22(22) $\begin{array}{rl}& {\phi }_{\tau }+\frac{1}{2}\nu {\sigma }_{\nu }^2{\phi }_{\mathrm{\nu \nu }}+\nu {\sigma }_{\nu }\rho {\phi }_{\mathrm{x\nu }}+\frac{1}{2}\nu {\phi }_{\mathrm{xx}}+({\sigma }_{\nu }^2-{\kappa }_{\nu }({\mu }_{\nu }-\nu )){\phi }_{\nu }\\ & +(q-r+\frac{\nu }{2}+{\sigma }_{\nu }\rho ){\phi }_x+({\kappa }_{\nu }-r)\phi =-K(K\tilde{V}-V_{\mathrm{data}})\mathrm{on}\mathrm{\Omega }\end{array}$(22) ) backwards in time. Then we have all the information available to compute the gradient and update the parameter set using a gradient step. The procedure is sketched in algorithm 1.

Since the parameter domain for ${\kappa }_{\nu }$, ${\mu }_{\nu }$, ${\sigma }_{\nu }$ and ρ is restricted, as well as the constraint that the Feller condition has to be fulfilled, we use the projected Armijo rule [18]. In the projected Armijo rule, we choose the maximum ${\sigma }_k\in \{1,1/2,1/4,\dots \}$ for which $f(\mathcal{P}({\xi }_k-{\sigma }_k\mathrm{\nabla f}({\xi }_k)))-f({\xi }_k)\le -\frac{\gamma }{{\sigma }_k}\Vert \mathcal{P}({\xi }_k-{\sigma }_k\mathrm{\nabla f}({\xi }_k))-{\xi }_k{\Vert }_2^2.$Here $\gamma \in (0,1)$ is a numerical constant, which is problem-dependent and typically chosen as $\gamma ={10}^{-4}$. We will use this value for the numerical results later on.

4. Discretization

In this section, we introduce a closure boundary condition at $\nu =0$ for the Heston and its adjoint and perform a domain truncation to discretize the problem. Following the Fichera theory and assuming that the Feller condition holds, no analytical boundary condition needs to be imposed at $\nu =0$, neither for the Heston model nor for its adjoint, since the PDEs have a pure outflow boundary at this point. Nevertheless, we need a closure condition for this boundary for the discretization. We suggest to follow Heston's approach, as discussed in Refs [3,10] and use the reduced hyperbolic formulation of the Heston PDE and its adjoint. At ${\mathrm{\Gamma }}_c$, we obtain (23) ${\tilde{V}}_{\tau }=r{\tilde{V}}_x+{\kappa }_{\nu }{\mu }_{\nu }{\tilde{V}}_{\nu }-r\tilde{V}$(23) for the log-transformed normalized PDE and similar for the adjoint (24) ${\phi }_{\tau }+({\sigma }_{\nu }^2-{\kappa }_{\nu }{\mu }_{\nu }){\phi }_{\nu }+(q-r+{\sigma }_{\nu }\rho ){\phi }_x+({\kappa }_{\nu }-r)\phi =-K(K\tilde{V}-V_{\mathrm{data}}).$(24) Now, we perform a domain truncation to obtain a rectangular grid, instead of a semi-unbounded domain and introduce the grid points. We consider uniform meshes in each direction and obtain for the spatial directions $x_i=x_{min}+i{\mathrm{\Delta }}_x$ for $i=0,\dots ,N_x$ with ${\mathrm{\Delta }}_x=(x_{max}-x_{min})/N_x$ and ${\nu }_j=j{\mathrm{\Delta }}_{\nu }$ for $j=0,\dots ,N_{\nu }$ with ${\mathrm{\Delta }}_{\nu }={\nu }_{max}/N_{\nu }$, as well as ${\tau }_k=k{\mathrm{\Delta }}_{\tau }$ for $k=0,\dots ,N_{\tau }$ with ${\mathrm{\Delta }}_{\tau }=T/N_{\tau }$ for the temporal direction.

Since this is a proof-of-concept, simple and well-known spatial and temporal discretization methods are used to illustrate our approach. For the time discretization, we use the well-known alternating direction implicit (ADI) method, in more detail the Hundsdorfer-Verwer scheme [9]. This scheme is of second order for any choice of θ, where θ is a measure of classification similar to the θ-method and is able to handle mixed derivative terms. To present the ADI method, we use a general second-order PDE formulation (25) $u_{\tau }+a_{11}{u}_{\mathrm{\nu \nu }}+2a_{12}{u}_{\mathrm{x\nu }}+a_{22}{u}_{\mathrm{xx}}+b_1{u}_{\nu }+b_2{u}_x+\mathrm{cu}=0.$(25) Since the log-transformed normalized Heston PDE and its adjoint are second-order PDEs and we further use the same discretization.

In the first step, we split the operator of the PDE into three operators (26) $\mathcal{F}(\tau )={\mathcal{F}}_0(\tau )+{\mathcal{F}}_1(\tau )+{\mathcal{F}}_2(\tau ),$(26) with $\begin{array}{rl}{\mathcal{F}}_0(\tau )& =2a_{12}{D}_{\mathrm{x\nu }}\\ {\mathcal{F}}_1(\tau )& =b_2{D}_x+a_{22}{D}_{\mathrm{xx}}-\frac{1}{2}\mathrm{cI}\\ {\mathcal{F}}_2(\tau )& =b_1{D}_{\nu }+a_{11}{D}_{\mathrm{\nu \nu }}-\frac{1}{2}\mathrm{cI},\end{array}$where $D_x$ describes the discretization matrix of the first derivative w.r.t. x and accordingly $D_{\mathrm{xx}}$ of the second derivative w.r.t. x, $D_{\nu }$ and $D_{\mathrm{\nu \nu }}$ of the first and second derivative w.r.t. ν, I denotes the identity matrix. The discretization matrices are derived using central finite differences. Let $u_k^{i,j}\approx u(x_i,{\nu }_j,{\tau }_k)$ and simplify $u_k\approx u(x,\nu ,{\tau }_k)$. In each time step, we have to solve the following system of equations (27) $\{\begin{array}{l}{Y}_0=u_k+\mathrm{\Delta }\tau \mathcal{F}({\tau }_k)u_k,\\ Y_1=Y_0+\theta \mathrm{\Delta \tau }({\mathcal{F}}_1({\tau }_{k+1})Y_1-{\mathcal{F}}_1({\tau }_k)u_k),\\ Y_2=Y_x+\theta \mathrm{\Delta \tau }({\mathcal{F}}_2({\tau }_{k+1},Y_2)-{\mathcal{F}}_2({\tau }_k)u_k),\\ {\tilde{Y}}_0=Y_0+\frac{1}{2}\mathrm{\Delta \tau }(\mathcal{F}({\tau }_{k+1})Y_2-\mathcal{F}({\tau }_k)u_k),\\ {\tilde{Y}}_1={\tilde{Y}}_0+\theta \mathrm{\Delta \tau }({\mathcal{F}}_1({\tau }_{k+1}){\tilde{Y}}_1-{\mathcal{F}}_1({\tau }_k)u_k),\\ {\tilde{Y}}_2={\tilde{Y}}_x+\theta \mathrm{\Delta \tau }({\mathcal{F}}_2({\tau }_{k+1}){\tilde{Y}}_2-{\mathcal{F}}_2({\tau }_k)u_k),\\ u_{k+1}={\tilde{Y}}_2.\end{array}$(27) We choose $\theta =3/4$ and improve the implementation as we used the approach in [16] by using a matrix based instead of a vector-based implementation of $u_k$.

At this point, we have to discuss the boundary conditions for the Heston model again. We start with the x dimension. We set $\tilde{V}(x_{min},\nu ,\tau )=\exp (-\mathrm{r\tau })$ and similar $\tilde{V}(x_{max},\nu ,\tau )=0$, as we suggest a sufficiently small $x_{min}$ and a sufficiently large $x_{max}$. For the variance boundaries, we follow the approach of Kùtik and Mikula [10]. Due to the Fichera theory, at $\nu =0$ we gain an outflow boundary and no information can enter the domain from the region $\nu <0$. Since the same holds for the adjoint, we propose the same approach there. In these, we extend the numerical domain for this case and introduce ghost cells $p_{i,0}=(x_i,{\nu }_0)$, $i=1,\dots ,N$ where ${\nu }_0=-{\mathrm{\Delta }}_y$ and determine the ghost cell values $u_{p_{i,0}}^n$, $i=1,\dots ,N$ by zero order extrapolation. We obtain a constant function (28) $u_{p_{i,0}}^n=u_{p_{i,1}}^n,i=1,\dots ,N.$(28) Now we consider the boundary, where $\nu \to \mathrm{\infty }$ and the truncation shrinks the domain from $(0,\mathrm{\infty })\to (0,{\nu }_{max})$. Dirichlet boundary condition from Heston would cause an unnatural jump in the solution. Therefore different strategies have been developed to overcome this issue, e.g. [10]. Again, we follow Kùtik and Mikula and impose artificial homogeneous Neumann boundary conditions. The choice is motivated by the variance independence of the original boundary condition from Heston, for sufficiently large ${\nu }_{max}=\mathcal{O}(1)$. Again, we use zero-order extrapolation to implement this condition. We add ghost cells $p_{i,N_{\nu }+1}$, $i=1,\dots ,N$ where ${\nu }_{N_{\nu }+1}={\nu }_{max}+{\mathrm{\Delta }}_y$ and determine their values $u_{p_{i,N_{\nu }+1}}^n$, $i=1,\dots ,N$ by a constant function (29) $u_{p_{i,N_{\nu }+1}}^n=u_{p_{i,N_{\nu }}}^n,i=1,\dots ,N.$(29) For a discussion about different numerical boundary conditions for the variance of the Heston model, we refer to [3]. The integrals appearing in the gradient are computed by the trapezoidal rule.

5. Numerical results

Following [3], we use (30) $K=1.0,r=0.1,{\xi }_{\mathrm{ref}}=(5.0,0.07,0.5,-0.5)$(30) for generating an artificial $V_{\mathrm{data}}$ for each time step ${\tau }_k$. For the discretization, we use $N_x=79,N_{\nu }=39,N_{\tau }=59.$ As bounds for the projected Armijo rule for $\xi =({\kappa }_{\nu },{\mu }_{\nu },{\sigma }_{\nu },\rho )$, we set $0<{\kappa }_{\nu }<8,0<{\mu }_{\nu }<1,0<{\sigma }_{\nu }<1,-1<\rho <1.$ Note that the projected Armijo-rule ensures that the Feller condition holds within each optimization step. Hence we are in the case of an outflow boundary. We set the maximal iteration value for the calibration to 20. For the initial guesses ${\xi }_{\mathrm{init}}=({\kappa }_{\nu }^{\mathrm{init}},{\mu }_{\nu }^{\mathrm{init}},{\sigma }_{\nu }^{\mathrm{init}},{\rho }^{\mathrm{init}})$ we used generated random numbers within a maximal percentage difference from ${\xi }_{\mathrm{ref}}$. We use four different percentages 10, 25, 50 and 75 and for each we generate five sets. The initial parameters as well as the calibrated parameters in the constant case are given in Table 1 for ${\kappa }_{\nu }$, Table 2 for ${\mu }_{\nu }$, Table 3 for ${\sigma }_{\nu }$ and Table 4 for ρ.

Table 1. Five initial values sets for ${\kappa }_{\nu }^{\mathrm{init}}$ for the different percentages and the corresponding calibrated values for C0.

${\kappa }_{\nu }^{\mathrm{init}}$	T1	T2	T3	T4	T5	${\kappa }_{\nu }^{\mathrm{cal}}$	T1	T2	T3	T4	T5
10	4.60	5.15	4.50	4.70	4.50	10	4.60	5.15	4.51	4.70	4.50
25	4.75	4.60	5.25	5.90	4.45	25	4.75	4.60	5.24	5.88	4.45
50	4.10	3.25	3.55	2.90	7.00	50	4.11	3.26	3.55	2.94	6.95
75	5.70	7.30	7.00	7.50	6.60	75	5.69	7.29	6.85	7.40	6.60

Table 2. Five initial values sets for ${\mu }_{\nu }^{\mathrm{init}}$ for the different percentages and the corresponding calibrated values for C0.

${\mu }_{\nu }^{\mathrm{init}}$	T1	T2	T3	T4	T5	${\mu }_{\nu }^{\mathrm{cal}}$	T1	T2	T3	T4	T5
10	0.0679	0.0637	0.0735	0.0714	0.0714	10	0.0598	0.0727	0.0575	0.0632	0.0577
25	0.0721	0.0777	0.0616	0.0630	0.0868	25	0.0655	0.0610	0.0720	0.0836	0.0615
50	0.1001	0.0406	0.0840	0.0448	0.0434	50	0.0476	0.0148	0.0484	0.0136	0.0882
75	0.0742	0.0532	0.0448	0.0399	0.0490	75	0.0847	0.1069	0.1005	0.1064	0.0928

Table 3. Five initial values sets for ${\sigma }_{\nu }^{\mathrm{init}}$ for the different percentages and the corresponding calibrated values for C0.

${\sigma }_{\nu }^{\mathrm{init}}$	T1	T2	T3	T4	T5	${\sigma }_{\nu }^{\mathrm{cal}}$	T1	T2	T3	T4	T5
10	0.465	0.545	0.510	0.515	0.480	10	0.451	0.559	0.484	0.502	0.455
25	0.575	0.570	0.505	0.415	0.480	25	0.565	0.546	0.523	0.455	0.429
50	0.580	0.360	0.630	0.345	0.655	50	0.501	0.292	0.581	0.261	0.713
75	0.705	0.470	0.235	0.350	0.565	75	0.718	0.562	0.390	0.489	0.630

Table 4. Five initial values sets for ${\rho }^{\mathrm{init}}$ for the different percentages and the corresponding calibrated values for C0.

${\rho }^{\mathrm{init}}$	T1	T2	T3	T4	T5	${\rho }^{\mathrm{cal}}$	T1	T2	T3	T4	T5
10	−0.550	−0.485	−0.520	−0.520	−0.495	10	−0.483	−0.547	−0.404	−0.460	−0.384
25	−0.410	−0.490	−0.605	−0.585	−0.380	25	−0.370	−0.390	−0.690	−0.829	−0.159
50	−0.555	−0.465	−0.350	−0.685	−0.725	50	−0.258	−0.242	−0.207	−0.441	−0.847
75	−0.695	−0.205	−0.775	−0.655	−0.150	75	−0.743	−0.972	−0.990	−0.990	−0.437

We observe that the calibration leaves ${\kappa }_{\nu }$ nearly untouched, whereas the other parameter values are significantly changed. This could be reasoned by the structure of the drift term ${\kappa }_{\nu }({\mu }_{\nu }-\nu )$, since ${\kappa }_{\nu }$ and ${\mu }_{\nu }$ are multiplied. As a optimization measure, we compute the relative reduction of the cost functional using (31) $r({\xi }_{\mathrm{init}})=100\cdot (1-\frac{J(\tilde{V},{\xi }_{\mathrm{opt}})}{J(\tilde{V},{\xi }_{\mathrm{init}})}).$(31) Table 5 shows the relative reduction of the cost functional for the different test cases for ${\xi }_{\mathrm{init}}$.

Table 5. Relative reduction of the cost functional computed with 31(31) $r({\xi }_{\mathrm{init}})=100\cdot (1-\frac{J(\tilde{V},{\xi }_{\mathrm{opt}})}{J(\tilde{V},{\xi }_{\mathrm{init}})}).$(31) using the different test cases for ${\xi }_{\mathrm{init}}$ and ${\xi }_{\mathrm{cal}}$ in the constant calibration setting.

C0	T1	T2	T3	T4	T5
10	43.12	77.38	65.22	63.22	52.75
25	58.46	78.38	62.14	66.06	64.53
50	83.09	20.84	52.76	30.98	72.62
75	16.23	72.89	82.89	82.09	80.43

As the calibrated values differ over the test cases, it is reasonable to assume, that we only find local minima. However the results are remarkable.

Since in the real market the parameter are not considered as constant, we improve the approach by considering different parameters, and some parameter sets as time-dependent. From the relative change in the constant calibration setting, we choose the following (additional) test cases listed in Table 6. Further the table includes the links to the cost functional reduction tables and figures of the cost function evolution with respect to the different cases.

Table 6. Different cases for calibration setting.

	${\kappa }_{\nu }$	${\mu }_{\nu }$	${\sigma }_{\nu }$	ρ
C0	constant	constant	constant	constant
C1	time-dependent	constant	constant	constant
C2	constant	time-dependent	constant	constant
C3	constant	constant	time-dependent	constant
C4	constant	constant	constant	time-dependent
C5	time-dependent	time-dependent	time-dependent	time-dependent
C6	constant	time-dependent	time-dependent	time-dependent
C7	constant	time-dependent	constant	time-dependent

For each case of the time-dependent test cases $V_{\mathrm{data}}$ is generated as before and ${\xi }^{\mathrm{init}}$ is assumed to be constant and thus also the same as before.

In general, the Figures 1–8 show the same cost functional evolution. The first test cases with the initial guess closest to the reference set have the smallest initial cost functional value. As the distance to the reference set increases, so does the initial cost function value. Accordingly, the final cost function value is higher for the sets with a more distant initial guess. But the $75\mathrm{\%}$ deviation within the initial guess shows the largest cost function reduction within the first steps. This is due to the existence of different minima, as these sets find the closest minima instead of the optimal minima. Overall, the reduction of the cost function per iteration decreases as the number of iterations increases. The number of iterations within the different test cases with a 10 % deviation is the same for all scenarios, except for T2 where it varies between 7 and 8. The scenarios with 75 % deviation give the same number of iterations in cases T1, T2 and T5, for cases T3 and T4 we observe iteration numbers with a total difference of one from each other. For the deviation of 25 %, we observe the same number of iterations for T1, a difference of one iteration for T2, T3 and T4 across all scenarios, and a difference of two for the last case T5. The largest difference within the number of iterations is found in the 50 % deviation. While T5 has the same number of iterations and T2 has a total difference of 2 for all cases. For the other test cases, Table 14 shows the number of iterations. These are the only cases where we see an increase in the number of iterations due to considering time-dependent parameters since the number of iterations increases when ${\mu }_{\nu }$ is considered time-dependent. These are highlighted in bold in the table. If we look at the corresponding Figures 3, 6, 7 as well as 8, we observe that the reduction of the cost function in the end is very small, s.t. the gain from iterations 10 forward is almost negligible. Thus, 10 iterations are sufficient for a sufficient cost function reduction in all cases.

Figure 1. Cost functional evolution per iteration for the test cases within the constant parameter calibration.

We have highlighted the largest relative cost functional reduction with bold numbers and the smallest relative cost functional value with italics over the Tables 5–13. Comparing the relative cost functional values for each case gives an absolute difference between the reductions, see Table 15. We observe that in 9 cases the difference is less than three, while in eight cases it is between three and five, and only in three cases it is greater than five. These are T1 with a deviation of 10 % with a difference of 5.91, a difference of 8.37 for T5 with a deviation of 25 %, and T4 with a deviation of 50 % with the highest difference of 12.93. Overall, we observe that the reductions between C0-C4 as well as C5-C7 are comparable, while a larger difference can be seen between these two case sets.

Table 7. Relative cost function reduction for the C1 calibration.

C1	T1	T2	T3	T4	T5
10	44.48	78.75	64.91	64.12	54.03
25	58.97	78.98	62.96	66.77	63.96
50	83.13	19.59	45.71	33.45	72.84
75	16.59	72.80	82.14	81.13	80.94

Table 8. Relative cost function reduction for the C2 calibration.

C2	T1	T2	T3	T4	T5
10	42.96	77.85	63.17	62.09	53.05
25	58.03	77.05	63.09	68.46	70.00
50	82.90	20.79	51.89	29.79	73.66
75	16.05	74.64	82.95	82.43	78.89

Table 9. Relative cost function reduction for the C3 calibration.

C3	T1	T2	T3	T4	T5
10	44.22	77.89	64.50	64.04	53.52
25	58.49	78.76	61.62	65.92	64.50
50	83.05	20.97	52.49	30.55	72.23
75	15.96	72.91	82.71	82.10	80.44

Table 10. Relative cost function reduction for the C4 calibration.

C4	T1	T2	T3	T4	T5
10	49.03	75.90	66.62	66.49	56.81
25	60.42	77.71	62.64	68.50	62.96
50	78.54	18.24	52.63	30.29	74.03
75	20.89	73.59	81.75	80.49	79.74

Table 11. Relative cost function reduction for the C5 calibration.

C5	T1	T2	T3	T4	T5
10	47.63	75.50	65.79	64.45	57.08
25	60.70	76.67	62.36	68.51	71.33
50	82.32	20.89	58.63	31.72	72.96
75	20.79	73.10	78.14	77.71	78.76

Table 12. Relative cost function reduction for the C6 calibration.

C6	T1	T2	T3	T4	T5
10	47.91	76.08	65.50	64.77	56.43
25	60.27	76.39	63.37	70.16	68.91
50	80.94	19.19	57.15	29.38	73.71
75	20.84	74.74	81.70	80.74	79.10

Table 13. Relative cost function reduction for the C7 calibration.

C7	T1	T2	T3	T4	T5
10	47.79	75.66	65.50	64.66	56.48
25	60.21	76.37	62.80	69.43	69.29
50	81.08	19.35	57.16	30.05	73.02
75	20.62	74.55	80.59	80.09	78.89

Table 14. Number of iterations for the different case scenarios for Tests 1, 3, and 5 within the $50\mathrm{\%}$ percentage difference in the initial parameter set.

Test cases	C0	C1	C2	C3	C4	C5	C6	C7
T1	12	12	13	12	12	17	16	16
T3	6	6	5	6	6	17	15	15
T4	9	9	19	9	9	19	19	19

Note: The bold cases illustrate the scenarios where ${\mu }_{\nu }$ is considered time dependent.

Table 15. Absolute difference between the largest and smallest relative cost functional reduction.

Difference	T1	T2	T3	T4	T5
10	5.91	2.39	2.12	3.27	4.33
25	2.67	2.61	1.75	4.24	8.37
50	4.59	2.73	12.92	4.07	1.41
75	4.93	1.94	4.81	4.29	2.85

To quantify the general behaviour for the different calibrations, we compute the average relative cost functional reduction by (32) $r_a({\xi }_{\mathrm{init}})=\frac{r({\xi }_{\mathrm{init}})}{20}$(32)

and summarize the results in Table 16. Note, that all cost function reduction averages are huge. We observe that a time-dependent calibration for one parameter alone (C1-C4) doesn't improve the cost functional reduction significantly. The first slightly improvement can be found by using at least two time-dependent parameters (C5-C7). Surprisingly C7 gives the best calibration results, where ${\kappa }_{\nu }$ is the only variable which is calibrated as a constant, even though in total it has only two times the smallest cost functional reduction and never the largest, as well as C5, where all parameters are calibrated as time-dependent gives the least improvement, when considering combinations of time-dependent parameter calibration. The fact that $V_{\mathrm{data}}$ is generated with constant parameters and the best case considers time-dependent parameters indicates that the time-dependency is a good way to overcome the local minima. Those results are inline with literature [15,17].

Table 16. Average cost functional reduction in percentage for the different cases over the 20 test cases, as well as a list for the corresponding relative cost functional reduction as well as the figures of the cost functional evolution.

Case Setup	C0	C1	C2	C3	C4	C5	C6	C7
$r_a({\xi }_{\mathrm{init}})$	61.30	61.31	61.49	61.34	61.86	62.25	62.36	66.12
Table	5	7	8	9	10	11	12	13
Figure	1	2	3	4	5	6	7	8

6. Conclusion and outlook

Our paper began with an introduction to the Heston model. The introduction was followed by the derivation of a gradient-based optimization, including the derivation of a gradient descent algorithm. In the next section, the discretization of the schemes and algorithms is presented and numerical results are illustrated. These show that the gradient descent algorithm is a feasible calibration technique. The relative cost function reduction is remarkable, even though we can expect to find only local minima. Thus, our approach follows recent research. Furthermore, the assumption of at least time-dependent parameters is better suited to the real market situation, since in the real market almost no parameter is constant. In future work, we plan to incorporate the gradient descent algorithm into the space mapping approach and to test the space mapping approach on real market data.

Figure 2. Cost functional evolution per iteration for the different test set in the case scenario C1.

Figure 3. Cost functional evolution per iteration for the different test set in the case scenario C2.

Figure 4. Cost functional evolution per iteration for the different test set in the case scenario C3.

Figure 5. Cost functional evolution per iteration for the different test set in the case scenario C4.

Figure 6. Cost functional evolution per iteration for the different test set in the case scenario C5.

Figure 7. Cost functional evolution per iteration for the different test set in the case scenario C6.

Figure 8. Cost functional evolution per iteration for the different test set in the case scenario C7.

Disclosure statement

No potential conflict of interest was reported by the author(s).