DIFFERENTIAL DEEP LEARNING AND
PHYSICS-INFORMED NEURAL NETWORKS FOR
FAST DERIVATIVE PRICING

Moez Mrad, Nazim Iassamen

► To cite this version:

Moez Mrad, Nazim Iassamen. DIFFERENTIAL DEEP LEARNING AND PHYSICS-INFORMED NEURAL NETWORKS FOR FAST DERIVATIVE PRICING. 2025. hal-05160171

HAL Id: hal-05160171

https://hal.science/hal-05160171v1

Preprint submitted on 13 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# DIFFERENTIAL DEEP LEARNING AND PHYSICS-INFORMED NEURAL NETWORKS FOR FAST DERIVATIVE PRICING

Moez Mrad*

Nazim Iassamen†

This version: 10 Jul 2025

## ABSTRACT

This article examines the Heston model with time-dependent parameters, known to be flexible and complex. Valuing derivative products within this model typically involves numerical methods, such as Monte Carlo (MC) simulations and partial differential equations (PDEs). However, these methods are often computationally intensive. Our approach leverages two advanced Deep Learning (DL) techniques in training an Artificial Neural Network (ANN) that approximates option prices in a Heston model with time-dependent parameters while reducing the number of labels needed for training. This methodology can also be extended to other dynamics, such as Bates or rough volatility models.

We use payoff realizations as labels for the accelerated offline generation of datasets and employ a combined Differential Deep Learning (DDL) and Physics-Informed Neural Networks (PINNs) regularization strategy to train the ANN. Our approach captures the dynamics of the underlying risk factors and handles a variety of payoffs using fewer and rapidly computed labels. It allows fast and accurate pricing of derivative products and precise sensitivities computation, which improves risk management.<fnref n="3" />

## 1 Introduction

Pricing models such as local volatility models (Cox and Ross [1976], Brigo and Mercurio [2001, 2002]) are used to value exotic options. However, they often oversimplify the forward smile. In contrast, stochastic volatility models (Heston [1993], Tompkins [2001]) offer a richer representation of smiles. Despite this advantage, valuing even vanilla European options under these models can be computationally demanding due to the lack of a straightforward closed-form solution. Instead, numerical methods such as MC simulations, PDEs, Fourier methods, or asymptotic approaches are necessary, resulting in lengthy calibration times to match implied volatility market data.

This paper examines the Heston model, which describes instantaneous variance ($\nu_t$) as a mean-reverting square-root process correlated with the stock price process ($S_t$). Its dynamics are given by the following stochastic differential equations (SDEs):

$$dS_t = r_t S_t dt + \sqrt{\nu_t} S_t dW_t^s, \quad (1)$$

$$d\nu_t = \kappa(\bar{\nu} - \nu_t)dt + \gamma\sqrt{\nu_t}dW_t^\nu \quad (2)$$

$(W_t^s)$ and $(W_t^\nu)$ are two standard Brownian motions with constant correlation $\rho$. The mean reversion rate is represented by $\kappa$, while $\gamma$ denotes the volatility of the variance $\nu_t$. The term $\bar{\nu}$ signifies the long-term mean

<footnote>*Crédit Agricole Corporate and Investment Bank - moez.mrad@ca-cib.com</footnote>

<footnote>†Crédit Agricole Corporate and Investment Bank - nazim.iassamen@ca-cib.com</footnote>

<footnote n="3">³The views expressed in this paper are the authors' own and do not reflect those of their employer.</footnote>
reversion level and ($r_t$) is the interest rate process. In the classic Heston model, all these parameters are assumed to be time-independent. When the Feller condition $2\kappa\bar{\nu} > \gamma^2$ is satisfied, the process $\nu_t$ remains strictly positive, thereby preventing it from becoming negative, which implies $P(\forall t > 0 : \nu_t > 0) = 1$. If this condition is not met, the behavior at the origin needs to be further examined (see Feller [1954] and Carr and Linetsky [2006] for more details).

In the Heston model, the bidimensional process that encompasses $\ln(S_t)$ and $\nu_t$ exhibits an affine structure, which means that its moment-generating function (MGF) is exponential-affine. This feature allows the use of semi-analytical formulas for pricing call and put options through the Fourier inversion of the characteristic function provided by Heston [1993], Lewis [2000], Lipton [2002]. Despite the increased complexity when introducing piecewise-constant parameters, these formulas are still tractable using PDE approaches (Mikhailov and Nogel [2003]) or applying Markov logic in tandem with affine models (Mrad and Triki [2014]).

When the Heston model parameters are time-dependent, an analytical solution is generally unavailable, making MC simulations indispensable. This conclusion is consistent with the findings of Benhamou et al. [2010]. Allowing for time-varying parameters is crucial for calibrating the model to capture the term structure observed in market data accurately. This requirement has driven advances in the design of MC simulations to enhance their efficiency, as reported in Andersen [2006] and Mrad and Triki [2014]. For an in-depth understanding of numerical PDEs, the readers are referred to Kluge [2002] and Itkin and Carr [2011]. A comprehensive review of the various strategies used to price options in a time-dependent Heston model can be found in Rouah [2015].

Analytical solutions for pricing exotic options, such as barrier or multidimensional options (e.g., basket options), are limited, even with constant coefficients. Traditional numerical methods, such as PDEs and MC simulations, for pricing derivatives without closed-form solutions can be computationally demanding, particularly in high dimensions. The situation becomes even more critical in high-frequency trading and risk management-type applications, where real-time decisions must be made. This work addresses the need for efficient solutions to these complex pricing problems.

Recent research has explored the application of machine learning (ML) techniques to accelerate option pricing functions. This acceleration improves the efficiency of model calibration and allows real-time valuation of option portfolios. Building upon the Universal Approximation Theorem, neural networks have emerged as promising candidates for approximating these pricing functions, as illustrated in Bayer and Stemper [2018] and Horvath et al. [2021a]. The main limitation is the extremely long time required to generate the training data, since prices (labels) are produced using resource-intensive numerical methods. This requirement hampers the pace of research and development. To mitigate this challenge, we use payoff realizations of options as labels. They offer an efficient solution when combined with their pathwise sensitivities, especially in time-varying stochastic volatility models or high-dimensional asset prices. This accelerates data generation and balances theoretical rigor and practical application.

In this paper, we use the Sobolev Training technique Czarnecki et al. [2017] also known as Deep Differential Learning (DDL) or Differential Machine Learning (DML) in finance, as in Huge and Savine [2020] in conjunction with Physics-Informed Learning (PIL), introduced by Raissi et al. [2019]. PIL uses Physics-Informed Neural Networks (PINNs), an advanced ML tool suitable for handling PDEs. We aim to use PINNs to solve the Backward Kolmogorov PDE, that plays a central role in pricing options in the Heston model framework.

Our approach combines two advanced DL techniques to minimize the number of required labels. We utilize payoff realizations as labels for accelerated offline dataset generation and employ DDL and PINNs to train an ANN that approximates option prices in a Heston model with time-dependent parameters. This combination captures the dynamics of the underlying risk factors and enables efficient handling of exotic payoffs using fewer and rapidly computed labels, thus allowing for fast and accurate derivative pricing. Furthermore, our method provides accurate sensitivity calculations for effective risk management.
The cornerstone of this work is to leverage advanced DL methods to address the pricing challenges posed by
the time-dependent Heston model. Its main contributions are summarized as follows:

* Introducing an innovative combination of DDL and PINNs regularization strategies for training ANNs to approximate option prices in the time-dependent Heston model. Empirical evaluations show significant improvements in computational efficiency and pricing accuracy.

* Using MC techniques to generate payoff realizations as labels for the accelerated offline production of datasets. This approach simplifies the training process and ensures robust approximations even with limited data.

* Developing methods for calculating the sensitivity of labels to inputs, particularly relevant when using option prices as labels in the Heston model with piecewise parameters. We also employ complex-step differentiation to determine the sensitivity of payoff realizations (used as labels) to inputs, where these realizations are simulated using the Quadratic-Exponential (QE) scheme.

* Evaluating DDL's effectiveness using prices and payoff realizations across various Heston model parameter time-dependencies. We demonstrate that DDL is robust to noisy labels, particularly with limited datasets and high-dimensional feature spaces. Similarly to data augmentation in computer vision, label differentials in DDL amplify the training data and guide the ANN in capturing the option value contour.

* Enabling the computation of the Hessian of a Feed Forward Neural Network (FFNN) w.r.t. its input, facilitating the application of PINNs, and reducing the training phase duration.

* Highlighting the capabilities of the proposed methodology, which extends beyond accurate pricing to provide comprehensive sensitivity analyses critical for effective risk management.

The paper is divided into several sections. Section 2 provides the foundation by discussing the computation of prices and the simulation of payoff realizations for a time-dependent Heston model. It also introduces the fundamental concepts of DDL and PINNs, along with the overall training loss specification that merges both techniques. Section 3 specifically delves into the computation of price differentials, particularly in the context of piecewise parameters. Section 4 explores the calculation of sensitivities for payoff realizations with respect to various parameters of the time-dependent Heston model. Section 5 analyzes the regularization capabilities of DDL and PINNs. It examines the performance with respect to different label types and feature space dimensionality. Finally, Section 6 summarizes the key findings and contributions of this work.

## 2 Background

We aim to predict the value of a derivative product based on its features (such as strike price, maturity, etc.), market conditions, and the time-dependent parameters of the Heston model. To achieve this, we will use both the payoff's realizations and prices as training labels. We plan to train ANNs using DDL and PINNs.

Using the option's price as a label is feasible for Vanilla European options under a Heston model with piece-
wise time-dependent parameters, as we will demonstrate in the upcoming section. When no semi-closed
formula is available for the option price, we can use payoff realizations as labels instead. We will explain the
discretization of the Heston model with time-varying parameters, enabling the generation of option payoff
realization labels through MC simulations based on a set of attributes that capture the features of both the
options and the Heston model.

2.1 Prices as labels: Semi-closed formulas

2.1.1 Vanilla Options Pricing

Let $V_t$ be the spot price of a European option at time $t$. We define its normalized forward price as
$\hat{V}_{t,T} := \frac{V_t}{KB(t,T)}$. Here, $B(t,T) := E^Q [\frac{\beta_T}{\beta_t} | \mathcal{F}_t]$ denotes the $t$-price of a zero-coupon bond maturing at time $T$. We assume, without loss of generality, that interest rates are deterministic, leading to $B(t,T) = \frac{\beta_T}{\beta_t}$. For the more general case involving stochastic rates, refer to Mrad and Triki [2014].
The pricing equations for the normalized European Call and Put options are in line with the formulas presented by Lewis [2000] and Lipton [2001, 2002].

$$ \hat{C}_t = e^{X_{t,T}} - \pi^{-1} \int_0^{+\infty} \mathcal{R} \left[ e^{z(y)X_{t,T}} \varphi_{t,T} \circ z(y) \right] \|z(y)\|^{-2} dy \quad (3) $$

and thanks to the call-put parity relationship

$$ \hat{P}_t = 1 - \pi^{-1} \int_0^{+\infty} \mathcal{R} \left[ e^{z(y)X_{t,T}} \varphi_{t,T} \circ z(y) \right] \|z(y)\|^{-2} dy \quad (4) $$

$\mathcal{R}$ denotes the real part of a complex number. The mapping $z$ is defined as $z : \mathbb{R} \to \mathbb{C}, y \mapsto 1/2 + iy$, $X_{t,T} := \ln \left( \frac{S_t}{KB(t,T)} \right)$ denotes the forward log-moneyness and $\varphi_{t,T}$ is the t-MGF of $X_{t,T}$ defined as follows:

$$ \varphi_{t,T}(u) := \mathbf{E}^{\mathcal{Q}} \left[ e^{uX_{t,T} | \mathcal{F}_t} \right], u \in \mathbb{C} \quad (5) $$

The above pricing formulas remain valid regardless of the specific dynamics of the underlying asset. However, the ease of numerical computation is primarily determined by the properties of the MGF. If the underlying asset price follows a Lévy process (examples include Black-Scholes, Heston, Bates, Merton's Jump Model, Variance Gamma, etc.), the MGF exhibits an affine structure due to the independent increments property of the Lévy process. Depending on the time dependence of the parameters, the affine character of the MGF in the Heston model can be leveraged to value European Call and Put options via Equations (3) and (4).

### 2.1.2 Time-dependent Heston model

Consider a fully filtered probability space given by $(\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \ge 0}, \mathcal{Q})$, satisfying the standard conditions. In this context, $\mathcal{Q}$ denotes the risk-neutral measure corresponding to the numéraire of the savings account $\beta_t := e^{\int_0^t r_s ds}$, where $(r_t)$ is a risk-free interest rate process. We turn our attention to a time-dependent Heston model. For the sake of notational simplicity, we omit the second index in the process $X_{t,T}$ defined previously. The two-dimensional process encompassing $(X_t, \nu_t)$ fulfills the following SDE:

$$ \begin{aligned} dX_t &= \frac{1}{2}\nu_t dt + \sqrt{\nu_t}dW_t^s, \\ d\nu_t &= \kappa(t)(\bar{\nu}(t) - \nu_t)dt + \gamma(t)\sqrt{\nu_t}dW_t^{\nu}, \\ d\langle W_t^s, W_t^{\nu} \rangle &= \rho(t)dt, \end{aligned} \quad (6) $$

where $\kappa, \gamma, \bar{\nu}$, and $\rho$ are deterministic functions that vary with time. From the dynamics presented above, we observe that the two-dimensional process $(X_t, \nu_t)$ is affine and that its t-MGF is:

