1 Introduction↩︎

Volatility in modern financial markets is neither constant nor exogenous. It evolves randomly across multiple time scales and interacts with other sources of uncertainty, including liquidity conditions, interest rates, and latent risk factors. Under the classical Black–Scholes framework and its stochastic-volatility extensions, derivative prices can be characterized as solutions of partial differential equations (PDEs)  [1]–[3].

When the market environment itself evolves randomly, the pricing operator becomes stochastic, leading naturally to stochastic partial differential equations (SPDEs). Such SPDE formulations arise in a variety of financial settings, including stochastic volatility surface models [4], stochastic term-structure models [5]–[7], and stochastic forward-curve models in energy markets [8]. In these applications, the solution of the SPDE typically represents the value of a contingent claim under a random market environment. Beyond computing option prices, practical applications require the estimation of risk sensitivities, commonly known as Greeks, including Delta, Vega, Gamma, and higher-order parameter sensitivities. These quantities play a central role in dynamic hedging, model calibration, risk management, and uncertainty quantification [9].

The numerical approximation of SPDEs has been studied extensively over the past several decades [10], [11]. A wide range of deterministic and stochastic discretization techniques have been developed, including finite-difference methods, finite-element methods, spectral methods, and stochastic Galerkin approaches  [12]–[14]. These methods have achieved considerable success in the simulation and analysis of stochastic systems arising in physics, engineering, and finance.

However, the direct numerical treatment of financial SPDEs remains challenging. In multi-factor stochastic-volatility and term-structure models, the effective state dimension can grow rapidly. As a consequence, grid-based discretizations often suffer from the curse of dimensionality, resulting in substantial computational costs for both spatial discretization and the solution of large-scale linear systems  [15], [16].

Moreover, many quantities of practical interest, including option prices, risk measures, and Greek sensitivities, are naturally expressed as expectations or nested expectations of the underlying stochastic system. Accurate estimation of such quantities often requires substantial sampling effort in addition to the numerical solution of the SPDE itself.

These challenges motivate the search for alternative formulations that avoid direct discretization of the underlying SPDE. An attractive approach is provided by probabilistic representations, which reformulate the problem in terms of stochastic differential equations and expectation estimation.

A major development in this direction was the introduction of backward stochastic differential equations (BSDEs) by Pardoux and Peng [17]. BSDEs provide a nonlinear extension of the classical Feynman–Kac formula and establish a probabilistic representation for broad classes of semilinear parabolic PDEs  [18], [19]. Since their introduction, BSDEs have become an important tool in stochastic control [20], [21], mathematical finance [22], and nonlinear expectation theory  [23], [24].

To extend the probabilistic correspondence between BSDEs and PDEs to stochastic partial differential equations, Pardoux and Peng introduced backward doubly stochastic differential equations (BDSDEs) [25]. They established a stochastic Feynman–Kac formula showing that solutions of a broad class of quasilinear SPDEs can be represented by solutions of BDSDEs driven simultaneously by a forward Brownian motion and a backward Brownian motion. For the class of stochastic option-pricing SPDEs considered in this work, let \(X^{t,x}\) denote the forward state process describing the underlying risk factors, initialized from \(x\) at time \(t\). Then the SPDE solution admits a BDSDE representation of the form \[\begin{align} u(s,X_s^{t,x})=Y_s^{t,x},\quad X_t^{t,x}=x, \quad s\in[t,T], \end{align}\] where \(Y_s^{t,x}\) solves the associated BDSDE. This representation transforms the original SPDE problem into the estimation of stochastic expectations and avoids direct spatial discretization of the underlying equation. Moreover, the coexistence of forward and backward sources of randomness naturally leads to conditional and nested expectation structures, providing a probabilistic foundation for Monte Carlo, multilevel Monte Carlo, and quantum mean estimation methodologies.

The numerical solution of BSDEs and BDSDEs has attracted considerable attention over the past several decades. For BSDEs, a variety of time-discretization, regression-based, and Monte Carlo methods have been developed [26]–[28]. These approaches provide practical algorithms for approximating BSDE solutions and have led to a rich literature on numerical methods for high-dimensional PDEs through probabilistic representations. More recently, BSDE formulations have inspired a variety of machine-learning approaches for high-dimensional PDEs, including Deep BSDE methods, deep backward dynamic programming schemes, and related neural-network-based algorithms  [29]–[31]. These methods substantially improve scalability in high-dimensional settings and have become an active research direction at the interface of scientific computing and machine learning.

For BDSDEs, several numerical approximation schemes have been proposed, including Euler-type and regression-based methods  [32], [33]. In particular, [34] developed a first-order scheme based on a two-sided Itô–Taylor expansion. Since many of these schemes rely on conditional expectations with respect to the forward Brownian motion, the forward randomness is integrated out at each time step. As a consequence, they do not directly provide the pathwise strong approximations required for multilevel Monte Carlo coupling and Greek estimation.

Many quantities arising from BDSDE representations, including option prices and Greek sensitivities, are naturally expressed as stochastic expectations. To estimate such quantities efficiently, multilevel Monte Carlo (MLMC), introduced by Giles [35], exploits strong couplings between successive discretization levels to achieve substantial computational savings over standard Monte Carlo methods.

MLMC has been successfully applied to a wide range of problems in computational finance, including Greek estimation, efficient risk measurement, and basket option pricing  [36]–[39]. It has also been extended to stochastic partial differential equations [40]–[42]. In particular, [43] developed a multilevel Monte Carlo framework based on a Milstein finite difference discretization for SPDEs and demonstrated its effectiveness in the pricing of basket credit derivatives. To the best of our knowledge, however, the combination of MLMC and BDSDE representations has received little attention in the existing literature.

Recent advances in quantum computing have opened new possibilities for further accelerating stochastic simulation. Quantum amplitude estimation and related quantum mean estimation algorithms reduce the sampling complexity of expectation estimation from \(\mathcal{O}(\epsilon^{-2})\) to \(\widetilde{\mathcal{O}}(\epsilon^{-1})\), thereby providing a quadratic speedup over classical Monte Carlo methods [44]–[47]. Building on these developments, quantum-accelerated multilevel Monte Carlo (QA-MLMC) methods have recently emerged as a powerful framework for expectation estimation. [48] combined quantum mean estimation with MLMC and established quantum speedups for stochastic differential equations. Subsequent developments further extended and refined this framework, including quadratic speedups for nonlinear and nested expectation problems and related stochastic simulation tasks  [49]–[52].

At the same time, quantum algorithms have also been investigated for financial PDEs and derivative pricing. Examples include quantum algorithms for option valuation and financial simulation [53]–[55], as well as recent end-to-end quantum PDE frameworks for option pricing under Black–Scholes and Heston-type models [56]. These works demonstrate the potential of quantum computation for high-dimensional problems arising in quantitative finance. On the other hand, quantum algorithms have also been explored for the numerical solution of BSDEs. For example, [57] introduced a quantum least-squares Monte Carlo approach for solving BSDEs.

However, most existing quantum approaches mainly focus on PDEs, SDEs, or BSDEs. To the best of our knowledge, although quantum simulation algorithms for stochastic differential equations are recently studied [58]–[62], quantum algorithms for SPDEs based on BDSDE representations have not yet been systematically studied. The extension of these ideas to SPDEs faces several challenges.

A key difficulty is that achieving the optimal complexity \(\widetilde{\mathcal{O}}(\epsilon^{-1})\) requires not only quantum acceleration, but also effective multilevel coupling and sufficiently accurate pathwise discretizations. In the standard multilevel complexity analysis, if the bias, variance, and cost exponents are denoted by \((\alpha,\beta,\gamma)\), then achieving the optimal quantum complexity typically requires \[\beta\ge 2\gamma,\] which is stronger than the classical MLMC requirement \(\beta\ge\gamma\). Equivalently, if a discretization has strong convergence order \(r\) so that the level-difference variance behaves like \(\mathcal{O}(h^{2r})\), then the optimal QA-MLMC regime requires \[r=\frac{\beta}{2}\ge \gamma.\] For direct discretizations of a \(d\)-dimensional SPDE, the cost exponent \(\gamma\) is often large because each sample involves the spatial degrees of freedom of the discretized random field. For instance, if a tensor-product grid with mesh size \(h_\ell\) is used in each spatial coordinate and the cost per grid point is uniformly bounded, then the number of spatial grid points is proportional to \(h_\ell^{-d}\). Hence \(\gamma_{\mathrm{SPDE}} = d\). Thus, achieving the optimal QA-MLMC complexity would require a strong approximation order comparable to the spatial dimension (\(r\ge \gamma_{\mathrm{SPDE}}= d\)), which is generally unrealistic in high-dimensional settings.

Even after passing to a BDSDE formulation, a standard strong-error order-\(1/2\) discretization (e.g., Euler discretization) is still insufficient for the full quantum speedup. When \(\gamma=1\) and \(r=1/2\), the variance exponent is only \(\beta=1\), and the resulting QA-MLMC complexity is typically \(\widetilde{\mathcal{O}}(\epsilon^{-3/2})\) rather than \(\widetilde{\mathcal{O}}(\epsilon^{-1})\). This is why the strong-error order one Forward–Backward Taylor discretization developed in this work is essential: it yields \(r=1\), reaching the critical regime \(\beta=2\gamma\) and enabling the full quadratic quantum speedup. The comparison is summarized in Table 1.

These considerations indicate that extending quantum multilevel methods to SPDEs requires more than a direct application of existing QA-MLMC theory.

Existing QA-MLMC frameworks are primarily designed for expectation estimation problems arising from SDE-type path simulations. In contrast, the BDSDE representation of SPDE solutions introduces an additional backward source of randomness, which naturally gives rise to conditional and nested expectation structures. These structures require quantum multilevel estimators that are adapted to the two sources of randomness.

Moreover, most existing multilevel and quantum approaches in this area focus primarily on pricing problems. By contrast, efficient estimators for first- and second-order Greeks under stochastic-environment models are much less developed. This leaves open the construction of multilevel and quantum estimators for sensitivity analysis in the SPDE and BDSDE setting.

1.1 Contributions↩︎

In this work, we develop a unified computational framework for stochastic derivative-pricing SPDEs based on their BDSDE representations. The framework combines multilevel Monte Carlo, conditional and nested quantum estimators, and new Forward–Backward Taylor discretization schemes with strong-error order one convergence, allowing derivative pricing and Greek estimation to be treated within a common probabilistic setting.

1.1.0.1 Modeling and algorithms.

1. BDSDE-based framework for SPDE pricing and sensitivity analysis.

   Section 2 introduces the SPDE models arising in stochastic-environment financial markets and develops a BDSDE-based probabilistic representation. This representation reformulates pricing and Greek estimation problems as conditional and nested expectation estimation tasks, providing a unified framework for direct pricing, first-order sensitivities, and second-order sensitivities. The resulting estimation framework serves as the foundation for all subsequent quantum algorithms developed in this work.

2. Conditional quantum-accelerated multilevel Monte Carlo for SPDEs.

   Building on the BDSDE-based estimation framework, we develop a conditional QA-MLMC methodology for approximating \[\begin{align} u(t,x;B) =\mathbb{E}_W[P(t,x;B)\mid B]. \end{align}\] The key ingredients are the conditional quantum level-\(\ell\) difference evaluator (Algorithm 1), which constructs quantum estimators for the level differences appearing in the multilevel decomposition, and the conditional quantum-accelerated MLMC estimator \(\widehat{u}(t,x;B)\) (Algorithm 2), which combines these level differences across discretization levels to produce a complete conditional estimator of \(u(t,x;B)\).

3. Nested quantum-accelerated multilevel Monte Carlo for SPDEs.

   Building on recent quantum algorithms for nonlinear and nested expectation estimation  [50], [63], we specialize the nested QA-MLMC framework to the BDSDE representations arising from SPDEs. This yields a nested quantum-accelerated estimator (Algorithm 4) for quantities of the form \[\begin{align} \mathbb{E}_B\!\left[\varphi(u(t,x;B))\right]. \end{align}\] The resulting methodology combines the coupled level difference evaluator (Algorithm 3) with a nested multilevel quantum estimation procedure, thereby extending quantum multilevel techniques to nested expectation structures arising naturally in stochastic-environment SPDE models.

4. Applications to derivative pricing and sensitivity analysis.

   In Section 4, we specialize the conditional and nested QA-MLMC frameworks to derivative pricing and sensitivity analysis. Using the probabilistic representations developed in this work, we construct quantum estimators for a broad class of first-order Greeks, including Delta, spot Vega, Rho, and general parameter Greeks; smooth second-order Greeks like Gamma, Vanna and Volga; extensions to nonsmooth payoffs via likelihood-ratio and Malliavin-weight representations; and Heston-type stochastic-volatility models.

5. Global strong-error order one analysis and numerical experiments.

   In Section 5, we establish a global strong-error framework for the BDSDE payoff functionals arising in pricing, first-order Greek estimation, and second-order Greek estimation. We then construct concrete Forward–Backward Taylor discretization operators and prove that they satisfy the first-order consistency and accumulated stability properties required by this framework. Consequently, the resulting estimators achieve global strong convergence of order one.

   We further validate this theory through numerical experiments in Section 6. For pricing, first-order Greek estimation, and second-order Greek estimation, the observed multilevel exponents satisfy \[(\alpha,\beta,\gamma)\approx(1,2,1),\] confirming the strong-error order one behavior predicted by the theory and the conditions required for the quadratic quantum speedup.

1.1.0.2 Theorems.

1. BDSDE representation and QA-MLMC foundations.

   Theorem 1, following [25], [64], provides the BDSDE representation for a class of linear backward SPDEs and serves as the probabilistic foundation of our framework. Theorem 2, following [48], recalls the general quantum-accelerated MLMC methodology that underlies our quantum complexity analysis.

2. Conditional quantum speedups for SPDE solution estimation.

   Theorem 3 shows that, conditioned on a realization of the environmental Brownian motion \(B\), the SPDE solution \[\begin{align} u(t,x;B)=\mathbb{E}_W[P(t,x;B)\mid B] \end{align}\] can be estimated with additive error \(\epsilon\) using computational cost \[\begin{align} \widetilde{\mathcal{O}}(\epsilon^{-1}), \end{align}\] yielding a quadratic quantum speedup over classical Monte Carlo methods.

3. Quantum speedups for nested SPDE expectations.

   Theorem 4 extends the conditional framework to nested quantities of the form \[\begin{align} \mathbb{E}_B[\varphi(u(t,x;B))]. \end{align}\]

   Under suitable regularity and moment assumptions, the resulting nested QA-MLMC estimator achieves additive error \(\epsilon\) with overall complexity \[\begin{align} \widetilde{\mathcal{O}}(\epsilon^{-1}), \end{align}\] thereby extending quantum speedups to SPDE quantities involving both forward and backward randomness.

4. Quantum estimation of first and second-order Greeks.

   Propositions 1 and 2 derive pathwise conditional representations for first-order and second-order Greeks. These representations express the Greeks as conditional expectations of explicitly constructed payoff functionals involving first-order and second-order variational processes.

   Building on these representations, Corollaries 1 and 2 show that both conditional and unconditional Greeks can be estimated with additive error \(\epsilon\) using quantum-accelerated MLMC methods with complexity \[\begin{align} \widetilde{\mathcal{O}}(\epsilon^{-1}). \end{align}\]

5. Global strong-error order one schemes via Forward–Backward Taylor discretization.

   Proposition 3, Proposition 6, and Proposition 9 establish a general strong-error framework showing that global strong-error order one follows whenever the constituent discretization operators satisfy suitable first-order consistency and stability properties.

   To realize this framework, we introduce a new family of Forward–Backward Taylor discretization schemes for the stochastic integrals arising in the BDSDE representations of pricing and Greek estimators. The construction explicitly captures mixed forward–backward iterated integrals and provides concrete discretization operators for pricing, first-order Greeks, and second-order Greeks.

   Propositions  4,  5,  7,  8,  10, and  11 verify that the Forward–Backward Taylor discretization satisfies the required strong-error order one consistency and stability conditions. Consequently, all resulting pricing and Greek estimators achieve global strong convergence of order one.

1.2 Organization↩︎

Section 1 introduces the background and motivation of the work, reviews the related literature, and summarizes the main contributions.

Section 2 presents the SPDE models considered in this paper and develops the BDSDE-based probabilistic representation and estimation framework that serves as the foundation of our methodology.

Section 3 develops the quantum-accelerated multilevel Monte Carlo framework. Subsection 3.1 introduces the conditional QA-MLMC estimator for approximating \(u(t,x;B)=\mathbb{E}_W[P(t,x;B)\mid B]\), while Subsection 3.2 develops the nested QA-MLMC estimator for quantities of the form \(\mathbb{E}_B[\varphi(u(t,x;B))]\).

Section 4 specializes the general framework to derivative pricing and sensitivity analysis. We derive conditional pricing representations, first- and second-order Greek representations, and the corresponding conditional and nested quantum estimators. Several extensions, including nonsmooth payoffs and Heston-type stochastic-volatility models, are also discussed.

Section 5 contains the main numerical analysis component of the paper. We develop strong-error order one numerical schemes for pricing and Greek estimators and establish the strong-convergence and stability results required for the quantum-accelerated multilevel complexity analysis.

Section 6 presents numerical experiments validating the proposed Forward–Backward Taylor discretization and the multilevel convergence properties across multiple realizations of the backward Brownian motion.

Finally, Section 7 concludes the paper by summarizing the main contributions and outlining possible directions for future work.

The appendices contain additional background material and technical proofs. In particular, Appendix 8 reviews the classical BDSDE representation of SPDEs, while Appendix 9 contains the proofs of the strong-error framework and the numerical analysis results developed in Section 5.

2 SPDE Formulation and Probabilistic Representation↩︎

