Quantum-accelerated multilevel Monte Carlo methods for stochastic differential equations in mathematical finance

Inspired by recent progress in quantum algorithms for ordinary and partial differential equations, we study quantum algorithms for stochastic differential equations (SDEs). Firstly we provide a quantum algorithm that gives a quadratic speed-up for multilevel Monte Carlo methods in a general setting. As applications, we apply it to compute expectation values determined by classical solutions of SDEs, with improved dependence on precision. We demonstrate the use of this algorithm in a variety of applications arising in mathematical finance, such as the Black-Scholes and Local Volatility models, and Greeks. We also provide a quantum algorithm based on sublinear binomial sampling for the binomial option pricing model with the same improvement.


Introduction
Differential equations are ubiquitous throughout mathematics, science, and engineering.Specifically, ordinary differential equations (ODEs) and partial differential equations (PDEs) characterize continuous processes of systems arising extensively in many fields, from solid mechanics, fluid dynamics, and electromagnetism to biology [24].Calculations of properties of such deterministic systems, typically require numerical schemes, that is one discretizes the differential equation in order to provide an approximate value of the quantity of interest [3].
For numerous systems arising in statistical physics, molecular dynamics, finance, and other real-world models, the dynamics is captured by a stochastic differential equation (SDE) [42,50].Given a typical SDE, a fundamental computational problem is to provide an expected value of a random variable Y , denoted E[Y ], which is a functional determined by the solution of the SDE.Such a computational problem has been widely studied in mathematical finance, where the quantity Y represents the payoff in option and derivative pricing.It is often computationally expensive to estimate E[Y ], since a scheme that approximates the SDE is necessarily run many times to average over the randomness.In this domain, Monte Carlo (MC) methods are basic tools with a provable complexity analysis.
Monte Carlo simulation, known for its flexibility and generality, refers to performing Monte Carlo methods in simulating SDEs [39,42].This method makes use of randomness to estimate E[Y ] of a SDE as introduced above.In general, Monte Carlo simulation generates k independent approximate samples from Y by performing a chosen scheme, and then outputs an average of those k outputs as an approximate expectation of Y .Assuming the variance of Y is bounded by σ 2 , according to Chebyshev's inequality, it suffices to utilize k = O(σ 2 / 2 ) samples to estimate E[Y ], where is the additive error [39].
Simulating SDEs using Monte Carlo is very useful when the cost of each sample is cheap.A typical example is the Geometric Brownian Motion (GBM) in the Black-Scholes (BS) model [6,7,37], which characterizes stock prices in the financial market.Classical algorithms are designed to directly sample from the solution of GBM, whose explicit form is known, rather than simulating the paths of the GBM itself.Assuming the cost of sampling the solution is ignorable, i.e.O (1), the computational complexity equals the number of samples, O(1/ 2 ).However, a challenge faced by Monte Carlo occurs when creating each sample is costly, so even quadratic dependence on may result in an computational complexity significantly larger than O(1/ 2 ) in practice.More concretely, we consider a general SDE without an explicit solution.To compute an approximate solution by numerical methods, we discretize the SDE on the time interval [0, T ] with the step size h, and perform a scheme of strong order r (which produces a random variable Y that is -close to Y , where = O(h r )) [42].Taking the Milstein scheme with strong order 1 as an example, i.e. r = 1, then T /h = Ω(1/ ) number of iterations is required to produce one sample.Under this circumstance, the computational complexity of classical Monte Carlo simulation is O(1/ 3 ) in total.Generally, when performing a scheme of strong order r, the complexity can be improved to O(1/ 2+1/r ).Usually, it is challenging to perform a scheme with large r, due to the higher smoothness requirement of the SDE, and, in practice, it is harder to implement explicit forms of higher strong order schemes [11,42,43].So the acceleration is moderate in practice.
To reduce such an expensive computational cost, multilevel Monte Carlo (MLMC) methods have attracted very considerable attention recently, and have successfully been applied to simulate SDEs with applications in finance [28,29].Recalling the goal of estimating E[Y ] of a general random variable Y , given a sequence of estimators P 0 , P l , . . ., P L that approximates Y with increasing accuracy and cost, multilevel Monte Carlo aims to estimate E[P L ] by simulating a sum of E[P l − P l−1 ], with different numbers of samples at each level l.The main improvement of multilevel Monte Carlo is from reducing the number of samples when the variance of P l − P l−1 is large.By balancing the sample numbers and variances for different P l − P l−1 and summing them together, multilevel Monte Carlo gives an optimal overall cost of estimating E[P L ] that approximates E[Y ] within the mean-squared error 2 .As described by Giles et al, multilevel Monte Carlo with a scheme of strong order r > 1, is capable of estimating E[Y ] of general SDEs with the overall cost O(1/ 2 ) ([29, Theorem 1]), where O neglects logarithmic factors.Compared to the standard Monte Carlo simulation with the same scheme, it removes a 1/ 1/r factor in the overall complexity.Moreover, this approach does not require the use of higher strong order schemes, and hence avoids the smoothness requirement and the implementation difficulty as well.We note that when the samples are allowed not to be random and independent, but based on certain lattice rules, it could lead to a O(1/ p ) cost with p < 2 under certain conditions [21,31].In our case we assume the samples are chosen randomly and independently, then O(1/ 2 ) is the best known complexity that classical algorithm can achieve.[29] However such a quadratic dependence on 1/ of the overall complexity that multilevel Monte Carlo can achieve is still far from ideal, particularly when applications to mathematical finance are considered.
Quantum computers are expected to outperform classical computers for solving a system of linear equations [1,2,13,32,34,45,58,59] and differential equations [4, 5, 12, 14-16, 18, 23, 46, 48, 52, 60, 61].Quantum algorithms for certain stochastic differential equations, such as simulating GBM of the Black-Scholes model as discussed above, have attracted increasing attention in quantum computational finance [8,26,33,53,55].Reference [33] claims an exponential speedup over classical algorithms to solve the Black-Scholes PDE, but does not include a detailed complexity analysis; it uses a very different approach to that used here, which does not seem to be easily extendible to general SDEs.For the payoff models, using quantum-accelerated Monte Carlo methods [51] based on sampling from the explicit solution of GBM, quantum algorithms can approximate the expected value of the price of portfolio within error in complexity O(1/ ), a quadratically improved dependence on compared to classical Monte Carlo simulation [41,53,55,57].However, previous quantum algorithms, that sample the explicit solution of GBM, cannot be extended to simulate general SDEs with no explicit formula for the solution.Reference [40] presented a practical quantum circuit for simulating the Local Volatility (LV) model, which generalizes the GBM and does not have an explicit solution [20,22].However, reference [40] did not provide a concrete complexity for simulating the LV model.In summary, quantum speedups for payoff models of general SDEs lacking explicit solutions are far from well established.
In this paper, we provide quantum algorithms for approximating classical outputs determined by general SDEs (i.e.ones with no explicit formulas for the solutions) with computational complexity O(1/ ).We will apply these to several payoff models of SDEs which arise in finance, and the outputs are the payoffs in pricing.To achieve such an improvement, we first propose a quantum-accelerated multilevel Monte Carlo (QA-MLMC) method in a general setting, and then apply QA-MLMC for general SDEs.Compared to the classical counterpart, QA-MLMC achieves a quadratic speedup in precision up to a logarithmic factor.The main ingredient to this acceleration is the quantum speedup of the Monte Carlo method [51].Roughly, to approximate E[P l − P l−1 ] in the MLMC approach as we discussed above, we only need to use O(1/ ) samples.We shall prove that this speedup is preserved in the telescoping sum L l=0 E[P l − P l−1 ] to approximate E[P L ].We remark that a somewhat similar idea was used in [51] for the special case of computing a partition function as a telescoping product of terms, each of which is approximated using QA-MC.
Instead of the mean-squared error considered in the classical case, QA-MLMC only returns an approximation of E[P L ] in the sense of the additive error.However, this will not incur any unfair comparison between MLMC and QA-MLMC since the two types of errors are almost equivalent (see Appendix A).Because of the different types of errors, MLMC and QA-MLMC have slight differences in their assumptions.This will become clear below when we apply QA-MLMC to solve practical problems in finance.
As discussed above, to solve a general SDE with no explicit solutions, usually we need a discretization scheme of strong order r.Classically, to apply MLMC, a numerical scheme of strong order r > 1 is usually sufficient to obtain a complexity of O(1/ 2 ).However, in the quantum case, to ensure a quadratic speedup, i.e. to obtain a complexity of O(1/ ), we have to apply a numerical scheme of strong order r > 2. Note that when an order r scheme is used, the Monte Carlo method can solve the SDE with a complexity of O(1/ 2+1/r ).And a direct corollary of [51] shows that in the quantum case this method can be improved to has a complexity of O(1/ 1+1/r ).However, for the reasons discussed above we cannot choose r as large as we want.So QA-MLMC still has advantages over QA-MC both in theory and in practice.