$$ \varphi(u, v, t, T) := \mathbf{E}^{\mathcal{Q}} \left[ e^{uX_T + v\nu_T |\mathcal{F}_t|} \right]; u, v \in \mathbb{C}, $$

satisfying:

$$ \varphi(u, v, t, T) = e^{A(u,v,t,T)+B(u,v,t,T)X_t+C(u,v,t,T)\nu_t}. $$

Here, $A, B$, and $C$ are the solutions of a Ricatti system of ODEs that are determined by recognizing the function $\varphi$ as satisfying the Backward Kolmogorov equation:

$$ \frac{\partial \varphi}{\partial t}(u, v, t, T) + \mathcal{L}\varphi(u, v, t, T) = 0; \quad \varphi(u, v, T, T) = e^{uX_T + v\nu_T}. \quad (7) $$

with, $\mathcal{L}$ is the infinitesimal generator of the process $(X_T, \nu_T)$:

$$ \mathcal{L} := -\frac{1}{2}\nu \frac{\partial}{\partial x} + \kappa(\bar{\nu} - \nu) \frac{\partial}{\partial \nu} + \gamma\rho\nu \frac{\partial^2}{\partial x \partial \nu} + \frac{1}{2}\nu \frac{\partial^2}{\partial x^2} + \frac{1}{2}\gamma^2\nu \frac{\partial^2}{\partial \nu^2}. \quad (8) $$

When the Heston parameters exhibit a general time dependence, the associated Ricatti system might not have an analytical solution. This lack of a solution complicates the Fourier-based approach to pricing European options. In the following paragraphs, we provide a detailed explanation of the analytical solution MGF for both constant and piecewise-constant parameters.
### 2.1.2.1 Constant parameters case

Solving Equation (7) for constant Heston parameters leads to the following Ricatti system:

$$ \frac{\partial A}{\partial t} = -\kappa \bar{\nu} C; \qquad A(u, v, T, T) = 0, \qquad (9) $$

$$ \frac{\partial B}{\partial t} = 0; \qquad B(u, v, T, T) = u, \qquad (10) $$

$$ \frac{\partial C}{\partial t} = -\gamma^2 C^2/2 + \kappa C - \rho\gamma BC + (B - B^2)/2; \qquad C(u, v, T, T) = v. \qquad (11) $$

The system above can be readily solved, yielding:

$$ A(u, v, t, T) = \kappa \bar{\nu} \left( \tau r_{-}(u) - \frac{2}{\gamma^2} \ln \left( \frac{1 - g(u, v)e^{-\delta(u)\tau}}{1 - g(u, v)} \right) \right); \qquad (t, T) \in \mathbb{R}^2, \qquad (12) $$

$$ B(u, v, t, T) = u; \qquad (t, T) \in \mathbb{R}^2, \qquad (13) $$

$$ C(u, v, t, T) = \frac{r_{+}(u)g(u, 0) - g(u, v)e^{-\delta(u)\tau}}{1 - g(u, v)e^{-\delta(u)\tau}}; \qquad (t, T) \in \mathbb{R}^2. \qquad (14) $$

where $\tau = T - t$ and:

$$ \delta(u) = \sqrt{(\kappa - \gamma \rho u)^2 - \gamma^2 u(u - 1)} $$

$$ r_{\pm}(u) = (\kappa - \gamma \rho u \pm \delta(u)) / \gamma^2 $$

$$ g(u, v) = (r_{-}(u) - v) / (r_{+}(u) - v) $$

### 2.1.2.2 Piecewise parameters case

We suppose the existence of at least one ($n \ge 1$) discontinuity times, which fall within the interval of $t$ and $T$ ($t \le T_1 < \dots < T_n \le T$). The functions $\kappa, \rho, \bar{\nu}$, and $\gamma$ remain constant within each interval $[T_j, T_{j+1}]$, where $j$ ranges from 0 to $n$ (with $T_0 = t$ and $T_{n+1} = T$). The value of $\varphi$ is calculated using an analytical backward recursion presented in Mrad and Triki [2014]. We utilize the outcomes from the previous section and apply the Markov chain rule. For any complex number $u \in \mathbb{C}$ and integer $k \in \langle 0, n \rangle$, the following holds:

**For $k = n$:**

$$ \begin{aligned} \varphi(u, T_n, T) &= \mathbf{E}^{\mathbf{Q}} [e^{-uX_T} | \mathcal{F}_{T_n}] \\ &= \varphi(u, 0, T_n, T) \\ &= e^{A_n(u, 0) + uX_{T_n} + C_n(u, 0)\nu_{T_n}} \end{aligned} $$

**For $k < n$:**

$$ \begin{aligned} \varphi(u, T_k, T) &= \mathbf{E}^{\mathbf{Q}} [\mathbf{E}^{\mathbf{Q}} [e^{-uX_T} | \mathcal{F}_{T_{k+1}}] | \mathcal{F}_{T_k}] \\ &= e^{A_{k+1}(u, 0)} \varphi(u, C_{k+1}(u, 0), T_k, T_{k+1}) \\ &= e^{A_k(u, 0) + uX_{T_k} + C_k(u, 0)\nu_{T_k}} \end{aligned} $$

The sequences $\{A_j\}_{j \in \langle 0, n \rangle}$ and $\{C_j\}_{j \in \langle 0, n \rangle}$ are computed backward as follows:

$$ \bullet \ A_n(u, v) = A(u, v, T_n, T) \text{ and } C_n(u, v) = C(u, v, T_n, T); \ u, v \in \mathbb{C} $$

* The functions $A_j$ and $C_j$ are determined for $j = n, \dots, 1$ via:

$$ C_{j-1}(u, v) = C(u, C_j(u, v), T_{j-1}, T_j); \quad u, v \in \mathbb{C} \qquad (15) $$

$$ A_{j-1}(u, v) = A_j(u, v) + A(u, C_j(u, v), T_{j-1}, T_j); \quad u, v \in \mathbb{C} \qquad (16) $$

As a result, we find:

$$ \varphi(u, t, T) = e^{A_0(u, 0) + uX_T + C_0(u, 0)\nu_T}; \quad u \in \mathbb{C} \qquad (17) $$

$\varphi(u, t, T)$ retains the same affine representation as proposed in Section 2.1.2.1. The distinction lies in the computation of the functions $A_0: u, v \mapsto A_0(u, v)$ and $C_0: u, v \mapsto C_0(u, v)$, which is now performed through a backward recursion.
## 2.2 Payoff realizations as labels: MC simulation

Consider a discrete time grid $T = \{t_i\}_{i \in \{1, N\}}$. Our objective is to generate a random trajectory of the process $(X_t, \nu_t)$ for a given set of Heston model parameters (attributes), starting from initial values $(X_0, \nu_0)$, which are derived from market data attributes. Utilizing the option's payoff, which may be path-dependent, we derive a label. We will also compute the sensitivities of this label with respect to its attributes.

To construct a random path of $(X_t, \nu_t)$, we sample $(X_{t_{i+1}}, \nu_{t_{i+1}})$ given $(X_{t_i}, \nu_{t_i})$ for an arbitrary time increment $\delta_i = t_{i+1} - t_i$. By iterative application of this single-period scheme, a complete trajectory $(X_t, \nu_t)_{t \in T}$ can be generated. When simulating a Heston model with piecewise-constant parameters, it is crucial to incorporate discontinuity dates into the discretization time grid. These parameters remain constant between consecutive discontinuity dates, which significantly simplifies the simulation procedure. In this work, we use the Quadratic Exponential (QE) scheme introduced by Andersen [2006] to discretize the Heston model. Its main steps are outlined below.

### 2.2.1 Simulation of $X$

We will use $X$ and $\nu$ to represent discrete-time approximations of the asset log-price and volatility respectively. Similarly to Broadie and Kaya [2006], we can compute $X_{t_{i+1}}$ from $X_{t_i}$ and $\nu_s$, $s \in [t_i, t_{i+1}]$:

$$X_{t_{i+1}} = X_{t_i} + \rho\gamma^{-1}(\nu_{t_{i+1}} - \nu_{t_i} - \kappa\bar{\nu}\delta_i) + (\kappa\rho\gamma^{-1} - 1/2) \int_{t_i}^{t_{i+1}} \nu_u du + \sqrt{1-\rho^2} \int_{t_i}^{t_{i+1}} \sqrt{\nu_u} dW_u^s \quad (18)$$

Given $\nu_{t_i}, \nu_{t_{i+1}}$ and $\int_{t_i}^{t_{i+1}} \nu_u du, X_{t_{i+1}} - X_{t_i}$ is Gaussian with variance $(1 - \rho^2) \int_{t_i}^{t_{i+1}} \nu_u du$. To simulate $X_{t_{i+1}}$, we need to:

1. Simulate $\nu_{t_{i+1}}$ given $\nu_{t_i}$.

2. Simulate $\int_{t_i}^{t_{i+1}} \nu_u du$ given $\nu_{t_i}$ and $\nu_{t_{i+1}}$.

Numerous studies, such as Chan and Joshi [2013], have addressed these steps comprehensively. Two main discretization methods are:

* **Exact:** Developed by Broadie and Kaya [2006] and improved by Glasserman and Kim [2011] and Chan and Joshi [2013], precisely samples $\nu_{t_{i+1}}$ and $\int_{t_i}^{t_{i+1}} \nu_u du$, but requires more time due to its reliance on a noncentral chi-square distribution.

* **Quasi exact:** These faster methods approximate $\nu_{t_{i+1}}$ by sampling with distributions that match its first two moments, with an adjustment term to correct for martingale deviations. Examples include the QE scheme proposed by Andersen [2006], its extension by Mrad and Triki [2014], and Alfonsi's second order scheme (Alfonsi [2010]).

### 2.2.2 Simulation of $\nu$

The instantaneous volatility process, $\nu$, is frequently referred to as the mean-reverting square root process. It possesses several well-known characteristics, some of which are listed below. More details can be found in the seminal works of Cox et al. [1985] and Dufresne [2001].

**Proposition 2.1.** Let $G(y, m, \alpha)$ represent the cumulative distribution function for a non-central chi-square random variable with $m$ degrees of freedom and non-centrality parameter $\alpha$:

$$G(y, m, \alpha) = e^{-\alpha/2} \sum_{k=0}^{\infty} \frac{(\alpha/2)^k}{k! 2^{m/2+k} \Gamma(m/2+k)} \int_{0}^{y} x^{m/2+k-1} e^{-x/2} dx$$

Let's define $d$ and $n(\tau)$ as follows:

$$d = \frac{4\kappa\bar{\nu}}{\gamma^2}, \quad n(\tau) = \frac{4\kappa e^{-\kappa\tau}}{\gamma^2(1 - e^{-\kappa\tau})}, \quad \tau > 0$$

Assuming $T = t + \tau$, conditional on $\nu_t, \nu_T$ is distributed as $e^{-\kappa\tau}/n(\tau)$ times a non-central chi-square random variable with $d$ degrees of freedom and non-centrality parameter $\nu_t n(\tau)$. That is:

$$Q(\nu_T < x | \nu_t) = G\left(\frac{x.n(\tau)}{e^{-\kappa\tau}}, d, \nu_t.n(\tau)\right)$$
In the subsequent corollary, we provide the conditional mean and variance of $\nu_T$ given $\nu_t$.

**Corollary 1.** For two sequential discretization times $t < T$ and $\tau = T - t$, given $\nu_t$, we have:

$$ \mathbf{E}(\nu_T|\nu_t) := \mu = \bar{\nu} + (\nu_t - \bar{\nu})e^{-\kappa\tau}, \qquad (19) $$

$$ \mathbf{V}(\nu_T|\nu_t) := \sigma^2 = \frac{\gamma^2 \nu_t e^{-\kappa \tau}}{\kappa} (1 - e^{-\kappa \tau}) + \frac{\gamma^2 \bar{\nu}}{2\kappa} (1 - e^{-\kappa \tau})^2. $$

For two successive discretization times $t < T$, $\nu_T|\nu_t$ can be generated by inverting the cumulative distribution $x \mapsto \mathbf{Q}(\nu_T < x|\nu_t)$, as detailed in Proposition 2.1. The QE scheme approximates $\nu_T$ using a distribution that matches the first two moments specified in Equation (19). The non-central chi-square distribution approaches a Gaussian distribution as the non-centrality parameter increases indefinitely. When $\nu_t$ is large enough, $\nu_T|\nu_t$ can be approximated as a Gaussian variable with the first two moments given in Equation (19). To maintain non-negativity, Andersen [2006] utilizes a squared Gaussian random variable:

$$ \hat{\nu}_T|\hat{\nu}_t = a(b + \varepsilon^\nu)^2 \qquad (20) $$

For smaller $\nu_t$, the non-centrality parameter tends towards zero, and the distribution of $\nu_T|\nu_t$ becomes proportional to a regular chi-square distribution with $\frac{4\kappa\nu}{\gamma^2}$ degrees of freedom. Practically, when $\frac{4\kappa\bar{\nu}}{\gamma^2} \ll 2$, this suggests that, for smaller $\nu_t$, the density of $\nu_T$ will concentrate significantly around 0. Using the asymptotic chi-square density as a reference, an approximate density is proposed:

$$ \mathbf{Q}(\hat{\nu}_T \in [x; x + dx] | \hat{\nu}_t) \approx (q\delta_0(x) + (1-q)\phi e^{-\phi x}) dx; x > 0 \qquad (21) $$

Here, $\delta_0$ denotes the standard Dirac function, while $q$ and $\phi$ represent constants that need to be determined. The occurrence of a probability mass at the origin is noteworthy. The cumulative distribution function can be defined as:

$$ F(x) = \mathbf{Q}(\hat{\nu}_T \le x | \hat{\nu}_t) = q + (1-q)(1 - e^{-\phi x}); \quad x > 0 $$

The inverse of this function is expressed as:

$$ F^{-1}(u) = \begin{cases} 0 & \text{if } 0 \le u \le q \\ \frac{1}{\phi} \ln \left( \frac{1-q}{1-u} \right) & \text{if } q < u \le 1 \end{cases} \qquad (22) $$

The coefficients $a, b, q$, and $\phi$ are defined in the following proposition:

**Proposition 2.2.** Let $\mu$ and $\sigma$ be as defined in (19) and let $\lambda = \sigma^2/\mu^2$.

* For $\lambda \le 2$, we have $b^2 = 2\lambda^{-1} - 1 + \sqrt{2\lambda^{-1}\sqrt{2\lambda^{-1} - 1}}$ and $a = \frac{\mu}{1+b^2}$.

* For $\lambda \ge 1$, we have $q = \frac{\lambda^{-1}}{\lambda+1} \in [0, 1]$, and $\phi = \frac{2}{\mu(\lambda+1)} > 0$.

The QE algorithm is summarized in the following steps, where $\lambda_c \in [1, 2]$ is an arbitrary level (e.g. $\lambda_c = 1.5$):

**Algorithm 1** Simulation using QE algorithm

1: Compute $\mu$ and $\sigma^2$ using Equation (19) given $\hat{\nu}_t$.

2: Compute $\lambda = \sigma^2/\mu^2$.

3: Generate a standard uniform random variable $U^\nu$.

4: **if** $\lambda \le \lambda_c$ **then**

5: Determine $a$ and $b$ using Proposition 2.2.

6: Calculate $\varepsilon^\nu = \mathcal{N}^{-1}(U^\nu)$, where $\mathcal{N}$ is the standard Gaussian cumulative function.

7: Employ Equation (20) and set $\hat{\nu}_T = a(b + \varepsilon^\nu)^2$.

8: **else**

9: Determine $\phi$ and $q$ using Proposition 2.2.

10: Employ Equation (22) and set $\hat{\nu}_T = F^{-1}(U^\nu)$.

11: **end if**
## 2.3 Artificial Neural Networks

In recent years, several works have leveraged standard DL techniques, training ANNs to approximate the valuation functions of European Vanilla Options when the dynamics of the underlying risk factors do not yield an explicit pricing formula. Using trained ANNs allows for ultra-fast calibration of risk factor dynamics to implied volatility surfaces, as first introduced by Hernandez [2017]. Further developments include Bayer et al. [2019] and Horvath et al. [2021b] for the Rough Bergomi model and Liu et al. [2019] for the Heston model with constant parameters.

The main limitation of the above-mentioned approaches is the requirement to generate large datasets of accurate option prices. Depending on the complexity of the dynamics of the risk factors involved, this process can take several days. In contrast, we use an approach that avoids the need to generate price datasets. Instead, our method favors the use of payoff realizations for a given set of attributes (risk factors dynamics parameters, initial market data, option strike price, and maturity).

We focus on FFNNs, known for approximating high-dimensional functions. An FFNN with $K$ layers maps an input vector $\xi \in \mathbb{R}^n$ to an output $\zeta \in \mathbb{R}$. At layer $k$, weights $w_k \in \mathbb{R}^{n_{k-1} \times n_k}$ and biases $b_k \in \mathbb{R}^{n_k}$ define the transformation:

$$h_0 = \xi, \quad h_k = \sigma_k(h_{k-1})w_k + b_k, \quad k \in \langle 1, K \rangle, \quad \zeta = h_k. \qquad (23)$$

Here, $h_k \in \mathbb{R}^{n_k}$ represents the $n_k$ neurons in layer $k$, and $\sigma_k: \mathbb{R} \to \mathbb{R}$ is the identity function for $k=0$ and a non-linear activation function otherwise. The first ($k=0$) and last ($k=K$) layers are, respectively, the input and output layers, while $k \in \langle 1, K-1 \rangle$ are the hidden layers.

### 2.3.1 Differential Deep Learning

In their seminal paper, Czarnecki et al. [2017] introduced the Sobolev Training method, a technique in which models are trained to approximate target functions by aligning both their output values and their derivatives to arbitrary degrees. This approach found application in computational finance under the name of Differential Machine Learning (DML), and was first implemented by Huge and Savine [2020]. We adopt the more specific term Differential Deep Learning (DDL) to describe this methodology. DDL is especially advantageous for complex stochastic models that lack closed-form solutions. It extends traditional deep learning applied to FFNNs by incorporating gradients of labels w.r.t. input into the training process. The following diagram illustrates how DDL is applied to a single-output FFNN. Beyond the traditional output node, represented by $\zeta$, DDL incorporates a vector output denoted by $\tilde{\xi}$. This new output is the gradient of $\zeta$ w.r.t. the input layer $\xi$. As detailed in Section 2.3.1.1, both $\zeta$ and $\tilde{\xi}$ are used during training.

<figure><img src="image_2.png" /><figcaption>Figure 1: Extended FFNN for DDL</figcaption></figure>

When labels are generated from payoff realizations, the training process displays a higher variance than when prices are used, leading to noisier optimizer updates. This challenge can be addressed by expanding the dataset and using larger batch sizes during training loss optimization. If data is scarce or the generation of additional data is time consuming, regularization strategies, such as adding a L1 (Lasso) or L2 (Ridge) penalty term, can be beneficial. DDL includes a penalty term that accounts for the aggregate of weighted
mean squared errors (MSE). This error is computed between the derivatives of the ANN output and the derivatives of the labels, with both derivatives taken w.r.t. the features. This method enhances learning and acts as an effective data augmentation technique, expanding the dataset without introducing bias.

### 2.3.1.1 Training loss

We illustrate the DDL training approach by considering a function $\pi$ as the learning target. Rather than using only the output values $\pi(\xi_i)$ for the training points $\xi_i$, we can also leverage the values of its first-order derivatives w.r.t. the input, denoted as $\nabla_{\xi} \pi(\xi_i)$. Hence, instead of using a standard training set of pairs $\{(\xi_i, \pi(\xi_i))\}_{i \in \langle 1, n \rangle}$, we use a triplet $\{(\xi_i, \pi(\xi_i), \nabla_{\xi} \pi(\xi_i))\}_{i \in \langle 1, n \rangle}$.

Given a neural network model $\mathcal{M}$ parameterized by $\omega$, our objective is to minimize the empirical error for $\pi$ following a specific loss function $l_0$: $\frac{1}{N} \sum_{i=1}^{N} l_0(\mathcal{M}(\xi_i|\omega), \pi(\xi_i))$. However, in DDL, this goal evolves into:

$$ \frac{1}{N} \sum_{i=1}^{N} [l_0(\mathcal{M}(\xi_i|\omega), \pi(\xi_i)) + l_1(\nabla_{\xi} m(\xi_i|\omega), \nabla_{\xi} \pi(\xi_i))] \quad (24) $$

Here, $l_1$ represents a loss function that quantifies the error in the first-order derivatives. This adjustment encourages the neural network to embed the derivatives of the target function into its own derivatives. In particular, such a model can be trained employing traditional backpropagation and customary optimizers. When using the MSE loss function, the DDL training loss is defined as:

$$ \mathbf{E} \left[ (\mathcal{M}(\xi|\omega) - \pi(\xi))^2 \right] + \sum_{j=1}^{J} \alpha_j E \left[ (\partial_j \mathcal{M}(\xi|\omega) - \partial_j \pi(\xi))^2 \right] $$

Where, $\forall j \in \langle 1, J \rangle$, $\alpha_j \in \mathbb{R}^{+,*}$ and $\partial_j$ is the derivative w.r.t. the j-th component of $\xi$ (that has dimension J).

### 2.3.1.2 Differentials of the Prediction

In this section, we detail the derivation of the first- and second-order differentials (gradient and Hessian) of the output value $\zeta = h_k$ from the model $\mathcal{M}$ w.r.t. the inputs $\xi = h_0$.

**First Order Differentials:** To compute the gradient of $\mathcal{M}$ w.r.t. inputs $h_0$, denoted as $\bar{h}_0$, we differentiate Equation (23) in reverse order for layers $k \in \{K, \dots, 1\}$:

$$ \bar{h}_k = \bar{\zeta} = \mathbf{1} \quad ; \quad \bar{h}_{k-1} = (\bar{h}_k w_k^T) \otimes \sigma_k'(h_{k-1}) \mathbf{1} \quad ; \quad \bar{\xi} = \bar{h}_0 \quad (25) $$

Here, $\bar{\zeta} = \frac{\partial \zeta}{\partial \zeta} = \mathbf{1}$, and $\mathbf{1}$ is a tensor of ones that match the shape of $h_K$. The symbol $\otimes$ denotes the element-wise product. We can construct an FFNN using backpropagation as described in Equation (25). This network, which shares weights with the original one defined by Equation (23), employs neurons that serve as adjoints to the original network.

**Second Order Differentials:** In this work, the DDL approach utilizes first-order derivatives, although Czarnecki et al. [2017] suggests that it is theoretically possible to include higher-order derivatives for labels, as long as they can be efficiently computed. Here, we provide the Hessian of the ANN with respect to its input, as this result will be crucial when implementing PINNs<fnref n="4" />.

The Hessian, denoted as $\bar{\bar{\xi}}$, provides the 2nd order differentials. It can be obtained by differentiating Equation (23) twice in reverse order. A straightforward application of differentiation rules yields the following:

$$ \bar{\bar{\xi}} = \sum_{j=1}^{K} \bar{h}_j w_0^T $$

<footnote n="4">⁴For continuity, we discussed this result in the background section, following the paragraph on the FFNN's first-order derivative with respect to its input. To our knowledge, this is a novel contribution of this paper.</footnote>
∀j ∈ (1, K), the term $\bar{\bar{h}}_j$ is determined backward. For the topmost layer l = K, we initialize:

$$ \bar{\bar{h}}_j = 1 $$

Subsequently, for layers below, we differentiate according to:

$$ \bar{\bar{h}}_j = \begin{cases} \left( \bar{\bar{h}}_j w_l^T \right) \otimes \sigma_l'(h_l) & \text{if } l \neq j \\ \left( \bar{\bar{h}}_j w_j^T \right) \otimes \frac{\partial h_j}{\partial \xi} \otimes \sigma_j''(h_j) & \text{otherwise} \end{cases} \quad (26) $$

Furthermore, the derivatives w.r.t. the input, $\frac{\partial h_j}{\partial \xi}$, are determined recursively by:

$$ \frac{\partial h_{j+1}}{\partial \xi} = \begin{cases} \frac{\partial h_j}{\partial \xi} \otimes \sigma'(h_j) w_n^T & \text{for } j > 1 \\ w_0^T & \text{if } j = 1 \end{cases} $$

In the above expressions:

* $\bar{\bar{\xi}}$ represents the Hessian or second-order differential of $\zeta$ relative to $\xi$.

* $w_l$ are the weights associated with the $l^{th}$ layer.

* $\sigma'_l$ and $\sigma''_l$ are the 1st and 2nd derivatives of the activation function applied to neurons in the layer $l$.

### 2.3.2 Physics-Inspired Neural Networks (PINNs)

PINNs are machine learning models that integrate established physical laws, typically expressed as PDEs, into their learning process. This crucial integration acts as a powerful regularizer, steering the network towards solutions that adhere to the underlying physics. This not only improves the accuracy and reliability of the model, but also allows it to perform well even with limited training data. As a consequence, PINNs have gained attention in various fields, including fluid dynamics and quantum physics, demonstrating their promise for modeling complex systems.

#### 2.3.2.1 Objective

PINNs seamlessly integrate differential equations, representing physical laws, into the structure of the neural network, as shown in Raissi et al. [2019]. To illustrate, consider the following canonical PDE:

$$ \begin{aligned} \mathcal{N}[u(\xi)] &= 0, & \xi \in D, \\ \mathcal{B}[u(\xi)] &= h(\xi), & \xi \in \partial D, \end{aligned} \quad (27) $$

where $\mathcal{N}$ is a differential operator. The function $u(x)$ is the latent solution of the PDE with $D$ representing the space-time domain. The operator $\mathcal{B}$ indicates the arbitrary initial or boundary conditions associated with the problem and $h$ means the boundary function. In particular, boundary conditions could manifest as Dirichlet, Neumann, Robin, or periodic. Such a formulation is capable of representing a wide range of problems in mathematical physics, encompassing conservation laws, diffusion processes, among others.

In the PINNs framework, the function $u$ is approximated using an ANN, parameterized by a set of parameters $p$. The neural network emulates the differential equations by finding the optimal parameters $p^*$. This is achieved by minimizing a composite loss function that includes the differential equation $L_{\mathcal{N}}$, the boundary conditions $L_{\mathcal{B}}$, and any available data $L_{\mathcal{D}}$. Each component of the loss function is appropriately weighted:

$$ p^* = \arg\min_{p} (\omega_{\mathcal{N}} L_{\mathcal{N}}(p) + \omega_{\mathcal{B}} L_{\mathcal{B}}(p) + \omega_{\mathcal{D}} L_{\mathcal{D}}(p)). \quad (28) $$

By integrating the residuals of the differential equations into the loss function, PINNs ensure that the model consistently abides by fundamental physical laws throughout its learning process.
<figure><img src="image_3.png" /><figcaption>Figure 2: Extended FFNN for PINNs</figcaption></figure>