2.1 From Stochastic Volatility Models to SPDEs↩︎

We begin with the classical Heston stochastic model [2] under the risk-neutral measure: \[\begin{align} \tag{1} \mathrm dS_t&=rS_t\,\mathrm dt+ \sqrt{V_t}\,S_t\,\mathrm dW_t^{(S)},\\ \tag{2} \mathrm dV_t&=\kappa(\theta-V_t)\,\mathrm dt +\xi\sqrt{V_t}\,\mathrm dW_t^{(V)}, \end{align}\] where \(S_t\) denotes the asset price, \(V_t\) is the instantaneous variance process, \(r\) is the risk-free interest rate, \(\kappa\) is the mean-reversion rate, \(\theta\) is the long-run variance level, \(\xi\) is the volatility-of-volatility parameter, and \[\begin{align} \label{eq:Heston32cor} \mathrm d\langle W^{(S)},W^{(V)}\rangle_t =\rho\,\mathrm dt. \end{align}\tag{3}\]

For a European payoff \(g(S_T)\), the option price is \[\begin{align} u(t,s,v) =\mathbb{E}^{\mathbb{Q}} \left[e^{-r(T-t)}g(S_T) \mid S_t=s,\,V_t=v \right]. \end{align}\]

Under standard regularity assumptions, the standard arbitrage-free pricing arguments [1], [65] imply that the option price \(u(t,s,v)\) satisfies the partial differential equation (PDE) \[\begin{align} \partial_t u+\mathcal{L} u-r u=0, \end{align}\] with terminal condition \[u(T,s,v)=g(s),\] where the generator \(\mathcal{L}\) associated with 1 –3 is given by \[\begin{align} \mathcal{L} =rs\partial_s +\kappa(\theta-v)\partial_v +\frac{1}{2} vs^2\partial_{ss} +\rho\xi sv\,\partial_{sv} +\frac{1}{2}\xi^2v\,\partial_{vv}. \end{align}\]

The above pricing equation is deterministic because all market factors are assumed to be either deterministic or fully represented by the finite-dimensional state variables \((S_t,V_t)\).

In modern financial markets, however, quantities such as volatility surfaces, liquidity conditions, interest-rate environments, and latent risk factors may themselves evolve randomly over time. As a consequence, the coefficients of the pricing equation become random, and the option value must be viewed as a random field.

Motivated by this observation, we consider stochastic option-pricing SPDEs of the form \[\begin{align} \label{eq:spde0} \mathrm du=\bigl(\mathcal{L} u -r u+F\bigr)\,\mathrm dt+G\,\mathrm dB_t, \end{align}\tag{4}\] with terminal condition \[\begin{align} u(T,s,v)=g(s). \end{align}\] Here, \(\mathcal{L}\) denotes the stochastic-volatility generator, \(r\) is the discount rate, \(F\) is a source term, and \(G\) characterizes the stochastic forcing induced by the random market environment.

Direct numerical discretization of 4 is computationally demanding, particularly in high-dimensional settings. To facilitate efficient numerical approximation, we employ a probabilistic representation based on backward doubly stochastic differential equations (BDSDEs). This representation transforms the SPDE problem into the estimation of conditional expectations over stochastic paths and serves as the foundation of the algorithms developed in this paper.

2.2 BDSDE Representation and Estimation Framework↩︎

A key advantage of stochastic partial differential equations is that, under suitable regularity conditions, their solutions admit probabilistic representations through backward doubly stochastic differential equations (BDSDEs). Introduced by Pardoux and Peng [25], BDSDEs extend the classical Feynman–Kac formula and provide a probabilistic representation for a broad class of quasilinear SPDEs.

In this work, the SPDE is treated through its associated BDSDE representation. Since the SPDE–BDSDE correspondence is classical, we only summarize the formulation required for the subsequent numerical analysis. Further details can be found in Appendix 8.

In our applications, we are primarily interested in the linear setting. Following [64], we consider the following class of linear backward SPDEs.

Theorem 1 reduces the solution of the SPDE to the evaluation of a path functional \(P(t,x;B)\). For a fixed realization of the backward Brownian motion \(B\), the only remaining randomness in ?? comes from the forward Brownian motion \(W\). Consequently, the SPDE solution admits the conditional expectation representation \[\begin{align} u(t,x;B) =\mathbb{E}_W \bigl[P(t,x;B)\,\big|\,B\bigr]. \end{align}\] This observation transforms the original SPDE problem into the estimation of a conditional expectation with respect to the forward diffusion.

In many applications, the quantity of ultimate interest is obtained after averaging over the random environment generated by \(B\). This leads to the nested expectation \[\begin{align} U(t,x) =\mathbb{E}_B \bigl[u(t,x;B)\bigr] =\mathbb{E}_B \Big[\mathbb{E}_W \bigl[P(t,x;B)\,\big|\,B\bigr]\Big]. \end{align}\]

Consequently, the probabilistic representation naturally gives rise to two estimation problems. The first is the conditional estimation problem \[\begin{align} u(t,x;B) =\mathbb{E}_W \bigl[P(t,x;B)\,\big|\,B\bigr], \end{align}\] where the realization of \(B\) is fixed. The second is the nested estimation problem \[\begin{align} U(t,x) =\mathbb{E}_B \bigl[u(t,x;B)\bigr], \end{align}\] which requires averaging over both sources of randomness.

These two estimation problems form the basis of the multilevel Monte Carlo and quantum algorithms developed in the subsequent sections.

3 Quantum-accelerated MLMC for SPDEs↩︎

The probabilistic representation developed in the previous section reduces the SPDE problem to the estimation of conditional and nested expectations associated with the path functional \(P(t,x;B)\). To estimate efficiently, we employ multilevel Monte Carlo (MLMC) together with quantum mean estimation.

Suppose that a sequence of approximations \(\{\Phi_\ell\}_{\ell\ge0}\) to a random variable \(\Phi\) is available and satisfies \[\begin{align} \left|\mathbb{E}[\Phi]-\mathbb{E}[\Phi_\ell]\right|=\mathcal{O}(2^{-\alpha\ell}), \quad \operatorname{Var}(\Phi_\ell-\Phi_{\ell-1}) =\mathcal{O}(2^{-\beta\ell}), \quad \operatorname{Cost}(\Phi_\ell) =\mathcal{O}(2^{\gamma\ell}), \end{align}\] for some \(\alpha,\beta,\gamma>0\). The MLMC estimator is based on the telescoping identity \[\begin{align} \mathbb{E}[\Phi_L]=\mathbb{E}[\Phi_0] +\sum_{\ell=1}^{L}\mathbb{E}[\Phi_\ell-\Phi_{\ell-1}]. \end{align}\] Under the standard MLMC framework [35], [66], if \(\beta>\gamma\), an estimator with mean-square error \(\mathcal{O}(\epsilon^2)\) can be constructed with overall computational complexity \(\mathcal{O}(\epsilon^{-2})\).

Recent advances in quantum computing provide the possibility of further accelerating expectation estimation. In particular, quantum mean estimation algorithms reduce the sampling complexity from \(\mathcal{O}(\epsilon^{-2})\) to \(\widetilde{\mathcal{O}}(\epsilon^{-1})\), yielding a quadratic speedup over classical Monte Carlo methods. The following result will be used as the basic quantum subroutine.

Using the standard powering argument [67], the success probability can be amplified from a constant to \(1-\delta\) with an additional \(\mathcal{O}(\log(1/\delta))\) overhead.

For notational simplicity, we write \[\operatorname{QME}(\Phi;\epsilon,\delta)\] for a quantum mean estimation subroutine which returns an estimate \(\widehat \mu\) of \(\mathbb{E}[\Phi]\) such that \[\Pr\left( |\widehat\mu-\mathbb{E}[\Phi]|\le \epsilon \right) \ge 1-\delta .\]

Replacing the classical sample average on each MLMC level by QME leads to quantum-accelerated MLMC (QA-MLMC). Our algorithm is mainly based on the following quantum-accelerated Monte Carlo methods proposed in [48].

However, the complexity requirements for QA-MLMC are substantially more restrictive. In particular, achieving the optimal quantum complexity \(\widetilde{\mathcal{O}}(\epsilon^{-1})\) requires \(\beta\ge 2\gamma\), which is stronger than the classical MLMC requirement \(\beta\ge\gamma\). For discretization-based SPDE solvers, the cost exponent \(\gamma\) is often large due to the curse of dimensionality. For example, under standard tensor-product spatial discretizations of a \(d\)-dimensional SPDE, one typically has \(\gamma\approx d.\) Since the MLMC variance exponent typically satisfies \(\beta=2r\), where \(r\) is the strong convergence order, the condition \(\beta\ge2\gamma\) would require \(r\ge d,\) namely, a strong convergence order at least equal to the spatial dimension, which is generally unrealistic in practice.

The BDSDE representation developed in the previous section avoids direct spatial discretization of the SPDE and converts the problem into the estimation of stochastic path functionals. This structure enables the construction of efficient multilevel estimators with favorable variance decay properties and provides a natural foundation for quantum acceleration.

In the remainder of this section, we first consider the conditional estimation problem and then extend the framework to the nested setting.

3.1 Conditional Quantum-accelerated MLMC Estimator↩︎

The probabilistic representation introduced in Subsection 2.2 reduces the pricing problem to the estimation of the conditional expectation \[\begin{align} u(t,x;B) =\mathbb{E}_W[P(t,x;B)\mid B]. \end{align}\] To construct the conditional quantum-accelerated estimators, we first develop suitable discretizations for the stochastic quantities appearing in this representation. Recalling that \[\begin{align} P(t,x;B) =\Phi(t,T)G(X_T^{t,x}) +\int_t^T\Phi(t,r)\left(F(r,X_r^{t,x})+d(r) H(r,X_r^{t,x})\right)\,\mathrm{d}r +\int_t^T\Phi(t,r)H(r,X_r^{t,x})\mathrm d \overleftarrow{B}_r, \end{align}\] we observe that there are three components that require discretization: the evolution of the SDE, the exponential weight \(\Phi\), and the integrals over the time interval \([t,T]\).

In classical numerical schemes, these components are often discretized using Euler-type methods, which already provide sufficient accuracy. However, in quantum algorithms, higher-order convergence is often required. Therefore, we introduce a unified notation \(\mathcal{S}\) to denote the discretization operators associated with each component, allowing us to systematically track and analyze their respective convergence properties.

We introduce three discretization operators. Let \[\begin{align} \Delta W_k:=W_{t_{k+1}}-W_{t_k}, \quad \Delta\overleftarrow B_k:=B_{t_k}-B_{t_{k+1}}, \quad t_k=t+kh . \end{align}\]

(1) State discretization \(\mathcal{S}_X\). Given the current state \(X_k\), the time \(t_k\), the stepsize \(h\), and the forward Brownian increment \(\Delta W_k\), we set \[\begin{align} \label{eq:def32of32S95X} X_{k+1}=\mathcal{S}_X(X_k,t_k,h;\Delta W_k). \end{align}\tag{5}\]

(2) Exponential weight discretization \(\mathcal{S}_\Phi\). Given the current weight \(\Phi_k\), the time \(t_k\), the stepsize \(h\), the forward Brownian increment \(\Delta W_k\), and the backward Brownian increment \(\Delta\overleftarrow B_k\), we set \[\begin{align} \label{eq:def32of32S95phi} \Phi_{k+1}=\mathcal{S}_\Phi(\Phi_k,t_k,h; \Delta W_k,\Delta\overleftarrow B_k). \end{align}\tag{6}\]

(3) Integral discretization \(\mathcal{S}_{\mathrm{int}}\). Given the current term \(Y_k\), weight \(\Phi_k\), state \(x_k\), the time \(t_k\), the stepsize \(h\), and the Brownian increments \((\Delta W_k,\Delta\overleftarrow B_k)\), we set \[\begin{align} \label{eq:def32of32S95int} Y_{k+1}=Y_k+ \mathcal{S}_{\mathrm{int}}\left( \Phi_k,x_k,t_k,h; \Delta W_k,\Delta\overleftarrow B_k \right). \end{align}\tag{7}\]

Assume that reversible unitary implementations of the discretization operators \(\mathcal{S}_X,\mathcal{S}_\Phi\) and \(\mathcal{S}_{\mathrm{int}}\) are available in the following standard form: \[\begin{align} O_X &: |x,s,h,\Delta W\rangle|0\rangle \mapsto |x,s,h,\Delta W \rangle\ket{\mathcal{S}_X(x,s,h;\Delta W)},\\ O_\Phi &: |\Phi,s,h,\Delta W,\Delta\overleftarrow B\rangle|0\rangle \mapsto |\Phi,s,h,\Delta W,\Delta\overleftarrow B\rangle \ket{\mathcal{S}_\Phi(\Phi,s,h;\Delta W,\Delta\overleftarrow B)},\\ O_\mathrm{int}&: \ket{Y,\Phi,x,s,h;\Delta W,\Delta\overleftarrow B}\ket{0} \mapsto \ket{Y,\Phi,x,s,h;\Delta W,\Delta\overleftarrow B} \ket{Y+\mathcal{S}_\mathrm{int}\left(\Phi,x,s,h;\Delta W,\Delta\overleftarrow B\right)}. \end{align}\]

For notational simplicity, we combine these elementary reversible updates into a single one-step oracle. Suppressing unchanged input registers, work registers, and the standard uncomputation of ancillas, we write, for \(\eta\in\{h,2h\}\), \[\mathcal{U}_{\mathcal{S},\eta}(s;\Delta W,\Delta\overleftarrow B): \ket{X,\Phi,Y} \mapsto \ket{ \mathcal{S}_X(X,s,\eta;\Delta W), \mathcal{S}_\Phi(\Phi,s,\eta;\Delta W,\Delta\overleftarrow B), Y+\mathcal{S}_{\rm int}(\Phi,X,s,\eta;\Delta W,\Delta\overleftarrow B) }.\]

The returned quantum state \(\ket{\Delta_{\mathcal{S}}(h;B)}\) encodes a coupled multilevel difference sample conditioned on \(B\). For notational uniformity, define \[\Delta_0(t,x;B):=P_0(t,x;B), \quad \Delta_\ell(t,x;B) := P_\ell^{f}(t,x;B)-P_{\ell}^{c}(t,x;B), \quad \ell\ge1.\] The fine and coarse trajectories are generated using the shared Brownian coupling \[\Delta W_k^c = \Delta W_{2k}^f+\Delta W_{2k+1}^f, \quad \Delta\overleftarrow B_k^c = \Delta\overleftarrow B_{2k}^f + \Delta\overleftarrow B_{2k+1}^f.\]

For a fixed realization of \(B\), the conditional multilevel identity gives \[u_L(t,x;B) = \mathbb{E}_W[P_L(t,x;B)\mid B] = \sum_{\ell=0}^{L} \mathbb{E}_W[\Delta_\ell(t,x;B)\mid B].\] By Lemma 1, each conditional level mean can be estimated by applying QME to the corresponding quantum encoding.

Algorithm 2 combines the conditional multilevel telescoping decomposition with quantum mean estimation applied independently on each level. The resulting estimator approximates the conditional quantity \[u_L(t,x;B) =\mathbb{E}_W[P_L(t,x;B)\mid B].\]

3.2 Nested Quantum-accelerated MLMC Estimator↩︎

Furthermore, let \(\varphi:\mathbb{R}\to\mathbb{R}\) be a globally Lipschitz function. In many applications, the quantity of interest is obtained by averaging over the random environment generated by the backward Brownian motion \(B\). This leads to the nested expectation \[\begin{align} \mathbb{E}_B\bigl[\varphi(u(t,x;B))\bigr]. \end{align}\] Following the framework developed in [50], we construct quantum-accelerated estimators for this quantity.

4 Applications to Derivative Pricing and Sensitivity Analysis↩︎

This section specializes the conditional and nested QA-MLMC framework of Sections 2–3 to derivative pricing and sensitivity analysis. The common structure is the following: after fixing the backward Brownian path \(B\), each target quantity is represented as \[\mathbb{E}_W[\mathcal{P}(t,x;B)\mid B],\] where \(\mathcal{P}(t,x;B)\) is a scalar path functional of the forward diffusion. Once such a representation is available, the same dyadic fine–coarse coupling and the same quantum mean estimation routine used for the price can be applied to the corresponding payoff register.

Throughout this section, unless otherwise stated, we work in the scalar case \(k=1\). Vector-valued claims are obtained componentwise. We write \(\mathbb{E}_W[\cdot\mid B]\) for expectation over the forward Brownian motion after conditioning on the realization of the backward Brownian motion \(B\). The terminal function in the BDSDE representation is denoted by \(G\); when the terminal condition is written as \(g\) elsewhere in the paper, one simply sets \(G=g\).

4.1 Conditional Pricing under Common Noise↩︎

By ?? –?? , the conditional value of the linear backward SPDE is \[\label{eq:app-price-representation} u(t,x;B) = \mathbb{E}_W\!\left[P(t,x;B)\mid B\right],\tag{8}\] where \[\begin{align} \label{eq:app-payoff-functional} P(t,x;B) &= \Phi(t,T)\,G\!\left(X_T^{t,x}\right) +\int_t^T \Phi(t,s) \Big(F(s,X_s^{t,x})+d(s)H(s,X_s^{t,x})\Big)\,\mathrm ds +\int_t^T \Phi(t,s)H(s,X_s^{t,x})\,\mathrm d \overleftarrow{B}_s. \end{align}\tag{9}\] Here \(X^{t,x}\) is the forward diffusion driven by \(W\), while the backward stochastic integral is evaluated along the fixed path \(B\).

4.1.0.1 Pure discounted pricing.