Algorithm Model Result
Classical MC with direct sampling [7,37] Black-Scholes model  As applications, we apply QA-MLMC to solving various payoff models, which satisfy our smoothness requirements and are of great interest in mathematical finance.Examples include the well-known Black-Scholes model that prices a variety of financial derivatives; the Local Volatility model that generalizes the Black-Scholes model by treating volatility as a function of the asset and the time; the Greeks that label the sensitivity of the price of the option for hedge portfolios [37]; and the binomial option pricing model as introduced above.
For the analytically solvable Black-Scholes model, in which QA-MC has been applied to reduce the complexity from O( −2 ) to O( −1 ) [55], we verify QA-MLMC is able to achieve the same quantum speedup.For the rest of the models, we establish the first quantum acceleration to achieve the complexity O( −1 ), by applying QA-MLMC to the Local Volatility model and Greeks.Table 1 compares the performance of our approaches to classical ones for various financial models with respect to the dependence on error tolerance .
Furthermore, we also study the Black-Scholes model with European and digital payoffs in numerical experiments.We provide concrete numerical implementations of the schemes of strong order up to 3 to test those parameters with different payoff functions, and numerical results are in good agreement with our theoretical estimates.This provides evidence that our complexities for MLMC and QA-MLMC are reasonable and sharp.
We also consider the binomial option pricing model (BOPM), also known as the binomial lattice model (BLM), which provides an alternative option pricing model different from the Black-Scholes model, by constructing a binomial tree with the same expectation and variance as the Geometric Brownian Motion [19,37,56].Inspired by the binomial structure, it is natural to simulate BOPM by performing sublinear binomial sampling, a recently developed technique that samples a binomial tree in sublinear time [9,25].This technique has been applied to construct a fast random walk, which is a specific binomial tree, to develop fast classical and quantum algorithms for heat equation [46].Observing that sublinear binomial sampling provides a numerical method to output an estimate of Y with ignorable cost per sample, we can estimate E[Y ] preserving the computational complexity the same as the sampling complexity, neglecting logarithmic factors.Therefore, we propose classical and quantum algorithms for BOPM, with the dependence on precision O(1/ 2 ) and O(1/ ), respectively.Complementing the multilevel Monte Carlo methods, this provides an alternative method for efficiently studying properties of struc-tured stochastic models in mathematical finance.
Overall, our theoretical and numerical results provide promise for potential applications of quantum computing in computational finance.
The rest of this paper is structured as follows.Section 2 introduces the payoff problem of SDEs that we study, and quantum-accelerated Monte Carlo for solving SDEs.Section 3 states the main theorem for the quantum-accelerated multilevel Monte Carlo method.Section 4 applies the quantum-accelerated multilevel Monte Carlo method to the SDE problem.Section 5 presents the Black-Scholes model as an application for estimating the price of the option.Section 6 generalizes the Black-Scholes model to the Local Volatility model.Section 7 introduces the Greeks as an application for estimating the sensitivity of the price.Section 8 covers the binomial lattice model as an alternative option pricing model.Section 9 includes a discussion and raises some open problems.Appendix A discusses the relationship between mean square error and additive error.Finally, Appendix B presents numerical results of several schemes for the Black-Scholes model with European and digital payoffs.
2 Payoff models of general SDEs for t ∈ [0, T ], where X t ∈ R is an Itô process, W t is a standard Brownian motion.The payoff problem we are concerned with is as follows.
Problem 1. Assuming there exists an oracle O I that samples from an initial distribution π 0 to produce X 0 as an initial condition, an oracle O W that samples from the Brownian motion W t , and an oracle O P that produces the payoff P(X) as a functional of a specific X, and oracles O µ and O σ that produce µ(X, t) and σ(X, t) as functions of X and t, respectively.Given an evolution time T > 0, we aim to compute within an error > 0, where X T is generated by (2.1).
Problem 1 is widely investigated in mathematical finance.Different kinds of financial models, such as option pricing and Greeks, can be formulated as Problem 1 with various assumptions, which are introduced in detail in Section 4. Note that if π 0 = δ X 0 , then X ∈ π 0 → X = X 0 , and Problem 1 is reduced to a deterministic initial value problem.
Before we proceed, it is helpful to clarify the meaning of the errors.By saying to compute E[P] within an error > 0 in Problem 1, there are two different scenarios (here the random variable Y is an estimator of E[P]): • the additive error |Y − E[P]| is bounded by with probability at least 0.99.
Most research on classical Monte Carlo simulation uses mean-squared error, whereas research on the quantum-accelerated Monte Carlo method typically uses additive error.To reduce technical difficulty and be consistent with existing literature, we will follow this convention, using mean-squared error in our classical algorithms and using additive error in our quantum algorithms.We note however that this still allows fair comparison between classical and quantum algorithms, because these two types of error bounds are indeed almost equivalent, which is elaborated in detail in Appendix A.