### 2.3.2.2 Heston Model PDE

In the context of the Heston model with piecewise parameters, we focus on a specific time interval, $[T_i, T_{i+1}]$. Here, the consecutive times $T_i$ and $T_{i+1}$ are part of the time grid discussed in Section 2.1.2.2.

The price $u(t, x, v)$ of a European option, with payoff $\pi(X_T)$ and maturing at $T$, satisfies the following PDE:

$$ \frac{\partial u}{\partial t} + \mathcal{L}_i[u] = 0; \quad t \in [T_i, T_{i+1}], \qquad (29) $$

with the operator $\mathcal{L}_i$ defined as:

$$ \mathcal{L}_i := -\frac{1}{2}\nu\frac{\partial}{\partial x} + \kappa_i(\bar{\nu}_i - \nu)\frac{\partial}{\partial \nu} + \gamma_i\rho_i\nu\frac{\partial^2}{\partial x\partial \nu} + \frac{1}{2}\nu\frac{\partial^2}{\partial x^2} + \frac{1}{2}\gamma_i^2\nu\frac{\partial^2}{\partial \nu^2}. \qquad (30) $$

Referring back to Section 2.1.2.2, $u(t, x, v)$ can be expressed as a function of $\tau_0, \dots, \tau_n$ where $\tau_j = T_{j+1} - T_j$, with $T_0 = t$ and $T_{n+1} = T$. This expression simplifies to:

$$ -\frac{\partial u}{\partial \tau_i} + \mathcal{L}_i[u] = 0. \qquad (31) $$

For the time-dependent Heston model, the residual PINN loss $L_N$ derived from Equation (28) can be represented as $L_N = \sum_i L_{N_i}$, where the specific loss term is:

$$ L_{N_i} = \text{MSE}(-\partial_{\tau_i} u + \mathcal{L}_i[u]). \qquad (32) $$

This PINN loss, which encapsulates the Heston model's parameters over time, acts as an important regularizer during the learning process.

### 2.3.2.3 Training/Collocation points

PINNs utilize the concept of collocation points, which significantly enforce PDEs as constraints. The key features of these collocation points are the following:

* Collocation points in PINNs are a set of points within the input space where the residual of the PDE is applied as a constraint. The concept of these points is derived from the spectral methods of numerical analysis used to solve differential equations.

* These points are advantageous as they do not necessitate labeled data, only the input parameters. The computational cost primarily stems from calculating the network's differentials to express the PDE loss. This cost-effective strategy differentiates PINNs from traditional supervised learning methods.
* Increasing the density of sampled points, particularly near the boundaries of the parameter space, proves beneficial. This tactic addresses areas with sparse training data and complex model behaviors, where the model exhibits high sensitivity to small changes in input parameters.

* Despite skewing the collocation points towards certain areas of the input space, the i.i.d. property is preserved as long as they are randomly drawn from the same probability distribution. This ensures the validity of the statistical assumptions that underlie the learning process.

* Latin Hypercube sampling is used to effectively span the search space. Subsequently, inverse sampling replaces the inner uniform distribution with a Beta distribution. This distribution, parameterized by α and β, allows the densification of constraint regions for experimental purposes.

* In high-dimensional spaces, resampling new collocation points at different training iterations can be an effective and inexpensive strategy to better cover the domain. This can ensure more comprehensive and accurate learning, particularly for complex physical systems modeled in high dimensions.

### 2.3.3 Overall Training Loss

We define our total training loss, $L_{TOT}$, as a composition of three terms:

$$L_{TOT} = \underbrace{\omega_{STD}L_{STD} + \omega_{DDL}L_{DDL}}_{\text{Trained with labeled data}} + \underbrace{\omega_{PINN}L_{PINN}}_{\text{Trained with unlabeled data}} \quad (33)$$

The weights $\omega_{STD}$, $\omega_{DDL}$, and $\omega_{PINN}$ (which belong to $\mathbb{R}^+$) are designed to adjust the contribution of each component during training, possibly depending on the training epoch. Given our context of a Heston model with piecewise parameters, the losses $L_{DDL}$ and $L_{PINN}$ are computed as accumulations of elementary losses: $L_{DDL} = \sum_{i} L_{i,DDL}$ and $L_{PINN} = \sum_{j} L_{j,PINN}$.

For $L_{DDL}$, we focus solely on the penalization term that accounts for the discrepancy between the label differentials and the ANN differentials w.r.t. the Heston model parameters. In contrast, $L_{PINN}$ employs only the collocation points, i.e. the unlabeled data, to address the problem in inverse mode. This involves exclusively the boundary loss $L_B$ and the residual loss $L_N$ from the PDE (both detailed in Equation (28)).

<figure><img src="image_4.png" /><figcaption>Figure 3: Extended FFNN for PINNs and DDL</figcaption></figure>
In the following, we detail each of the main training loss components:

* $L_{STD}$ represents the MSE between the labels - obtained from the option's payoff $\pi$ for every underlying asset MC path and each instance of Heston model parameters - and the ANN, which uses the Heston model parameters as input.

* $L_{i,DDL}$ constitutes the MSE between the differential of the labels and the differential of the ANN w.r.t. the $i$-th feature, i.e., the $i$-th Heston model parameter.

* $L_{j,PINN}$ is the MSE between zero and the Heston differential operator linked with the $j$-th time interval $[T_j, T_{j+1}]$, as elaborated in Equation (31). Importantly, this loss component only involves the ANN and is independent of labeled data.

As mentioned in the previous paragraph, there might be situations where it is beneficial to position collocation points both within the feature domain and close to its borders. By employing a Beta distribution, we can sample i.i.d. collocation points near the boundaries of the feature domain.

Consequently, when necessary, we split the PINN loss into two distinct components. The first penalization, denoted by $L_{inner_PDE}$, targets unlabeled collocation points located within the interior of the feature domain. In contrast, the second penalization, represented by $L_{brdrs_PDE}$, focuses on points near the domain boundaries. These penalties are weighted by $\omega_{inner_PDE}$ and $\omega_{brdrs_PDE}$, respectively.

The following figure visually provides the positioning of these collocation points within the feature domain.

<figure><img src="image_5.png" /><figcaption>Figure 4: Collocation points distributed in the feature domain's interior and its periphery.</figcaption></figure>

# 3 Prices as Labels

## 3.1 Differentials w.r.t. $\bar{\nu}, \kappa, \gamma$ and $\rho$

Given the European option prices, as described by Equations (3) and (4), we seek to compute the sensitivities w.r.t. the parameters $\bar{\nu}, \kappa, \gamma$ and $\rho$. These parameters can be piecewise. To do this, we must differentiate the MGF w.r.t. $\theta \in \{\bar{\nu}, \kappa, \gamma, \rho\}$, yielding:

$$ \partial_{\theta} \hat{C}_t = \partial_{\theta} \hat{P}_t = -\pi^{-1} \int_{0}^{+\infty} \mathcal{R} \left[ e^{z(y)X_{t,T}} \partial_{\theta} \varphi_{t,T} \circ z(y) \right] \|z(y)\|^{-2} dy \quad (34) $$
### 3.1.1 Case with Constant Parameters

To streamline the notation, we omit the dependence on $u, v$, and $\tau = T-t$, letting the MGF $\varphi = e^{A+BX_t+C\nu_t}$, where $A, B$, and $C$ are defined in Equation (12). We can write:

$$
\begin{align*}
\beta &= \kappa - u\gamma\rho \\
\delta &= \sqrt{\beta^2 - \gamma^2 u(u-1)} \\
r_{\pm} &= (\beta \pm \delta)\gamma^{-2} \\
g &= (r_{-} - v)/(r_{+} - v) \\
z &= (1 - ge^{-\delta\tau})/(1 - g) \\
A &= \kappa\bar{\nu}(\tau r_{-} - 2\gamma^{-2}\log(z)) \\
C &= (r_{-} - ge^{-\delta\tau})/(1 - ge^{-\delta\tau})
\end{align*}
$$

Our current interest lies in evaluating $\partial_{\theta} A$, where $\theta \in \{\bar{\nu}, \kappa, \gamma, \rho\}$. For $\theta = \bar{\nu}$, it is straightforward to show that $\partial_{\bar{\nu}} A = A/\bar{\nu}$ and $\partial_{\bar{\nu}} C = 0$. For $\theta \in \{\kappa, \gamma, \rho\}$, we calculate $\partial_{\theta}$ as follows:

$$
\begin{align*}
\partial_{\theta}\beta &= 1_{\theta=\kappa} - u(\rho1_{\theta=\gamma} - \gamma1_{\theta=\rho}) \\
\partial_{\theta}\delta &= \delta^{-1}(\partial_{\theta}\beta \cdot \beta - \gamma u(u-1)1_{\theta=\gamma}) \\
\partial_{\theta}r_{\pm} &= -2r_{\pm} \cdot \gamma^{-1} \cdot 1_{\theta=\gamma} + \gamma^{-2} (\partial_{\theta}\beta \pm \partial_{\theta}\delta) \\
\partial_{\theta}g &= (\partial_{\theta}r_{-} - g \cdot \partial_{\theta}r_{+})/(r_{+} - v) \\
\partial_{\theta}z &= \frac{g}{1-g} \left( \frac{\partial_{\theta}g}{g^2} (z-1) + \tau \partial_{\theta}\delta e^{-\delta\tau} \right) \\
\partial_{\theta}A &= \kappa\bar{\nu} (\tau \partial_{\theta}r_{-} - 2\gamma^{-2}\partial_{\theta}z/z) + 2\gamma^{-1} (\kappa\bar{\nu}\tau r_{-} - A) \cdot 1_{\theta=\gamma} + \kappa^{-1} A \cdot 1_{\theta=\kappa} \\
\partial_{\theta}C &= (C - r_{+})\delta^{-1}\partial_{\theta}\delta - (\gamma^2/2)\delta^{-1}(C - r_{+})^2 (\partial_{\theta}g + \tau\partial_{\theta}\delta) e^{-\delta\tau} + \partial_{\theta}r_{+}
\end{align*}
$$

### 3.1.2 Case with piecewise Parameters

We suppose at least one discontinuity ($n \ge 1$) between $t$ and $T$, the parameters $\kappa, \rho, \bar{\nu}$, and $\gamma$ remain constant within each interval $[T_j, T_{j+1}]$ ($T_0 = t$ and $T_{n+1} = T$). The $\varphi$ value at $T_k$ is defined as:

$$ \varphi(u, T_k, T) = e^{A_k(u,0)+uX_{T_k}+C_k(u,0)\nu_{T_k}} $$

The sequences $\{A_j\}_{j \in \langle 0, n \rangle}$ and $\{C_j\}_{j \in \langle 0, n \rangle}$ are defined in a backward manner as explained in Section 2.1.2.2. We simplify the notation by suppressing the dependence on $u, v, \tau_j = T_{j+1} - T_j$, and $\theta_j$, where $\theta_j$ is a representative of any one of the four Heston parameters for the interval $[t_j, t_{j+1}]$.

* For $j = n$, we set $A_n = A(u, v, \tau_n, \theta_n)$ and $C_n = C(u, v, \tau_n, \theta_n)$.

* For $j = n, \dots, 1$, we calculate $A_j$ and $C_j$ as:

$$ C_{j-1} = C(u, C_j, \tau_{j-1}, \theta_{j-1}), \qquad (35) $$

$$ A_{j-1} = A_j + A(u, C_j, \tau_{j-1}, \theta_{j-1}). \qquad (36) $$

Using the chain rule of differentiation, we compute the intermediary derivatives of the sequences $\{A_j\}_{j \in \langle 0, n \rangle}$ and $\{C_j\}_{j \in \langle 0, n \rangle}$ w.r.t. $\theta_k$, for $j \in \langle 1, k+1 \rangle$. This results in:

$$ \partial_{\theta_k} C_{j-1} = \partial_v C(u, C_j, \tau_{j-1}, \theta_{j-1}) \cdot \partial_{\theta_k} C_j, \qquad (37) $$

$$ \partial_{\theta_k} A_{j-1} = \partial_{\theta_k} A_j + \partial_v A(u, C_j, \tau_{j-1}, \theta_{j-1}) \cdot \partial_{\theta_k} C_j, \qquad (38) $$

Here, $\partial_v C(u, C_j, \tau_{j-1}, \theta_{j-1})$ represents the derivative of the function $C(u, C_j, \tau_{j-1}, \theta_{j-1})$ w.r.t. its second input, and similarly for $\partial_v A(u, C_j, \tau_{j-1}, \theta_{j-1})$. They are respectively given by:

$$ \partial_v A(u, v, \tau, \theta) = \frac{2\kappa\bar{\nu}}{\gamma^2} \frac{\partial_v g}{(1-g^2)} e^{-\delta\tau} \quad \text{and} \quad \partial_v C(u, v, \tau, \theta) = \frac{(C-1)^2}{r_- - 1} \partial_v g e^{-\delta\tau} $$

with

$$ \partial_v g = \frac{2\beta}{\gamma^2(r_+ - v)^2} $$
The full differentiation is given by:

$$ \partial_{\theta_k} C_0 = \prod_{j=1}^{k} \partial_v C(u, C_j, \tau_{j-1}, \theta_{j-1}) \cdot \partial_{\theta_k} C_k, \quad (39) $$

$$ \partial_{\theta_k} A_0 = \partial_{\theta_k} A_k + \sum_{j=1}^{k} \partial_v A(u, C_j, \tau_{j-1}, \theta_{j-1}) \cdot \left( \prod_{l=j+1}^{k} \partial_v C(u, C_l, \tau_{l-1}, \theta_{l-1}) \right) \cdot \partial_{\theta_k} C_k, \quad (40) $$