For a European claim with discounted terminal payoff only, take \[\label{eq:app-pure-pricing} G(x)=g(x),\quad F\equiv0,\quad H\equiv0,\quad c(s)=-r(s),\quad d(s)=\widetilde{c}(s)=0 .\tag{10}\] Then \[\label{eq:app-european-price} u(t,x;B) = \mathbb{E}_W\!\left[ \exp\!\left(-\int_t^T r(s)\,\mathrm ds\right) g\!\left(X_T^{t,x}\right) \;\middle|\; B \right].\tag{11}\] In this pure terminal-payoff case, the dependence on \(B\) disappears unless the forward model itself contains common-noise coefficients. Nonzero \(F\) covers running cashflows, and nonzero \(H\) encodes exposure to the backward common-noise factor.

4.1.0.2 Compatibility with the QA-MLMC algorithms.

For each level \(\ell\), let \(P_\ell(t,x;B)\) be the time-discretized version of 9 generated by the operators \[(\mathcal{S}_X,\mathcal{S}_\Phi,\mathcal{S}_{\mathrm{int}})\] and by the same dyadic Brownian coupling as in Algorithm 1. The conditional telescoping identity is \[\label{eq:app-price-telescoping} \mathbb{E}_W[P_L(t,x;B)\mid B] = \mathbb{E}_W[P_0(t,x;B)\mid B] + \sum_{\ell=1}^L \mathbb{E}_W[P_\ell(t,x;B)-P_{\ell-1}(t,x;B)\mid B].\tag{12}\] Therefore Algorithm 2 estimates \(u(t,x;B)\). The unconditional price \[\label{eq:app-unconditional-price} U(t,x):=\mathbb{E}_B[u(t,x;B)]\tag{13}\] is obtained by Algorithm 4 with \(\varphi(z)=z\).

4.2 First-order Greeks for Smooth Payoffs↩︎

We first treat pathwise first-order Greeks. In this subsection the coefficients \(c,\widetilde{c},d\) are deterministic functions of time, as in the linear BDSDE representation. Hence the exponential weight \(\Phi(t,s)\) depends on \(W\) and \(B\), but not on the initial state \(x\). If one allows state-dependent discounting or state-dependent coefficients in \(\Phi\), then the corresponding spatial derivatives of \(\Phi\) must be added to the formulas below.

4.2.1 Delta↩︎