Monte Carlo method
For simplicity we assume that the costs of querying O I , O W , O P , O µ and O σ are O (1).For our quantum algorithms, we assume that the oracles can produce coherent superpositions corresponding to these distributions.This assumption can usually be satisfied, for example given a classical algorithm that generates samples based on uniformly random bits; see [51] for a discussion.In particular, note that we will not actually need to put the initial distribution (which may be quite complicated) into superposition -the algorithm will compute E[P(X T )] for a given starting position X 0 , and sampling X 0 can be performed classically.If a sample of X T generated by (2.1) can be obtained with cost O(1), the computational complexity of solving Problem 1 equals the total number of samples required to estimate E[P(X T )] within .While for general SDEs, the cost of simulating (2.1) to calculate X T each time should be taken into account.
For instance, we consider the widely used Milstein discretization with the step size h, giving where k ∈ [n] 0 , where n = T /h.This has computational cost O(1/h) to produce a sample X n that approximates X T .Based on Monte Carlo methods [39,42], a simple estimation for E[P(X T )] is to generate N samples, each independently outputs P(X T ), and then to produce an average,

.4)
Note that the mean-squared error can be decomposed as which can be bounded by O(N −1 + h 2 ).Here we let V[Y ] denote the variance of a random variable Y .If we seek to bound the mean-squared error by 2 , then we can take This result is first stated in [28].In general, for high order schemes, the complexity of classical Monte Carlo method can be bounded as follows.
Proposition 1.Let Y be an estimator of E[P(X T )] with bounded variance, based on a numerical discretization with time step size h using a numerical scheme with strong order r (defined in (4.2)).Then in order to achieve the accuracy of mean-squared error 2 , the computational complexity of the Monte Carlo method is O( −2−1/r ).
Proof.For simplicity we assume P is globally Lipschitz continuous (we refer to Section 4 for a more general analysis).To bound the mean-squared error, we only need to bound the right hand side of the equation (2.5).Note that and In order to bound the mean-squared error by 2 , it suffices to choose N ∼ −2 and h ∼ 1/r , thus the complexity becomes O(N/h) = O( −2−1/r ).
Proposition 1 tells us that in the classical Monte Carlo method, the total complexity can be indeed improved by using a higher order scheme, but it is bottlenecked by O( −2 ) due to the variance of the estimator.To minimize the overall computational cost without needing higher order schemes, we will introduce an advanced Monte Carlo approach, known as multilevel Monte Carlo [28], in Section 3.
However before doing so, we show that quantum-accelerated Monte Carlo, as described in [51] can give a speed up of non-multilevel methods.

Quantum-accelerated Monte Carlo method
In [51], Montanaro showed that using a quantum computer, the number of samples used in Monte Carlo can be reduced quadratically.
Lemma 1 (Theorem 5 of [51]).Let A be a (classical or quantum) algorithm.Let v(A) be the random variable corresponding to v(x) when the outcome of A is x.Assume that then there is a quantum algorithm that estimates E[v(A)] up to additive error with success probability at least 2/3 by using samples.
The powering lemma stated below can increase the success probability of Lemma 1 to 1 − δ for any arbitrarily small δ.
Lemma 2 (Lemma 1 of [51]).Let A be a (classical or quantum) algorithm which aims to estimate some quantity µ, and whose output μ satisfies |µ − μ| ≤ except with probability γ, for some fixed γ < 1/2.Then, for any δ > 0, it suffices to repeat A O(log 1/δ) times and take the median to obtain an estimate which is accurate to within with probability at least 1 − δ.
In Lemma 1 and Lemma 2, the randomized or quantum algorithm A is used to produce a random variable X and then compute the payoff function P(X) ∈ R. For the detailed implementation of A, we refer to the start of Section 2 of [51] for more details.

Theorem 1.
Let A be an algorithm that generates a sample of numerical solution X n of the SDE using a numerical discretization with time step size h using r-th order scheme in the sense that E| X n −X T | = O(h r ).Assume that P( X n ) has bounded variance independent of h.Then there exists a quantum algorithm that achieves the accuracy of additive error with probability at least 0.99, with computational complexity O( −1−1/r ).
Proof.Similarly to Proposition 1, for technical simplicity we assume P is globally Lipschitz continuous (we refer to Section 4 for analysis on more general payoff functions).By Lemma 1 and Lemma 2, there exists a quantum algorithm that generates an estimator Y such that with probability at least 0.99 using O( −1 ) queries to A. Furthermore, , and the complexity of each query to A becomes O( −1/r ).Combining the above two estimates, in order to bound the additive error |Y − E[P(X T )]| by with probability at least 0.99, the total complexity is O( −1−1/r ).

Multilevel Monte Carlo method
Heinrich [35] developed the first work on multilevel Monte Carlo (MLMC) methods for parametric integration, then Giles [28] introduced MLMC to simulate SDEs.In the following, we first briefly introduce this method.Then we show how to accelerate this method using a quantum-accelerated Monte Carlo method [51].For more about MLMC, we refer to the survey paper [29].
Instead of focusing on SDEs, we review the idea of MLMC in the general setting.Now let P be a random variable, our goal is to estimate E[P ], given a sequence P 0 , P 1 , . . ., P L that approximates P with increasing accuracy, but also increasing cost.For instance, for the SDE (2.1), P is the payoff, P l = P( X n l ) with n l = T /h l = 2 l T .Now we have the following telescoping sum where P −1 = 0. We can estimate E[P L ] by using the Monte Carlo method to approximate each term E[P l − P l−1 ].So we obtain the following approximation of The superindex l means that the samples are generated independently.Note that for SDE, P l − P l−1 comes from two discrete approximations with different timesteps but the same Brownian path.To generate random samples [28] is as follows.First constructing the Brownian increments for the simulation of the discrete path leading to the evaluation of P (l,i) l . Then summing them in groups of size 2 to give the discrete Brownian increments for the evaluation of To determine the cost, the MLMC approach considers the mean-squared error, which can be decomposed as follows For l ≥ 0, let C l , V l be the cost and variance of one sample of P l − P l−1 .
Then the overall cost and variance of Y given in (3.2) is L l=0 N l C l and L l=0 N −1 l V l .To minimize the cost with a fixed variance 2 /2, we introduce the Lagrange multiplier λ 2 and minimize This leads to We restate the rigorous argument in [29] as follows.
Lemma 3 (Theorem 1 of [29]).Let P be a random variable and P l be the corresponding level l numerical approximation.Let Y l be an approximation of E[P l − P l−1 ] based on Monte Carlo method such that the expected cost and variance of one sample is C l and V l respectively.If there exist positive constants α, β, γ such that α then for any < 1/e there exists an L such that Y = L l=0 Y l has a mean-squared error with bound E(Y − E[P ]) 2 ≤ 2 .Moreover, the total computational cost is