In terms of implementation, Equations (37) and (38) are more practical than Equations (39) and (40). The former pair allows for an efficient computation by recursion, whereas the latter requires storing intermediate values of the derivatives in memory, hence increasing computational complexity. It is also worth mentioning that the derivative $\partial_{\theta_k} C_k$ is computed using the results obtained in the case of constant parameters.

## 3.2 Differentials w.r.t. $X_t$ and $\nu_t$

We use (3) and (4) and compute the following derivatives w.r.t. $X_t$:

$$ \partial_{X_t} \hat{P}_t = \partial_{X_t} \hat{C}_t - e^{X_t} = -2\pi^{-1} \int_{0}^{+\infty} \mathcal{R} \left[ z(y) e^{z(y)X_t} \varphi_{t,T} \circ z(y) \right] \|z(y)\|^{-2} dy \quad (41) $$

The derivatives w.r.t. $\nu_t$:

$$ \partial_{\nu_t} \hat{P}_t = \partial_{\nu_t} \hat{C}_t = -\pi^{-1} \int_{0}^{+\infty} \mathcal{R} \left[ e^{z(y)X_t} C_0(z(y), 0) \varphi_{t,T} \circ z(y) \right] \|z(y)\|^{-2} dy \quad (42) $$

## 3.3 Differentials w.r.t. durations $\tau_i$

Since the forward log-moneyness $x_{t,T}$ depends on the residual maturity $T = \sum_{i=0}^{n} \tau_i$, we can determine the sensitivity of the option price to $\tau_i$ using the Backward Kolmogorov PDE equation:

$$ \partial_{\tau_i} \hat{P}_t = \mathcal{L}_i \hat{P}_t, \quad (43) $$

where the infinitesimal generator $\mathcal{L}_i$ of the process $(X_t, \nu_t)$ over the interval $[t_i, t_{i+1}]$ is defined as:

$$ \mathcal{L}_i := -\frac{1}{2}\nu\frac{\partial}{\partial x} + \kappa_i(\bar{\nu}_i - \nu)\frac{\partial}{\partial \nu} + \gamma_i\rho_i\nu\frac{\partial^2}{\partial x\partial \nu} + \frac{1}{2}\nu\frac{\partial^2}{\partial x^2} + \frac{1}{2}\gamma_i^2\nu\frac{\partial^2}{\partial \nu^2}. $$

To compute $\partial_{\tau_i} \hat{P}_t$, we need the log-price derivatives $\partial_{X_t}^2, \partial_{X_t} \partial_{\nu_t}, \partial_{\nu_t}^2$. We obtain:

$$ \partial_{X_t}^2 \hat{P}_t = -4\pi^{-1} \int_{0}^{+\infty} \mathcal{R} \left[ z(y)^2 e^{z(y)X_t} \varphi_{t,T} \circ z(y) \right] \|z(y)\|^{-2} dy, \quad (44) $$

$$ \partial_{X_t} \partial_{\nu_t} \hat{P}_t = -2\pi^{-1} \int_{0}^{+\infty} \mathcal{R} \left[ z(y) C_0(z(y), 0) e^{z(y)X_t} \varphi_{t,T} \circ z(y) \right] \|z(y)\|^{-2} dy, \quad (45) $$

$$ \partial_{\nu_t}^2 \hat{P}_t = -\pi^{-1} \int_{0}^{+\infty} \mathcal{R} \left[ C_0(z(y), 0)^2 e^{z(y)X_t} \varphi_{t,T} \circ z(y) \right] \|z(y)\|^{-2} dy. \quad (46) $$

# 4 Payoff Realizations as Labels

## 4.1 Training with Payoff realizations

In this section, we present how to use payoff realizations as labels instead of prices. Let $\hat{V}_t(\xi)$ denote an option's price, where $\xi$ is the features combination (including Heston model parameters, option's features, Initial market data) used to evaluate the price $\hat{V}_t(\xi)$. This price can be expressed as:

$$ \hat{V}_t(\xi) = \mathbf{E} \left[ \pi \left( X_{T,T}^{\xi} \right) | \mathcal{F}_t \right] \quad (47) $$
where $\pi$ is the option's payoff, and $X_{T,T}^{\xi}$ is derived from the Heston Model dynamics with parameters $\xi$. When using the standard MSE loss function with prices as labels, the training loss function is defined as:

$$ \mathbf{E} \left[ \left( \mathcal{M}(\xi|\omega) - \hat{V}_t(\xi) \right)^2 \right] $$

with $\mathcal{M}$ an ANN parameterized by $\omega$. Applying the tower rule for conditional expectations leads to:

$$ \mathbf{E} \left[ \left( \mathcal{M}(\xi|\omega) - \hat{V}_t(\xi) \right)^2 \right] = \mathbf{E} \left[ \left( \mathcal{M}(\xi|\omega) - \pi(X_{T,T}^{\xi}) \right)^2 \right] - \mathbf{VAR} \left[ \pi(X_{T,T}^{\xi}) \right] $$

Since the variance term in the Equation above is independent of the model to be trained $\mathcal{M}$, minimizing the MSE using the prices as labels is equivalent to using the payoff realizations as labels. In the DDL case, the regularization term in the training loss leverages the sensitivity of the payoff's realization to the features $\xi$.

## 4.2 Payoff Realizations Differentials

In this section, we temporarily set aside the sensitivity of payoff realizations to the durations $\tau_i$. We focus on developing methods to compute the differentials of labels, particularly for cases involving non-differentiable functions. To this end, we use the complex-step differentiation technique, which is well suited for functions that are not smoothly differentiable at certain points, such as the square root function near zero.

### 4.2.1 Complex-Step Differentiation

This approach extends a real-valued function $f(x)$ in the complex plane to evaluate its derivative at a particular point using a complex perturbation. According to Lai and Crassidis [2008], the first and second-order derivatives can be approximated using the following expressions:

$$ f'(x) = \frac{\mathcal{I}\left(f\left(x + i\frac{2}{3}h\right) - f\left(x - i\frac{2}{3}h\right)\right)}{\sqrt{3}h} + O(h^4), \qquad (48) $$

$$ f''(x) = \frac{\mathcal{I}\left(f\left(x + i\frac{1}{2}h\right) + f\left(x - i\frac{1}{2}h\right)\right)}{h^2} + O(h^4). \qquad (49) $$

Here, $\mathcal{I}$ denotes the imaginary part, and $O(h^4)$ represents the error term<fnref n="5" />. The Complex-Step method has several advantages over traditional numerical differentiation:

* **Accuracy:** It provides highly accurate derivatives, even for small step sizes, without suffering from subtractive cancellation errors.

* **Robustness:** It doesn't rely on finite differences, making it more stable and reliable.

* **Speed:** It enables efficient computation of gradient projections on random vectors (see Czarnecki et al. [2017]), particularly useful for high-dimensional input, as shown in Saadeddine [2022].

Implementing Complex-Step differentiation is advantageous but requires careful consideration, particularly for non-smooth payoff functions. A key requirement is to ensure that the computational pipeline supports complex numbers. To enhance differentiability, non-smooth functions are smoothed. For example, the Softplus function approximates the ReLU function used in European option pay-offs, while the sigmoid function approximates the Heaviside function in QE discretization.

Let $\eta$ be a smoothing factor in $\mathbb{R}^+$. The holomorphic extension of the sigmoid function associated with the parameter $\eta$ is relatively straightforward: $\text{Sigmoid}_{\eta}(z) := 1/(1+e^{-\eta z})$, $\forall z \in \mathbb{C}$. It retains many features of the real-valued sigmoid function while introducing more intricate behavior in the complex plane. In contrast, the holomorphic extension of the Softplus function is slightly more involved. It is given by:

$$ \text{Softplus}_{\eta}(z) = z - \eta^{-1} \left( \ln |\text{Sigmoid}_{\eta}(z)| + i \cdot \text{Arg}(\text{Sigmoid}_{\eta}(z)) \right) $$

where $\ln$ denotes the natural logarithm, $|\cdot|$ represents the complex modulus, and Arg is the principal argument defined on $(-\pi, \pi]$. This formulation provides a smooth and analytical extension of the Softplus function to the complex plane, maintaining holomorphicity wherever $\text{Sigmoid}_{\eta}(z) \neq 0$.

<footnote n="5">⁵We use Equations (48) and (49), which have an error order of $h^4$. The parameter $h$ is chosen between $10^{-4}$ and $10^{-3}$ aligning with the smoothing factors used to smooth the option's payoff and other functions involved in the MC simulation. This enables the computation of samples using second-order derivatives in the PINNs approach.</footnote>
### 4.2.2 Sensitivities w.r.t. durations $\tau_i$

Similarly to the approach taken when using prices as labels, the Backward Kolmogorov equation can be used. This technique is based on the connection between the partial derivative of a certain price function $\hat{\pi}(x_{T,T}^\xi)$ w.r.t. a given time duration $\tau_i$ and the infinitesimal generator $\mathcal{L}_i$ of the stochastic process $(X_{t,T}, \nu_t)$ over the interval $[t_i, t_{i+1}]$. Mathematically, this relationship is expressed as

$$ \partial_{\tau_i} \hat{\pi}(x_{T,T}^\xi) = \mathcal{L}_i \hat{\pi}(x_{T,T}^\xi), \qquad (50) $$

To facilitate the computation of the second partial derivative $\partial_x^2 \hat{\pi}(x_{T,T}^\xi)$, it is essential to work with a smooth approximation of $\pi$. If $\pi$ is not smoothed, the second partial derivative becomes a Dirac function, making MC estimation intractable. In the case of European Call and Put options, the Softplus function can serve as an effective smooth approximation, enabling the computation of the required derivatives.

# 5 Numerical Implementation and Results

## 5.1 Generating Data

This section outlines the methodology for creating datasets to train, validate, and test our ANNs. We apply machine learning to approximate a time-dependent Heston model featuring numerous parameters. This method is advantageous for managing large portfolios or complex financial metrics such as Value-at-Risk and Risk Weighted Assets, critical for market and counterparty risk management.

Given practical constraints in generating extensive offline labels, even with advanced technologies, we deliberately selected datasets containing 2048 and 8192 data points. Such scales are standard practice for evaluating complex metrics and examining their sensitivities to market changes. Selecting this scale allows us to effectively showcase the advantages of both the DDL and the combined DDL-PINN approaches, especially when dealing with limited and noisy payoff realizations rather than numerically computed prices.

### 5.1.1 Attributes generation via LHS

We sample the attributes composed of the Heston model's piecewise parameters, denoted by $\{\bar{\nu}_i, \kappa_i, \gamma_i, \rho_i, \tau_i\}_{i \in \langle 0, n \rangle}$. In addition, we consider the option and market parameters $X_0$, $\nu_0$, and $T = \sum_{i=0}^n \tau_i$, all of which are defined within their respective ranges.

Assuming that the features and labels $\{(\xi_m, \zeta_m)\}_{m \in \langle 1, M \rangle}$ are i.i.d., we use MC techniques for data generation. To avoid sparsity issues, especially with extensive input ranges, such as time-to-maturities from 1 day to 5 years, we consider Latin hypercube sampling (LHS). It is a stratified method that generates random samples over equidistant intervals in a $N$-dimensional space, ensuring that each sample occupies a unique position on $(N-1)$-dimensional hyperplanes. This ensures a non-collapsing, uniformly spread dataset, adaptable to data density and position. Thus, we use LHS to generate the feature domain, as summarized in Table 1.

<table><caption>Table 1: Features range</caption><thead><tr><td>Parameter</td><td>Range</td></tr></thead><tbody><tr><td>Forward log-moneyness: X<sub>0</sub></td><td>[ln 0.3, ln 2]</td></tr><tr><td>Initial variance: ν<sub>0</sub></td><td>[0.01,0.04]</td></tr><tr><td>Option's maturity: T</td><td>[1 day, 5 years]</td></tr><tr><td>Mean reversion's speeds: κ<sub>i</sub></td><td>[0.1, 5]</td></tr><tr><td>Variance's long-term averages: ν&#773;<sub>i</sub></td><td>[0.1,0.3]</td></tr><tr><td>Variance's volatilities: γ<sub>i</sub></td><td>[0.01,1]</td></tr><tr><td>Asset-Variance correlations: ρ<sub>i</sub></td><td>[-0.99, 0]</td></tr></tbody></table>

### 5.1.2 Labels generation

Given a specific set of attributes, including the Heston model's piecewise parameters, option features, and market data, we utilize the pricing formulas presented in Sections 2.1.1 and Section 3 to compute the option price label and its associated differentials. Using the same attributes, we generate an MC path as detailed in Section 4, resulting in both the realization of a payoff and its differentials.
### 5.1.2.1 Constant parameters

In the case of a Heston model with constant parameters, there are seven attributes: the five Heston parameters $κ, \bar{ν}, γ, ρ, τ$, along with log-moneyness and the initial stochastic variance. We begin by generating a training dataset consisting of 2048 samples. Labels are computed as both option prices and payoff realizations using the QE scheme with daily time steps. Calculating prices and their sensitivities to input parameters required approximately 6 hours on a single A100 GPU with 40GB of RAM, while generating 2048 noisy payoff realizations took less than 15 minutes. This results in two distinct training datasets: 1_Step_2048_Payoff, containing 2048 price labels along with their sensitivities to model attributes, and 1_Step_2048_Payoff, containing 2048 samples labeled with payoff realizations from the QE scheme and their corresponding sensitivities.

