Bermudan option pricing by quantum amplitude estimation and Chebyshev interpolation

Pricing of financial derivatives, in particular early exercisable options such as Bermudan options, is an important but heavy numerical task in financial institutions, and its speed-up will provide a large business impact. Recently, applications of quantum computing to financial problems have been started to be investigated. In this paper, we first propose a quantum algorithm for Bermudan option pricing. This method performs the approximation of the continuation value, which is a crucial part of Bermudan option pricing, by Chebyshev interpolation, using the values at interpolation nodes estimated by quantum amplitude estimation. In this method, the number of calls to the oracle to generate underlying asset price paths scales as O˜(ϵ−1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\widetilde{O}(\epsilon ^{-1})$\end{document}, where ϵ is the error tolerance of the option price. This means the quadratic speed-up compared with classical Monte Carlo-based methods such as least-squares Monte Carlo, in which the oracle call number is O˜(ϵ−2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\widetilde{O}(\epsilon ^{-2})$\end{document}.

In this paper, we focus on Bermudan option pricing and consider how to speed it up by quantum algorithms.Let us briefly describe the problem.An option is a financial contract between two parties, the option buyer and seller, which conveys the option buyer the right to buy some underlying assets such as stocks and bonds from or sell them to the option seller, at some specified price on some date.Or, more generally, it can be regarded as a contract, in which the option buyer receives some amount of money (payoff) determined in reference to the underlying asset price, from the option seller.There are some kinds of options with respect to timing of exercise of right.In an European option, the option buyer can exercise the right at one predetermined date, which is called the maturity.On the other hand, there are early-exercisable options, in which the option buyer can choose the exercise date.In an American option, the option buyer can exercise the right at any time before the final maturity T .In a Bermudan option, exercise of right is possible on any of finite predetermined dates including T .We hereafter call such dates exercise dates.
Major financial firms hold large portfolios of a wide variety of options, and therefore pricing them is an important task for business.However, it is also a difficult task.Basically, the option price is the expected value of the payoff under some stochastic model describing random time evolution of the underlying assets.Although European options can be sometimes priced easily, for example, by some analytic formulas, pricing Bermudan and American options typically involves heavy numerical calculations.The difficulty partly stems from the nature of the problem as dynamic programming.That is, pricing early-exercisable options contains determining the optimal exercise time as a crucial part.Although there are some kinds of methods which aim to reflect early exercise to the option price, each one has pros and cons.
One major category of pricing methods is the Monte Carlobased method 3 , in which we generate many sample paths of evolution of underlying asset prices, and estimate the expected payoff as an averaged payoff over the paths.This approach has an advantage in the case of multiple underlying assets.Namely, in this approach, the estimation error on the option price decays as O(N −1/2 ) 4 when the sample number N increases, regardless of the number of the underlying assets d.In other words, it suffices to take O(ǫ −2 ) samples in order to achieve the error tolerance ǫ on the option price.This contrasts to other methods, for example approaches based on solving partial differential equations [27,28], whose complexity is O((1/ǫ) poly(d) ).On the other hand, in the Monte Carlobased methods, it is difficult to precisely determine the optimal exercise time, and we have to approximate this in some way.In many cases, this is done through approximation of the continuation value, which is the option price at each exercise date in the case that the option buyer forgoes the exercise.The option should be exercised if the payoff is larger than the continuation value, and should not be exercised otherwise.
In this category, the least-squares Monte Carlo (LSM) [29] is widely used.LSM estimates the continuation value at each exercise date by linear regression using the generated sample paths as training data, and then, going backward from the final maturity to the present, finds the present option price.
Note that this method can also price American options approximately, replacing exercisability at any point in the con-tinuous time period with that at discrete dates with sufficiently small intervals.
In this paper, we propose a new method for Bermudan option pricing, combining Chebyshev interpolation and quantum algorithm for Monte Carlo integration [30][31][32], which is based on quantum amplitude estimation (QAE) [31,[33][34][35][36][37][38][39][40][41][42].As far as the author knows, this is the first proposal on the quantum method for Bermudan option pricing.Chebyshev interpolation is a widely used method for function approximation 5 , and has already been used in some (classical) methods for Bermudan option pricing [44][45][46][47][48][49][50].In the proposed method, given the access to the quantum circuit (or, the oracle) for time evolution of underlying asset prices, we calculate the continuation values at the interpolation nodes by the quantum algorithm, and find Chebyshev interpolation using these values.Importantly, this method outputs an estimation of the option price with the error at most ǫ, calling the oracle only O(ǫ −1 ) times.Thus, as we commonly observed in applications of QAE to various kinds of problems, we obtain the quadratic speed-up compared with the classical Monte Carlo-based methods such as LSM and the Chebyshev interpolation-based methods.
The rest of this paper is organized as follows.In Section II, we briefly explain Chebyshev polynomials and function approximation by them.We present how to calculate the coefficients of Chebyshev expansion in general and the upper bound for the approximation error.In Section III, we present the general formulation of Bermudan option pricing and explain LSM as a typical classical solutions.In Section IV, we explain QAE and QAE-based Monte Carlo integration.Then, in Section V, we present the new algorithm for Bermudan option pricing based on Chebyshev interpolation and QAE.We also present an upper bound on the price error in the method, and that on the complexity sufficient to achieve the given error tolerance.Section VI summarizes this paper.All proofs are presented in the appendix.

A. Notations
We here explain the notations used in this paper.N denotes the set of all positive integers, and For a, b ∈ N 0 , δ a,b denotes the Kronecker delta, which is 1 if a = b or 0 otherwise.For d ∈ N and l 1 , l 2 ∈ N d 0 , we also define δ l 1 , l 2 , which is 1 if l 1 = l 2 and 0 otherwise.
For a measure space (Ω, F , µ) and p ∈ R ≥1 , L p (Ω, µ) denotes the L p space on it.
The indicator function 1 C takes 1 if the condition C is satisfied, and 0 otherwise.
In this paper, we consider quantum states of systems consisting of some quantum registers with some qubits.For x ∈ R, |x denotes one of the computational basis states on some register, whose bit string corresponds to the binary representation of x with truncation at some digit.For d ∈ N and x = (x 1 , ..., x d ) T ∈ R d , | x := |x 1 ... |x d is the state on the d-register system.

II. APPROXIMATION OF FUNCTIONS BY CHEBYSHEV INTERPOLATION
For l ∈ N 0 , the l-th Chebyshev polynomial (of the first kind) is defined as where x ∈ [−1, 1].One of its important properties is the discrete orthogonality: for any m ∈ N 0 and l 1 , Here, x m, j is the Chebyshev node defined as for j ∈ [m] 0 .x m,0 , ..., x m,m are the zeros of T m+1 (x).We also define the tensorized Chebyshev polynomials on a general hyperrectangle in R d , where d ∈ N.That is, given with for every l = (l 1 , ..., l d ) T ∈ N d 0 and S = (S 1 , ..., S d ) T ∈ D, where For the above polynomials, the orthogonality relation is now for every m ∈ N and l 1 , with j = ( j 1 , ..., j d ) T ∈ [m] d 0 .We can use the above polynomials for function approximation.Given D as (4) and m ∈ N, we define the Chebyshev interpolation of a function f : D → R as for every S ∈ D, where the coefficient a f, l is calculated by for every l ∈ [m] d 0 .This is in fact an interpolation, since D .The error in the above approximation has been investigated in [47,51].They gave the error bound, making an assumption on analyticity of the interpolated function f .We here present the theorem on such a error in the case where we are given the values of f at the Chebyshev nodes with some errors [47].However, let us make some definitions prior to the theorem.For ρ ∈ R >1 , the Bernstein ellipse B ρ is defined as the open region in the complex plane bounded by the ellipse 1  2 u + u −1 u ∈ C, |u| = ρ .We also define the generalized Bernstein ellipse as Furthermore, we define the multivariate version of the Lebesgue constant of the Chebyshev nodes: for every m ∈ N, where for every j ∈ [m] 0 .As [47] showed, holds, which is derived from the well-known upper bound Λ 1,m ≤ 2 π log(m + 1) + 1 [43].Then, the theorem is as follows 6 .
Theorem 1.Let d and m be positive integers.Let D be a hyper-rectangle like (4).Let f : D → R be a function that has an analytic extension to B D,ρ for some ρ ∈ R >1 .Besides, assume that sup s∈B D,ρ | f ( s)| ≤ B for some B ∈ R.Moreover, suppose that we are given a real number f j for every j ∈ [m] d 0 , and that there exists holds.Here, for every S ∈ D, f ( S ) is defined as with the coefficients ã l calculated by for every l ∈ [m] d 0 , and

A. General formulation
In this paper, we consider pricing a Bermudan option with d ∈ N underlying assets and K ∈ N exercise dates t 1 , ..., t K , which satisfy t 0 < t 1 < ... < t K with t 0 := 0 being the present and t K := T ∈ R + being the final maturity.This is formulated as follows.Under some filtered probability space (Ω, F , (F t ) t≤0 , P), consider the S-valued Markov process S (t) := (S 1 (t), ..., S d (t)) T , where S is a subset of R d equipped with its Borel σ-algebra inherited from R d , and S 0 := S (0) is deterministic.S (t) corresponds to the values of the underlying asset prices at time t, or transformations of them by some function (for example, the logarithms of the asset prices).We are mainly interested in its values at t 1 , ..., t K , that is, the discrete-time process S k = (S 1,k , ..., S d,k ) := S (t k ), k ∈ [K] 0 .We hereafter denote an instance of this process, which is a (K + 1)-tuple of elements of S, as S = ( S 0 , S 1 , ..., S K ).Besides, suppose that we are given the function f , where ρ k is the image probability measure on S induced by S k .This corresponds to the payoff which arises by the exercise at t k .Although we assume that the risk-free rate is 0 for simplicity in this paper, we can consider that f pay k is the discounted payoff, that is, the product of the payoff and the discount factor.
Then, the price of the Bermudan option at t k with S k = s ∈ S is given as where E[•] denotes the (conditional or unconditional) expected value with respect to P, and T k , k ∈ [K] is the set of all {k, ..., K}-valued stopping times.In particular, the present option price is where T = T 1 .
The problem to find V 0 can be written as a kind of dynamic programming, that is, for every s ∈ S, and Here, for every k ∈ [K − 1] and s ∈ S, is called the continuation value.This corresponds to the option price at t k in the case that the option has not been exercised at that time and that S k = s.
Note that this problem can be seen as finding the optimal exercise date τ op ∈ T , which maximizes (20).This can be recursively determined as and τ op = τ 1 .Also note that for every k ∈ [K − 1], which means that Q k is the expected value of the payoff under the exercise strategy τ k+1 .
Omitting some technical details, we describe the outline of LSM as follows.As a preparation, for every k ∈ [K − 1], we determine the set of functions , the set of all real-coefficient polynomials on R d of degree at most m ∈ N, and we hereafter consider this.
Next, we generate N samp ∈ N ≥2 sample paths of underlying asset prices, which are denoted as . Then, we determine the stopping time, which approximate the optimal one τ op , by the following procedure.First, we set τ or, in other words, best fits to the realized payoffs under the stopping time τ k+1,i on the sample paths.It is guaranteed by statistical leaning theory that fitting to the sample values of the payoff, which distribute around the continuation value, yields the approximation of the continuation value [55,57,[59][60][61].
Then, we set for every i ∈ [N samp ].By repeating this until we reach k = 1, we get τ 1,i , and finally Let us make some comments on the procedure.First, note that it is assumed that we can generate sample paths S i .In the usual situation where the stochastic differential equation (SDE) for S (t) is given, we can use some method for numerical simulation of SDEs such as Euler-Maruyama method.Second, we mention how to find g k minimizing (26).Note that this is just least-squares linear regression, since R d,m is a vector space.Therefore, we can solve this by various methods, for example, solving the normal equation of linear regression, some numerical optimization, and so on.
Then, let us mention the relationship between the error and the sample number in LSM.According to [61], under some technical assumptions, taking appropriately large m, we obtain the error bound on the option price which scales as Here, E samp [•] denotes the expectation with respect to randomness of samples, and n and p are the quantities which characterize smoothness of Q 1 ( S ), ..., Q K−1 ( S ) and boundedness of the norms of f pay 1 ( S ), ..., f pay K ( S ), respectively (see [61] for more details).For larger p and n, the RHS of (28) decreases faster against the increase of N samp .In the limit of n, p → ∞, which means that Q k 's are highly smooth and the norms of f pay k 's are well-bounded, the RHS of (28) becomes O(N −1/2 samp ), which coincides with the well-known error decay rate in Monte Carlo integration.Conversely, in this limit, it is sufficient to take O(ǫ −2 ) samples in order to achieve the error tolerance ǫ.

IV. QUANTUM AMPLITUDE ESTIMATION AND QUANTUM ALGORITHM FOR MONTE CARLO INTEGRATION
A. Quantum amplitude estimation (QAE) We here briefly review QAE.Consider the system consisting of a quantum register R 1 and a qubit R 2 .Suppose that we are given the oracle A, which transforms |0 |0 , the state in which all qubits in R 1 and R 2 are set to |0 , into with some a ∈ (0, 1).Here, the first and second kets correspond to R 1 and R 2 respectively, and |ψ 0 , |ψ 1 are some quantum states.Then, our goal is estimating a, which is the probability to obtain 1 in R 2 when we measure |Ψ , with the error tolerance ǫ.There exist some algorithms which output such an estimation with O(ǫ −1 ) calls to A and its inverse A † in total [31,[33][34][35][36].Although these QAE algorithms output a number close to a not with certainty but with some probability, we can enhance the success probability by running QAE many times and take the median of the outputs [30,62].Let us define the (N QAE , N rep )-QAE as the method for estimating a which runs N rep rounds of QAE and makes N QAE calls to A and A † in total in each round.Obviously, in the (N QAE , N rep )-QAE, A and A † are called N QAE N rep times in total.Then, combining Theorem 12 in [33] and Lemma 6.1 in [62], we obtain the following theorem.
Theorem 2. Suppose that we are given the accesses to A in (29) and its inverse A † .Then, for any γ ∈ (0, 1) and ǫ ∈ (0, 0.1), a (N QAE , N rep )-QAE, where the positive integers N QAE and N rep satisfy and where a is defined as (29), with probability higher than 1 − γ.

B. Quantum algorithm for Monte Carlo integration
One application of QAE is the algorithm for Monte Carlo integration, that is, the method to calculate expected values.Suppose that we want to calculate E[F( X)], the expected value of F( X), where X is some real vector-valued stochastic variable and F is a real-valued function acting on X.We also assume that the range of F is in [0, 1], and, if not, we make this hold by adding and/or multiplying some constants to F. Furthermore, suppose that we can use the following oracles O X and O F .O X is the oracle to generate the state in which the distribution of X is encoded.That is, O X operates on a quantum register and transform the state with all qubits set to |0 into where Here, we assume that the set of all the values that X can take is finite.If X is continuous, we need some discretization.How to create states corresponding to widely used distributions such as normal distribution has been investigated [13,63].The second oracle O F operates on a two-register system, and, using the first register as the input x, outputs F( x) into the second register.That is, for any x in the domain of F, By these oracles, the following computation is possible.
Preparing two registers R 1 , R 2 and a qubit R 3 , and initializing all of them to |0 , we perform where the first, second and third kets correspond to R 1 , R 2 and R 3 , respectively.We use O X and O F at the first and second arrows, respectively.The transformation at the third arrow is done by arithmetic circuits [64] and controlled rotation gates.Note that the probability to obtain 1 in R 3 when we measure the final state in (36) is Therefore, using the whole operation in (36) as the oracle A, we can estimate E[F( X)] by QAE.

V. BERMUDAN OPTION PRICING BY CHEBYSHEV INTERPOLATION AND QAE
Now, let us present the method for Bermudan option pricing by Chebyshev interpolation and QAE.

A. Assumptions
We begin with making some assumptions necessary to execute the proposed method.The first one is as follows.
Assumption 1.We are given the access to the oracle O step,k , which generates the state corresponding to the probability distribution of S k+1 conditional on S k .That is, for every k ∈ [K − 1] 0 and S ∈ S, where S k+1 ( S ) is the set of possible values of S k+1 under the condition that S k = S , and p k+1 ( s; S ) : We here make comments on how to implement O step,k .As mentioned in Section III, usually, following some SDE and some numerical method such as Euler-Maruyama, we can generate random sample values of S k+1 with the given value of S k as the initial condition.Implementations of such a calculation on quantum circuits have been discussed in the previous papers [7,9,13].That is, we can prepare the states corresponding to some (discretely approximated) random variables (e.g. standard normal) on the other registers, and, using them at discretized time steps, generate the path of S (t) from t k to t k+1 .This yields the state like (37).We should also note that, in Assumption 1, it is assumed that S k+1 can take only a finite number of values for the fixed S k .This is not the case in the most models of S (t), in which it takes continuous values.However, under the aforementioned implementations for time evolution of S (t), in which both time and random variables are discretely approximated, the number of the possible values of S k necessarily becomes finite.
Hereafter, we are mainly interested in the number of calls to O step,k in calculating the option price as a measure of complexity, since calculation for time evolution of underlying asset prices is typically the most time-consuming part in option pricing.
The second assumption is as follows.Here, I A denotes the set of all real-valued functions on a given subset A ⊆ R d .Assumption 2. For every k ∈ [K − 1], we are given the following such that the following (i) and (ii) are satisfied.
(i) There exists ǫ OB k ∈ R + such that either or is satisfied for any s ∈ S \ D k .Here, F k [•] is the 'flat extrapolation operator' defined as for any F ∈ I D k and s = (s 1 , ..., s d ) T ∈ S.

(ii) If, for some G ∈ I D k , we have the access to the oracle
for any s ∈ D k , we also have the access to the oracle O G , which acts as Here, G k [•] is defined as for any H ∈ I D k and s ∈ S, where A k is a subset of S \ D k such that ( 39) and ( 40) hold for any s ∈ A k and any s ∈ (S \ D k ) \ A k , respectively.
We also define G K [H]( s) := H( s) for any H ∈ I S and s ∈ S Roughly speaking, this assumption means that, when some of underlying asset prices are extremely large or small, we can approximate the option value V k by some known and easily computable function V OB k or the flat extrapolation of V k from moderate underlying asset prices.Postponing explanation on why this assumption is necessary to Section V B, we here see that it is actually satisfied in some typical settings in option pricing.For example, let us consider a basket put option, whose payoff function is under some model in which S 1 (t), ..., S d (t) are unbounded from above but bounded from below, say, by 0, as the Black-Scholes model.Then, in each of the following situations, (39) or (40) holds.
• If some of S 1,k , ..., S d,k are extremely large, the option is far out-of-money, and therefore its price is almost 0.
• If some of S 1,k , ..., S d,k are smaller than the sufficiently small thresholds L 1,k , ..., L d,k ∈ R + respectively, but the others are not, setting the former to the thresholds hardly affects the option price.
• If all of S 1,k , ..., S d,k are sufficiently close to 0, the option is exercised at t k , and therefore Thirdly, we make the following assumption, which is necessary for bounding the interpolation error in the proposed method.
Assumption 3.For every k ∈ [K − 1], Q k ( S ) has an analytic extension to B D k ,ρ k , where D k is given in Assumption 2 and ρ k is some real number greater than 1, and holds, where B k is some positive real number.

B. The proposed method
Under these assumptions, we can construct the procedure for Bermudan option pricing based on QAE and Chebyshev interpolation.This is also a backward calculation similarly to LSM; we sequentially calculate the approximate continuation value Q k and option price V k at t k , going from the final maturity to the present.Roughly, the outline is as follows.As preparation, for every k ∈ [K − 1], we set m k ∈ N, the degree of Chebyshev polynomials used for the approximation, and the hyper-rectangle We begin the iterative calculation by setting V K ( S ) := f .Using these, we construct Q k , the Chebyshev interpolation of the approximate continuation value, and set V k ( S ) = max{ Q k ( S ), f pay k ( S )} for every S ∈ S. We repeat these steps until we reach k = 1.Finally, we estimate the expected value of V 1 ( S 1 ) by QAE again, and let the result be an approximation of V 0 .
The fully detailed procedure is shown in Algorithm 1.

Algorithm 1 The method for Bermudan option pricing based on Chebyshev interpolation and QAE
Input: , the degree of Chebyshev polynomials.
, the hyper-rectangle for Chebyshev interpolation.Here, -QAE, obtain an estimation P k, j of the probability P k, j to observe 1 on the last qubit in measuring |Ψ k, j in (49), and let (2 for every S ∈ D k , with ãk, l calculated as (48) for l ∈ J k .7: Set V k ( S ) := max f pay k ( S ), Q k ( S ) for any S ∈ D k .8: end for 9: Using (N QAE 0 , N rep 0 )-QAE, obtain an estimation P 0 of the probability P 0 to observe 1 on the last qubit in measuring |Ψ 0 in (52), and output (2 Some additional explanations should be made.The first one is on |Ψ k, j in Step 4. For every k ∈ [K − 1] and j ∈ J k , given the approximation V k+1 ∈ I D k+1 of V k+1 , we generate the state |Ψ k, l on the appropriate multi-register system with the last one being single-qubit, by the following operation: where O step,k in Assumption 1 and O V k+1 in Assumption 2 are used at the second and third arrows, respectively.Note that the probability to obtain 1 on the last qubit in measuring |Ψ k, l is where Therefore, as long as . This is why we can obtain approximations of the continuation values at Chebyshev nodes by Step 4, with the errors from QAEs being also small.
Second, let us explain the state |Ψ 0 in Step 9. Given V 1 , we can generate |Ψ 0 similarly to |Ψ k, j as where the last ket corresponds to a single-qubit register.Since the probability P 0 to obtain 1 on the last qubit in measuring |Ψ 0 satisfies where we can obtain an approximation of V 0 by Step 9, as long as and the QAE error is small.Lastly, let us comment on the reason why Assumption 2 is necessary.This is because we have to handle underlying asset prices out of D k+1 in Step 4 and 9, or, more specifically, in generating |Ψ k, j and |Ψ 0 .In fact, when we generate |Ψ k, j , S k+1 can be out of D k+1 with some probability.In particular, when |Ψ k, j corresponds to a Chebyshev node S D k ,m k j close to the boundary of D k , or, in other words, the condition that S k is close to the boundary of D k is imposed, such a probability becomes non-negligible.

C. Evaluation of the error
Then, let us consider the error on the present option price in the proposed method.First, we have the following theorem.

Theorem 3. Under Assumptions 1 to 3, consider Algorithm 1. Suppose that, for every k ∈
is satisfied, where ǫ QAE k is some positive real number.Moreover, suppose that holds, where, for k and The proof is presented in Appendix A.

D. Complexity
Based on Theorem 3, we can evaluate the complexity, that is, the number of calls to O step,k sufficient to achieve the desired level of the error on the present option price.
with probability higher than 0.99.
The proof is presented in Appendix B.
We here explain why the parameters are set as above.As we see in the proof in Appendix B, m 1 , ..., m K−1 satisfying (61) make the first term in the RHS in (57) that is, the total number N tot of calls to {O step,k } k=0,...,K−1 divided by the QAE repetition number N rep , is minimized under the constraint that, if all the QAEs in Algorithm 1 succeed, the third term in the RHS in (57) is smaller than ǫ V max 1 /2.Finally, {N rep k } k=0,...,K−1 are determined so that the probability that these QAEs all succeed becomes higher than 0.99 = 1 − 0.01.In total, Algorithm 1 with the setting in Corollary 1 gives an approximation of V 0 with an error at most ǫ V max 1 with probability higher than 0.99.
Note that, in reality, it is difficult to set m k to m th k , since ρ k and B k are usually unknown.In practice, we might set them to some conservatively large values, based on, for example, the calculation results of some benchmark pricing problems for various {m k } k=1,...,K−1 .Besides, note that, in the above setting, the half of the error tolerance ǫ is assigned to the interpolation error and another half is assigned to the QAE error.Although we can of course change this assignment ratio, it affects the complexity only logarithmically, since the sufficient levels of {m k } k=1,...,K−1 are logarithmically affected by such a change and so are {N Let us consider the dependency of the total complexity on the error tolerance ǫ.We see that for every k ∈ [K − 1], and that where polyloglog(•) means polylog log(•) .Combining these with (65), we obtain which eventually beats LSM's complexity O(ǫ −2 ) for small ǫ.

E. Comparison with existing Chebyshev interpolation-based methods
In fact, the idea that we approximate the continuation value by Chebyshev interpolation is not novel.There are some classical methods for Bermudan option pricing based on Chebyshev interpolation [44][45][46][47][48][49][50].However, in addition to whether we use QAE or other classical methods for calculating the nodal continuation values, there are the following differences between the above proposed method and the existing methods.
First, we note that we do not have to use Monte Carlo for calculating the continuation value, and [44,45] actually used other methods.These papers considered the situation where the transition probability of the underlying asset prices can be easily calculated, and computed the continuation value by the numerical integration of the product of the transition probability and the option value at the next exercise date.Note that this way is possible only for some simple models for underlying asset evolution such as the Black-Scholes model.On the other hand, in more complicated settings where, for example, we price a multi-asset option under the stochastic local volatility model, Monte Carlo can be the sole solution, and such a time-consuming situation is a meaningful target for quantum speed-up.However, combining Chebyshev interpolation with various methods might be an interesting possibility also in the quantum setup, and worth to be investigated as a future work.
Let us also mention the differences from [48].The major difference is that, in the method in [48], the continuation value is not the target of either Monte Carlo integration or Chebyshev interpolation.Instead, the method calculates the conditional expectations of Chebyshev polynomials by Monte Carlo or other methods, and find the Chebyshev interpolation of not the continuation value but the option price at each exercise date.This approach saves the computational time when we price many options under a same model, since we can reuse the conditional expectations for interpolations in pricing different options.Considering the quantum version of this approach might be interesting too.
F. Exponential factor with respect to the number of exercise dates in the error bound Now, let us make a comment on the factor Λ 1,K−1 in (57), which is reflected into the polyloglog factors in (69).This exponentially depends on the number of exercise dates K. Therefore, it seems that the error on the option price exponentially grows as K increase, and so does the complexity sufficient to achieve a given error tolerance.Similar situations arose in the error analyses for LSM [53,55,57,[59][60][61] and classical Chebyshev interpolation-based methods [48].
However, we should note that ( 57) is an upper bound on the error, and that the actual error might not necessarily grows exponentially against K.In fact, in the numerical experiment in [48], where American options were approximately priced as Bermudan options with a small exercise date interval, the error was suppressed even if hundreds or thousands of exercise dates were set.
Let us now consider why the factor exponentially depending on K appears in (57).In the derivation of (57) described in Appendix A, we make Assumption 3 on the analyticity and boundedness of the continuation values Q k , and apply Theorem 1 to Chebyshev interpolation of Q k in Algorithm 1.Since with some errors, in interpolation, the term like the second term in the RHS in (1) arises in the upper bound of the difference between Q k and the interpolant Q k , and is amplified at every later interpolation.
On the other hand, we can consider that the actual target of Chebyshev interpolation is not as not an estimate on . Then, we can make an assumption not on analyticity and boundedness of Q k but those of Q k .This leads to a different error bound than Theorem 3. Actually, if we make Assumption 4 instead of Assumption 3, we have an error bound with no exponential factor as shown in Theorem 4.
holds, where B k is some positive real number.
The proof is presented in Appendix C. Note that a similar point has been made for LSM in [59] (Theorem 3.1).Of course, Q k is defined as (51) with V k+1 , which is the intermediate output in Algorithm 1, and making assumptions on such a thing does not lead to self-contained discussion.It is more desirable to derive the error bound under assumptions on Q k and/or other quantities determined independently from pricing algorithms.We leave considering whether we can obtain an error bound similar to Theorem 4 under such assumptions or not as a future work.

G. Quantization of LSM
Lastly, we make a comment on whether we can consider the quantum algorithm for LSM.Since there are some quantum algorithms for linear regression [65][66][67][68][69][70][71], we naturally wonder that we can apply these to LSM and then obtain speed-up.However, this is not so straightforward, since most of these algorithms output the regression result as a quantum state, in which the regression coefficients are amplitude-encoded.Fortunately, some algorithms [67,71] output the regression coefficient as classical data.In particular, the algorithm in [71] has the complexity of O(D 7/2 κ 4 /ǫ) with the tolerance ǫ, the explanatory variable number D, and the condition number κ of the design matrix, from which we expect the quadratic speedup of LSM with respect to ǫ.Nevertheless, applying this algorithm to LSM is not immediate either, because of some points to be considered.For example, the complexity has strong dependence on D and κ, which might make the algorithm disadvantageous.Therefore, it will be crucial to evaluate these, especially κ, in addition to finding the basis function set which makes κ as small as possible.We will consider this direction in the future work.

VI. SUMMARY
In this paper, we considered application of quantum algorithms to Bermudan option pricing.Since there are QAEbased algorithms for Monte Carlo integration, which provide quadratic speed-up compared with the classical counterparts, and applications of them to some option pricing problems have been investigated, it is natural to consider to apply them to Bermudan option pricing.One crucial issue in this problem is how to approximate the continuation value Q k , which determines the optimal exercise date.In order to cope with this, we considered combination of QAE and Chebyshev interpolation.That is, the proposed method estimates the values of Q k on the interpolation nodes by QAE, and, using such estimates, find a Chebyshev interpolation as an approximation of Q k .We presented the calculation procedure in detail, along with the error bound and the complexity, which corresponds to the number of calls to the oracle for underlying asset price evolution, sufficient to achieve the desired error tolerance ǫ.As expected, this method has the complexity depending on ǫ as O(ǫ −1 ), which means the quadratic speed-up compared with LSM, the typical classical algorithm for Bermudan option pricing.
As a future work, it is interesting to consider the quantum version of LSM, as mentioned in Section V G. Besides, it is also meaningful to extend the proposed method to other types of dynamic programming, which is ubiquitous in many fields of science and industry.
holds.We prove this by induction.For k = K −2, (A1) implies that Similarly, if (A13) hold for k ∈ {2, ..., K − 2}, (A1) implies that Therefore, (A13) is proved for every k Finally, the claim is proved as follows.We see that where we used (A7) and (A13) at the second and last inequalities, respectively.Combining this and (56), we obtain On the other hand, Theorem 2 implies that, for every k ∈ [K−1] and j ∈ J k , using N QAE k set as (62), Step 4 in Algorithm 1 gives us P k, j , an estimate on P k, j in (50), satisfying P k, j − P k, j ≤ ǭk , (B2) and then with some probability.Similarly, it is implied that, with N QAE 0 set as (62), Step 9 gives us P 0 , an estimation on P 0 , satisfying and then V 0 satisfying with some probability.Therefore, when all of these succeed, holds, according to Theorem 3. Here, we used (B1), (B3), (B5), and the assumption that ǫ OB 1 , ..., ǫ OB K−1 are 0, along with simple algebra.
The remaining task is proving that the probability P all that all of the estimations in Steps 4 and 9 succeed is larger than 0.99 under the setting of N rep k as (64).Note that, according to Theorem 2, with the setting as (64), the probability that Step 4 for a set of k ∈ [K − 1] and j ∈ J k outputs P k, j satisfying (B2) is higher than Similarly, the probability that Step 9 outputs P 0 satisfying (B4) is also higher than (B7 Proof.Because of (55), Theorem 1 implies that holds with ǫ K = 0 and ǫ OB K = 0, which leads to similarly to (A12).Therefore, for every k ∈ [K − 1], This implies Finally, combining (C5) with (56) and |V 0 − V 0 | ≤ ǫ 1 + ǫ OB 1 , which we can see as (A16), we obtain )

payK
( S ) for every S ∈ S.Then, for k ∈ [K − 1], given V k+1 , we estimate the expected value of V k+1 ( S k+1 ) under the condition that S k = S D k ,m k j for every Chebyshev node S D k ,m k j by QAE, and denote the estimation as Q QAE k, j
k ∈ [K−1] and S ∈ D k .Besides, for every k ∈ [K−2] and S ∈ D k , Q k ( S ) − Q k ( S ) ≤ ǫ k+1 + ǫ OB k+1 holds as (A8), where ǫ k is defined as (A2) for every k ∈ [K − 1], whereas Q K−1 ( S ) = Q K−1 ( S ) by definition of G K [•].Combining these, we see that, for every k ∈ [K − 1] and S ∈ D k , ). Besides, the total number of these estimations is N est .Combining these, we obtain a lower bound of P all as Appendix C: Proof of Theorem 4