Quantum-accelerated multilevel Monte Carlo method
Recall that in MLMC, in the telescoping sum (3.1), the Monte Carlo method is utilized to approximate each mean value E[P l − P l−1 ], in which we can obtain a quadratic speedup using a quantum computer.To approximate the telescoping sum (3.1) via (3.2), we have a classical algorithm to do the sampling, and thus the assumption of the existence of the classical algorithm to do the sampling is satisfied for MLMC in Lemma 1.The proof of the following theorem is similar to that of Lemma 3.
Theorem 2. Let P denote a random variable, and let P l (l = 0, 1, . . ., L) denote a sequence of random variables such that P l approximates P at level l.Further define P −1 = 0. Let C l be the cost of sampling from P l , and let V l be the variance of P l − P l−1 .If there exist positive constants α, β = 2 β, γ such that α ≥ min( β, γ) and then for any < 1/e there is a quantum algorithm that estimates E[P ] up to additive error with probability at least 0.99, and with cost Proof.In MLMC, we use the telescoping sum obtained by the quantumaccelerated Monte Carlo method (see Lemma 1).By Lemma 1 and Lemma 2, for any l ≥ 0, to make sure As a result, the error in approximating E[P ] satisfies the bound so that 2 −αL ≤ /2.This ensures the first error term of (3.9) is bounded by /2.As for the second term, we choose different N l based on the values l , β, γ.Choosing δ = 1/(100(L + 1)), then we can make sure that the success probability is at least (1 − δ) L+1 ≥ e −1/100 > 0.99.We now split into cases to analyse and optimise the complexity of this approach.In the following analysis, for notational convenience we just write N l = 2 − βl / l .But in the end, the cost should be multiplied by (log 1/δ)(log 1/ ) 3/2 (log log 1/ ) for all cases. (a).
. The second error term of (3.9) is bounded by Also we have where the "+1" term is caused by the ceiling function.The cost is bounded by In the above, we used L l=0 2 γl = O(2 γL ) = O( −γ/α ) by our choice (3.10) and the fact that α ≥ γ.

Quantum-accelerated MLMC for solving SDEs
Let us discuss how to apply MLMC to solve Problem 1 with stochastic differential equation (2.1).

Preliminary
Throughout the paper we make the following assumptions on the coefficients of the SDE and the payoff function.
Assumption 1.We assume µ and σ are globally Lipschitz continuous, i.e., there exists a constant L such that We remark that Assumption 1 implies at most linear growth of µ and σ, and there exists a unique strong solution of SDE (2.1) [42].
We say a numerical approximation X k with time step size h = T /n is of strong order r, if for any m ≥ 1, there exists a constant C m such that One class of general high order schemes is the Taylor-Itô scheme of the general form [42] where f α 's are the coefficient functions (depending on µ and σ) and I α are multiple Itô integrals over the time interval [kh, (k + 1)h].For instance, we may consider the Euler-Maruyama scheme (of strong order 1/2) for k ∈ [n] 0 , or the Milstein scheme (of strong order 1) for k ∈ [n] 0 , where ∆W k are i.i.d.normal random variables with expected value zero and variance h.We remark that there exists another kind of general high order schemes called Taylor-Stratonovich schemes [42], which is easier to implement.We will discuss it in Appendix B.

Assumption 2. The coefficient functions f α are globally Lipschitz continuous with respect to x.
There exist a Taylor-Itô scheme and a Taylor-Stratonovich scheme satisfying Assumption 2, which can achieve strong order r = k/2 for all k ≥ 1.We refer to [42, Section 10] for more details.
We are also given an assumption for the final payoffs: Assumption 3. We assume the payoff function is piecewise Lipschitz continuous, i.e., there exist constants In the reminder of this paper, we consider the SDEs, stochastic schemes, and payoffs satisfying Assumption 1, Assumption 2, Assumption 3, respectively.We remark that our results also hold true for high-dimensional systems of SDEs, given that the payoff function is "piecewise Lipschitz continuous" in some sense (e.g.all the discontinuous points are jump discontinuous points and form several separable hyperplanes).For technical simplicity, we will only focus on the analysis of SDE in one dimension in this section.

Method and theory
To solve a SDE problem, we apply the standard multilevel Monte Carlo method and regard the Taylor-Itô scheme as the discretization subroutine [28,29].
At the high level, we estimate the discretized path X k (k ∈ [n] 0 ) for different numbers of iterations n, and perform the quantum oracle to evaluate P(x) for any x.Setting n l = 2 l for l = 1, . . ., L, we apply the quantumaccelerated multilevel Monte Carlo to estimate E[P(X T )].At the lower level, we divide [0, T ] by a uniform partition 0 = t 0 < t 1 < . . .< t n l = T with h = T /n l = T /2 l on the l-level discretization of (2.1), and perform stochastic numerical schemes to approximate X T by X n l .
To estimate the complexity of QA-MLMC, we need to figure out the parameters α, β, γ in Theorem 2.Here o(1) refers to an arbitrarily small real positive number.Furthermore, if the payoff function P is globally Lipschitz continuous, then the estimates on the parameters can be improved to α = r, β = 2r, γ = 1.
Proof.The estimate on γ comes from the construction of the algorithm.At the level l, we use a time step size T /2 l to discretize the path, and compute the expectation at the final time.The dominant computational cost comes from simulating the path, which requires 2 l time steps.At the level l +1, we halve the time step size, and the number of the time steps is doubled.Since we are using the same numerical scheme at each level, the computational cost of propagating a single step remains the same, thus the total computational cost at the level l + 1 is doubled.This indicates that γ = 1.We now focus on the estimate of α and β.We first consider the general payoff function satisfying Assumption 3.
The proof is inspired by [30].For a sample of X T and a sample of numerical approximation X n with time step size h, we define a linear path Λ q} be the number of the discontinuity points along the path.Note that P is at least piecewise Lipschitz continuous, all the discontinuity points are of jump discontinuity.We hereby define the maximum size of the jump to be J, i.e.J = max 1≤j≤q | lim The first part is bounded by O(h r ).The second part can be bounded as, for a large integer m, In the second inequality, the first term O(h α ) follows from that X T is a continuouslydistributed random variable with a bounded density due to the Picard iteration used to establish existence and uniqueness [50] under the global Lipschitz continuous assumption.The second term follows from the Markov inequality.Therefore we have holds for all m, which implies α = rm/(m + 1) = r − o (1).
The estimate of β is similar to that of α.Note that The first part is bounded by O(h 2r ).The second part is bounded by O(h β ) + O(h m(r−β) ) for an arbitrarily large integer m, for the same reason in estimating α.It follows that which implies β = rm/(m + 1) = r − o(1).Finally, if further the payoff function P is globally Lipschitz everywhere, then we have P(M ( X n , X T ) ≥ 1) = 0.It is straightforward to conclude from the previous analysis that α = r, β = 2r.
We remark that the estimates of α and β are possibly not sharp for some of the Taylor-Itô schemes.For example, if the payoff function P is a linear function, then α will be exactly the weak convergence order, which is, for many numerical schemes, larger than the strong convergence order.Nevertheless, our Proposition 2 holds true for more general payoff functions and general high order schemes, and it suffices for QA-MLMC to achieve speedup over classical algorithms.In Appendix B, we perform careful numerical tests of the values of α and β under different smoothness assumptions.We observe that our estimate for β is sharp for both Lipschitz continuous payoff function and discontinuous payoff function, and our estimate for α is sharp for discontinuous payoff functions, while larger α is observed for smoother payoff functions.
Proposition 2 is sufficient to determine the complexity of both classical and quantumaccelerated MLMC for solving SDEs.We start with the classical case.A discussion about certain numerical schemes for typical payoffs has been proposed in Section 5 of [29].For high-order schemes and general payoffs, we state the result as follows.