Using payoff realizations as labels introduces higher variance during training compared to prices, which can result in more erratic optimizer updates. However, generating MC payoff realizations, including their sensitivities, is significantly more efficient than computing prices. This computational advantage enables the use of a larger dataset of 8192 samples, denoted 1_Step_8192_Payoff, which helps reduce overfitting and improve generalization.

### 5.1.2.2 Yearly piecewise parameters

In the yearly piecewise parameter scenario (comprising 5 steps), the model takes in 27 attributes. To detail, each of the 5 steps is governed by five Heston parameters: $κ, \bar{ν}, γ, ρ, τ$, summing up to $5 \times 5$. There are also two features derived from market data and option specifics, namely the log-moneyness and the initial instantaneous variance. Thus, the aggregate parameter count is $5 \times 5 + 2$.

As described in Section 3, computing prices for time-dependent Heston parameters is noticeably more time consuming compared to the constant parameter scenario. This emphasizes the advantage of leveraging payoff realizations. As in the previous case, we generate 3 training datasets: 5_Steps_2048_Payoff, 5_Steps_2048_Payoff, and 5_Steps_8192_Payoff.

Additionally, we constructed a validation dataset using prices; it is noteworthy that there was no need to consider label sensitivities with respect to attributes for this validation set, which accelerated its generation. Consequently, its creation was faster than that of 5_Steps_2048_Payoff. In total, the data generation process spanned nearly 24 hours. It is remarkable that out of this time frame, generating the MC labels— for 5_Steps_2048_Payoff and 5_Steps_8192_Payoff—consumed less than 30 minutes.

### 5.1.2.3 Monthly piecewise parameters

In the case of a Heston model with monthly piecewise parameters, relying on prices as labels becomes computationally impractical. For example, generating just 2048 price labels could require several weeks of computation. As a result, we construct two training datasets, 60_Steps_2048_Payoff and 60_Steps_8192_Payoff, containing 2048 and 8192 data points, respectively. The labels for these datasets are generated using MC simulations combined with a daily time-step QE scheme. In this setting, the number of input attributes increases to 302, corresponding to $60 \times 5 + 2$: five Heston parameters per monthly interval, along with log-moneyness and the initial instantaneous variance.

## 5.2 Normalizing Data

In DL, data normalization is an essential preprocessing step that eliminates the disparities in the variables scales. This technique helps prevent larger-valued variables from unduly influencing the model without losing or distorting the differences between values. A common method for normalizing the data involves computing the mean $\bar{\xi}_j$, and the standard deviation $s_j$, for each variable $\xi_j$. Each observation $\xi_{ij}$ is then normalized using the expression:

$$ \xi_{ij} = (\xi_{ij} - \bar{\xi}_j) / s_j $$
This transformation results in a distribution for each variable with a mean of zero and a variance of one, a process often referred to as reduced centered normalization.

For the training data, normalization is applied to both inputs and labels. In the DDL approach, their derivatives also undergo a specific normalization process. This includes applying a reduced-centred transformation to the inputs and labels. The differentials of the labels with respect to the features are normalized by adjusting them by multiplying by the standard deviations of the inputs and division by the standard deviation of the labels, as described by Huge and Savine [2020]. This particular adjustment ensures that the differentials are consistent with the transformed space, which facilitates the learning process.

The validation set is normalized using the same mean and standard deviation derived from the training dataset. To compare predictions with exact prices, the predictions are scaled back to the original units by reversing the normalization transformation.

## 5.3 Performance Assessment

This section rigorously evaluates the methodologies and algorithms described previously, focusing on their key features and applications. Particular attention is paid to the three Heston parameter time-dependency scenarios introduced in Section 5.1.2. Our analysis indicates that DDL induces a regularization effect, largely attributable to DDL loss. We also examine the impact of incorporating the PINN loss as an additional regularization technique. We view PINNs as part of a broader regularization strategy, where the combined use of DDL and PINN losses synergistically enhances accuracy.

Importantly, as noted in Section 2.3.2.1, the PINNs approach can be interpreted as a method to solve inverse problems. This is achieved using only the PDE residual loss, $L_N$, and the boundary loss, $L_B$, without relying on the labeled data. In this setting, only unlabeled collocation points are utilized. We consider this to be a specific weight configuration within the broader range of weight combinations explored in the combined DDL-PINN strategy.

### 5.3.1 DDL regularization

In this section, we analyze the impact of DDL regularization on the three time-dependent Heston parameter scenarios described in Section 5.1.2. Using prices as labels is feasible only in cases where the Heston parameters were either constant or yearly piecewise. When valuing a vanilla European option, it is important to note that the computational speed of the Lipton formula is slower than that of closed-form solutions like the Black-Scholes model. However, the constant-parameter case remains more efficient than the Heston model with piecewise parameters. In particular, computations for the yearly piecewise parameters are significantly more resource-intensive.

To compute the price (label) sensitivities to Heston parameters (attributes), we compare a direct computation via Automatic Adjoint Differentiation (AAD)<fnref n="6" /> to the method detailed in Section 3. The latter approach reduced the computation time of the differentials of the labels by 25%. Despite this efficiency, producing 2048 labels and their differentials for the monthly piecewise Heston parameters remains impractical, as the process can take weeks. Consequently, our DDL regularization analysis using prices as labels is limited to the constant-parameter and yearly piecewise parameter scenarios.

When using payoff realizations as labels, we construct all three Heston parameter time-dependency scenarios and train models on datasets of size 2048 and 8192.

Furthermore, recall that the DDL loss is a weighted sum of MSEs measuring the discrepancy between the ANN model's sensitivity to specific features (i.e., Heston model parameters) and the sensitivity of the label to the same feature. For a time-dependent Heston model, the cumulative DDL loss is given by $L_{DDL} = \sum_i L_{i,DDL}$. Each component $L_{i,DDL}$ is defined by:

$$L_{i,DDL} = \mathbf{E} \left[ \left( \partial_i \mathcal{M}(\xi|\omega) - \partial_i \pi(\xi) \right)^2 \right]. \qquad (51)$$

<footnote n="6">⁶We leveraged the `autograd` function from PyTorch</footnote>
Since DDL loss acts as a regularizer, a weight $\omega_{DDL}$ is incorporated into the total loss to help prevent overfitting. Huge and Savine [2020] suggests that a commonly chosen value for $\omega_{DDL}$ is $1/N$, where N denotes the number of features in the DDL loss. However, our experiments indicate that this heuristic does not hold universally. In particular, the initial magnitudes of the losses $L_{STD}$ and $L_{DDL}$ significantly influence the choice of penalization weight, especially when dealing with noisy labels (e.g., payoff realizations). We found that when the training set is small (as in our case) and the number of attributes increases, assigning a higher weight to $L_{DDL}$ in the overall training loss leads to better results.

### 5.3.1.1 Prices as Labels

In this case, the training uses 2048 data points, and the ANN validation uses 512 data points. We utilize the following training hyperparameters:

<table><caption>**Table 2:** Training hyperparameters when using prices as labels</caption><thead><tr><td>Hyperparameter</td><td>Value</td></tr></thead><tbody><tr><td>Epochs</td><td>100</td></tr><tr><td>Batch size</td><td>64</td></tr><tr><td>Learning rate x-axis</td><td>[0.0, 0.2, 0.6, 0.9, 1.0]</td></tr><tr><td>Learning rate y-axis</td><td>[1e-8, 1e-1, 1e-2, 1e-5, 1e-6]</td></tr></tbody></table>

The learning rate schedule provided in the table above follows a piecewise linear function. The values on the x-axis represent percentages of the total number of epochs, and the values on the y-axis correspond to the learning rate values at those points. The following table summarizes the results, presenting weight parameters across different sets along with their corresponding MSE values:

<table><caption>**Table 3:** Summary of $\omega_{DDL}$ and corresponding MSE values for different sets and parameters' time-dependency.</caption><thead><tr><td>Set</td><td>ω<sub>DDL</sub></td><td>Constant</td><td>Yearly</td></tr></thead><tbody><tr><td>Set 1</td><td>0.00</td><td>10.77 × 10<sup>-6</sup></td><td>50.07 × 10<sup>-6</sup></td></tr><tr><td>Set 2</td><td>0.04</td><td>5.30 × 10<sup>-6</sup></td><td>18.43 × 10<sup>-6</sup></td></tr><tr><td>Set 3</td><td>0.10</td><td>5.46 × 10<sup>-6</sup></td><td>19.37 × 10<sup>-6</sup></td></tr><tr><td>Set 4</td><td>0.20</td><td>7.17 × 10<sup>-6</sup></td><td>22.57 × 10<sup>-6</sup></td></tr><tr><td>Set 5</td><td>0.50</td><td>8.32 × 10<sup>-6</sup></td><td>25.90 × 10<sup>-6</sup></td></tr><tr><td>Set 6</td><td>1.00</td><td>14.60 × 10<sup>-6</sup></td><td>29.70 × 10<sup>-6</sup></td></tr></tbody></table>

The table above shows that when using the datasets 1_Step_2048.Price and 5_Steps_2048.Price, which include both prices and their respective differentials with respect to input features, the ANN achieves remarkable accuracy without requiring DDL penalization. In the Heston model with constant parameters, there are 7 attributes, whereas in the yearly piecewise parameter case (with 5 steps), there are 27 attributes. This explains the observed difference in MSE magnitudes between the constant and yearly cases.

In particular, when the labels are numerically computed prices, they are less prone to noise. Consequently, the task reduces to an interpolation problem for the ANN. In this case, the primary role of DDL is to mitigate the numerical imprecisions that may arise during price computation, as detailed in Section 3.

DDL improves prediction accuracy even in this case. For instance, in the constant parameter scenario, the lowest validation MSE is achieved with the weightings of Set 2, and a similar trend is observed for yearly parameters. DDL penalization improves alignment with labels by enforcing penalties on incorrect differentials. Note that this regularization approach does not introduce bias. Even when trained solely on label differentials, the model converges to an approximation that, while potentially offset by a constant, remains appropriately aligned. The standard loss term minimizes this constant offset to zero.
The figure below shows the precision of the price predictions for selected sets in the preceding table:

<figure><img src="image_6.png" /><figcaption>a: Set 1: Constant</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>b: Set 2: Constant</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>c: Set 1: Yearly</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>d: Set 2: Yearly</figcaption></figure>

<caption>**Figure 5:** Predicted vs ground truth option prices for different weight sets and parameters time-dependency cases.</caption>

The following figure shows the progression of the training loss components over epochs. For both constant and yearly Heston parameters, Set 1 shows significantly higher DDL, inner, and border PINN losses compared to other sets, as it only optimizes the standard loss. In contrast, Set 2 achieves the lowest training loss. Interestingly, using only DDL penalization performs well on inner- and border-PINN losses, even though these are not included in the training loss. This can be attributed to the fact that they involve several first-order derivatives.

<figure><img src="image_10.png" /><figcaption>**Figure 6:** Training Losses for: DDL - Price - 2048</figcaption></figure>

Next, we evaluate the ANN's performance on the test set. The two figures below present the MSE and MAE, which quantify the deviation between the ground-truth prices and the ANN's predictions.
<figure><img src="image_11.png" /><figcaption>Figure 7: Evaluation metrics: DDL - Price - 2048</figcaption></figure>

The trained ANNs demonstrate broadly comparable performance in predicting option prices; however, Set 2 consistently emerges as the optimal configuration. As shown in the previous table, it achieves the lowest training and evaluation MSE, indicating strong generalization without overfitting. Notably, Set 2 also maintains minimal DDL and PINN losses, In practical applications, it is crucial to select an ANN that accurately captures sensitivities and adheres to the PDE that governs the underlying asset dynamics. Such an ANN inherently possesses strong dynamic hedging capabilities.

### 5.3.1.2 Payoff Realizations as Labels

Using payoff realizations instead of prices for training enables faster data generation. However, this approach introduces specific difficulties, notably when the training process experiences greater variance compared to the use of prices. This increased variance can produce noisier optimizer updates due to significant fluctuations in the gradient estimators.

To mitigate high variance during the optimization of the training loss, a sufficiently large dataset and larger batch sizes can be used. If data is limited or generating additional data is time-consuming, traditional methods such as regularization of L1 (Lasso) or L2 (Ridge) can be beneficial. The subsequent section will illustrate the effectiveness of DDL penalization, especially when working with limited training data (e.g., 2048 samples) and a large number of features (e.g., monthly piecewise Heston model).

In addition to the primary training datasets that contain 2048 and 8192 data points, we employ two validation sets for the ANN with 512 and 1024 data points, respectively. Despite using payoff realizations as labels, model accuracy is evaluated based on the MSE between the ANN's predicted prices and benchmark prices, computed as described in Section 3. The following training hyperparameters are used:

<table><caption>Table 4: Training hyperparameters when payoff realizations are labels</caption><thead><tr><td>Hyperparameter</td><td>Value</td></tr></thead><tbody><tr><td>Epochs</td><td>300</td></tr><tr><td>Batch size</td><td>256</td></tr><tr><td>Learning rate x-axis</td><td>[0.0, 0.2, 0.6, 0.9, 1.0]</td></tr><tr><td>Learning rate y-axis</td><td>[1e-3, 1e-4, 1e-5, 1e-5, 1e-6]</td></tr><tr><td>number of hidden layers</td><td>3</td></tr><tr><td>number of neurons per hidden layer</td><td>64</td></tr></tbody></table>
The results, showcasing various weight parameter sets and their corresponding MSE values, are presented in the following table:

<table><caption>**Table 5:** Summary of ω<sub>DDL</sub> and corresponding MSE values for different sets and parameters' time-dependency.</caption><thead><tr><th>Set</th><th>ω<sub>DDL</sub></th><th>Constant</th><th>Yearly</th><th>Monthly</th></tr></thead><tbody><tr><td colspan="5" style="text-align: center; font-weight: bold;">2048 payoff realizations</td></tr><tr><td>Set 1</td><td>0.00</td><td>18.30 × 10<sup>-4</sup></td><td>21.33 × 10<sup>-4</sup></td><td>275.25 × 10<sup>-4</sup></td></tr><tr><td>Set 2</td><td>0.04</td><td>18.29 × 10<sup>-4</sup></td><td>15.03 × 10<sup>-4</sup></td><td>86.17 × 10<sup>-4</sup></td></tr><tr><td>Set 3</td><td>0.10</td><td>18.83 × 10<sup>-4</sup></td><td>12.34 × 10<sup>-4</sup></td><td>53.33 × 10<sup>-4</sup></td></tr><tr><td>Set 4</td><td>0.20</td><td>19.94 × 10<sup>-4</sup></td><td>10.53 × 10<sup>-4</sup></td><td>32.51 × 10<sup>-4</sup></td></tr><tr><td>Set 5</td><td>0.50</td><td>21.43 × 10<sup>-4</sup></td><td>8.90 × 10<sup>-4</sup></td><td>14.78 × 10<sup>-4</sup></td></tr><tr><td>Set 6</td><td>1.00</td><td>22.87 × 10<sup>-4</sup></td><td>8.17 × 10<sup>-4</sup></td><td>7.37 × 10<sup>-4</sup></td></tr><tr><td colspan="5" style="text-align: center; font-weight: bold;">8192 payoff realizations</td></tr><tr><td>Set 1</td><td>0.00</td><td>7.66 × 10<sup>-4</sup></td><td>13.71 × 10<sup>-4</sup></td><td>542.77 × 10<sup>-4</sup></td></tr><tr><td>Set 2</td><td>0.04</td><td>6.87 × 10<sup>-4</sup></td><td>8.50 × 10<sup>-4</sup></td><td>108.88 × 10<sup>-4</sup></td></tr><tr><td>Set 3</td><td>0.10</td><td>6.66 × 10<sup>-4</sup></td><td>6.74 × 10<sup>-4</sup></td><td>58.12 × 10<sup>-4</sup></td></tr><tr><td>Set 4</td><td>0.20</td><td>6.53 × 10<sup>-4</sup></td><td>5.66 × 10<sup>-4</sup></td><td>29.98 × 10<sup>-4</sup></td></tr><tr><td>Set 5</td><td>0.50</td><td>6.48 × 10<sup>-4</sup></td><td>4.63 × 10<sup>-4</sup></td><td>10.13 × 10<sup>-4</sup></td></tr><tr><td>Set 6</td><td>1.00</td><td>6.66 × 10<sup>-4</sup></td><td>4.07 × 10<sup>-4</sup></td><td>5.04 × 10<sup>-4</sup></td></tr></tbody></table>

The table above provides several key insights:

* The effectiveness of DDL penalization becomes more pronounced as the number of attributes increases, as seen in the monthly piecewise case, which results in 302 attributes.

* Looking at Sets 4 to 6, we observe that increasing the size of the training dataset reduces the MSE when applying DDL regularization with weights between 0.2 and 1.0.

* In the monthly case (e.g. Set 1 and even Sets 2–3), increasing the training size (2048 to 8192) does not reduce the MSE. The high-variance loss from noisy payoff labels creates an irregular optimization landscape, making the network prone to overfitting noise, particularly when regularization is weak.

* The regularization effect of DDL intensifies as the dimension of the feature space dimension increases. For instance, for a training dataset of size 2048, the MSE decreases from 275.25 × 10⁻⁴ with standard loss to 7.37 × 10⁻⁴. A similar improvement is observed with 8192 data points.

* When the training dataset size increases (e.g., 8192) while the feature dimension remains relatively small (e.g., 7 in the constant Heston parameters case), selecting the appropriate DDL penalization weight is crucial. This is similar to scenarios where prices are used as labels.

Next, we present a figure that illustrates the evolution of the accuracy of the ANN model for the monthly piecewise parameter case using 2048 training data points<fnref n="7" />.

<figure><img src="image_12.png" /><figcaption>**Figure 8:** ANN predicted prices compared to ground truth prices for different weight sets</figcaption></figure>

<footnote n="7">⁷For conciseness, we have chosen not to include additional visualizations that lead to the same conclusions.</footnote>
The following figure illustrates the evolution of various training loss components over epochs for the monthly piecewise Heston parameter case. In Set 1, the DDL loss, along with several PINN losses, stabilizes at a significantly higher value than in other sets. This behavior is expected since Set 1 optimizes only the standard loss component. Incorporating DDL regularization leads to a notable reduction in validation loss and PINN-related training loss components, although they are not explicitly activated in this context. This highlights the beneficial impact of DDL regularization, yielding significant improvements despite the limited number of training samples.

<figure><img src="image_13.png" /><figcaption>Figure 9: Training Losses for: DDL - Price - 2048 - Monthly piecewise case</figcaption></figure>

Next, we present the evaluation loss for the monthly piecewise parameter case with 2048 training data points. The following figures illustrate the MSE and MAE, which quantify the deviation between the ground-truth prices and the ANN's predictions.

<figure><img src="image_14.png" /><figcaption>Figure 10: Evaluation metrics: DDL - Price - 2048 - Monthly piecewise case</figcaption></figure>

Increasing the weight of DDL loss in the training objective improves the prediction accuracy of the ANN. This improvement is attributed to the regularization effect of the DDL loss, which mitigates the overfitting observed in Set 1.
### 5.3.2 Joint DDL & PINN Regularization

This section examines the impact of combining the DDL and PINN regularization methods. This joint approach strengthens regularization, particularly when training on small labeled datasets (2048 or 8192) with a substantial number of Heston parameters. We consider the monthly piecewise configuration, which involves 302 parameters, and a new more granular weekly configuration, where the parameter count rises to 1262 (= 252 × 5 + 2).

The labeled data is used to train $L_{STD}$ and $L_{DDL}$, while the collocation points are employed for $L_{PINN}$. We subsequently demonstrate that integrating PINN regularization with DDL enhances model accuracy. The intuition is that DDL regularization, with a well-chosen $\omega_{DDL}$, already guides the ANN parameters towards a favorable region of the parameter space. However, some overfitting may still occur, particularly when the training set is limited to 2048 labeled samples. The PINN term provides an additional layer of regularization, enforcing that the trained ANN also satisfies the following constraints:

* The time-dependent Heston PDE: This is enforced by minimizing the inner (resp. boundary) PDE loss, which corresponds to the PDE residual at collocation points located within (resp. near the boundaries of) the feature domain. The weights $\omega_{in.pde}$ and $\omega_{out.pde}$ control their respective contributions to the overall loss.

* Boundary conditions: Near the option's maturity, the trained ANN must converge to the intrinsic value of the option. We enforce a Neumann-type boundary condition that ensures that the ANN matches this value at $\tau_i = 0$, while also maintaining consistent sensitivities with respect to the input features. The contributions of these boundary losses are weighted by $\omega_{exp}$ and $\omega_{exp.s}^8$

We consider two strategies to combine the losses of DDL and PINN. In the first, $L_{DDL}$ and $L_{PINN}$ coexist throughout the entire training process. In the second, we dynamically modulate the influence of DDL penalization during the training period by adjusting the weights assigned to $L_{STD}$ and $L_{DDL}$ as a function of the training epoch. To do this, we introduce two parameters: $s_{PRM} \in \mathbb{R}^+$ and $i_{PRM} \in \langle 1, I \rangle$, where $I$ is the total number of training epochs. These parameters define a multiplier $\mu_{PRM}$ for $L_{STD}$ and $L_{DDL}$, which decays exponentially with training epoch:

$$ \mu_{PRM}(i) := \exp(-s_{PRM} (i - i_{PRM})^{+} / I), \quad i \in \langle 1, I \rangle $$

The overall training loss at epoch $i \in \langle 1, I \rangle$ is then expressed as:

$$ L_{TOT}(i) = \mu_{PRM}(i) \times \underbrace{(\omega_{STD} L_{STD} + \omega_{DDL} L_{DDL})}_{\text{Trained with labeled data}} + \underbrace{(\omega_{in.p} L_{in.p} + \omega_{bd.p} L_{bd.p} + \omega_{exp} L_{exp} + \omega_{exp.s} L_{exp.s})}_{\text{Trained with unlabeled data}} $$

#### 5.3.2.1 Monthly piecewise parameters

The training procedure uses 2048 unlabeled data points, consisting of 1024 collocation points within the feature domain and an additional 1024 points near the domain border. For validation, 1024 data points are used. We use the following training hyperparameters:

<table><caption>**Table 6:** Hyperparameters using DDL-PINN penalization for 60 steps</caption><thead><tr><th>Labeled data size</th><th>2048</th><th>8192</th></tr></thead><tbody><tr><td>Epochs</td><td>200</td><td>150</td></tr><tr><td>Batch size</td><td>512</td><td>512</td></tr><tr><td>Learning rate x-axis</td><td>[0.0, 0.3, 1.0]</td><td>[0.0, 0.4, 1.0]</td></tr><tr><td>Learning rate y-axis</td><td>[8e-3, 1e-3, 2e-4]</td><td>[1e-2, 1e-3, 1e-4]</td></tr><tr><td>number of hidden layers</td><td>3</td><td>3</td></tr><tr><td>number of neurons per hidden layer</td><td>49</td><td>64</td></tr></tbody></table>

The following tables present the results for various sets of hyperparameters, along with their corresponding minimal and final MSE values. Note that the hyperparameter sets yielding the lowest MSE values were refined through hyperparameter tuning.

<footnote n="8">⁸This work only considers boundary conditions related to the underlying asset S. Note that additional boundary conditions for stochastic volatility $\nu$ may be included during training, depending on its domain.</footnote>
<table><caption>**Table 7:** Performance of various DDL and PINN configurations in the monthly case.</caption><thead><tr><th>Set</th><th>ω<sub>STD</sub></th><th>ω<sub>DDL</sub></th><th>ω<sub>in.PDE</sub></th><th>ω<sub>bd.p</sub></th><th>ω<sub>exp</sub></th><th>ω<sub>exp.s</sub></th><th>s<sub>PRM</sub></th><th>i<sub>PRM</sub></th><th>MIN MSE</th><th>LAST MSE</th></tr></thead><tbody><tr><td colspan="11" style="text-align: center; font-weight: bold;">2048 payoff realizations</td></tr><tr><td>Set 1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>-</td><td>-</td><td>9.60 × 10<sup>-4</sup></td><td>17.42 × 10<sup>-4</sup></td></tr><tr><td>Set 2</td><td>1</td><td>1</td><td>0.2</td><td>0.1</td><td>0.1</td><td>0.1</td><td>-</td><td>-</td><td>8.98 × 10<sup>-4</sup></td><td>14.08 × 10<sup>-4</sup></td></tr><tr><td>Set 3</td><td>1</td><td>50</td><td>0</td><td>0</td><td>0</td><td>0</td><td>-</td><td>-</td><td>2.92 × 10<sup>-4</sup></td><td>3.17 × 10<sup>-4</sup></td></tr><tr><td>Set 4</td><td>1</td><td>50</td><td>0.2</td><td>0.1</td><td>0.1</td><td>0.1</td><td>-</td><td>-</td><td>2.69 × 10<sup>-4</sup></td><td>3.07 × 10<sup>-4</sup></td></tr><tr><td>Set 5</td><td>1</td><td>50</td><td>0.2</td><td>0.1</td><td>0.1</td><td>0.1</td><td>3</td><td>60</td><td>2.68 × 10<sup>-4</sup></td><td>2.98 × 10<sup>-4</sup></td></tr><tr><td>Set 6</td><td>1</td><td>50</td><td>0.2</td><td>0.1</td><td>0.1</td><td>0.1</td><td>3</td><td>100</td><td>2.69 × 10<sup>-4</sup></td><td>2.90 × 10<sup>-4</sup></td></tr><tr><td colspan="11" style="text-align: center; font-weight: bold;">8192 payoff realizations</td></tr><tr><td>Set 1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>-</td><td>-</td><td>5.84 × 10<sup>-4</sup></td><td>12.94 × 10<sup>-4</sup></td></tr><tr><td>Set 2</td><td>1</td><td>20</td><td>0</td><td>0</td><td>0</td><td>0</td><td>-</td><td>-</td><td>3.15 × 10<sup>-4</sup></td><td>3.32 × 10<sup>-4</sup></td></tr><tr><td>Set 3</td><td>1</td><td>12.5</td><td>0</td><td>0</td><td>0</td><td>0</td><td>-</td><td>-</td><td>3.11 × 10<sup>-4</sup></td><td>3.26 × 10<sup>-4</sup></td></tr><tr><td>Set 4</td><td>1</td><td>12.5</td><td>1.0</td><td>1.0</td><td>0.5</td><td>0.2</td><td>-</td><td>-</td><td>1.56 × 10<sup>-4</sup></td><td>1.68 × 10<sup>-4</sup></td></tr><tr><td>Set 5</td><td>1</td><td>12.5</td><td>1.0</td><td>1.0</td><td>0.5</td><td>0.2</td><td>3</td><td>37</td><td>1.48 × 10<sup>-4</sup></td><td>1.73 × 10<sup>-4</sup></td></tr><tr><td>Set 6</td><td>1</td><td>12.5</td><td>1.0</td><td>1.0</td><td>0.5</td><td>0.2</td><td>3</td><td>112</td><td>1.56 × 10<sup>-4</sup></td><td>1.64 × 10<sup>-4</sup></td></tr></tbody></table>