Assume that the first state component is the spot, so that \(x=(S,\chi)\), where \(\chi\) denotes all remaining factors. The conditional Delta is \[\label{eq:app-delta} \Delta(t,x;B) := \partial_S u(t,x;B) = \mathbb{E}_W\!\left[P^{(S)}(t,x;B)\mid B\right].\tag{14}\] In the pure pricing case 10 , if \(G(x)=g(S)\) depends only on the terminal spot, then \[\label{eq:app-delta-pure} P^{(S)}(t,x;B) = \Phi(t,T)\, g'(S_T^{t,x})\,\partial_S S_T^{t,x}.\tag{15}\] Thus Delta is obtained by propagating only the tangent direction \(J_s^{t,x}e_S\).

4.2.2 Spot Vega↩︎

If the state contains an instantaneous variance or volatility factor \(V\), the spot Vega is the sensitivity with respect to that state variable: \[\label{eq:app-spot-vega} \mathcal{V}_{\mathrm{spot}}(t,x;B) := \partial_V u(t,x;B) = \mathbb{E}_W\!\left[P^{(V)}(t,x;B)\mid B\right].\tag{16}\] If \(G(x)=g(S)\) depends only on the terminal spot, then \[\label{eq:app-spot-vega-pure} P^{(V)}(t,x;B) = \Phi(t,T)\, g'(S_T^{t,x})\,\partial_V S_T^{t,x}.\tag{17}\] The dependence on the initial variance or volatility is transmitted through the tangent flow.

4.2.3 Rho and general parameter Greeks↩︎

State sensitivities are only one class of Greeks. Let \(\vartheta\) be a scalar parameter entering the forward coefficients, the terminal payoff, the running terms, or the exponential weight. We write \[X^{\vartheta,t,x},\quad \Phi^\vartheta,\quad P^\vartheta(t,x;B),\quad u^\vartheta(t,x;B) = \mathbb{E}_W[P^\vartheta(t,x;B)\mid B].\] Assume differentiability in \(\vartheta\) and sufficient moment bounds for differentiation under \(\mathbb{E}_W[\cdot\mid B]\).

Define the parameter tangent process \[\label{eq:app-parameter-tangent} U_s^\vartheta := \partial_\vartheta X_s^{\vartheta,t,x}, \quad U_t^\vartheta=0 .\tag{18}\] It solves \[\begin{align} \label{eq:app-parameter-tangent-sde} dU_s^\vartheta &= \Big[ \partial_\vartheta b^\vartheta(X_s^{\vartheta,t,x}) + \nabla_x b^\vartheta(X_s^{\vartheta,t,x})U_s^\vartheta \Big]\,ds + \sum_{a=1}^{d} \Big[ \partial_\vartheta\sigma_{\cdot a}^\vartheta(X_s^{\vartheta,t,x}) + \nabla_x\sigma_{\cdot a}^\vartheta(X_s^{\vartheta,t,x})U_s^\vartheta \Big]\,dW_s^a . \end{align}\tag{19}\]

When the coefficients \(c^\vartheta,d^\vartheta,\widetilde{c}^\vartheta\) in \(\Phi^\vartheta\) depend on \(\vartheta\), define \[\begin{align} \label{eq:app-logPhi-derivative} \Lambda_\vartheta(t,s) &:= \int_t^s \partial_\vartheta c^\vartheta(r)\,\mathrm dr + \int_t^s \partial_\vartheta d^\vartheta(r)\, \mathrm d \overleftarrow{B}_r + \int_t^s \partial_\vartheta\widetilde{c}^\vartheta(r)\,\mathrm dW_r - \int_t^s \Big( \widetilde{c}^\vartheta(r)\partial_\vartheta\widetilde{c}^\vartheta(r) - d^\vartheta(r)\partial_\vartheta d^\vartheta(r) \Big)\,\mathrm dr . \end{align}\tag{20}\] Then \[\label{eq:app-Phi-derivative} \partial_\vartheta\Phi^\vartheta(t,s) = \Phi^\vartheta(t,s)\Lambda_\vartheta(t,s).\tag{21}\] For notational compactness, write \[A^\vartheta(s,x) := F^\vartheta(s,x)+d^\vartheta(s)H^\vartheta(s,x).\] Then \[\label{eq:app-general-parameter-greek} \partial_\vartheta u^\vartheta(t,x;B) = \mathbb{E}_W\!\left[Q^{(\vartheta)}(t,x;B)\mid B\right],\tag{22}\] where \[\begin{align} \label{eq:app-general-parameter-greek-payoff} Q^{(\vartheta)}(t,x;B) &= \partial_\vartheta\Phi^\vartheta(t,T)\, G^\vartheta(X_T^{\vartheta,t,x}) + \Phi^\vartheta(t,T) \Big( \partial_\vartheta G^\vartheta(X_T^{\vartheta,t,x}) + \nabla_x G^\vartheta(X_T^{\vartheta,t,x})^\top U_T^\vartheta \Big) \nonumber\\ &\quad + \int_t^T \Big[ \partial_\vartheta\Phi^\vartheta(t,s)\, A^\vartheta(s,X_s^{\vartheta,t,x}) + \Phi^\vartheta(t,s) \Big( \partial_\vartheta A^\vartheta(s,X_s^{\vartheta,t,x}) + \nabla_x A^\vartheta(s,X_s^{\vartheta,t,x})^\top U_s^\vartheta \Big) \Big]\,\mathrm ds\notag\\ &\quad+ \int_t^T \Big[ \partial_\vartheta\Phi^\vartheta(t,s)\, H^\vartheta(s,X_s^{\vartheta,t,x}) + \Phi^\vartheta(t,s) \Big( \partial_\vartheta H^\vartheta(s,X_s^{\vartheta,t,x}) + \nabla_x H^\vartheta(s,X_s^{\vartheta,t,x})^\top U_s^\vartheta \Big) \Big]\, \mathrm d \overleftarrow{B}_s. \end{align}\tag{23}\]

4.2.3.1 Discount-only Rho.

If the short rate enters only through the discount coefficient \(c^\eta(s)=-r^\eta(s)\), with \[r^\eta(s)=r(s)+\eta\psi(s),\] and all other coefficients, including the forward dynamics, are kept fixed, then in the pure pricing setting 10 , \[\label{eq:app-rho-discount-only} \mathrm{Rho}^{\mathrm{disc}}_{\psi}(t,x;B) := \left.\partial_\eta u^\eta(t,x;B)\right|_{\eta=0} = - \mathbb{E}_W\!\left[ \left(\int_t^T\psi(s)\,\mathrm ds\right) \Phi(t,T)G(X_T^{t,x}) \;\middle|\;B \right].\tag{24}\] For a constant short-rate shift, \(\psi\equiv1\), \[\label{eq:app-rho-discount-only-constant} \mathrm{Rho}^{\mathrm{disc}}(t,x;B) = -(T-t)u(t,x;B).\tag{25}\] This identity is a discount-only identity. In risk-neutral equity models where \(r\) also enters the drift of \(X\), the full Rho must instead be computed from 23 , including the parameter tangent \(U^\vartheta\).

4.2.3.2 Discrete parameter Greeks.

For the augmented Milstein dynamics of Section 5, let \[\widetilde{X}_n=(S_n,V_n,D_n,Y_n,Z_n)\] and write one step as \[\widetilde{X}_{n+1} = \mathcal{S}(\widetilde{X}_n;\xi_n),\] where \(\xi_n\) collects the Brownian increments and any iterated stochastic integrals required by the scheme. A discrete parameter Greek is obtained by propagating \[\label{eq:app-discrete-parameter-tangent} \dot{\widetilde{X}}_{n+1} = \nabla_{\widetilde{X}}\mathcal{S}(\widetilde{X}_n;\xi_n) \dot{\widetilde{X}}_n + \partial_\vartheta\mathcal{S}(\widetilde{X}_n;\xi_n).\tag{26}\] The terminal Greek payoff is obtained by differentiating the terminal payoff register, for example \[D_NG(S_N)+Y_N+Z_N .\] This covers full Rho as well as model-parameter sensitivities such as \(\partial_\kappa u\), \(\partial_\xi u\), \(\partial_\rho u\), and \(\partial_\theta u\).

4.3 Second-order Greeks for Smooth Payoffs↩︎

Second-order Greeks are obtained by differentiating the pathwise first-order representation once more.

4.3.1 Gamma↩︎

When \(x=(S,\chi)\) and the first component is the spot, the conditional Gamma is \[\label{eq:app-gamma} \Gamma(t,x;B) := \partial_{SS}u(t,x;B) = \mathbb{E}_W\!\left[P^{(SS)}(t,x;B)\mid B\right].\tag{27}\] In the pure pricing case 10 , if \(G(x)=g(S)\), then \[\label{eq:app-gamma-pure} P^{(SS)}(t,x;B) = \Phi(t,T) \left[ g''(S_T^{t,x})(\partial_SS_T^{t,x})^2 + g'(S_T^{t,x})\partial_{SS}^2S_T^{t,x} \right].\tag{28}\]

4.3.2 Vanna and Volga↩︎

If the state contains a variance or volatility factor \(V\), then \[\label{eq:app-vanna-volga} \mathrm{Vanna}(t,x;B) := \partial_{SV}u(t,x;B) = \mathbb{E}_W\!\left[P^{(SV)}(t,x;B)\mid B\right], \quad \mathrm{Volga}(t,x;B) := \partial_{VV}u(t,x;B) = \mathbb{E}_W\!\left[P^{(VV)}(t,x;B)\mid B\right].\tag{29}\] If \(G(x)=g(S)\) depends only on the terminal spot, then \[\begin{align} P^{(SV)}(t,x;B) &= \Phi(t,T) \left[ g''(S_T^{t,x}) \partial_SS_T^{t,x}\, \partial_VS_T^{t,x} + g'(S_T^{t,x})\partial_{SV}^2S_T^{t,x} \right], \tag{30} \\ P^{(VV)}(t,x;B) &= \Phi(t,T) \left[ g''(S_T^{t,x}) (\partial_VS_T^{t,x})^2 + g'(S_T^{t,x})\partial_{VV}^2S_T^{t,x} \right]. \tag{31} \end{align}\]

4.4 Greek Estimators in the MLMC and QA-MLMC Pipeline↩︎

The formulas above identify the exact scalar random variables whose conditional expectations are the Greeks. To use them in the algorithms, one augments the one-step discretization by tangent variables.

4.4.0.1 First-order Greeks.

Let \[\xi_k\] collect all random inputs required by the one-step scheme, including \(\Delta W_k\), \(\Delta\overleftarrow B_k\), and, when used, mixed iterated integrals such as \(J_{t_k,h}^{WB}\). A first-order Greek discretization has the form \[\label{eq:app-discrete-jacobian} J_{k+1} = \mathcal{S}_J(J_k,X_k,t_k,h;\Delta W_{t_k},\Delta\overleftarrow B_{t_k}),\tag{32}\] together with a Greek-payoff update \[\label{eq:app-discrete-first-order-payoff} Y_{k+1}^{(i)} = Y_k^{(i)} + \mathcal{S}_{\mathrm{int}}^{(i)} \big( \Phi_k,X_k,J_k,t_k,h;\Delta W_{t_k},\Delta\overleftarrow B_{t_k} \big),\tag{33}\] whose terminal value approximates ?? .

4.4.0.2 Second-order Greeks.

For second-order sensitivities, add \[\label{eq:app-discrete-second-tangent} K_{k+1}^{(ij)} = \mathcal{S}_K\left(K_k^{(ij)},J_k,X_k,t_k,h;\Delta W_{t_k},\Delta\overleftarrow B_{t_k}\right),\tag{34}\] and \[\label{eq:app-discrete-second-order-payoff} Y_{k+1}^{(ij)} = Y_k^{(ij)} + \mathcal{S}_{\mathrm{int}}^{(ij)} \big( \Phi_k,X_k,J_k,K_k^{(ij)},t_k,h;\Delta W_{t_k},\Delta\overleftarrow B_{t_k} \big),\tag{35}\] whose terminal value approximates ?? .

Using the same dyadic Brownian coupling as in Algorithm 1, define the level differences \[\label{eq:app-greek-level-differences} \Delta_\ell^{(i)} := P_\ell^{(i,f)}-P_{\ell-1}^{(i,c)}, \quad \Delta_\ell^{(ij)} := P_\ell^{(ij,f)}-P_{\ell-1}^{(ij,c)} .\tag{36}\] These level differences are fed into Algorithm 2 in exactly the same way as the price level differences.

4.5 Nonsmooth Payoffs: Likelihood-ratio and Malliavin-weight Greeks↩︎

The pathwise Greek formulas above were stated under differentiability assumptions on the terminal function \(G\). Hence they do not directly cover nonsmooth payoffs. For Lipschitz payoffs with isolated kinks, such as vanilla calls, first-order pathwise estimators may still be justified under suitable non-atomicity or density assumptions on \(X_T^{t,x}\). However, discontinuous payoffs, barrier-type path functionals, and higher-order Greeks generally require additional smoothing, likelihood-ratio, or Malliavin integration-by-parts arguments [68]–[70].

For illustration, consider the pure terminal-value case 10 , so that \(F=H=0\), \(c=-r\), and \(d=\widetilde{c}=0\). In this case \(\Phi(t,T)\) is independent of the forward Brownian motion after conditioning on \(B\). Assume that the forward diffusion is sufficiently smooth and uniformly elliptic, with \[a(s,x):=\sigma(s,x)\sigma(s,x)^\top \succeq \lambda I_d\] for some \(\lambda>0\).

If the transition density \(p(t,x;T,y)\) exists and is differentiable with respect to the initial condition \(x\), then the likelihood-ratio representation has the form [68] \[\label{eq:app-lr-greek} \partial_{x_i}u(t,x;B) = \mathbb{E}_W\!\left[ \Phi(t,T)G(X_T^{t,x}) \,\partial_{x_i}\log p(t,x;T,X_T^{t,x}) \;\middle|\;B \right].\tag{37}\]

Alternatively, one may use a Malliavin-weight representation. Define \[\label{eq:app-bel-weight-general} \mathcal{W}_i(t,T) := \frac{1}{T-t} \int_t^T \left( \sigma(s,X_s^{t,x})^\top a(s,X_s^{t,x})^{-1} J_s^i \right)^\top \mathrm dW_s .\tag{38}\] When \(\sigma\) is square and invertible, this reduces to \[\label{eq:app-bel-weight-square} \mathcal{W}_i(t,T) = \frac{1}{T-t} \int_t^T \left( \sigma(s,X_s^{t,x})^{-1}J_s^i \right)^\top \mathrm dW_s .\tag{39}\] The Bismut–Elworthy–Li formula then gives [69], [71], [72] \[\label{eq:app-bel-directional} \partial_{x_i}u(t,x;B) = \mathbb{E}_W\!\left[ \Phi(t,T)G(X_T^{t,x})\,\mathcal{W}_i(t,T) \;\middle|\;B \right].\tag{40}\] Thus the derivative is again represented as the conditional expectation of a scalar path functional, with the payoff derivative replaced by a stochastic weight.

The simple form 40 uses the fact that, in the pure terminal-value setting, \(\Phi(t,T)\) is independent of \(W\). If \(\widetilde{c}\neq0\), or more generally if the exponential factor or running payoff contains \(W\)-dependent terms, the Malliavin integration-by-parts formula must also account for the Malliavin derivative of those terms. Such general weighted representations can still be incorporated into the same conditional-expectation framework, but the weight and payoff register must be modified accordingly.

The conditional and nested QA-MLMC architecture can be applied to these weighted-payoff estimators once the corresponding discretized weighted payoff registers satisfy the same bias, variance-decay, cost, and moment assumptions used in Corollaries 1 and 2. In particular, one needs \[\left| \mathbb{E}_W[P_{\ell}^{\mathrm{w}}-P^{\mathrm{w}}\mid B] \right| = O(2^{-\ell}), \quad \operatorname{Var}_W(P_{\ell}^{\mathrm{w}}-P_{\ell-1}^{\mathrm{w}}\mid B) = O(2^{-2\ell}), \quad C_\ell=O(2^\ell),\] together with sufficient moment or tail bounds for the likelihood-ratio or Malliavin weight. If the weight is unbounded, clipping or truncation may be used only after controlling the induced bias at the target accuracy. Barrier and other path-dependent discontinuous payoffs require path-dependent Malliavin weights or conditional-smoothing arguments and are not covered by the terminal-value formula above [69], [73], [74].

4.6 Heston-type Stochastic-volatility Models↩︎

As a canonical two-factor example, consider the Heston stochastic-volatility model [2]. To match the independent Brownian-input convention used in the algorithms, write \[\label{eq:app-heston} \begin{align} dS_s &= rS_s\,ds+\sqrt{V_s}S_s\,dW_s^{(1)},\\ dV_s &= \kappa(\theta-V_s)\,ds + \xi\sqrt{V_s} \left( \rho\,dW_s^{(1)} + \sqrt{1-\rho^2}\,dW_s^{(2)} \right), \end{align}\tag{41}\] where \(W^{(1)}\) and \(W^{(2)}\) are independent Brownian motions. Equivalently, one may use correlated Brownian motions with \(d\langle W^{(1)},W^{(2)}\rangle_s=\rho\,ds\).

With state vector \(x=(S,V)\), the conditional sensitivities above become \[\Delta=\partial_Su,\quad \mathcal{V}_{\mathrm{spot}}=\partial_Vu,\quad \Gamma=\partial_{SS}u,\quad \mathrm{Vanna}=\partial_{SV}u,\quad \mathrm{Volga}=\partial_{VV}u.\] Rho and the model-parameter Greeks \[\partial_\kappa u,\quad \partial_\theta u,\quad \partial_\xi u,\quad \partial_\rho u\] are computed using the parameter-tangent formulas 19 or, at the discrete level, 26 .

Once the chosen model and discretization satisfy the bias, variance, stability, and moment assumptions stated in Corollaries 1 and 2, the same conditional and nested QA-MLMC complexity bounds apply to Delta, spot Vega, Gamma, Vanna, Volga, full Rho, and the Heston parameter Greeks.

5 Strong-Error Order One Schemes for Pricing and Greek Estimators↩︎

[34] developed a first-order scheme for BDSDEs based on a two-sided Itô–Taylor expansion. However, since the scheme relies on conditional expectations with respect to the forward Brownian motion at each time step, it integrates out the forward randomness and therefore does not directly provide the pathwise strong approximations required by our conditional-on-\(B\) MLMC framework. Inspired by the underlying Itô–Taylor expansion, we develop new discretization schemes and prove conditional strong-error order one convergence. Consequently, the resulting multilevel differences exhibit second-order variance decay, which is essential for retaining the quadratic quantum speedup in the QA-MLMC framework.

The BDSDE representations derived in the previous section provide probabilistic formulas for option prices and Greeks. To enable efficient multilevel Monte Carlo simulation, we construct strong-error order one discretization schemes for direct pricing estimators, first-order Greek estimators, and second-order Greek estimators.

The presentation of this section follows a common pattern for direct pricing estimators, first-order Greek estimators, and second-order Greek estimators.

For each estimator, we first derive a sufficient condition under which a generic discretization operator yields global strong convergence of order one. We then introduce a concrete Forward–Backward Taylor discretization operator and prove that it satisfies the required conditions. As a result, all three estimators admit globally first-order accurate approximations.

5.1 Forward–Backward Taylor Discretization Framework↩︎

A central contribution of our numerical scheme is a family of Forward–Backward Taylor discretizations for the path functionals appearing in the BDSDE representations of prices and Greeks.

The key challenge is that the backward stochastic integrals \[\begin{align} \int Q(r,X_r)\,\mathrm d \overleftarrow{B}_r, \end{align}\] where \(Q\) denotes a generic integrand and may represent different functions in the pricing and Greek representations. Since the integrand depends on the forward diffusion \(X\), which is itself driven by an independent Brownian motion \(W\), pathwise approximations compatible with the conditional-on-\(B\) MLMC framework must retain the dependence on both sources of randomness throughout the discretization.

Our approach is based on Taylor expansions of the stochastic integrands. The resulting discretizations naturally involve the mixed forward–backward iterated integrals \[\begin{align} J^{WB}_{s,h} &= \Big(J_{ij}^{WB}(s,h)\Big)_{1\le i\le \ell,\;1\le j\le d}, \quad J_{ij}^{WB}(s,h) := \int_s^{s+h}(W_r^j-W_s^j)\, \mathrm d \overleftarrow{B}_r^i, \tag{42} \\ J^{BB}_{s,h} &= \Big(J_{ij}^{BB}(s,h)\Big)_{1\le i,j\le \ell}, \quad J_{ij}^{BB}(s,h) := \int_s^{s+h} \bigl(B_{s}^j-B_r^j\bigr)\, \mathrm d \overleftarrow{B}_r^i, \tag{43} \end{align}\] together with the usual Brownian increments. These terms capture the interaction between the forward Brownian motion \(W\) and the backward Brownian motion \(B\), and play a crucial role in achieving strong-error order one.

Throughout this section, we consider a uniform partition \[t_k=t+kh,\quad k=0,\ldots,N,\] with stepsize \(h=(T-t)/N\). We denote the forward and backward Brownian increments by \[\Delta W_{t_k} = W_{t_{k+1}}-W_{t_k}, \quad \Delta\overleftarrow B_{t_k} = B_{t_k}-B_{t_{k+1}},\] respectively.

Unless stated otherwise, all \(L^q\) norms are taken with respect to the joint law of \((W,B)\): \[\|\cdot\|_{L^q} := \|\cdot\|_{L^q_{W,B}}.\] When conditioning on a fixed realization of the backward Brownian motion \(B\), we write the conditional norm as \(\|\cdot\|_{L_W^q}\).

5.2 Direct Pricing Estimator↩︎

Recall that the discretization operators \(\mathcal{S}_X\), \(\mathcal{S}_\Phi\), and \(\mathcal{S}_{\mathrm{int}}\) were introduced in 5 , 6 , and 7 , respectively. The goal of this subsection is to establish the strong approximation properties of these operators and to derive the accumulated stability estimates required for the global strong-error order one analysis.

After defining the strong-error order associated with each individual discretization step, we now turn to the study of the overall strong-error of the full discretized scheme. Before carrying out the detailed analysis, we first introduce the following accumulated stability property.

A detailed proof is provided in Appendix 9, Paragraph 9.1.0.1. Consequently, on a dyadic grid with \(h_\ell=2^{-\ell}\) and \(N_\ell=(T-t)2^\ell\), Lemma 2 implies the corresponding fixed-\(B\) conditional \(L_W^2\) strong-error estimate, up to the logarithmic factor \((1+\ell)^{a/2}\) for any \(a>1\).

For the forward diffusion process \(X\), we employ the Milstein discretization operator \[\begin{align} \mathcal{S}_X^{\mathrm{mil}} (x,s,h;\Delta W_s) = x +b(x)h +\sigma(x)\Delta W_s + \sum_{i=1}^{d}\sum_{j=1}^{d} L_i\sigma_j(x) \int_s^{s+h} (W_r^i-W_s^i)\,d\mathrm W_r^j, \label{eq:SX-Mil} \end{align}\tag{44}\] where \(L_i = \sum_{n=1}^{d} \sigma_{ni}(x)\partial_{x_n}.\)

For the scalar weight process \(\Phi\), we use the exponential-type discretization \[\begin{align} \mathcal{S}_\Phi (\Phi,s,h; \Delta W_s, \Delta\overleftarrow B_s) = \Phi \exp\!\Big( C_{s,h} +\widetilde{c}(s+\tfrac h2)\cdot\Delta W_s +d(s+\tfrac h2)\cdot\Delta\overleftarrow B_s -\tfrac12 Q_{s,h} \Big), \label{eq:Sphi} \end{align}\tag{45}\] where \[C_{s,h} = \frac{h}{6} \bigl( c(s)+4c(s+\tfrac h2)+c(s+h) \bigr),\] and \[Q_{s,h} = \frac{h}{6} \bigl( q(s)+4q(s+\tfrac h2)+q(s+h) \bigr), \quad q(r)=\|\widetilde{c}(r)\|^2-\|d(r)\|^2.\]

These discretizations are standard and achieve strong-error order one for the corresponding stochastic differential equations. The main additional difficulty in the present setting arises from the accumulated payoff term \[Y_{t_k} = \int_t^{t_k} \Phi(t,r) \bigl( F(r,X_r^{t,x}) + d(r)H(r,X_r^{t,x}) \bigr)\,\mathrm dr + \int_t^{t_k} \Phi(t,r)H(r,X_r^{t,x}) \mathrm d \overleftarrow{B}_r,\] whose discretization involves a backward stochastic integral. A naive Euler-type approximation yields only global strong-error order \(1/2\) and is therefore insufficient for the multilevel framework developed in this work. To overcome this difficulty, we introduce a Forward–Backward Taylor discretization operator \(\mathcal{S}_{\mathrm{int}}^{\mathrm{FBT}}\), which incorporates suitable higher-order corrections for the backward stochastic integral.

Consider the integral discretization operator \(\mathcal{S}_{\mathrm{int}}^{\mathrm{FBT}}\) defined by \[\begin{align} \mathcal{S}_{\mathrm{int}}^{\mathrm{FBT}} (\Phi,X,s,h;\Delta W_s,\Delta \overleftarrow B_s) =&\, \Phi h\bigl(F(s,X)+d(s) H(s,X)\bigr) +\Phi H(s,X)\cdot \Delta\overleftarrow B_s \notag\\ &+\Phi \sum_{i=1}^{\ell}\sum_{j=1}^{d} \bigl(H_x(s,X)\sigma(X)+\widetilde{c}(s)H(s,X)\bigr)_{ij} J_{ij}^{WB}(s,h) +\Phi \sum_{i=1}^{\ell}\sum_{j=1}^{\ell} d_j(s)H_i(s,X) J_{ij}^{BB}(s,h), \label{eq:Sint-FBT} \end{align}\tag{46}\] where \(J^{WB}_{s,h}\) and \(J^{BB}_{s,h}\) are defined in (42 )–(43 ).

Note that the first two terms in 46 coincide with the direct discretization approximation, while two additional terms \(J_{ij}^{WB}(s,h)\) and \(J_{ij}^{BB}(s,h)\) appears. These extra terms serve as higher-order corrections to the backward stochastic integral \[\int_s^{s+h}\Phi(t,r)H(r,X_r)\cdot \mathrm d \overleftarrow{B}_r,\] and are obtained from the taylor expansion of the integrand.

The following two properties explain why the integral discretization operator is defined in this way. Detailed proofs are provided in Appendix 9, Paragraphs 9.1.0.2 and 9.1.0.3.

5.3 First-Order Greek Estimators↩︎

The first-order Greek representation involves, in addition to the forward diffusion \(X\) and the weight process \(\Phi\), the Jacobian flow \(J^{t,x}\) associated with the forward SDE and a corresponding accumulated payoff term.

The purpose of this subsection is twofold. We first develop a general strong-error framework showing that the convergence rate of the first-order Greek estimator is determined by the approximation properties of the discretization operators for \(J^{t,x}\) and the accumulated payoff term. We then construct suitable discretization operators satisfying the requirements of the framework, leading to a global strong-error order one approximation for the first-order Greek estimator.

The following proposition establishes a general strong-error estimate for the first-order Greek payoff estimator. It shows that the convergence rate of the payoff approximation is determined by the approximation properties of the underlying discretization operators together with a suitable stability property of \(\mathcal{S}_{\mathrm{int}}^{(i)}\). A detailed proof is provided in Appendix 9, Paragraph 9.2.0.1.

Applying Lemma 2 with \[\Pi_{\ell,k} = P_{t_k^\ell}^{(i)} - P_{\ell,k}^{(i,h_\ell)}\] yields the fixed-\(B\) conditional first-order Greek estimate with logarithmic loss.

We now construct discretization operators satisfying the assumptions of Proposition 6.

For the Jacobian flow \(J^{t,x}\), we use the standard Milstein discretization. Since \[dJ_s^{t,x} = \nabla b(X_s^{t,x})J_s^{t,x}\,ds + \sum_{a=1}^{d} \nabla \sigma_{\cdot a}(X_s^{t,x})J_s^{t,x}\,dW_s^a , \quad J_t^{t,x}=I_d,\] we define \[\begin{align} \mathcal{S}_J^{\mathrm{mil}} (J,X,s,h;\Delta W_s,\Delta\overleftarrow B_s) =J +\nabla b(X)J\,h +\sum_{a=1}^{d} \nabla\sigma_{\cdot a}(X)J\,\Delta W_s^a +\sum_{a=1}^{d}\sum_{b=1}^{d} L_a\!\big(\nabla\sigma_{\cdot b}J\big)(X,J) \int_s^{s+h} (W_r^a-W_s^a)\,\mathrm dW_r^b . \end{align}\] Under the regularity assumptions imposed on \(b\) and \(\sigma\), this discretization has strong-error order one. The remaining task is the construction of a suitable integral discretization operator \(\mathcal{S}_{\mathrm{int}}^{(i),\mathrm{FBT}}\) for the accumulated payoff term in the first-order Greek representation.

As in the direct pricing case, the construction is based on a first-order stochastic Taylor expansion of the corresponding integrands. The resulting forward–backward Taylor correction terms involve both mixed forward–backward iterated integrals and purely backward iterated integrals.

Fix \(1\le i\le d\). Recall that the first-order Greek integral is given by \[\begin{align} Y_{t_k}^{(i)} =&\int_t^{t_k}\Phi(t,s) \nabla_x\!\Big(F(s,X_s^{t,x})+d(s)H(s,X_s^{t,x})\Big)^\top J_s^{t,x}e_i\,\mathrm ds +\int_t^{t_k}\Phi(t,s)\nabla_xH(s,X_s^{t,x})^\top J_s^{t,x}e_i\,\mathrm d \overleftarrow{B}_s. \label{eq:first-order-greek-integral} \end{align}\tag{47}\] Here \(J_s^{t,x}=\nabla_x X_s^{t,x}\) is the Jacobian flow and \(e_i\) is the \(i\)-th unit vector in \(\mathbb{R}^d\).

Consider the first-order Greek integral discretization operator \[\begin{align} &\mathcal{S}_{\mathrm{int}}^{(i),\mathrm{FBT}} \left(\Phi,X,J,s,h;\Delta W_s,\Delta\overleftarrow B_s\right)\notag\\ =& \Phi h\,\nabla_x\!\Big(F(s,X)+d(s)H(s,X)\Big)^\top Je_i +\Phi\left(\nabla_xH(s,X)^\top Je_i\right)\cdot\Delta\overleftarrow B_s +\Phi\sum_{j=1}^{\ell}\sum_{\alpha=1}^{\ell}d_\alpha(s) \left(\nabla_xH_j(s,X)^\top Je_i\right) J_{j\alpha}^{BB}(s,h)\notag\\ &+\Phi\sum_{j=1}^{\ell}\sum_{a=1}^{d} \bigg[\nabla_x^2H_j(s,X)[\sigma_{\cdot a}(X),Je_i] +\nabla_xH_j(s,X)^\top\big(\nabla_x\sigma_{\cdot a}(X)Je_i\big) +\widetilde{c}_a(s)\nabla_xH_j(s,X)^\top Je_i\bigg]J_{ja}^{WB}(s,h). \label{eq:Sint-first-greek-mil} \end{align}\tag{48}\]

The form of the operator is motivated by a first-order expansion of the integrands with respect to the forward diffusion, the Jacobian flow, and the weight process. A detailed derivation is provided in Appendix 9, Paragraph 9.2.0.2.

The following propositions show that the above Forward–Backward Taylor discretization achieves strong-error order one and satisfies the accumulated stability estimate required in Proposition 6. Detailed proofs are provided in Appendix 9, Paragraphs 9.2.0.3 and 9.2.0.4.

5.4 Second-Order Greek Estimators↩︎

The second-order Greek representation involves, in addition to the forward diffusion \(X^{t,x}\) and the Jacobian flow \(J^{t,x}\), the second-order variational process \(K^{(ij),t,x}\) associated with the forward SDE. For fixed \(1\le i,j\le d\), the process \(K^{(ij),t,x}\) satisfies \[\begin{align} dK_s^{(ij),t,x} =\Big[\nabla_x b(X_s^{t,x})K_s^{(ij),t,x} +\nabla_x^2 b(X_s^{t,x}) [J_s^i,J_s^j]\Big]\,ds +\sum_{a=1}^{d} \Big[\nabla_x\sigma_{\cdot a}(X_s^{t,x}) K_s^{(ij),t,x} +\nabla_x^2\sigma_{\cdot a}(X_s^{t,x}) [J_s^i,J_s^j] \Big]\,dW_s^a, \end{align}\] where \(J_s^i=J_s^{t,x}e_i.\)

The purpose of this subsection is twofold. We first develop a general strong-error framework showing that the convergence rate of the second-order Greek estimator is determined by the approximation properties of the discretization operators associated with \(K^{(ij),t,x}\) and the corresponding accumulated payoff term. We then construct suitable discretization operators satisfying the requirements of the framework, leading to a global strong-error order one approximation for the second-order Greek estimator.

The following proposition establishes a general strong-error estimate for the second-order Greek payoff estimator. In particular, if the underlying discretization operators satisfy suitable strong-error and stability properties, then the resulting payoff approximation inherits the same convergence rate. The proof is deferred to Appendix 9, Paragraph 9.3.0.1.

Applying Lemma 2 with \[\Pi_{\ell,k} = P_{t_k^\ell}^{(ij)} - P_{\ell,k}^{(ij,h_\ell)}\] yields the fixed-\(B\) conditional second-order Greek estimate with logarithmic loss.

We now construct discretization operators satisfying the assumptions of Proposition 9.

For the second-order variational process \(K^{(ij),t,x}\), we use the Milstein discretization \[\begin{align} \mathcal{S}_K^{\mathrm{mil}} \left(K,J^i,J^j,X,s,h;\Delta W_s,\Delta\overleftarrow B_s\right) = K &+A_K(X,J^i,J^j,K)\,h +\sum_{a=1}^{d} B_{K,a}(X,J^i,J^j,K)\,\Delta W_s^a\notag\\ &+\sum_{a=1}^{d}\sum_{b=1}^{d} \mathcal{L}_a^{X,J,K}B_{K,b}(X,J^i,J^j,K) \int_s^{s+h} (W_r^a-W_s^a)\,\mathrm dW_r^b . \end{align}\] Here \[\begin{align} A_K(X,J^i,J^j,K) =\nabla_x b(X)K +\nabla_x^2 b(X)[J^i,J^j], \quad B_{K,a}(X,J^i,J^j,K) =\nabla_x\sigma_{\cdot a}(X)K +\nabla_x^2\sigma_{\cdot a}(X)[J^i,J^j]. \end{align}\] The differential operator \(\mathcal{L}_a^{X,J,K}\) is the diffusion vector field of the extended process \((X,J^i,J^j,K)\) associated with the \(a\)-th Brownian component, namely \[\begin{align} \mathcal{L}_a^{X,J,K} =\sigma_{\cdot a}(X)\cdot\nabla_X +\big(\nabla_x\sigma_{\cdot a}(X)J^i\big)\cdot\nabla_{J^i} +\big(\nabla_x\sigma_{\cdot a}(X)J^j\big)\cdot\nabla_{J^j} +\Big(\nabla_x\sigma_{\cdot a}(X)K +\nabla_x^2\sigma_{\cdot a}(X)[J^i,J^j] \Big)\cdot\nabla_K . \end{align}\] Thus the approximation is generated by \[K_{k+1}^{(ij,h)} = \mathcal{S}_K^{\mathrm{mil}} \left( K_k^{(ij,h)},J_k^{(i,h)},J_k^{(j,h)},X_k^{(h)},t_k,h; \Delta W_{t_k},\Delta\overleftarrow B_{t_k} \right), \quad K_0^{(ij,h)}=0 .\] Under the regularity assumptions imposed on the coefficients, this Milstein approximation achieves strong-error order one.

Similarly to the first order case, the integral discretization operator \(\mathcal{S}_{\mathrm{int}}^{(ij),\mathrm{FBT}}\) is obtained from a first order stochastic Taylor expansion of the integrands appearing in the second-order Greek representation. resulting discretization. A detailed derivation is provided in Appendix 9, Paragraph 9.3.0.2.

Fix \(1\le i,j\le d\). For simplicity, we write \[\begin{align} X_s=X_s^{t,x}, \quad J_s=J_s^{t,x}, \quad K_s^{(ij)}=K_s^{(ij),t,x}, \quad \Phi_s=\Phi(t,s), \quad J_s^j=J_s^{t,x}e_j, \end{align}\] For the second-order Greek integral \[\begin{align} Y_{t_k}^{(ij)} =&\int_t^{t_k}\Phi_s \Big(\nabla_x^2(F+dH)(s,X_s)[J_s^j,J_s^i] +\nabla_x(F+dH)(s,X_s)^\top K_s^{(ij)}\Big)\,\mathrm ds\notag\\ &+\int_t^{t_k}\Phi_s\Big( \nabla_x^2H(s,X_s)[J_s^j,J_s^i] +\nabla_xH(s,X_s)^\top K_s^{(ij)}\Big) \mathrm d \overleftarrow{B}_s, \end{align}\] we define the discretization operator by \[\begin{align} &\mathcal{S}_{\mathrm{int}}^{(ij),\mathrm{FBT}} \left(\Phi,X,J,K,s,h; \Delta W_s,\Delta\overleftarrow B_s\right)\notag\\ =&\Phi h \Big(\nabla_x^2(F+dH)(s,X)[J^j,J^i] +\nabla_x(F+dH)(s,X)^\top K\Big) +\Phi\Big(\nabla_x^2H(s,X)[J^j,J^i] +\nabla_xH(s,X)^\top K\Big) \cdot\Delta\overleftarrow B_s\notag\\ &+\Phi\sum_{\nu=1}^{\ell}\sum_{a=1}^{d} \bigg[\nabla_x^3H_\nu(s,X) [\sigma_{\cdot a}(X),J^i,J^j] +\nabla_x^2H_\nu(s,X)[(\nabla_x\sigma_{\cdot a}(X)J)^j,J^i] +\nabla_x^2H_\nu(s,X) [J^j,(\nabla_x\sigma_{\cdot a}(X)J)^i]\notag\\ & +\sigma_{\cdot a}(X)^\top\nabla_x^2H_\nu(s,X)K +\nabla_xH_\nu(s,X)^\top \big(\nabla_x\sigma_{\cdot a}(X)K +\nabla_x^2\sigma_{\cdot a}(X)[J^i,J^j] \big)\notag\\ & +\widetilde{c}_a(s) \Big(\nabla_x^2H_\nu(s,X)[J^j,J^i] +\nabla_xH_\nu(s,X)^\top K\Big)\bigg] J_{\nu a}^{WB}(s,h)\notag\\ &+\Phi\sum_{\nu=1}^{\ell}\sum_{\alpha=1}^{\ell} d_\alpha(s)\Big((J^j)^\top \nabla_x^2H_\nu(s,X)J^i +\nabla_xH_\nu(s,X)^\top K\Big) J_{\nu\alpha}^{BB}(s,h), \label{eq:Sint-second-greek-mil-full} \end{align}\tag{49}\] where \(J^i=Je_i\), \(J^j:=Je_j\), and \((\nabla_x\sigma_{\cdot a}(X)J)^i =\nabla_x\sigma_{\cdot a}(X)J^i.\)

The following two propositions verify the assumptions of Proposition 9. More precisely, we show that \(\mathcal{S}_{\mathrm{int}}^{(ij),\mathrm{FBT}}\) has strong-error order one and satisfies the required accumulated stability estimate. Detailed proofs are provided in Appendix 9, Paragraphs 9.3.0.3 and 9.3.0.4.

6 Experiments on Conditional MLMC↩︎

In this section, we perform classical Monte Carlo experiments to validate the discretization operators developed in Section 5.

More specifically, we numerically investigate the multilevel convergence rates \(\alpha\), \(\beta\), \(\gamma\) defined through \[\begin{align} |\mathbb{E}[P_\ell-P]| =\mathcal{O}(h_\ell^\alpha), \quad \mathrm{Var}(P_\ell-P_{\ell-1}) =\mathcal{O}(h_\ell^\beta), \quad C_\ell=\mathcal{O}(h_\ell^{-\gamma}), \end{align}\] where \(P_\ell\) denotes the level-\(\ell\) approximation and \(C_\ell\) is the average cost of generating one sample on level \(\ell\). According to Theorem 3 and Corollary 1, the proposed numerical scheme is expected to yield \[\begin{align} \alpha=1,\quad\beta=2,\quad\gamma=1. \end{align}\] The purpose of the experiments is therefore to verify these predicted rates for the pricing estimator, the first-order Greek estimator, and the second-order Greek estimator.

We consider the SPDE 4 with coefficients \[\begin{align} c(s)&=-0.30+0.3s, \quad d(s)=0.15+0.2s, \quad \tilde{c}(s)=0.20+0.2s,\\ F(s,x)&=0.05e^{-s}\sin x,\quad H(s,x)=0.08e^{-s}\cos x. \end{align}\] The forward diffusion is given by \[\begin{align} dX_s=b(X_s)\,ds+\sigma(X_s)\,dW_s, \end{align}\] where \[\begin{align} b(x)=\kappa(\mu-x), \quad \sigma(x)=\nu(1+0.5x), \end{align}\] with parameters \[\begin{align} \kappa=1.2,\quad \mu=0,\quad \nu=0.35. \end{align}\] The terminal payoff is chosen as \[\begin{align} G(x)=\sin(x). \end{align}\] Unless otherwise specified, we take \[\begin{align} t=0,\quad T=1, \quad x_0=0.3. \end{align}\]

For each experiment, we fix the maximal refinement level \[\begin{align} L_{\max}=12 \end{align}\] and use \[\begin{align} N=2\times 10^4 \end{align}\] Monte Carlo samples on each level.

For the bias estimate, the exact target \(P\) is approximated by the finest-level approximation \(P_{L_{\max}}\). Thus, we estimate \[\begin{align} \left|\mathbb{E}[P_\ell-P_{L_{\max}}]\right| \end{align}\] for \(\ell<L_{\max}\), and obtain \(\alpha\) from the regression \[\begin{align} \left| \mathbb{E}[P_\ell-P_{L_{\max}}] \right| =\mathcal{O}(h_\ell^\alpha). \end{align}\]

For the variance estimate, we use the coupled level difference \[\begin{align} \Delta_\ell=P_\ell-P_{\ell-1}, \end{align}\] and estimate \[\begin{align} \mathrm{Var}(\Delta_\ell) =\mathcal{O}(h_\ell^\beta). \end{align}\]

Finally, \(C_\ell\) denotes the average computational cost of generating one level-\(\ell\) sample, and \(\gamma\) is estimated from \[\begin{align} C_\ell =\mathcal{O}(h_\ell^{-\gamma}). \end{align}\]

For each realization of the backward Brownian path \(B\), we estimate the exponents \(\alpha\), \(\beta\), and \(\gamma\). The experiment is repeated for ten independent realizations of \(B\), and the reported values correspond to the empirical mean of the resulting estimates.

In particular, the mixed forward–backward iterated integral \[\begin{align} J^{WB}_{s,h} =\int_s^{s+h} (W_r-W_s)\,\mathrm d \overleftarrow{B}_r \end{align}\] is sampled using the Gaussian approximation \[\begin{align} J^{WB}_{s,h} \approx\frac{1}{2}\,\Delta W\,\Delta B +\sqrt{\frac{h\Delta W^2+h \Delta B^2+h^2}{12}}\,Z, \quad Z\sim N(0,1), \end{align}\] where \(Z\) is independent of \(\Delta W\) and \(\Delta B\). This approximation matches the conditional mean and variance of \(J^{WB}_{s,h}\) given \((\Delta W,\Delta B)\).

The empirical rates shown in Figure 5 and Tables 2–4 report the estimated multilevel exponents for the direct pricing estimator, the first-order Greek estimator, and the second-order Greek estimator.

For the direct pricing estimator, the averaged exponents over ten independent realizations of the backward Brownian motion are \[\begin{align} \alpha=1.145,\quad \beta=2.045,\quad \gamma=0.987. \end{align}\] The first-order Greek estimator gives \[\begin{align} \alpha=1.118,\quad \beta=2.000,\quad \gamma=0.998, \end{align}\] while the second-order Greek estimator gives \[\begin{align} \alpha=1.001,\quad \beta=2.038,\quad \gamma=0.996. \end{align}\]

Across all three estimators, the variance and cost exponents are particularly stable with respect to the realization of the backward Brownian motion. Moreover, the estimated exponents remain close to \[\begin{align} (\alpha,\beta,\gamma)=(1,2,1), \end{align}\] which provides strong numerical evidence for the expected multilevel scaling relations \[\begin{align} \alpha=1, \quad \beta=2, \quad \gamma=1. \end{align}\] These results are in excellent agreement with the theoretical predictions obtained from the global strong-error framework and the associated Forward–Backward Taylor schemes developed in Section 5. Consequently, the assumptions required by Theorem 3 and Corollary 1 are numerically supported.

7 Discussion and Future Work↩︎

In this paper, we developed a quantum-accelerated multilevel Monte Carlo (QA-MLMC) framework for stochastic partial differential equations (SPDEs) arising from stochastic-environment financial models. The main idea is to exploit the BDSDE representation of the SPDE solution to reformulate pricing and sensitivity estimation as conditional and nested expectation estimation problems. These expectation problems can then be estimated by QA-MLMC, leading to quadratic quantum speedups for derivative pricing and Greek estimation in stochastic environments.

The proposed framework contains two main components. First, for a fixed realization of the backward Brownian motion, we construct a conditional QA-MLMC estimator for the SPDE quantity of interest. Second, for quantities involving an additional average over the random environment, we develop a nested QA-MLMC estimator. Under the multilevel assumptions on bias, variance, and cost, these estimators achieve quantum sampling complexity of order \[\widetilde{\mathcal{O}}(\epsilon^{-1}),\] for an additive error tolerance \(\epsilon\). We also show that the same algorithmic structure applies not only to prices, but also to first-order and second-order Greeks.

A key numerical ingredient is the Forward–Backward Taylor scheme developed in this work. Unlike discretization methods that eliminate the forward Brownian randomness through conditional expectations at each time step, the Forward–Backward Taylor discretization keeps the joint pathwise dependence on the forward and backward Brownian motions. This feature is essential for constructing coupled level differences in the conditional multilevel estimator. The strong-error order one convergence results for pricing and Greek estimators provide the numerical foundation for the complexity analysis of the proposed quantum algorithms.

There are several natural directions for future work. First, although the present paper focuses on SPDEs arising from stochastic-environment financial models, the conditional and nested estimation structure is more general. It would be interesting to extend the proposed QA-MLMC framework to broader classes of SPDEs and stochastic systems with common noise, random coefficients, or random media, where BDSDE-type representations can be used to connect SPDE solutions with expectation estimation problems.

Second, the Forward–Backward Taylor discretization introduced in this work suggests several further numerical developments. One direction is to construct higher-order or adaptive forward–backward schemes that remain compatible with conditional multilevel couplings. Another is to analyze such schemes under weaker regularity assumptions, more general noise structures, or path-dependent functionals. A better understanding of the relation between strong convergence, level variance, and sample cost would also help optimize the resulting MLMC and QA-MLMC complexity.

Finally, the nested estimator developed here could be extended to more general nested or nonlinear quantities, such as risk measures and multi-layer conditional expectations. These problems arise naturally in financial applications and in stochastic systems with random environments, and provide a promising setting for further applications of quantum computing.

8 Terminal-value SPDE and Its BDSDE Representation↩︎

It is well known that the Feynman-Kac formula provides a classical probabilistic representation for solutions of linear parabolic PDEs [76]. Given a terminal-value problem of the form \[\partial_t u(t,x) + \mathcal{L} u(t,x) - r(t,x)u(t,x) + F(t,x) = 0, \quad u(T,x)=G(x),\] the Feynman–Kac formula represents the solution \(u(t,x)\) as the conditional expectation of a functional of the solution to an associated SDE. More precisely, \(u(t,x)\) can be written as \[u(t,x)=\mathbb{E}\left[\Phi\left(X^{t,x}_{T}\right)\right],\] where \(X^{t,x}\) solves the SDE whose infinitesimal generator coincides with the differential operator \(\mathcal{L}\).

An important extension of the Feynman-Kac framework was achieved through the theory of backward stochastic differential equations (BSDEs). The seminal work of Pardoux and Peng established that BSDEs provide probabilistic representations for semilinear terminal-value PDEs [17], [18]. In particular, [22] discusses applications of BSDEs in mathematical finance, including the pricing of European options.

However, When the evolution equation contains an additional noise term of the form \(G(t,x)\mathrm{d}W_t\), the solution becomes a stochastic process in both time and space, leading to a stochastic partial differential equation. In this setting, neither the classical Feynman–Kac formula nor the BSDE framework is applicable.

Backward doubly stochastic differential equations (BDSDEs)  [25], [64] were introduced precisely to overcome this limitation. By incorporating an additional backward Itô integral with respect to a time-reversed Brownian motion, BDSDEs extend the BSDE framework and provide a rigorous probabilistic representation for terminal-value SPDEs.

In contrast to classical BSDEs, which are associated with deterministic or random-coefficient partial differential equations, BDSDEs involve two sources of randomness: a forward Brownian motion and a backward stochastic integral. This additional backward component allows BDSDEs to faithfully capture the intrinsic randomness of SPDE solutions while preserving the backward-in-time structure imposed by the terminal condition.

Our goal is to give a probabilistic representation for the solution of the quasilinear backward SPDEs: \[\label{eq:SPDE40with32BDSDE41} \begin{cases} \partial_t u(t,x) = \left[\mathcal{L}u(t,x) + f\left(t,x,u(t,x),(\sigma^\top\nabla u)(t,x)\right)\right]\,\mathrm{d}t + h(t,x,u(t,x),(\sigma^\top\nabla u)(t,x))\,\mathrm{d}B_t,\quad t\in[0,T],\\ u(T,x)=G(x). \end{cases}\tag{50}\] There \(u:\mathbb{R_+}\times \mathbb{R}^d\rightarrow \mathbb{R}^k\) and \(\mathcal{L}u=\left( Lu_1,\cdots ,Lu_k \right)^\top\) with \(L=\frac{1}{2}\sum\limits_{i,j=1}^d(\sigma\sigma^\top)_{ij}\frac{\partial^2}{\partial x_i\partial x_j}+\sum\limits_{i=1}^db_i\frac{\partial}{\partial x_i}.\) Alternatively, we can write the SPDE 50 in the integral form \[\begin{align} \label{eq:SPDE40with32BDSDE41-integral} u(t,x)=G(x)+\int_t^T \left[\mathcal{L}u(s,x) + f\left(s,x,u(s,x),(\sigma^\top\nabla u)(s,x)\right)\right]\,\mathrm{d}s + \int_t^T h(s,x,u(s,x),(\sigma^\top\nabla u)(s,x))\mathrm d \overleftarrow{B}_s. \end{align}\tag{51}\] The backward Itô integral with respect to \(\mathrm d \overleftarrow{B}_s\) is defined as \[\begin{align} \int_t^T \phi(s) \mathrm d \overleftarrow{B}_s := \int_0^{T-t} \phi_{T-s}\,\mathrm{d}B_s', \end{align}\] where \(B_s' = B_T - B_{T-s}\) is the time-reversed Brownian motion. Equivalently, for any partition \(t=t_0<\cdots<t_n=T\), \[\begin{align} \int_t^T \phi(s)\mathrm d \overleftarrow{B}_s = \lim_{|\Pi|\to 0}\sum_{i=0}^{n-1} \phi_{t_{i+1}}\left(B_{t_i}-B_{t_{i+1}}\right), \end{align}\] with convergence in \(L^2(\Omega)\).

The connection with BDSDEs is obtained as follows. For each \((t,x)\in\mathbb{R}_+\times \mathbb{R}^d\), let \(\{X_s^{t,x};t\leq s\leq T\}\) be the solution of the SDE: \[\label{eq:SDE40with32BDSDE41} \begin{cases} \mathrm{d}X_s^{t,x}=b(X_s^{t,x})\,\mathrm{d}s +\sigma(X_s^{t,x})\,\mathrm{d}W_s,\quad s\in[t,T],\\ X_t^{t,x}=x. \end{cases}\tag{52}\] Assume that the SPDE 50 has a classical solution. Then the couple \(\left(Y_s^{t,x},Z_s^{t,x}\right)\) where \[\begin{align} Y_s^{t,x}=u(s,X_s^{t,x}),\quad Z_s^{t,x}=(\sigma^\top\nabla u)(s,X_s^{t,x}) \end{align}\] verify the following BDSDE: \[\begin{align} \label{eq:BDSDE} Y_s^{t,x}=G(X_T^{t,x})+\int_s^T f(r,X_r^{t,x},Y_r^{t,x},Z_r^{t,x})\,\mathrm{d}r +\int_s^T h(r,X_r^{t,x},Y_r^{t,x},Z_r^{t,x})\mathrm d \overleftarrow{B}_r-\int_s^T Z_r^{t,x}\,\mathrm{d}W_r,\quad t\leq s\leq T, \end{align}\tag{53}\] or alternatively \[\label{eq:BDSDE-d} \begin{cases} \mathrm{d}Y_s^{t,x}= - f(s,X_s^{t,x},Y_s^{t,x},Z_s^{t,x})\,\mathrm{d}s - h(s,X_s^{t,x},Y_s^{t,x},Z_s^{t,x})\mathrm d \overleftarrow{B}_s + Z_s^{t,x}\,\mathrm{d}W_s, \quad s\in[t,T],\\ Y_T^{t,x}=G(X_T^{t,x}). \end{cases}\tag{54}\]

In the following, based on [25], [64], we provide a more detailed exposition of results related to BDSDEs.

8.0.0.1 Fundamentals of BDSDEs

Let \((\Omega,\mathcal{F},\mathbb{P})\) be a probability space and \(T>0\) be a fixed terminal time. Let \(\{W_{t}, 0\leq t\leq T\}\) and \(\{B_{t},0\leq t\leq T\}\) be two mutually independent standard Brownian motion processes with values in \(\mathbb{R}^{d}\) and in \(\mathbb{R}^{l}\), respectively. For each \(t\in [0,T]\), we define \[\begin{align} \mathcal{F}_{t} := \mathcal{F}_{0,t}^{W} \vee \mathcal{F}_{t,T}^{B}, \end{align}\] where for any process \(\{\eta_{t}\}, \;\mathcal{F}_{s,t}^{\eta}=\sigma\{\eta_{r}-\eta_{s};s\leq r\leq t\}\vee \mathcal{N}\) and \(\mathcal{N}\) is the class of \(\mathbb{P}\)-null sets of \(\mathcal{F}\). Note that the collection \(\{ \mathcal{F}_{t}; t\in[0,T] \}\) is neither increasing nor decreasing, and it does not constitute a filtration.

We recall some notation from Pardoux and Peng [25]. For \(k\in\mathbb{N}\), we denote by \(C^k(\mathbb{R}^p;\mathbb{R}^q)\) the space of \(C^k\) functions from \(\mathbb{R}^p\) to \(\mathbb{R}^q\), by \(C^k_b(\mathbb{R}^p;\mathbb{R}^q)\) the subspace of functions whose partial derivatives up to order \(k\) are bounded, and by \(C^k_p(\mathbb{R}^p;\mathbb{R}^q)\) the space of functions whose partial derivatives up to order \(k\) have at most polynomial growth at infinity.

Let \(M^{2}([0,T]; \mathbb{R}^{n})\) be the set of \(n\)-dimensional jointly measurable stochastic processes \(\{ \varphi_{t}; t\in[0,T] \}\) which satisfy:
(i) \(\|\varphi \|_{M^{2}}^{2} =\mathbb{E}\left[\int_0^T |\varphi_{t}|^{2} \mathrm dt\right]<\infty\);
(ii) \(\varphi_{t}\) is \(\mathcal{F}_{t}\)-measurable for a.e. \(t\in[0,T].\)

Similarly, let \(S^{2}([0,T]; \mathbb{R}^{n})\) be the set of \(n\)-dimensional continuous stochastic processes, which satisfy:
(i) \(\|\varphi \|_{S^{2}}^{2} =\mathbb{E}\left[\sup\limits_{0\leq t\leq T} |\varphi_{t}|^{2}\right]<\infty\);
(ii) \(\varphi_{t}\) is \(\mathcal{F}_{t}\)-measurable for a.e. \(t\in[0,T].\)

The following two theorems establish the connection between solutions of BDSDEs and SPDEs in the general setting.

8.0.0.2 Linear SPDEs and BDSDEs

In our applications of SPDEs, we are mainly interested in the linear case. Therefore, following [64], we present the linear setting that will be used in this work.

In the linear setting, consider \[\begin{align} f(t,x,y,z)=F(t,x)+c(t)y+\widetilde{c}(t)z,\quad h(t,x,y,z)=H(t,x)+d(t)y, \end{align}\] where \(c,\widetilde{c}\) are bounded deterministic functions. Then the BDSDE 53 becomes \[\begin{align} \label{eq:linear32BDSDE} Y_s^{t,x}=G(X_T^{t,x})+\int_s^T \left[F(r,X_r^{t,x})+c(r)Y_r^{t,x}+\widetilde{c}(r)Z_r^{t,x}\right]\,\mathrm{d}r +\int_s^T \left[H(r,X_r^{t,x})+d(r)Y_r^{t,x}\right]\mathrm d \overleftarrow{B}_r- \int_s^T Z_r^{t,x}\,\mathrm{d}W_r,\quad t\leq s\leq T, \end{align}\tag{55}\] and the SPDE 50 becomes \[\label{eq:linear32SPDE40with32BDSDE41} \begin{cases} \mathrm{d} u(t,x) = \Bigl[\mathcal{L}u(t,x) + F(t,x) + c(t)u(t,x) + \widetilde{c}(t)(\sigma^\top\nabla u)(t,x)\Bigr]\,\mathrm{d}t + \Bigl[H(t,x)+d(t)u(t,x)\Bigr]\,\mathrm{d}B_t,\quad t\in[0,T],\\ u(T,x)=G(x). \end{cases}\tag{56}\]

9 Proofs for the Strong-error Framework↩︎

This appendix collects the technical proofs underlying the strong-error analysis developed in Section 5. The corresponding discretization operators are introduced in the main text, and here we establish the local consistency, strong-error, and accumulated stability estimates required for the global strong-error order one convergence results.

9.0.0.1 Notation and conventions.

Throughout this appendix, \(C\) and \(C_q\) denote generic constants, independent of the time step \(h\), whose values may change from line to line. Whenever a uniform grid is used, we write \[t_k=t_0+kh,\quad k=0,\ldots,N,\quad t_N=T,\] with \(h=(T-t_0)/N\). For simplicity, denote \[A(r,x):=F(r,x)+d(r)H(r,x).\]

Unless explicitly stated otherwise, all \(L^p\) norms are taken with respect to the joint law of the Brownian motions \((W,B)\). That is, \[\|\cdot\|_{L^p}:=\|\cdot\|_{L^p_{W,B}} .\] When a fixed realization of the backward Brownian motion \(B\) is considered, we write the conditional norm explicitly as \(L_W^p\).

9.0.0.2 Forward–backward information and martingale estimates.

In the BDSDE setting, the natural information at a grid point \(t_k\) consists of the forward \(W\)-information up to \(t_k\) and the backward \(B\)-information from \(t_k\) to \(T\). At the grid level we write this information as \[\mathcal{G}_k := \sigma(W_{t_0},\ldots,W_{t_k}) \vee \sigma(B_{t_k}-B_{t_m}: k\le m\le N) \vee\mathcal{N} .\] The family \((\mathcal{G}_k)_{k=0}^N\) is not a filtration in the ordinary increasing-time sense: the \(W\)-part is increasing in \(k\), while the \(B\)-part is decreasing in \(k\). Thus, when estimating accumulated local fluctuations below, we do not regard \((\mathcal{G}_k)\) itself as a filtration.

Instead, each martingale-type contribution is treated by splitting it into its forward and backward parts. The \(W\)-terms are estimated as ordinary forward martingale differences and the \(B\)-terms are estimated as backward martingale differences, or equivalently as ordinary martingale differences after introducing the reversed Brownian motion \[\widehat B_u:=B_T-B_{T-u},\quad 0\le u\le T .\] The Burkholder–Davis–Gundy estimates used below are always applied in this split sense, and the resulting forward and reversed-time estimates are combined by the triangle and Minkowski inequalities.

9.0.0.3 A joint-to-conditional estimate.

Most of the strong-error estimates in this appendix are first proved under the joint law of the two Brownian motions \((W,B)\). However, in conditional setting, we need to interpret the estimate after fixing a realization of the backward Brownian motion \(B\), so that the remaining randomness comes only from the forward Brownian motion \(W\). The following elementary lemma provides this passage from joint \(L^p_{W,B}\)-bounds to conditional \(L^p_W\)-bounds for \(\mathbb{P}_B\)-almost every realization of \(B\).

The point of the lemma is that a deterministic joint strong-error estimate on a sequence of dyadic time steps can be converted into an almost-sure-in-\(B\) conditional estimate, at the cost of a harmless logarithmic-type factor. More precisely, if an error family satisfies \[\mathbb{E}_{W,B} \left[ \sup_{0\le k\le N_\ell} |\Pi_{\ell,k}|^p \right] \le C_p2^{-p\ell},\] then for almost every fixed \(B\), the conditional \(W\)-error has the same dyadic decay rate up to the factor \((1+\ell)^a\), with a finite random constant depending on \(B\). This loss is mild and is sufficient for the pathwise-in-\(B\) estimates used below.

9.1 Direct Pricing Estimator↩︎

9.1.0.1 Proof of Proposition 3

9.1.0.2 Proof of Proposition 4

9.1.0.3 Proof of Proposition 5

9.2 First-order Greek Estimators↩︎

9.2.0.1 Proof of Proposition 6

9.2.0.2 Construction of \(\mathcal{S}_{\mathrm{int}}^{(i),\mathrm{FBT}}\)

Fix \(1\le i\le d\). Recall that the first-order Greek integral is given by \[\begin{align} Y_{t_k}^{(i)} =&\int_t^{t_k}\Phi(t,s) \nabla_x\!\Big(F(s,X_s^{t,x})+d(s)H(s,X_s^{t,x})\Big)^\top J_s^{t,x}e_i\,\mathrm ds +\int_t^{t_k}\Phi(t,s)\nabla_xH(s,X_s^{t,x})^\top J_s^{t,x}e_i\,\mathrm d \overleftarrow{B}_s. \label{kdjlsair} \end{align}\tag{65}\] Here \(J_s^{t,x}=\nabla_x X_s^{t,x}\) is the Jacobian flow and \(e_i\) is the \(i\)-th unit vector in \(\mathbb{R}^d\).

For notational simplicity, throughout this subsection we write \[X_s=X_s^{t,x},\quad J_s=J_s^{t,x},\quad \Phi_s=\Phi(t,s).\]

Recall the mixed iterated integral \[\begin{align} J^{WB}_{s,h} = \Big(J_{ij}^{WB}(s,h)\Big)_{1\le i\le \ell,\;1\le j\le d}, \quad J_{ij}^{WB}(s,h) := \int_s^{s+h} (W_r^j-W_s^j)\, \mathrm d \overleftarrow{B}_r^i. \label{eq:def-JWB-first-greek} \end{align}\tag{66}\] Similarly, define the backward iterated integral \[\begin{align} J^{BB}_{s,h} = \Big(J_{ij}^{BB}(s,h)\Big)_{1\le i,j\le \ell}, \quad J_{ij}^{BB}(s,h) := \int_s^{s+h} (B_{s}^j-B_r^j)\, \mathrm d \overleftarrow{B}_r^i. \label{eq:def-JBB-first-greek} \end{align}\tag{67}\]

We define the first-order Greek integral discretization operator by \[\begin{align} &\mathcal{S}_{\mathrm{int}}^{(i),\mathrm{FBT}} \left(\Phi,X,J,s,h;\Delta W_s,\Delta\overleftarrow B_s\right)\notag\\ =& \Phi h\,\nabla_x\!\Big(F(s,X)+d(s)H(s,X)\Big)^\top Je_i +\Phi\left(\nabla_xH(s,X)^\top Je_i\right)\cdot\Delta\overleftarrow B_s +\Phi\sum_{j=1}^{\ell}\sum_{\alpha=1}^{\ell}d_\alpha(s) \left(\nabla_xH_j(s,X)^\top Je_i\right) J_{j\alpha}^{BB}(s,h)\notag\\ &+\Phi\sum_{j=1}^{\ell}\sum_{a=1}^{d} \bigg[\nabla_x^2H_j(s,X)\big[\sigma_{\cdot a}(X),Je_i\big] +\nabla_xH_j(s,X)^\top\big(\nabla_x\sigma_{\cdot a}(X)Je_i\big) +\widetilde{c}_a(s)\nabla_xH_j(s,X)^\top Je_i\bigg]J_{ja}^{WB}(s,h). \label{axogtvsf} \end{align}\tag{68}\]

The construction of this operator is natural, as it is obtained from first-order Taylor expansions of the integrands. For \(r\in[s,s+h]\), the local expansions are \[\begin{align} X_r-X_s =\sigma(X_s)(W_r-W_s)+R^X_{s,r}, \label{eq:X-expansion-first-greek} \end{align}\tag{69}\] and \[\begin{align} J_r-J_s =\sum_{a=1}^d\nabla_x\sigma_{\cdot a}(X_s)J_s (W_r^a-W_s^a)+R^J_{s,r}. \label{eq:J-expansion-first-greek} \end{align}\tag{70}\] Moreover, \[\begin{align} \Phi(t,r) =\Phi(t,s)+\Phi(t,s)\widetilde{c}(s)\cdot(W_r-W_s) +\Phi(t,s)d(s)\cdot(B_{s}-B_r) +R^\Phi_{s,r}. \label{eq:Phi-expansion-first-greek} \end{align}\tag{71}\] Under the standard smoothness and moment assumptions, for every \(q\ge2\), \[\begin{align} \|R^X_{s,r}\|_{L^q} +\|R^J_{s,r}\|_{L^q} +\|R^\Phi_{s,r}\|_{L^q} \le C(r-s). \label{eq:basic-remainder-first-greek} \end{align}\tag{72}\]

For each \(1\le j\le \ell\), expanding with respect to \(X_r\), we obtain \[\begin{align} &\nabla_xH_j(r,X_r)^\top J_se_i =\nabla_xH_j(s,X_s)^\top J_se_i +(X_r-X_s)^\top\nabla_x^2H_j(s,X_s) J_se_i+R^{X,H,i}_{j,s,r}. \label{eq:X-part-taylor-first-greek} \end{align}\tag{73}\] For the Taylor remainder in the spatial expansion of \(\nabla_xH_j(r,X_r)^\top J_se_i\), we have \[\begin{align} |R^{X,H,i}_{j,s,r}| \le C\left(|r-s|+|X_r-X_s|^2\right)|J_se_i|. \label{eq:RXH-pointwise-bound-first-greek} \end{align}\tag{74}\] Hence, using \(\|X_r-X_s\|_{L^{2q}}\le C(r-s)^{1/2}\) and the uniform moment bounds of \(J_s\), we obtain \[\begin{align} \|R^{X,H,i}_{j,s,r}\|_{L^q} \le C(r-s). \label{eq:RXH-bound-first-greek} \end{align}\tag{75}\] Using 69 , we deduce \[\begin{align} &(X_r-X_s)^\top\nabla_x^2H_j(s,X_s)J_se_i =\sum_{a=1}^d\sigma_{\cdot a}(X_s)^\top \nabla_x^2H_j(s,X_s)J_se_i(W_r^a-W_s^a) +\widetilde{R}^{X,H,i}_{j,s,r}. \label{eq:X-H-expanded-no-Psi-first-greek} \end{align}\tag{76}\]

Next, expanding with respect to the Jacobian flow \(J_r\). Using 70 we obtain \[\begin{align} &\nabla_xH_j(s,X_s)^\top (J_r-J_s)e_i =\sum_{a=1}^d\nabla_xH_j(s,X_s)^\top \Big(\nabla_x\sigma_{\cdot a}(X_s)J_se_i\Big) (W_r^a-W_s^a)+R^{J,H,i}_{j,s,r}. \label{eq:J-H-expanded-no-Psi-first-greek} \end{align}\tag{77}\] There \(R^{J,H,i}_{j,s,r} =\nabla_xH_j(s,X_s)^\top R^J_{s,r}e_i,\) and hence \(\|R^{J,H,i}_{j,s,r}\|_{L^q} \le C(r-s).\)

Combining 73 and 77 , we finally obtain \[\begin{align} &\nabla_xH_j(r,X_r)^\top J_re_i\notag\\ =&\nabla_xH_j(s,X_s)^\top J_se_i +\sum_{a=1}^d\bigg[\sigma_{\cdot a}(X_s)^\top \nabla_x^2H_j(s,X_s)J_se_i +\nabla_xH_j(s,X_s)^\top \Big(\nabla_x\sigma_{\cdot a}(X_s)J_se_i\Big)\bigg](W_r^a-W_s^a) +R^{H,i}_{j,s,r}, \label{eq:H-gradient-expansion-no-Psi-first-greek} \end{align}\tag{78}\] where \(\|R^{H,i}_{j,s,r}\|_{L^q}\le C(r-s)\).

Combining 71 and 78 , we obtain \[\begin{align} &\Phi(t,r)\nabla_xH_j(r,X_r)^\top J_re_i\notag\\ =&\Phi(t,s)\nabla_xH_j(s,X_s)^\top J_se_i\notag\\ &+\Phi(t,s)\sum_{a=1}^d \bigg[\sigma_{\cdot a}(X_s)^\top \nabla_x^2H_j(s,X_s)J_se_i +\nabla_xH_j(s,X_s)^\top \Big(\nabla_x\sigma_{\cdot a}(X_s)J_se_i\Big) +\widetilde{c}_a(s) \nabla_xH_j(s,X_s)^\top J_se_i \bigg](W_r^a-W_s^a)\notag\\ &+\Phi(t,s)\sum_{\alpha=1}^{\ell} d_\alpha(s) \left(\nabla_xH_j(s,X_s)^\top J_se_i\right) (B_{s}^{\alpha}-B_r^\alpha) +R^{\Phi H,i}_{j,s,r}. \label{eq:Phi-H-gradient-expansion-first-greek} \end{align}\tag{79}\] There \[\begin{align} |R^{\Phi H,i}_{j,s,r}| \le C|R^\Phi_{s,r}| +C|R^{H,i}_{j,s,r}| +C\Big(|W_r-W_s|+|B_{s}-B_r|\Big)|W_r-W_s|. \label{eq:RPhiH-pointwise-bound-first-greek} \end{align}\tag{80}\] Using \[\|W_r-W_s\|_{L^{2q}}\le C(r-s)^{1/2}, \quad\|B_{s}-B_r\|_{L^{2q}}\le C(r-s)^{1/2}\le Ch^{1/2},\] we obtain \[\begin{align} \|R^{\Phi H,i}_{j,s,r}\|_{L^q} \le C(r-s). \label{eq:RPhiH-bound-sharp-first-greek} \end{align}\tag{81}\] In particular, \[\begin{align} \int_s^{s+h} \|R^{\Phi H,i}_{j,s,r}\|_{L^q}^2\,\mathrm dr \le Ch^3. \label{eq:RPhiH-integrated-bound-first-greek} \end{align}\tag{82}\]

Finally, substituting 79 into the backward stochastic integral gives \[\begin{align} &\int_s^{s+h}\Phi(t,r)\nabla_xH(r,X_r)^\top J_re_i\, \mathrm d \overleftarrow{B}_r\notag\\ =&\Phi(t,s) \left(\nabla_xH(s,X_s)^\top J_se_i\right) \cdot\Delta\overleftarrow B_s\notag\\ &+\Phi(t,s) \sum_{j=1}^{\ell}\sum_{a=1}^{d} \bigg[\sigma_{\cdot a}(X_s)^\top \nabla_x^2H_j(s,X_s)J_se_i +\nabla_xH_j(s,X_s)^\top \big(\nabla_x\sigma_{\cdot a}(X_s)J_se_i\big) +\widetilde{c}_a(s) \nabla_xH_j(s,X_s)^\top J_se_i \bigg] J_{ja}^{WB}(s,h) \notag\\ &+ \Phi(t,s) \sum_{j=1}^{\ell}\sum_{\alpha=1}^{\ell} d_\alpha(s) \left( \nabla_xH_j(s,X_s)^\top J_se_i \right) J_{ j\alpha}^{BB}(s,h) +R^{B,i}_{s,h}. \label{eq:backward-expansion-first-greek} \end{align}\tag{83}\] Here the remainder is given by \[\begin{align} R^{B,i}_{s,h}:=\sum_{j=1}^{\ell} \int_s^{s+h}R^{\Phi H,i}_{j,s,r} \mathrm d \overleftarrow{B}_r^j. \label{eq:RB-definition-first-greek} \end{align}\tag{84}\]

Similarly, for the time integral part, we can expand the integrand as \[\begin{align} &\Phi(t,r) \nabla_x\!\Big(F(r,X_r)+d(r)H(r,X_r)\Big)^\top J_re_i\notag\\ =&\Phi(t,s)\nabla_x\!\Big(F(s,X_s)+d(s)H(s,X_s)\Big)^\top J_se_i\notag\\ &+\Phi(t,s)\sum_{a=1}^{d} \bigg[\sigma_{\cdot a}(X_s)^\top \nabla_x^2\!\Big(F(s,X_s)+d(s)H(s,X_s)\Big)J_se_i +\nabla_x\!\Big(F(s,X_s)+d(s)H(s,X_s)\Big)^\top \big(\nabla_x\sigma_{\cdot a}(X_s)J_se_i\big)\notag\\ & +\widetilde{c}_a(s) \nabla_x\!\Big(F(s,X_s)+d(s)H(s,X_s)\Big)^\top J_se_i\bigg] (W_r^a-W_s^a)\notag\\ &+\Phi(t,s)\sum_{\alpha=1}^{\ell} d_\alpha(s)\nabla_x\!\Big(F(s,X_s)+d(s)H(s,X_s)\Big)^\top J_se_i\, (B_{s}^\alpha-B_r^\alpha)+R^{D,i}_{s,r} \label{eq:drift-integrand-expansion-first-greek} \end{align}\tag{85}\] with \(\|R^{D,i}_{s,r}\|_{L^q} \le C_q(r-s).\)

Integrating 85 over \([s,s+h]\), we obtain \[\begin{align} &\int_s^{s+h}\Phi(t,r) \nabla_x\!\Big(F(r,X_r)+d(r)H(r,X_r)\Big)^\top J_re_i \,\mathrm dr =\Phi(t,s)h \nabla_x\!\Big(F(s,X_s)+d(s)H(s,X_s)\Big)^\top J_se_i +A^{D,i}_{s,h}+M^{D,i}_{s,h}, \label{eq:drift-expansion-first-greek} \end{align}\tag{86}\] where \[\begin{align} A^{D,i}_{s,h}=\int_s^{s+h}R^{D,i}_{s,r}\,\mathrm dr, \label{eq:AD-explicit-first-greek} \end{align}\tag{87}\] and \[\begin{align} M^{D,i}_{s,h} =&\Phi(t,s) \sum_{a=1}^{d} \bigg[ \sigma_{\cdot a}(X_s)^\top \nabla_x^2\!\Big(F(s,X_s)+d(s)H(s,X_s)\Big) J_se_i +\nabla_x\!\Big(F(s,X_s)+d(s)H(s,X_s)\Big)^\top \big(\nabla_x\sigma_{\cdot a}(X_s)J_se_i\big)\notag\\ &+\widetilde{c}_a(s) \nabla_x\!\Big(F(s,X_s)+d(s)H(s,X_s)\Big)^\top J_se_i\bigg] \int_s^{s+h}(W_r^a-W_s^a)\,\mathrm dr\notag\\ &+\Phi(t,s)\sum_{\alpha=1}^{\ell} d_\alpha(s)\nabla_x\!\Big(F(s,X_s)+d(s)H(s,X_s)\Big)^\top J_se_i \int_s^{s+h}(B_{s}^\alpha-B_r^\alpha)\,\mathrm dr. \label{eq:MD-explicit-first-greek} \end{align}\tag{88}\] Here \(A^{D,i}_{s,h}\) is the accumulated Taylor remainder in the time integral and is of finite-variation type. The term \(M^{D,i}_{s,h}\) collects the first-order stochastic fluctuations of \(X\), \(J\), and \(\Phi\) inside the time integral. It is centered conditionally on the information at time \(s\) and is therefore treated as the local martingale-type contribution. We have \[\begin{align} \|A^{D,i}_{s,h}\|_{L^q} \le\int_s^{s+h}\|R^{D,i}_{s,r}\|_{L^q}\,\mathrm dr \le C\int_s^{s+h}(r-s)\,\mathrm dr \le Ch^2. \label{eq:AD-bound-first-greek} \end{align}\tag{89}\] Moreover, since \[\begin{align} \left\|\int_s^{s+h}(W_r^a-W_s^a)\,\mathrm dr\right\|_{L^q} \le Ch^{3/2}, \quad \left\| \int_s^{s+h}(B_{s}^\alpha-B_r^\alpha)\,\mathrm dr \right\|_{L^q} \le Ch^{3/2}, \end{align}\] and the coefficients have uniformly bounded moments, we also have \[\begin{align} \|M^{D,i}_{s,h}\|_{L^q} \le Ch^{3/2}. \label{eq:MD-bound-first-greek} \end{align}\tag{90}\]

Combining 83 and 86 , we have the one-step consistency relation \[\begin{align} &\int_s^{s+h}\Phi(t,r) \nabla_x\!\Big(F(r,X_r)+d(r)H(r,X_r)\Big)^\top J_re_i\,\mathrm dr+ \int_s^{s+h}\Phi(t,r)\nabla_xH(r,X_r)^\top J_re_i\, \mathrm d \overleftarrow{B}_r\notag\\ =&\mathcal{S}_{\mathrm{int}}^{(i),\mathrm{FBT}} \left(\Phi(t,s),X_s,J_s,s,h; \Delta W_s,\Delta\overleftarrow B_s\right) +R^{B,i}_{s,h}+A^{D,i}_{s,h}+M^{D,i}_{s,h}. \label{eq:one-step-consistency-first-greek} \end{align}\tag{91}\]

9.2.0.3 Proof of Proposition 7

9.2.0.4 Proof of Proposition 8

9.3 Second-order Greek Estimators↩︎

9.3.0.1 Proof of Proposition 9

9.3.0.2 Construction of \(\mathcal{S}_{\mathrm{int}}^{(ij),\mathrm{FBT}}\)

For the discretization operator \[\begin{align} &\mathcal{S}_{\mathrm{int}}^{(ij),\mathrm{FBT}} \left(\Phi,X,J,K,s,h; \Delta W_s,\Delta\overleftarrow B_s\right)\notag\\ =&\Phi h \Big((J^j)^\top\nabla_x^2(F+dH)(s,X)J^i +\nabla_x(F+dH)(s,X)^\top K\Big) +\Phi\Big((J^j)^\top\nabla_x^2H(s,X)J^i +\nabla_xH(s,X)^\top K\Big) \cdot\Delta\overleftarrow B_s\notag\\ &+\Phi\sum_{\nu=1}^{\ell}\sum_{a=1}^{d} \bigg[\nabla_x^3H_\nu(s,X) [\sigma_{\cdot a}(X),J^i,J^j] +\big((\nabla_x\sigma_{\cdot a}(X)J)^j\big)^\top \nabla_x^2H_\nu(s,X)J^i +(J^j)^\top\nabla_x^2H_\nu(s,X) (\nabla_x\sigma_{\cdot a}(X)J)^i\notag\\ & +\sigma_{\cdot a}(X)^\top\nabla_x^2H_\nu(s,X)K +\nabla_xH_\nu(s,X)^\top \big(\nabla_x\sigma_{\cdot a}(X)K +(J^i)^\top\nabla_x^2\sigma_{\cdot a}(X)J^j \big)\notag\\ & +\widetilde{c}_a(s) \Big((J^j)^\top\nabla_x^2H_\nu(s,X)J^i +\nabla_xH_\nu(s,X)^\top K\Big)\bigg] J_{\nu a}^{WB}(s,h)\notag\\ &+\Phi\sum_{\nu=1}^{\ell}\sum_{\alpha=1}^{\ell} d_\alpha(s)\Big((J^j)^\top \nabla_x^2H_\nu(s,X)J^i +\nabla_xH_\nu(s,X)^\top K\Big) J_{\nu\alpha}^{BB}(s,h), \label{rtqbisdj} \end{align}\tag{102}\] where \(J^i=Je_i\), \(J^j:=Je_j\), and \((\nabla_x\sigma_{\cdot a}(X)J)^i =\nabla_x\sigma_{\cdot a}(X)J^i,\) the construction is natural, as it is obtained from first-order Taylor expansions of the integrands. For \(r\in[s,s+h]\), we will repeatedly use the following first-order expansions: \[\begin{align} X_r-X_s =&\sum_{a=1}^d\sigma_{\cdot a}(X_s)(W_r^a-W_s^a) +R^X_{s,r}, \tag{103}\\ J_r-J_s =&\sum_{a=1}^d \nabla_x\sigma_{\cdot a}(X_s)J_s (W_r^a-W_s^a) +R^J_{s,r}, \tag{104}\\ K_r^{(ij)}-K_s^{(ij)} =&\sum_{a=1}^d \Big(\nabla_x\sigma_{\cdot a}(X_s)K_s^{(ij)} +\left(J_s^i\right)^\top\nabla_x^2\sigma_{\cdot a}(X_s)J_s^j \Big)(W_r^a-W_s^a) +R^K_{s,r}, \tag{105}\\ \Phi(t,r) =&\Phi(t,s) +\Phi(t,s)\widetilde{c}(s)\cdot(W_r-W_s) +\Phi(t,s)d(s)\cdot(B_{s}-B_r) +R^\Phi_{s,r}. \tag{106} \end{align}\] Moreover, under the standard smoothness and moment assumptions, for every \(q\ge2\), \[\begin{align} \|R^X_{s,r}\|_{L^q} +\|R^J_{s,r}\|_{L^q} +\|R^K_{s,r}\|_{L^q} +\|R^\Phi_{s,r}\|_{L^q} \le C(r-s). \end{align}\]

For each \(1\le \nu\le \ell\), we expand the backward integrand \[\begin{align} (J_r^j)^\top\nabla_x^2H_\nu(r,X_r)J_r^i +\nabla_xH_\nu(r,X_r)^\top K_r^{(ij)} \end{align}\] around \((s,X_s,J_s,K_s^{(ij)})\).

First, by Taylor expansion, \[\begin{align} \nabla_x^2H_\nu(r,X_r) =\nabla_x^2H_\nu(s,X_s) +\nabla_x^3H_\nu(s,X_s)[X_r-X_s] +R^{\nabla^2H_\nu}_{s,r}, \end{align}\] where \(\nabla_x^3H_\nu(s,X_s)[X_r-X_s]\) denotes the derivative of the Hessian matrix in the direction \(X_r-X_s\). Since \[\begin{align} \|X_r-X_s\|_{L^{2q}}\le C(r-s)^{1/2} \end{align}\] and the fourth spatial derivatives of \(H_\nu\) are bounded, the remainder satisfies \[\begin{align} \left\|R^{\nabla^2H_\nu}_{s,r}\right\|_{L^q} \le C(r-s). \end{align}\] Using the first-order expansions of \(X_r\) and \(J_r\), we have \[\begin{align} (J_r^j)^\top\nabla_x^2H_\nu(r,X_r)J_r^i =&(J_s^j)^\top\nabla_x^2H_\nu(s,X_s)J_s^i +(J_s^j)^\top \big(\nabla_x^3H_\nu(s,X_s)[X_r-X_s]\big) J_s^i\notag\\ &+(J_r^j-J_s^j)^\top \nabla_x^2H_\nu(s,X_s)J_s^i +(J_s^j)^\top \nabla_x^2H_\nu(s,X_s)(J_r^i-J_s^i) +R^{(1),H,(ij)}_{\nu,s,r}. \end{align}\] Here the terms containing products such as \((J_r^j-J_s^j)^\top\nabla_x^2H_\nu(s,X_s)(J_r^i-J_s^i),\) as well as the products involving \(R^{\nabla^2H_\nu}_{s,r}\), are absorbed into \(R^{(1),H,(ij)}_{\nu,s,r}\). Substituting the expression for \(X_r-X_s\) given in 103 , we get \[\begin{align} &(J_s^j)^\top \big(\nabla_x^3H_\nu(s,X_s)[X_r-X_s]\big)J_s^i =\sum_{a=1}^{d} \nabla_x^3H_\nu(s,X_s) [\sigma_{\cdot a}(X_s),J_s^i,J_s^j] (W_r^a-W_s^a) +\widetilde{R}^{(1),H,(ij)}_{\nu,s,r}. \end{align}\] Similarly, using 104 , we obtain \[\begin{align} &(J_r^j)^\top\nabla_x^2H_\nu(r,X_r)J_r^i\notag\\ =&(J_s^j)^\top\nabla_x^2H_\nu(s,X_s)J_s^i +\sum_{a=1}^{d} \bigg[\nabla_x^3H_\nu(s,X_s)[\sigma_{\cdot a}(X_s),J_s^i,J_s^j] +\big((\nabla_x\sigma_{\cdot a}(X_s)J_s)^j\big)^\top \nabla_x^2H_\nu(s,X_s)J_s^i\notag\\ & +(J_s^j)^\top\nabla_x^2H_\nu(s,X_s) (\nabla_x\sigma_{\cdot a}(X_s)J_s)^i \bigg](W_r^a-W_s^a) +R^{(2),H,(ij)}_{\nu,s,r}. \label{eq:Hessian-part-expanded-second-greek} \end{align}\tag{107}\] The remainder satisfies \[\begin{align} \|R^{(2),H,(ij)}_{\nu,s,r}\|_{L^q}\le C(r-s). \end{align}\]

Next, consider the term \(\nabla_xH_\nu(r,X_r)^\top K_r^{(ij)}.\) A first-order Taylor expansion gives \[\begin{align} \nabla_xH_\nu(r,X_r)^\top K_r^{(ij)} =&\nabla_xH_\nu(s,X_s)^\top K_s^{(ij)} +(X_r-X_s)^\top\nabla_x^2H_\nu(s,X_s)K_s^{(ij)} +\nabla_xH_\nu(s,X_s)^\top(K_r^{(ij)}-K_s^{(ij)}) +R^{(3),H,(ij)}_{\nu,s,r}. \label{eq:gradient-K-term-Taylor-second-greek} \end{align}\tag{108}\] Using 103 , we get \[\begin{align} (X_r-X_s)^\top\nabla_x^2H_\nu(s,X_s)K_s^{(ij)} =\sum_{a=1}^d \sigma_{\cdot a}(X_s)^\top \nabla_x^2H_\nu(s,X_s)K_s^{(ij)} (W_r^a-W_s^a) +\widetilde{R}^{(3),H,(ij)}_{\nu,s,r}. \label{eq:X-K-term-second-greek} \end{align}\tag{109}\] By the expansion of \(K_r^{(ij)}-K_s^{(ij)}\) given in 105 , \[\begin{align} \nabla_xH_\nu(s,X_s)^\top(K_r^{(ij)}-K_s^{(ij)}) =&\sum_{a=1}^d \nabla_xH_\nu(s,X_s)^\top \Big(\nabla_x\sigma_{\cdot a}(X_s)K_s^{(ij)} +(J_s^i)^\top\nabla_x^2\sigma_{\cdot a}(X_s)J_s^j\Big) (W_r^a-W_s^a) +\widetilde{R}^{(4),H,(ij)}_{\nu,s,r}. \label{eq:K-H-term-second-greek} \end{align}\tag{110}\] Combining 108 , 109 , and 110 , we obtain \[\begin{align} &\nabla_xH_\nu(r,X_r)^\top K_r^{(ij)} =\nabla_xH_\nu(s,X_s)^\top K_s^{(ij)}\notag\\ &+\sum_{a=1}^d \bigg[\sigma_{\cdot a}(X_s)^\top \nabla_x^2H_\nu(s,X_s)K_s^{(ij)} +\nabla_xH_\nu(s,X_s)^\top \Big(\nabla_x\sigma_{\cdot a}(X_s)K_s^{(ij)} +(J_s^i)^\top\nabla_x^2\sigma_{\cdot a}(X_s)J_s^j\Big)\bigg] (W_r^a-W_s^a) +R^{(4),H,(ij)}_{\nu,s,r}. \label{eq:gradient-K-expanded-second-greek} \end{align}\tag{111}\]

Adding 107 and 111 , we finally get \[\begin{align} &(J_r^j)^\top\nabla_x^2H_\nu(r,X_r)J_r^i +\nabla_xH_\nu(r,X_r)^\top K_r^{(ij)}\notag\\ =&(J_s^j)^\top\nabla_x^2H_\nu(s,X_s)J_s^i +\nabla_xH_\nu(s,X_s)^\top K_s^{(ij)}\notag\\ &+\sum_{a=1}^{d}\bigg[\nabla_x^3H_\nu(s,X_s) [\sigma_{\cdot a}(X_s),J_s^i,J_s^j] +\big((\nabla_x\sigma_{\cdot a}(X_s)J_s)^j\big)^\top \nabla_x^2H_\nu(s,X_s)J_s^i +(J_s^j)^\top \nabla_x^2H_\nu(s,X_s) (\nabla_x\sigma_{\cdot a}(X_s)J_s)^i\notag\\ & +\sigma_{\cdot a}(X_s)^\top\nabla_x^2H_\nu(s,X_s)K_s^{(ij)} +\nabla_xH_\nu(s,X_s)^\top \Big(\nabla_x\sigma_{\cdot a}(X_s)K_s^{(ij)} +(J_s^i)^\top\nabla_x^2\sigma_{\cdot a}(X_s)J_s^j \Big)\bigg] (W_r^a-W_s^a)+R^{H,(ij)}_{\nu,s,r}. \label{eq:H-second-integrand-expansion} \end{align}\tag{112}\]

It remains to estimate the remainder. The terms collected in \(R^{H,(ij)}_{\nu,s,r}\) are of the following types: \[\begin{align} &|r-s|, \quad|X_r-X_s|^2, \quad|J_r-J_s|^2, \quad|K_r^{(ij)}-K_s^{(ij)}|^2,\\ &|X_r-X_s|\,|J_r-J_s|, \quad|X_r-X_s|\,|K_r^{(ij)}-K_s^{(ij)}|, \quad|J_r-J_s|\,|K_r^{(ij)}-K_s^{(ij)}|,\\ &|R^X_{s,r}|, \quad|R^J_{s,r}|, \quad|R^K_{s,r}|. \end{align}\] Using \[\begin{align} \|X_r-X_s\|_{L^{2q}} +\|J_r-J_s\|_{L^{2q}} +\|K_r^{(ij)}-K_s^{(ij)}\|_{L^{2q}} \le C(r-s)^{1/2}, \end{align}\] and \[\begin{align} \|R^X_{s,r}\|_{L^q} +\|R^J_{s,r}\|_{L^q} +\|R^K_{s,r}\|_{L^q} \le C(r-s), \end{align}\] we obtain \[\begin{align} \|R^{H,(ij)}_{\nu,s,r}\|_{L^q} \le C(r-s). \end{align}\]

Combining 106 and 112 , we obtain \[\begin{align} &\Phi(t,r) \Big((J_r^j)^\top\nabla_x^2H_\nu(r,X_r)J_r^i +\nabla_xH_\nu(r,X_r)^\top K_r^{(ij)}\Big)\notag\\ =&\Phi(t,s) \Big((J_s^j)^\top\nabla_x^2H_\nu(s,X_s)J_s^i +\nabla_xH_\nu(s,X_s)^\top K_s^{(ij)}\Big)\notag\\ &+\Phi(t,s) \sum_{a=1}^{d} \bigg[\nabla_x^3H_\nu(s,X_s) [\sigma_{\cdot a}(X_s),J_s^i,J_s^j] +\big((\nabla_x\sigma_{\cdot a}(X_s)J_s)^j\big)^\top \nabla_x^2H_\nu(s,X_s)J_s^i +(J_s^j)^\top\nabla_x^2H_\nu(s,X_s) (\nabla_x\sigma_{\cdot a}(X_s)J_s)^i\notag\\ & +\sigma_{\cdot a}(X_s)^\top \nabla_x^2H_\nu(s,X_s)K_s^{(ij)} +\nabla_xH_\nu(s,X_s)^\top \Big(\nabla_x\sigma_{\cdot a}(X_s)K_s^{(ij)} +\nabla_x^2\sigma_{\cdot a}(X_s)[J_s^i,J_s^j]\Big)\notag\\ & +\widetilde{c}_a(s) \Big((J_s^j)^\top\nabla_x^2H_\nu(s,X_s)J_s^i +\nabla_xH_\nu(s,X_s)^\top K_s^{(ij)}\Big)\bigg](W_r^a-W_s^a)\notag\\ &+\Phi(t,s)\sum_{\alpha=1}^{\ell}d_\alpha(s) \Big((J_s^j)^\top\nabla_x^2H_\nu(s,X_s)J_s^i +\nabla_xH_\nu(s,X_s)^\top K_s^{(ij)}\Big) (B_{s}^{\alpha}-B_r^\alpha) +R^{\Phi H,(ij)}_{\nu,s,r}. \label{eq:PhiH-second-expansion} \end{align}\tag{113}\] The product remainder satisfies \[\begin{align} \|R^{\Phi H,(ij)}_{\nu,s,r}\|_{L^q} \le C(r-s). \end{align}\] Consequently, \[\begin{align} \sum_{\nu=1}^{\ell} \int_s^{s+h} \|R^{\Phi H,(ij)}_{\nu,s,r}\|_{L^q}^2 \,\mathrm dr \le Ch^3. \end{align}\]

Substituting 113 into the backward stochastic integral gives \[\begin{align} &\int_s^{s+h}\Phi(t,r) \Big((J_r^j)^\top\nabla_x^2H(r,X_r)J_r^i +\nabla_xH(r,X_r)^\top K_r^{(ij)}\Big) \mathrm d \overleftarrow{B}_r\notag\\ =&\Phi(t,s) \Big((J_s^j)^\top\nabla_x^2H(s,X_s)J_s^i +\nabla_xH(s,X_s)^\top K_s^{(ij)}\Big) \cdot\Delta\overleftarrow B_s\notag\\ &+\Phi(t,s)\sum_{\nu=1}^{\ell}\sum_{a=1}^{d} \bigg[\nabla_x^3H_\nu(s,X_s) [\sigma_{\cdot a}(X_s),J_s^i,J_s^j] +\big((\nabla_x\sigma_{\cdot a}(X_s)J_s)^j\big)^\top \nabla_x^2H_\nu(s,X_s)J_s^i\notag\\ & +(J_s^j)^\top \nabla_x^2H_\nu(s,X_s) (\nabla_x\sigma_{\cdot a}(X_s)J_s)^i +\sigma_{\cdot a}(X_s)^\top \nabla_x^2H_\nu(s,X_s)K_s^{(ij)}\notag\\ & +\nabla_xH_\nu(s,X_s)^\top \Big(\nabla_x\sigma_{\cdot a}(X_s)K_s^{(ij)} +(J_s^i)^\top\nabla_x^2\sigma_{\cdot a}(X_s)J_s^j\Big)\notag\\ & +\widetilde{c}_a(s) \Big((J_s^j)^\top\nabla_x^2H_\nu(s,X_s)J_s^i +\nabla_xH_\nu(s,X_s)^\top K_s^{(ij)}\Big)\bigg] J_{\nu a}^{WB}(s,h)\notag\\ &+\Phi(t,s) \sum_{\nu=1}^{\ell}\sum_{\alpha=1}^{\ell} d_\alpha(s) \Big((J_s^j)^\top\nabla_x^2H_\nu(s,X_s)J_s^i +\nabla_xH_\nu(s,X_s)^\top K_s^{(ij)}\Big) J_{\nu\alpha}^{BB}(s,h) +R^{B,(ij)}_{s,h}, \label{eq:backward-second-expansion} \end{align}\tag{114}\] where \[\begin{align} R^{B,(ij)}_{s,h} =\sum_{\nu=1}^{\ell}\int_s^{s+h} R^{\Phi H,(ij)}_{\nu,s,r} \mathrm d \overleftarrow{B}_r^\nu. \end{align}\]

For the time integral part, set \[\begin{align} \Lambda^{(ij)}(r) :=(J_r^j)^\top\nabla_x^2(F+dH)(r,X_r)J_r^i +\nabla_x(F+dH)(r,X_r)^\top K_r^{(ij)} . \end{align}\] By the same Taylor expansion as for the backward integrand, we have \[\begin{align} \Lambda^{(ij)}(r) =\Lambda^{(ij)}(s) +\sum_{a=1}^d \lambda_a^{(ij)}(s)(W_r^a-W_s^a) +R^{\Lambda,(ij)}_{s,r}, \label{eq:Lambda-expansion-second-greek} \end{align}\tag{115}\] where \[\begin{align} \lambda_a^{(ij)}(s) :=&\nabla_x^3(F+dH)(s,X_s) [\sigma_{\cdot a}(X_s),J_s^i,J_s^j] +\big((\nabla_x\sigma_{\cdot a}(X_s)J_s)^j\big)^\top \nabla_x^2(F+dH)(s,X_s)J_s^i\notag\\ &+(J_s^j)^\top\nabla_x^2(F+dH)(s,X_s) (\nabla_x\sigma_{\cdot a}(X_s)J_s)^i +\sigma_{\cdot a}(X_s)^\top \nabla_x^2(F+dH)(s,X_s)K_s^{(ij)}\notag\\ &+\nabla_x(F+dH)(s,X_s)^\top \Big(\nabla_x\sigma_{\cdot a}(X_s)K_s^{(ij)} +(J_s^i)^\top\nabla_x^2\sigma_{\cdot a}(X_s)J_s^j\Big). \end{align}\] The remainder satisfies \[\begin{align} \|R^{\Lambda,(ij)}_{s,r}\|_{L^q} \le C(r-s). \end{align}\] Combining 115 with 106 , we obtain \[\begin{align} \Phi(t,r)\Lambda^{(ij)}(r) =&\Phi(t,s)\Lambda^{(ij)}(s) +\Phi(t,s)\sum_{a=1}^d \Big(\lambda_a^{(ij)}(s) +\widetilde{c}_a(s)\Lambda^{(ij)}(s)\Big) (W_r^a-W_s^a)\notag\\ &+\Phi(t,s)\sum_{\alpha=1}^{\ell} d_\alpha(s)\Lambda^{(ij)}(s)(B_{s}^\alpha-B_r^\alpha) +R^{D,(ij)}_{s,r}. \label{eq:PhiLambda-expansion-second-greek} \end{align}\tag{116}\] The product remainder satisfies \[\begin{align} \|R^{D,(ij)}_{s,r}\|_{L^q} \le C(r-s)+C(r-s)^{1/2}(s+h-r)^{1/2}. \end{align}\] In particular, \[\begin{align} \int_s^{s+h} \|R^{D,(ij)}_{s,r}\|_{L^q}\,\mathrm dr \le Ch^2. \label{eq:PhiLambda-integrated-remainder-second-greek} \end{align}\tag{117}\] Integrating 116 over \([s,s+h]\) gives \[\begin{align} &\int_s^{s+h}\Phi(t,r) \Big((J_r^j)^\top\nabla_x^2(F+dH)(r,X_r)J_r^i +\nabla_x(F+dH)(r,X_r)^\top K_r^{(ij)}\Big) \,\mathrm dr\notag\\ =&\Phi(t,s)h \Big((J_s^j)^\top\nabla_x^2(F+dH)(s,X_s)J_s^i +\nabla_x(F+dH)(s,X_s)^\top K_s^{(ij)}\Big) +M^{D,(ij)}_{s,h}+A^{D,(ij)}_{s,h}, \end{align}\] where \[\begin{align} A^{D,(ij)}_{s,h} :=\int_s^{s+h} R^{D,(ij)}_{s,r}\,\mathrm dr, \end{align}\] and \[\begin{align} M^{D,(ij)}_{s,h}:=&\Phi(t,s) \sum_{a=1}^d\Big(\lambda_a^{(ij)}(s) +\widetilde{c}_a(s)\Lambda^{(ij)}(s)\Big) \int_s^{s+h}(W_r^a-W_s^a)\,\mathrm dr +\Phi(t,s)\sum_{\alpha=1}^{\ell} d_\alpha(s)\Lambda^{(ij)}(s) \int_s^{s+h}(B_{s}^\alpha-B_r^\alpha)\,\mathrm dr. \end{align}\] Here \(A^{D,(ij)}_{s,h}\) is the finite-variation remainder, while \(M^{D,(ij)}_{s,h}\) collects the first-order Brownian fluctuations of \(X\), \(J\), \(K\), and \(\Phi\) inside the time integral. By 117 , \[\begin{align} \|A^{D,(ij)}_{s,h}\|_{L^q} \le \int_s^{s+h} \|R^{D,(ij)}_{s,r}\|_{L^q}\,\mathrm dr \le Ch^2. \end{align}\] Moreover, using the uniform moment bounds of the coefficients and \[\begin{align} \left\| \int_s^{s+h}(W_r^a-W_s^a)\,\mathrm dr \right\|_{L^q} \le Ch^{3/2}, \quad \left\| \int_s^{s+h}(B_{s}^\alpha-B_r^\alpha)\,\mathrm dr \right\|_{L^q} \le Ch^{3/2}, \end{align}\] we obtain \[\begin{align} \|M^{D,(ij)}_{s,h}\|_{L^q} \le Ch^{3/2}. \end{align}\]

Combining the time integral expansion and the backward stochastic integral expansion, we obtain the one-step consistency relation \[\begin{align} &\int_s^{s+h}\Phi(t,r) \Big((J_r^j)^\top\nabla_x^2(F+dH)(r,X_r)J_r^i +\nabla_x(F+dH)(r,X_r)^\top K_r^{(ij)}\Big) \,\mathrm dr\notag\\ &+\int_s^{s+h}\Phi(t,r) \Big((J_r^j)^\top\nabla_x^2H(r,X_r)J_r^i +\nabla_xH(r,X_r)^\top K_r^{(ij)}\Big) \mathrm d\overleftarrow{B}_r\notag\\ &=\mathcal{S}_{\mathrm{int}}^{(ij),\mathrm{FBT}} \left(\Phi(t,s),X_s,J_s,K_s^{(ij)},s,h; \Delta W_s,\Delta\overleftarrow B_s\right) +A^{D,(ij)}_{s,h} +M^{D,(ij)}_{s,h} +R^{B,(ij)}_{s,h}. \label{eq:one-step-consistency-second-greek} \end{align}\tag{118}\]

9.3.0.3 Proof of Proposition 10

9.3.0.4 Proof of Proposition 11

References↩︎