Proposition 3. Consider Problem 1 for the stochastic differential equation (2.1) under Assumption 1, Assumption 2 and Assumption 3. Then MLMC with a numerical scheme for SDE of strong order r estimates E[P] up to additive error with probability at least
(4.12) Furthermore, if the payoff function P is globally Lipschitz continuous everywhere, then the cost can be improved to Proof.This is a straightforward result from Lemma 3 and Proposition 2.
Proposition 3 tells that, for general SDE (2.1) and payoff function satisfying Assumption 1, Assumption 2, and Assumption 3, it suffices for MLMC to use a numerical scheme of strong order r > 1 to obtain the complexity O( −2 ) , and using a numerical scheme of strong order 1, e.g., Milstein scheme, will lead to the complexity O( −2−o(1) ), i.e. almost quadratic dependence on 1/ .We note that if the payoff function is globally Lipschitz continuous everywhere, then it suffices to use a numerical scheme of strong order 1/2, e.g., Euler-Maruyama scheme, to achieve the complexity O( −2 ).
As a counterpart, we are ready to state our main theorem regarding the complexity of QA-MLMC for solving SDE.

Theorem 3. Consider Problem 1 for the stochastic differential equation (2.1) under Assumption 1, Assumption 2 and Assumption 3. Then QA-MLMC with a numerical scheme for SDE of strong order r estimates E[P] up to additive error with probability at least 0.99 in cost
Furthermore, if the payoff function P is globally Lipschitz continuous everywhere, then the cost can be improved to Proof.This is a straightforward result from Theorem 2 and Proposition 2.
Theorem 3 tells that, for general SDE (2.1) and payoff function satisfying Assumption 1, Assumption 2, and Assumption 3, it suffices for QA-MLMC to use a numerical scheme of strong order r > 2 to obtain the complexity O( −1 ), and using a numerical scheme of strong order 2 will lead to the complexity O( −1−o(1) ), i.e. almost linear dependence on 1/ .We note that if the payoff function is globally Lipschitz continuous everywhere, then it suffices to use a numerical scheme of strong order 1, e.g., Milstein scheme, to achieve the complexity O( −1 ).
Compared with the classical MLMC, the convergence order required to achieve possibly optimal complexity is higher.This is due to the tighter requirement on the parameters α, β and γ.Nevertheless, once the optimal complexity is reached, QA-MLMC achieves a quadratic speedup over the classical MLMC in terms of .

Generalization to the entire path dependence case
So far we have confined ourselves to the case that the payoff function P only depends on the final value X T of the stochastic process.However, the payoff functions of some options of widely practical interest, including Asian call and put options, are in general dependent on the stochastic integral along the entire trajectory of X t for any 0 ≤ t ≤ T .Furthermore, as being discussed later in Section 7, one of the most efficient ways to estimate the sensitivity of the portfolio involves calculating the expectation of a stochastic integral.Fortunately, QA-MLMC can be straightforwardly generalized to this situation, and we will elaborate it in this subsection.
We first allow the payoff function P to be possibly dependent on the stochastic integral of the entire trajectory through the following assumption: Assumption 4. We assume the payoff function has the form where P(X, y, z) is piecewise Lipschitz continuous with respect to X, y, z, and f and g are two globally Lipschitz continuous functions.
Such a payoff function can be transformed to only depend on the final value of another stochastic process by extending the system to higher dimension.Precisely, let us define then Y t and Z t , together with X t , satisfy the system of stochastic differential equations This can be formally written as a system of stochastic differential equations in higher dimension, i.e.
and thus estimating the expectation of P is equivalent to solving Problem 1 in terms of X t and P. By Assumption 4, we can check that X t and P satisfy the Assumption 1, Assumption 2 and Assumption 3. Therefore, according to Proposition 3 and Theorem 3, we have the following result.

Corollary 1. Consider the payoff function satisfying Assumption 4. Then, further under Assumption 1 and Assumption 2, MLMC with a numerical scheme for SDE of strong order r estimates E[P] up to additive error with probability at least 0.99 in cost
Furthermore, if the payoff function P is globally Lipschitz continuous everywhere, then the cost can be improved to Furthermore, if the payoff function P is globally Lipschitz continuous everywhere, then the cost can be improved to We next discuss several applications arising in mathematical finance.

Black-Scholes option pricing model
We consider the Black-Scholes model for option pricing [6,7,37], which is at the core of quantitative finance, as the first application of Problem 1.

Black-Scholes equation
The Black-Scholes model is used to price a variety of financial derivatives via a simple and analytical solvable model using a small number of input parameters.In Black-Scholes model, we are mainly interested in the following Geometric Brownian Motion where S t is the asset price, σ the velocity of the asset and µ is the risk-free interest rate.The rate of return on the asset µ is assumed to be a constant.There are some other assumptions on the model, which Fischer Black and Myron Scholes have thoroughly discussed in [6].For path-independent option, the price of option V (s, t) follows the Black-Scholes equation for s > 0, t ≤ T , with a terminal condition V (s, T ) = ψ(s), where ψ is the final payoff.The link between (5.2) and (5.1) is established by the Feyman-Kac formula.The solution of (5.2) can be represented as the expectation of the solution of (5.1).
Lemma 4. Consider a PDE defined in (5.2) subject to a terminal condition V (s, T ) = ψ(s).The solution V (s, t) can be written as a conditional expectation where S t is an Itô process driven by (5.1).
Using Itô lemma, (5.1) under the condition S t = s can be solved as (5.4) For European options, one can analytically solve for V (s, t), which has a deterministic expression.However, in the case of more complex payoff functions, one needs to resort to the Monte Carlo method, which gives an approximation of V (s, t) by averaging the payoff over samples W t .The quantum version of the Monte Carlo method can be applied to reduce the complexity from O( −2 ) to O( −1 ), and the examples of Black-Scholes model for European options and Asian options are thoroughly discussed in [55].Nevertheless, the analytical solution (5.4) can be used for benchmarking numerical simulations to understand the scaling behavior of discretization schemes.In Appendix B, we present some numerical tests on the Black-Scholes equation to demonstrate our theory results.In our method, we apply discretization schemes to (5.1) to simulate random paths of S t , such as Euler-Maruyama scheme and Milstein scheme (5.6) Here, (5.5) and (5.6) are reduced from (4.4) and (4.5) respectively.Then we estimate E[P (S T )] with P (S t ) = e −µ(T −t) ψ(S T ), where the quantum oracle (4.8) can be reduced to It is worthwhile to mention that our method is suitable for general SDEs and payoff functions.