When using 2048 or 8192 labeled training samples, we observe that Set 5 achieves the lowest minimal MSE (optimal when using early stopping), while Set 6 yields the lowest final MSE. For Sets 4 to 6, the combination of DDL and PINN regularization leads to better ANN performance compared to Set 3, which uses only DDL regularization. Furthermore, comparing the performance of Sets 5 and 6 to Set 4 suggests that applying DDL regularization that gradually vanishes after a certain training epoch helps reduce overfitting (as Set 4's last MSE is higher than Set 6's). The following figures illustrate the evolution of the valuation MSE with respect to the number of epochs:

<figure><img src="image_15.png" /><figcaption>**Figure 11:** Evaluation metrics for DDL and PINN combinations in the monthly piecewise case.</figcaption></figure>

### 5.3.2.2 Weekly piecewise parameters

The training procedure utilizes 256 unlabeled data points, consisting of 128 collocation points within the feature domain and an additional 128 points near the domain border. For validation, 1024 data points are used. We use the following training hyperparameters:
<table><caption>**Table 8:** Hyperparameters using DDL-PINN penalization for 252 steps</caption><thead><tr><th>Labeled data size</th><th>2048</th><th>8192</th></tr></thead><tbody><tr><td>Epochs</td><td>150</td><td>150</td></tr><tr><td>Batch size</td><td>512</td><td>512</td></tr><tr><td>Learning rate x-axis</td><td>[0.0, 0.2, 1.0]</td><td>[0.0, 0.2, 1.0]</td></tr><tr><td>Learning rate y-axis</td><td>[1e-2, 1e-3, 2e-4]</td><td>[1e-2, 1e-3, 1e-4]</td></tr><tr><td>number of hidden layers</td><td>3</td><td>3</td></tr><tr><td>number of neurons per hidden layer</td><td>64</td><td>64</td></tr></tbody></table>

Below, we provide the final and minimal MSE values for various parameters' sets:

<table><caption>**Table 9:** Performance of various DDL and PINN configurations in the weekly case.</caption><thead><tr><th>Set</th><th>ω<sub>STD</sub></th><th>ω<sub>DDL</sub></th><th>ω<sub>in_PDE</sub></th><th>ω<sub>bd.p</sub></th><th>ω<sub>exp</sub></th><th>ω<sub>exp.s</sub></th><th>s<sub>PRM</sub></th><th>i<sub>PRM</sub></th><th>MIN MSE</th><th>LAST MSE</th></tr></thead><tbody><tr><td colspan="11" style="text-align: center; font-weight: bold;">2048 payoff realizations</td></tr><tr><td>Set 1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>-</td><td>-</td><td>18.97 × 10<sup>-4</sup></td><td>18.97 × 10<sup>-4</sup></td></tr><tr><td>Set 2</td><td>1</td><td>50</td><td>0</td><td>0</td><td>0</td><td>0</td><td>-</td><td>-</td><td>2.44 × 10<sup>-4</sup></td><td>2.44 × 10<sup>-4</sup></td></tr><tr><td>Set 3</td><td>1</td><td>1500</td><td>0</td><td>0</td><td>0</td><td>0</td><td>-</td><td>-</td><td>3.17 × 10<sup>-4</sup></td><td>3.17 × 10<sup>-4</sup></td></tr><tr><td>Set 4</td><td>1</td><td>1500</td><td>1.0</td><td>1.0</td><td>0.5</td><td>0</td><td>-</td><td>-</td><td>1.76 × 10<sup>-4</sup></td><td>1.78 × 10<sup>-4</sup></td></tr><tr><td>Set 5</td><td>1</td><td>1500</td><td>1.0</td><td>1.0</td><td>0.5</td><td>0.2</td><td>-</td><td>-</td><td>1.76 × 10<sup>-4</sup></td><td>1.79 × 10<sup>-4</sup></td></tr><tr><td>Set 6</td><td>1</td><td>1500</td><td>1.0</td><td>1.0</td><td>0.5</td><td>0.2</td><td>3</td><td>112</td><td>1.74 × 10<sup>-4</sup></td><td>1.75 × 10<sup>-4</sup></td></tr><tr><td colspan="11" style="text-align: center; font-weight: bold;">8192 payoff realizations</td></tr><tr><td>Set 1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>-</td><td>-</td><td>7.45 × 10<sup>-4</sup></td><td>8.51 × 10<sup>-4</sup></td></tr><tr><td>Set 2</td><td>1</td><td>1000</td><td>0</td><td>0</td><td>0</td><td>0</td><td>-</td><td>-</td><td>2.25 × 10<sup>-4</sup></td><td>2.29 × 10<sup>-4</sup></td></tr><tr><td>Set 3</td><td>1</td><td>500</td><td>0</td><td>0</td><td>0</td><td>0</td><td>-</td><td>-</td><td>2.18 × 10<sup>-4</sup></td><td>2.23 × 10<sup>-4</sup></td></tr><tr><td>Set 4</td><td>1</td><td>500</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.2</td><td>3</td><td>112</td><td>1.68 × 10<sup>-4</sup></td><td>1.82 × 10<sup>-4</sup></td></tr><tr><td>Set 5</td><td>1</td><td>500</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.2</td><td>-</td><td>-</td><td>1.63 × 10<sup>-4</sup></td><td>1.63 × 10<sup>-4</sup></td></tr><tr><td>Set 6</td><td>1</td><td>1000</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.2</td><td>3</td><td>112</td><td>1.62 × 10<sup>-4</sup></td><td>1.66 × 10<sup>-4</sup></td></tr></tbody></table>

With 2048 or 8192 labeled samples, Set 6 achieves the lowest minimal MSE under early stopping. For final MSE, Set 6 and Set 5 perform best with 2048 and 8192 samples, respectively. For Sets 4 to 6, combining DDL and PINN regularization improves ANN performance over Set 3, which uses only DDL. Furthermore, the comparison between Sets 5 and 6 suggests that gradually reducing DDL regularization after a certain epoch enhances predictive accuracy under early stopping, as Set 6 achieves a lower minimum MSE than Set 5. The following figures show the evolution of valuation MSE over training epochs:

<figure><img src="image_16.png" /><figcaption>**Figure 12:** Evaluation metrics for DDL and PINN combinations in the weekly piecewise case.</figcaption></figure>
# 6 Conclusion

The primary contribution of this work is the combination of two advanced deep learning techniques to reduce the number of required labels. using payoff realizations as labels, we accelerate the offline generation of datasets. This approach enables us to train an ANN that approximates option prices within a Heston model featuring time-dependent parameters.

To implement PINNs and DDLs, we developed algorithms to compute label sensitivities to inputs, particularly when labels are interpreted as prices under piecewise Heston parameters. Additionally, we utilized complex step differentiation to evaluate the sensitivity of payoff realizations based on the QE scheme.

We showed that the DDL is robust to noisy labels and can be used in situations with limited datasets and high-dimensional feature spaces. DDL's label differentials, which are similar to data augmentation techniques in computer vision, can be used to augment the training dataset and help the ANN learn the intrinsic structure of the option price as a function of the features (time-dependent Heston model parameters, option features, and initial market data).

Additionally, we derived a formulation to compute the Hessian of an FFNN with respect to its input. This result enables the seamless integration of PINNs, further improving training when combined with DDL.

Finally, we highlight the versatility of the proposed deep learning framework. Beyond delivering accurate option prices, the trained ANN also provides sensitivities with respect to its input features, namely, the Heston model parameters, log-moneyness, and initial instantaneous variance. In addition, this ANN can be easily used for complex sensitivity analyzes, which are crucial for effective derivative portfolios risk management.
## References

A. Alfonsi. High order discretization schemes for the cir process: Application to affine term structure and heston models. *Mathematics of computation*, 79(269):209–237, 2010. pages 6

L. Andersen. Efficient simulation of the heston stochastic volatility model, 2006. URL https://ssrn.com/abstract=946405. pages 2, 6, 7

C. Bayer and B. Stemper. Deep calibration of rough stochastic volatility models. *arXiv preprint*, 2018. pages 2

C. Bayer, B. Horvath, A. Muguruza, B. Stemper, and M. Tomas. On deep calibration of (rough) stochastic volatility models, 2019. URL https://arxiv.org/abs/1908.08806. arXiv preprint: https://arxiv.org/abs/1908.08806. pages 8

E. Benhamou, E. Gobet, and M. Miri. Time dependent heston model. *SIAM Journal on Financial Mathematics*, 1(1):289–325, 2010. pages 2

D. Brigo and F. Mercurio. *Interest Rate Models, Theory and Practice*. Springer-Verlag, 2001. pages 1

D. Brigo and F. Mercurio. Lognormal-mixture dynamics and calibration to market volatility smiles. *International Journal of Theoretical and Applied Finance*, 5(4):427–446, 2002. pages 1

M. Broadie and O. Kaya. Exact simulation of stochastic volatility and other affine jump diffusion models. *Operations Research*, 54(2):217–231, 2006. pages 6

P. Carr and V. Linetsky. A jump to default extended cev model: an application of bessel processes. *Finance and Stochastics*, 10(3):303–330, 2006. pages 2

J. H. Chan and M. Joshi. Fast and accurate long-stepping simulation of the heston stochastic volatility model. *Journal of Computational Finance*, 16(3):47–97, 2013. pages 6

J. Cox and S. Ross. The valuation of options for alternative stochastic processes. *Journal of Financial Economics*, 3:145–166, 1976. pages 1

J.C. Cox, J.E. Ingersoll, and S.A. Ross. A theory of the term structure of interest rates. *Econometrica*, 53(2): 385–407, 1985. pages 6

W.M. Czarnecki, S. Osindero, M. Jaderberg, G. Świrszcz, and R. Pascanu. Sobolev training for neural networks, 2017. pages 2, 8, 9, 16

D. Dufresne. The integrated square-root process. Technical report, SSRN, 2001. pages 6

W. Feller. Diffusion processes in one dimension. *Transactions of the American Mathematical Society*, 77(1): 1–31, 1954. pages 2

P. Glasserman and K. K. Kim. Gamma expansion of the heston stochastic volatility model. *Finance and Stochastics*, 15(2):267–296, 2011. pages 6

A. Hernandez. Model calibration with neural networks. *Risk*, 2017. pages 8

S. Heston. A closed form solution for options with stochastic volatility with applications to bond and currency options. *The Review of Financial Studies*, 6(2):327–343, 1993. pages 1, 2

B. Horvath, A. Muguruza, and M. Tomas. Deep learning volatility: a deep neural network perspective on pricing and calibration in (rough) volatility models. *Quantitative Finance*, 21(1):11–27, 2021a. pages 2

B. Horvath, A. Muguruza, and M. Tomas. Deep learning volatility: a deep neural network perspective on pricing and calibration in (rough) volatility models. *Quantitative Finance*, 21(1):11–27, 2021b. pages 8

B.N. Huge and A. Savine. Differential machine learning, 2020. URL https://ssrn.com/abstract=3591734. pages 2, 8, 19, 20

A. Itkin and P. Carr. Jumps without tears: A new splitting technology for barrier options. *International Journal of Numerical Analysis and Modeling*, 8(4):667–704, 2011. pages 2

T. Kluge. *Pricing derivatives in stochastic volatility models using the Finite Difference method*. PhD thesis, Technical University, Chemnitz, 2002. pages 2

K.-L. Lai and J.L. Crassidis. Extensions of the first and second complex-step derivative approximations. *Journal of Computational and Applied Mathematics*, 219(1):276–293, 2008. ISSN 0377-0427. doi: https://doi.org/10.1016/j.cam.2007.07.026. URL https://www.sciencedirect.com/science/article/pii/S0377042707004086. pages 16
A.L. Lewis. *Option Valuation under Stochastic Volatility*. Finance Press, 2000. pages 2, 4

A. Lipton. *Mathematical Methods for Foreign Exchange: A Financial Engineer's Approach*. World Scientific Publishing Company, Incorporated, 2001. ISBN 9789810246150. pages 4

A. Lipton. The vol smile problem. *Risk Magazine*, pages 61–65, 2002. pages 2, 4

S. Liu, A. Borovykh, L.A. Grzelak, and C. Oosterlee. A neural network-based framework for financial model calibration, 2019. URL https://arxiv.org/abs/1904.10523. pages 8

S. Mikhailov and U. Nogel. Heston's stochastic volatility model: implementation, calibration and some extensions. *Wilmott Magazine*, 14(7):74–79, 2003. pages 2

M. Mrad and R. Triki. Hybrid equity-interest rate model with stochastic interest rate and volatility, 2014. URL https://ssrn.com/abstract=2434563. pages 2, 3, 5, 6

M. Raissi, P. Perdikaris, and G.E. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. *Journal of Computational Physics*, 378(C), 2019. doi: 10.1016/j.jcp.2018.10.045. pages 2, 10

F.D. Rouah. *The Heston Model and its Extensions in VBA*. John Wiley & Sons, 2015. ISBN 978-1-119-00330-4. pages 2

B. Saadeddine. Fast calibration using complex-step sobolev training, 2022. URL https://ssrn.com/abstract=4275942. pages 16

R.G. Tompkins. Stock index futures markets: stochastic volatility models and smiles. *The Journal of Future Markets*, 21(1):43–78, 2001. pages 1