Option pricing
Black-Scholes model provides a useful tool for option pricing, which was a longstanding problem in finance.Briefly speaking, an option is a contract that allows the holder to buy or sell a financial asset at a fixed price in the future.A call option is an option to buy an asset and a put option is an option to sell it.In this subsection, we would like to present some well-known call options as examples.Put options can be developed in a similar way.We first consider Lipschitz continuous options.One of the famous options is European option, in which the final payoff is given by where K > 0 is the strike price of the option.That means that if S T ≤ K, the option is worthless, and if S T > K, the holder can buy the asset for K dollars and sell it at market price, making a profit of S T − K.Note that the European option is path independent and only relies on the terminal price S T , without the consideration of the whole path.
There also exist path-dependent options relying on the path {S t }.One example is Asian option, in which the final payoff is considered to be with the strike K > 0. The payoff of Asian option is determined by the average of the asset price over [0, T ].
It follows from a theorem that the classical multilevel Monte Carlo method with Euler-Maruyama scheme achieves the complexity O( −2 ) for globally Lipschitz continuous options, which improves the complexity of standard Monte Carlo method O( −3 ).[28,29].
We then consider piecewise Lipschitz continuous options.A typical option is Digital option, in which the digital option (Cash-or-nothing option) gives the final payoff as the form ψ(S T ) = H(S T − K), (5.10) with the strike K > 0, where H is the Heaviside function.Clearly it is a non-Lipschitz continuous option.
For piecewise Lipschitz continuous options, classical multilevel Monte Carlo method with Euler-Maruyama scheme achieves the complexity O( −2.5 ), and it can be further improved to O( −2 ) by Milstein scheme and extreme paths [29].
By Theorem 3 and Corollary 2, the quantum-accelerated multilevel Monte Carlo method can be used to obtain the following result: Corollary 3. We consider Problem 1 given by a Geometric Brownian Motion (5.1) with an initial condition S t = s 0 .We are given the zeroth-order quantum oracle (5.7) for a final piecewise Lipschitz continuous payoff P as defined in (4.6) or (4.16).Then QA-MLMC with a numerical scheme for SDE of strong order r estimates E[P ] up to additive error with probability at least 0.99 in cost (5.11) Furthermore, if the payoff function P is globally Lipschitz continuous everywhere, then the cost can be improved to O( −1 ), r ≥ 1, O( −1/r ), r < 1. (5.12) For globally Lipschitz continuous P , it suffices to perform the Milstein scheme (5.6) with strong order r = 1 to achieve the complexity Õ( −1 ).While for piecewise Lipschitz continuous P , it requires to use a numerical scheme of strong order 5/2 to obtain the complexity Õ( −1 ).

Local Volatility model
The Local Volatility model generalizes the Black-Scholes model (5.1) by treating volatility σ as a function of the asset S t and the time t [20,22].The Local Volatility model is a kind of simplification of the stochastic volatility model, which assumes σ has a randomness of its own.Differing from the analytically solvable Black-Scholes model, we are required to simulate the SDE, as there is a lack of explicit solutions to estimate the price of the Local Volatility model.
The Local Volatility model characterizes S t by a generalized Geometric Brownian Motion as the form dS t = µS t dt + σ(S t , t)S t dW t (6.1) for t ≤ T , with an initial condition S t = s 0 .Here µ is the instantaneous risk-free interest rate, σ(S t , t) is the instantaneous volatility of the risky asset, which is sufficiently smooth with respect to S t and t, and W t is a standard Brownian motion.We further require (6.1) satisfies the assumptions in Section 4. By Feynman-Kac Formula, the SDE (6.1) corresponds to the Black-Scholes equation that describes the price of the option V = V (s, t) by for s > 0, t ≤ T , with a terminal condition V (s, T ) = ψ(s), where ψ is the final payoff.We establish the link between (6.1) and (6.2) by Feynman-Kac Formula.
Lemma 5. Consider a PDE defined in (6.2) subject to a terminal condition V (s, T ) = ψ(s).The solution V (s, t) can be written as a conditional expectation where S t is an Itô process driven by (6.1).
Thus, we can estimate E[P (S T )] with P (S t ) = e −µ(T −t) ψ(S T ) to obtain V (s, t) for general options.
We consider the same setting for payoffs in Section 5.By Theorem 3 and Corollary 2, the quantum-accelerated multilevel Monte Carlo method can be used to obtain the following results: Corollary 4. We consider Problem 1 given by a SDE (6.1) with an initial condition S t = s 0 .We are given the zeroth-order quantum oracle (5.7) for a final piecewise Lipschitz continuous payoff P as defined in (4.6) or (4.16).Then QA-MLMC with a numerical scheme for SDE of strong order r estimates E[P ] up to additive error with probability at least 0.99 in cost Furthermore, if the payoff function P is globally Lipschitz continuous everywhere, then the cost can be improved to (6.5)

Sensitivity of price
Besides the price, another important factor of financial products is the risk.If two financial products have similar expected return, the one with lower risk is usually preferable because people would like to ensure a positive rate of return earnings under most situations and are not willing to take a risk.One common way to reduce the risk is hedging, that is, forming a portfolio by buying several financial products with different reaction and sensitivity to the change of the market.For example, if there are two stocks, and one is positively correlated with the global economy and the other one is negatively correlated, then one may purchase both of them with a specific ratio to guarantee the return no matter whether the global economy flourishes or declines.Correctly identifying and measuring the sensitivity is crucial to hedging portfolios."The Greeks" label the sensitivity of the price of the option [37].They are partial derivatives with respect to the parameters such as the initial price, the starting time, the risk-free rate and the volatility.Computing these is essential to hedge portfolios, and therefore it is even more important than pricing the option itself.Different from the Black-Scholes model, there's no closed form for Greeks with general payoff functions and coefficients in SDEs.Let u(s, t) denote E(P(X T )) where X T is the solution of SDE with starting time t and initial condition X t = s.Different types of Greeks are as follow: • Gamma: • Vega: ∂u(s, t) ∂σ , ( • Theta: ∂u(s, t) ∂t , ( • Rho: ∂u(s, t) ∂r . (7.6) Numerical treatments of Greeks have been widely studied in e.g.[10,27,47,49,54].A natural way to compute the Greeks is finite difference, that is, to compute the expectation of the payoff function with different parameters and use finite difference to approximate the derivatives.Although finite difference approximation of derivatives suffers from numerical instability, using common stochastic paths for both estimators in finite difference can reduce the classical complexity.However, if the payoff function is not smooth enough, finite difference might not be effective.Here we consider an alternative approach using Malliavin calculus [27] to compute the Greeks.Although rigorous derivation of using Malliavin calculus to compute Greeks is technically involved, the basic idea is simply using an analog of calculus of variations and integration by parts formula in the stochastic calculus setting, and the outcome formula to compute Greeks is amazingly concise thus easy to implement.We will show how to combine it with QA-MLMC to achieve quadratic speedup.
For expository purpose let us consider the example of computing Delta in one dimension and omit a few technical details.We refer interested readers to [27] for elaborations and other types of the Greeks.Assume the process X t is given by (2.1), and for simplicity we assume that the system is autonomous, i.e. there is no explict time dependence in µ and σ.Define another stochastic process Y t to be Under Assumption 1, Assumption 2, Assumption 3, and several further technical assumptions that µ, σ are C 1 functions, σ is uniformly bounded away from 0 and the payoff function P has uniformly bounded second moment, [27, Proposition 3.2] tells that Delta can be represented via the following formula: Notice that here we only assume the payoff function P to be piecewise Lipschitz continuous.Therefore, we can compute Delta by sampling (X t , Y t ) according to (2.1) and (7.7), and then use QA-MLMC to compute the expectation (7.8), which is in the form of Assumption 4. By Corollary 2, we have the following complexity estimate: Corollary 5. Assume we are given the zeroth-order quantum oracle (5.7) for a final piecewise Lipschitz continuous payoff P as defined in (4.16).Then by sampling (X t , Y t ) according to (2.1) and (7.7) with a numerical scheme of strong order r and estimating Delta via (7.8) using QA-MLMC, the Delta can be approximated up to an additive error with probability at least 0.99 in cost (7.9) Furthermore, if the payoff function P is globally Lipschitz continuous, then the cost can be improved to (7.10)

Binomial option pricing model
The binomial option pricing model, a.k.a. the binomial lattice model, provides a discrete time (lattice-based) approximation to the Black-Scholes model [19,37,56].The bino-mial model assumes that movements in the price follow a binomial distribution, which approaches the log-norm distribution of the Geometric Brownian Motion.Given a risk-neutral measure and an initial condition S 0 , we divide [0, T ] by a uniform partition 0 = t 0 < t 1 < . . .< t n = T with h = T /n, and model the price of a stock in discrete time by a Markov chain where {Y k } k ∈ [n] 0 are i.i.d. with a two points "up" U and "down" D distribution Under this setting, the probability of the value S N is given by Thus, we create a binomial tree with distribution B(n, p) that describes the prices over time.
We require the conditional expectation of (8.1) matches the Geometric Brownian Motion (5.1), giving which is Practially, we could use E[S k+1 | S k ] ≈ S k (1 + rh) instead, which in fact corresponds with the expectation of the Euler-Maruyama scheme (5.5).Similarly, we require the conditional variance of (8.1) matches the geometry Brownian motion (5.1), .6) There are various binomial lattice models satisfying (8.4) and (8.6).The first and the most famous model is the CRR model [19].If we choose step size h ≤ σ 2 /r 2 , and assume Alternatively, we can also consider an equal probabilities model, the JR model [38].If we choose p = 1/2 to determine U and D, we would obtain Similar as before, we can use first-order approximations of exponential factors as values of U and D in practice.Now we perform BOPM to estimate expectation of option pricing E[P (S T )] The error of approximating general Black-Scholes option prices by CCR model (8.7) is O(1/ √ n) [36], while it can be further improved to O(1/n) for European option [44].
In general, BOPM works for various options with complexity O(n), since we need to calculate S n = U k D n−k S 0 .Even if U = 1/D, we should still multiphy U in total 2k − n times.In general, we can apply multilevel MC to reduce the cost of BOPM.Since we are required to choose n = O(1/ 2 ), its complexity is generally no better than Euler-Maruyama scheme.
But for the piecewise constant payoff, such as the digital option (5.10), it allows us to develop classical and quantum algorithms that only require O(log N ) time complexity.Given S 0 < • • • < S m , we define a piecewise constant payoff by where ψ j , ψ L , ψ R are constant.Noting that (8.9) is a special case of Assumption 3. We are also given a zeroth-order classical oracle which outputs the final payoff (8.9).The procedure is described as follows: given S 0 , U , D, N , we can determine the criteria k * j corresponding with S j , by When we input n, we just compare n with n * j , and determine the payoff (8.9).Under this setting, the complexity of (8.10) is independent on n, by avoiding n times multiplication for calculating We follow Theorem 12 of [46], in which we replace the fast random walk by BOPM.
Proposition 4. We consider Problem 1 given by a Geometric Brownian Motion (5.1) with an initial condition S 0 .We are given the zeroth-order classical oracle (8.10) for a final piecewise constant payoff as defined in (8.9), which has a bounded variance independent of h.There exists a classical algorithm that estimates E[P ]up to additive error with probability at least 0.99 in cost O( −2 ).(8.12) Proof.We aim to perform n = O(1/ 2 ) steps BOPM (8.1) to obtain S T within error , and then calculate P (S T ) as one sample of Monte Carlo simulation.Given U , D, S 0 sampled from π 0 in O(1), and the oracle (8.10) for (8.9), we do the sampling procedure as follows.
Inspired by Lemma 11 and Theorem 12 of [46], we determine n by sampling a binomial distribution B(n, p), and then compare to n * j to determine the payoff by (8.10).The sampling from a binomial distribution requires O(log n) expected samples and O(log n) expected time [9,25].Thus, the cost of each iteration is O(log n).
By Monte Carlo simulation, we repeat the above sampling procedure O( −2 ) times.Thus, we can estimate E[P ] with the final complexity O( −2 log n) = O( −2 log 1/ ) in time.
Similarly, we are given a zeroth-order quantum oracle by modifying (5.7) to be where ψ is the final payoff (8.9).The same as the procedure of (8.10), the complexity of (8.13) is independent on N .We follow Theorem 22 of [46], in which we replace the quantum walk by BOPM.
Theorem 4. We consider Problem 1 given by a Geometric Brownian Motion (5.1) with an initial condition S 0 .We are given the zeroth-order quantum oracle (8.13) for a final piecewise constant payoff as defined in (8.9), which has a bounded variance independent of h.There exists a quantum algorithm that estimates E[P ] up to additive error with probability at least 0.99 in cost O( −1 ).(8.14) Proof.Inspired by Theorem 22 of [46], we apply QA-MC to the random seed used as input to a procedure for sampling from the binomial distributions.According to Proposition 4, the cost of each iteration is O(log n) with n = O(1/ 2 ), and amplitude estimation gives the final complexity O( −1 (log 1/ ) 3/2 (log log 1/ ) 2 ) in time.

Discussion
We have presented quantum-accelerated multilevel Monte Carlo methods for stochastic processes.We apply our algorithm to several applications arising in mathematical finance, in which we classify different financial models corresponding with different payoffs in detail.
We have shown a quadratically improved dependence on precision can be achieved by our algorithm.This work raises several natural open problems.First, from the PDE perspective, we only deal with parabolic PDEs as an application of simulating SDEs.However, Poisson's equation and elliptic PDEs can be solved by classical multilevel Monte Carlo methods [17].Can we apply our algorithm to Poisson's equation, elliptic PDEs, or more general PDEs?
Second, we consider several types of financial models that can be presented as a SDE model Problem 1.However, we only consider time-independent payoffs P(X T ) that only rely on X T .For more general time-dependent payoffs P(X t ) where X t is the stochastic path in time t, such as the Lookback option [29], can we achieve such a quadratic speed-up for this general model?Furthermore, there are some financial models that can be solved by variational inequalities, such as American option.Can we develop corresponding quantum algorithms for such a generalization?Finally, we aim to output a classical value for estimating the mean of a payoff (Problem 1).Can we provide other meaningful characteristics of stochastic processes, or some processes beyond the SDE modelling by quantum computer?And can we find more practical quantum input-output models for potential applications in finance or other fields?

A Different types of errors
Let a random variable Y be an estimator of some unknown quantity a.In this paper we have considered two different types of errors, mean-squared error E[(Y − a) 2 ] and additive error |Y −a|.Here we will show that these two types of errors are indeed almost equivalent.More precisely, we will show that, up to some absolute pre-constants in the errors and logarithmic factors in the cost, the situation that mean-squared error is on the level of 2 indicates that the additive error is on the level of with probability at least 0.99, and vice versa.

Proposition 5.
Let A be a (classical or quantum) algorithm that generates a random variable Y to estimate some unknown quantity a.Assuming that E[|Y | 2+δ ] < ∞ for some δ > 0, we have: Together with Chebyshev's inequality, we have This indicates that a single sample of Y can estimate a up to additive error 3 with probability at least 3/4.Then by Lemma 2, it suffices to repeat A a constant number of times to boost the success probability to 0.99.This completes the proof of the first part.
is bounded by 2 2 .This completes the proof of the second part.
That is also why we prefer high order numerical schemes derived from Taylor-Stratonovich expansion, such as the strong order 2 scheme and strong order 3 scheme in the following.Since we consider Black-Scholes model, the strong order 2 scheme can be simplified to where μ = µ − 1 2 σ 2 and we use a nice property of multiple Stratonovich integral that The strong order 3 schemes for general SDEs are thoroughly discussed in [43].We choose the strong order 3.0 scheme based on Taylor-Stratonovich expansion and the corresponding discretization scheme for (B.8) is as follow.
where J (1,1,1,1,1) and J (1,1,1,1,1,1) are multiple Stratonovich integrals with As for options, we consider European option and Digital option pricing.For European option, the payoff function is For Digital option, the payoff function is where H is Heaviside step function.The numerical scalings of α and β are estimated using linear regression based on the last four points, in which the time step size h is small enough such that the corresponding numerical scheme has already well converged (as shown in the figures) and the numerical errors are dominated by the leading order term.Table 2 and Table 3 shows the estimate of β and α of these schemes respectively.In the case of European option, the estimates of β for five schemes all approximately equal twice the corresponding strong order, which implies that our theoretical estimate for β in Proposition 2 is sharp for Lipschitz continuous functions.The estimates of α with higher order schemes are approximately equal to the corresponding strong order, which also agrees well with our theoretical estimate for α.In the case of Digital option, other than the strong order 3 scheme (which will be discussed later), the estimates of β are roughly equal to the corresponding strong order of the schemes, which again verifies our theoretical results for non-Lipschitz payoff functions.We remark that, similarly to the European option case, some of the estimates for α are larger than the strong order.This is because our theoretical estimate for α is quite conservative that we only employ properties of strong order, while the convergence of payoff function is more related to the weak order, which describes the convergence order of the moments of the stochastic process.This also implies that our theoretical estimate for α might be improvable for particular schemes and payoff functions.
We notice that for Digital option, strong order 3 scheme cannot output the estimate of α and β.That is because strong order 3 scheme is of pretty high accuracy.When we compute the option price, since the expectation of S T under our parameter choices is larger than K, the Heaviside step function H outputs 1 with extremely high probability and the option price is nearly deterministic.The outcome of numerical experiment has such a small variance that the influence of the increase in level l on the mean and variance is hard to distinguish, which explains the results shown in Figure 5.
So we make some changes to the choice of parameters.We increase σ to 1.5, which increase the randomness, and choose K to be exactly the theoretical expectation 100e −0.05 .The numerical results are shown in Figure 6.The estimate of β is 2.957982, which agrees well with theoretical result, and the estimate of α is 2.438328, which might be estimated more accurately by using more samplings as well as smaller time step sizes.

− 2 MC 2 MC
with scheme of strong order r (Proposition 1) payoff models of general SDEs −2−1/r MLMC with scheme of strong order r > 1 (Proposition 3) payoff models of general SDEs −with scheme of strong order r (Theorem 1) payoff models of general SDEs −1−1/r QA-MLMC with scheme of strong order r > 2 (Theorem 3) payoff models of general SDEs −1 QA-MC with binomial sampling (Theorem 4) binomial option pricing model −1

1 Figure 1 : 1 Figure 2 :
Figure 1: estimates of β and α when using Euler Maruyama Scheme for European option (left) and Digital option (right).The top two plots show the variance versus level l, and the slope gives an estimate of β.The bottom two plots show the mean versus level l and the slope gives an estimate of α.

1 Figure 3 : 1 Figure 4 :
Figure 3: estimates of β and α when using Strong order 1.5 Scheme for European option (left) and Digital option (right).The top two plots show the variance versus level l, and the slope gives an estimate of β.The bottom two plots show the mean versus level l and the slope gives an estimate of α.

1 Figure 5 :
Figure 5: estimates of β and α when using Strong order 3 Scheme for European option (left) and Digital option (right).The top two plots show the variance versus level l, and the slope gives an estimate of β.The bottom two plots show the mean versus level l and the slope gives an estimate of α.

Figure 6 :
Figure 6: estimates of β and alpha when using Strong order 3 Scheme for Digital option with new parameters σ = 1.5 and K = 100e −0.05 .The left plot shows the variance versus level l, and the slope gives an estimate of β.The right plot shows the mean versus level l and the slope gives an estimate of α.

Table 1 :
Summary of the time complexities of classical and quantum algorithms for financial models with the additive error , in which logarithmic factors are omitted.
then there exists an algorithm which repeats A a constant number of times and outputs Y such that | Y − a| ≤ 3 with probability at least 0.99.2. If |Y − a| ≤ holds with probability at least 0.99, then there exists an algorithm which repeats A O(log(1/ )) times and outputs Y such that E[( Y − a) 2 ] ≤ 2 2 .Proof. 1.By E[(Y − a) 2 ] ≤ 2 and the bias-variance decomposition

Table 2 ,
3 and Figure1, 2, 3, 4, 5 show our numerical results on estimating α and β as well as the convergence performance of the errors at each layer of MLMC.The five schemes mentioned above are tested by running respectivly 10 6 , 10 6 , 10 7 , 10 8 , 10 9 independent simulations with the following choice of parameters:

Table 2 :
numerical estimates of β based on linear regression for five schemes

Table 3 :
numerical estimates of α based on linear regression for five schemes