Quantum-enhanced analysis of discrete stochastic processes

Discrete stochastic processes (DSP) are instrumental for modeling the dynamics of probabilistic systems and have a wide spectrum of applications in science and engineering. DSPs are usually analyzed via Monte-Carlo methods since the number of realizations increases exponentially with the number of time steps, and importance sampling is often required to reduce the variance. We propose a quantum algorithm for calculating the characteristic function of a DSP, which completely defines its probability distribution, using the number of quantum circuit elements that grows only linearly with the number of time steps. The quantum algorithm reduces the Monte-Carlo sampling to a Bernoulli trial while taking all stochastic trajectories into account. This approach guarantees the optimal variance without the need for importance sampling. The algorithm can be further furnished with the quantum amplitude estimation algorithm to provide quadratic speed-up in sampling. The Fourier approximation can be used to estimate an expectation value of any integrable function of the random variable. Applications in finance and correlated random walks are presented. Proof-of-principle experiments are performed using the IBM quantum cloud platform.


I. INTRODUCTION
Simulation of physical processes on quantum computers [1,2] has many facets, with recent developments improving the usage of this technology [3][4][5][6][7][8][9][10][11][12][13][14].While one obvious task for a quantum computer is to simulate quantum mechanical behaviors of nature [1], quantum simulation of probability distributions and stochastic processes has gained attention recently [15][16][17][18][19][20].Since quantum mechanics can be viewed as a mathematical generalization of probability theory, where non-negative real-valued probabilities are replaced by complex-valued probability amplitudes, quantum computing appears to be a natural tool for simulating classical probabilistic processes.Intuitively, some quantum advantage is expected since the probability amplitudes can interfere, unlike in the classical probabilistic computing.Indeed, a quadratic quantum speed-up has been reported for solving financial problems when compared to Monte Carlo simulations [16].
Although the main focus of this work is to achieve quantum sampling advantage for discrete stochastic processes (DSPs), we note in passing that quantum memory advantage for stochastic processes with causal structures has been shown elsewhere [18][19][20].
Classically, Monte Carlo methods are essential for estimating expected values of random variables in DSPs, since the number of realizations increases exponentially with the number of time steps.When the Monte Carlo sampling is repeated N times, the expectation value to be found converges with O(1/ √ N ) regardless of the number of realizations.The central limit theorem ensures this, but it is important to note that this convergence is only attained in the limit of N → ∞.When this limit is not nearly attained, it is often crucial to have sampling strategies to reduce the variance of the estimate.The reason for this is that less likely events are also less likely to be sampled from, while such events can be of high impact.This can create a bias in the computation as the estimation can be dominated by more probable, but less important values.To correctly sample with a given N that is not near the limit of the large number, it is beneficial to modify the probability of the process in a way that balances the importance of the events.This is called importance sampling [21][22][23][24][25][26][27].This makes apparent two problems with Monte Carlo algorithms.First, it may be difficult to calculate the probability of a particular realization of a random variable.Second, importance sampling strategy relies on sophisticated understanding of the probability distribution.
In this manuscript we show that the characteristic function of a DSPs can be efficiently calculated on a quantum computer, and in so doing we introduce an effect that can be called quantum brute-force.The random variables of the DSP of interest do not need to be identically and independently distributed, and non-Markovian processes can also be studied.The quantum state space that grows exponentially with the number of qubits is used to move along all paths of the discrete stochastic process simultaneously in quantum superposition, and hence the term quantum brute-force is adequate.
With a Pauli measurement scheme on a single qubit, essentially the probability of a Bernoulli trial needs to be arXiv:2008.06443v1[quant-ph] 14 Aug 2020 estimated which exhibits the optimal variance that can be achieved according to the central limit theorem.This leads to a crucial result that no sampling strategies are necessary.Moreover, we connect recent developments in quantum finance [15,16,28] to our method and show that sampling convergence can be improved by means of quantum amplitude estimation (AE) by a power of two.The methods introduced here point therefore to an exciting and promising use of quantum computers: less variance and faster convergence for Monte-Carlo sampling.
Applications of our method span extensively across many fields in science and engineering; any physical behaviour that can be mathematically modelled as a DSP can be studied in principle.In particular, we show how the found simulation procedure can be applied to option pricing theory and to correlated random walks which leads to various applications in biology, ecology and finance.For each example, we experimentally demonstrate the proof-of-principle using the IBM quantum cloud platform.

II. RESULTS
A discrete stochastic process can be described with n discrete random variables X l : Ω l → R, l = 1, . . ., n for some n ∈ N + , each having at most k non-zero realizations, i.e. given a sample space Ω l , there exist at most k elements x l,0 , . . ., x l,k−1 ∈ R with X l (Ω l ) = {x l,0 , . . ., x l,k−1 }.Hereinafter, we use the following notations.Each realization of the stochastic process is identified with an index vector j = (j 1 , . . ., j n ) ∈ K n with K := {0, . . ., k − 1} and is denoted by x(j) = (x 1,j1 , . . ., x n,jn ) .Moreover, sum {x(j)} = n l=1 x l,j l and x (m) (j) = (x 1,j1 , . . ., x m,jm ) with m ≤ n.The first and the second subscripts of the random variables and probabilities label the time step and the event, respectively.A quantity of interest for such processes is the expectation value of an integrable function f : R → R of the random variable Now, we explain how to encode the described stochastic process in a quantum state, and evaluate equation ( 1) by making a measurement on the quantum state.The quantum state consists of an index system and a data system defined in a Hilbert space where I (D) indicates the index (data) system and d is determined by the problem of interest.Each realization of the stochastic process is represented as a unitary operator U (•) : R n −→ B(H I ⊗ H D ) parametrized by some n-dimensional vector (B(H) is the space of linear operators on a Hilbert space H).Then a DSP can be represented as The factor of each part of the sum are denoted by p(j) with j p 2 (j) = 1 and the state of the index system is denoted by where |j l for j l = 0, . . ., k − 1 is the orthogonal basis.As described in Ref. [29], measuring an expectation value of an observable M on the data system of the final state in equation ( 2) yields the convex sum of independent expectation values measured from all k n trajectories as where M (j) = U † (j)M U (j).The coefficient p 2 (j) can be identified with the joint probability and the expectation value with the evaluation of f , i.e.
for a function f : R → R that we will specify below.
In the worst case, evaluation of equation ( 3) requires two expensive procedures as follows.First, k n probabilities need to be encoded as the amplitudes of k n computational basis given by n qudits of dimension k.This can be done with various quantum state preparation techniques with substantial amount of computational overhead [30,31].Next, k n unitary operators, conditioned on all possible index states, need to be applied on an input state |ψ .Such operators can be expressed as with V : R → H D .On the other hand, for many interesting DSPs, the number of necessary unitary operators can be reduced to O(nk).
Before explaining such exponential-reduction in detail, we introduce two results in the next two propositions (with proofs provided in Supplementary Information) to establish the grounds for the measurement scheme.
Proposition 1 (Pauli X and Y Measurement).Let √ 2, we find As the true value of E [cos(S n )] and E [sin(S n )] must be estimated, in general the convergence behaves according to the central limit theorem, which guarantees that the measurement statistics approaches to a normal distribution around a mean that corresponds to E [f (S n )] as the number of experiments N S → ∞.The speed of convergence is moreover given by O(1/ √ N S ).Taking into account that the Pauli measurements have two eigenvalues, the task is essentially estimating a probability of a Bernoulli trial, which specifies how the central limit theorem is realized.Given a confidence of 1 − α, the number of experiments to be within a margin of error > 0 is ) being the quantile function.In contrast, with classical Monte Carlo sampling the convergence rate by the central limit theorem is achieved with the caveat that the Monte Carlo simulation samples from an usually unknown stochastic process and concise estimates on the number of experiments given a margin of error are in general not easily accessible.As we see, the property that equation ( 2) encompasses all possible paths with the correct probability with a quantum measurement leads to the seemingly small advantage of knowing the convergence before hand, irrespective of the distribution of the underlying DSP.Contemplation on this fact reveals that this is no small feat: the quantum advantage lies in the fact that no sampling strategies are necessary.
The quantum amplitude estimation algorithm [32] can provide further speedup compared to the convergence rate of classical Monte Carlo method given by the central limit theorem, as suggested in Refs.[15,16,28,33].The algorithm uses m ancilla qubits in addition to data and index qubits and O(poly(m)) number of Grover-like iterators followed by O(m 2 ) Hadamard and controlled phaseshift gates for quantum fourier transform (QFT) [34] to estimate an amplitude with an error of O(1/2 m ).The number of repeated measurements needed to reach confidence that the estimation succeeded is independent of m.The quantum AE algorithm can be adapted to our method to construct an even more powerful strategy by formulating an AE problem as follows.Given the final state in equation ( 2) written as a linear combination |Ψ f = |Ψ 0 + |Ψ 1 by separating the full Hilbert space into two orthogonal subspaces, we estimate the amplitude defined as a = Ψ 1 |Ψ 1 .Then, the following proposition connects AE with the DSP simulation. with With AE the value ã will be estimated, hence , then there exists a similar decomposition The above result opens up an exciting avenue towards a fast Monte-Carlo alternative without the need of sampling strategies.These propositions in conjunction with theorem 1 which is stated in the following section show that a quantum computer can simulate the quantities E [cos(S n )] and E [sin(S n )] efficiently.As a result, one can calculate a random variable's characteristic function In order to extend the above ideas to estimate expectation values of a range of integrable functions f : R → R, a Fourier-series is used.If f is P -periodic, then is the Fourier-approximation of order L for f (x).By linearity of the expectation value, this approximation carries over to As a consequence, it is possible to approximate any such expectation value in O(LN ) experiments, where N is the number of shots per experiment.Convergence on Fourier-series is a rich and mature field [35][36][37] which establishes basic results about convergence and the rate of convergence of each coefficient for given properties of the function f .The underlying idea of simulating DSPs on a quantum computer has been laid out.Now we show that the simulation can be performed efficiently with a quantum circuit.

A. Independent increments
To deliver the underlying idea with a simple example, we start by presenting the case for DSPs with independent increments.To this end, the strategy taken is as follows.We first construct the state j∈K n p(j) |j ∈ H I , and systematically entangle the data system D with the index system.When each of the increments are independent of each other, i.e.P [X i = x, X j = y] = P [X i = x] P [X j = y] for pairwise different i = j, equation (1) can be written as < l a t e x i t s h a 1 _ b a s e 6 4 = " B c T 9 X 4 e d 9 5 j v z q X x p C D / L P N 2 + X 8 = " < l a t e x i t s h a 1 _ b a s e 6 4 = " 5 U / U 6 W p A O g W w q z g a 5 < l a t e x i t s h a 1 _ b a s e 6 4 = " r m E q u + Z v R Z P < l a t e x i t s h a 1 _ b a s e 6 4 = " 5 U / U 6 W p A O g W w q z g a 5 FIG. 1. Circuit Design of DSPs with Independent Increments.(a) A quantum circuit including preparation of the index system (operators in blue denoted by A l ) and the realization of c-U (j) (operators in green denoted by U l (x) with a verticle line connected to black squares).A black square is placed on a control register, and the green box denoted by U l (x)is placed on a target register.Unlike the conventional symbol for a controlled-NOT gate between qubits, the black square indicates that one of the k different unitary transformation is performed conditioned on the state of the control register of dimension k.(b) Schematic circuit representation of the applications of c-U l ∈ B(HI l ⊗ HD) for l = 1, . . ., n as used in (a).We use the notation ).The open circle with a number j means that the unitary operator V l j is applied if the control register state is |j .
index-state.Then it can be described as a product state Each evolution from |α(l) → |α(l + 1) is done by an unitary operator A l+1 ∈ B(H I l ) as given by A l+1 |0 = k−1 j=0 p l+1,j |j .After n levels, i.e. the application of 14) is created.Since each of the operators A l only operates on a separate subspace H I l , they commute and can be applied in parallel.This is a consequence of the independence of the increments X l .For example, an n-step DSP with k = 2 possible paths at each time step can be realized with n index qubits each prepared by a single-qubit unitary operation A l |0 = R y (θ l )|0 = cos(θ l /2)|0 + sin(θ l /2)|1 , where θ l is chosen to satisfy cos 2 (θ l /2) = p l,0 and sin 2 (θ l /2) = p l,1 .Now, for each step of the DSP, k unitary operators applied to the data system controlled by an index qudit split the data space into k spaces, each attached to an orthogonal subspace of the index system.In other words, in each time step, the data system undergoes k independent trajectories.Thus n steps of k controlled unitary operations allows for the encoding of k n independent realizations of a DSP to the index-data quantum state.To this end, we identify to each realization x l,j l of the random variable X l to the application of the operator V (x l,j l ) ∈ B(H D ) to the data system.In fact, the operator will be defined in such a way, that the lth index-level state's branches |j l l with the amplitudes p l,j l (for j l ∈ K) will be controlling the operator, thereby identifying the probability p 2 l,j l = P [X l = x l,j l ] with the occurrence of V (x l,j l ).We denote a projection to the lth index subsystem I l as where 1l k is the k-by-k unit matrix and its perpendicular pendant as Π ⊥ lj .Then the following theorem establishes the desired result (see Supplementary Information for the proof).
with c-V lj l (x) ∈ B(H I l ⊗ H D ).Furthermore, for l = 1, . . ., n and x ∈ R k , define the operator The application of nk controlled operations c-V lj l (x l,j l ) is thus necessary.See Fig. 1 for a depiction of equations ( 16) and (17).
Propositions 1, 2 and theorem 1 establish that for given expectation values of equations ( 6) and ( 7), respectively, with O(nk) controlled gates.The above algorithm can be visualized quite intuitively.Consider a set of unitary operators V = {V lj l = V (x l,j l ) : l = 1, . . ., n, j l ∈ K} as above.These operators span an ordered tree of height n and a maximal degree of k.The label of the nodes themselves are of secondary importance, but each edge is labelled by one of the operators of the family V in such a way that, given the level l (l = 1 is the root), each node of the lth level has k edges, each of them are labelled by V lj l , j l ∈ K.A path between any two nodes, i.e. (V lj l , . . ., V l+ij l+i ), can be interpreted as a concatenation of operators, i.e.V (j l , . . ., After l level operations, paths from the root to every nodes at (l + 1)th level in the operator tree are travelled simultaneously.Hence the data system is evolved by the operator

B. Path-dependent increments
For a DSP with path-dependent increments, the probability of a particular realization can be expressed as This is the first-order Markov chain that the probability to take a particular path in a given time step is determined by which path was taken in the previous step.Discussions thus far can easily be extended to such DSPs.The only difference is in the preparation of the index system, i.e. in the construction of the state j∈K n p(j) |j .Unlike in the case of the independent increments, one needs to introduce entangling operations directly within the index system.Without loss of generality, we use the index system consisting of qubits, i.e. two paths at each level, to describe the procedure.First, the index qubit for the first level is prepared in A 0 |0 = cos(θ 0 /2)|0 + sin(θ 1 /2)|1 as before, while the rest of the qubits are in |0 .Then the second index qubit is prepared by using the controlled operation |0 0| ⊗ R y (θ 2 ), controlled by the first index qubit, where the superscript indicates the path index of the previous step.Then the second index qubit is used as the control to prepare the third index qubit and so on.Thus, the index system state preparation can be done in n steps and the lth index qubit is prepared by applying a controlled operator ) l ).An example quantum circuit for preparing the index quantum state for a DSP of n = 3 time steps (levels), each branches to two | 0⟩ Step 0 Step 1 Step 2 FIG. 3. Quantum circuit to prepare the index system for simulating a path-dependent DSP.The figure depicts an example of the first-order Markov chain of two steps each consisting of two realizations.
possibilities, is shown in Fig. 3.
It is also straight-forward to generalize the above procedure to implement path-dependence between events that are more than one time-steps away, by placing the controlled rotation on any two index qubits.Moreover, more complicated path-dependence, such as higher-order Markov chains, can be realized by designing multi-qubit controlled rotations among index qubits.
Besides the index state preparation procedure, the analysis of the DSP follows exactly the same procedure as described in the previous section.

C. State-dependent increments
The increments in a given DSP can also vary in each step, determined by the state in the preceding step.In this case, the probability of a particular realization can be expressed as Intuitively, in order to control the increment given a particular realization of the current step, one needs to apply controlled operation to the index register, controlled by the data register in between each successive step.We leave the explicit details on how to realize the above process with quantum circuits as future work.

The Delta for European call option
An exciting potential application of quantum computing is financial analysis [15,16,28].In particular, the framework developed in this work can be employed to compute the Delta of an European call option.Let S t be stochastic process of an underlying asset then where Φ : R → [0, 1] is the cumulative distribution function (CDF) of the standard normal distribution, K > 0 is the strike price, r is the risk-free interest rate, T − t is the time to maturity and σ is the volatility of the underlying asset.We are interested in the expectation value of the Delta, E [f (S t )] when the underlying asset is described by a geometric Brownian motion, i.e.
S t = S 0 exp µ − σ 2 /2 t + σW t where µ ∈ R and S 0 denote the drift and the starting value, respectively.Indeed, equation ( 20) considers a value reminiscent of the log-return, which is modelled by the Brownian motion and can be approximated as a random walk by Donsker's invariance principle [38,39].The following result establishes the way to implement this on a quantum computer.
Proposition 3. Given n independent and identically distributed random variables X l with P and an initial starting point of with Sn = x 0 + n l=1 X l , we find that when Φ is limited on the interval [−P/2, P/2].Note that we define the sum with prime as the sum over all summands except for l = 0 (see Supplementary Information).
When using the Pauli measurement scheme as explained in Proposition 1, one needs to use V (vx) = R z (2vx) as operators with constants v = 2πl/P for l = −L, . . ., L to evaluate the characteristic function at those points for the Fourier-series approximation.This makes a total of 2L experiments (L measurements for each Pauli observable) with n Hadamard gates on n index qubits for encoding the probability information, n bit-flip gates on the index qubits for implementing the controlled operations that operate if the control qubit is 0, and 2n controlled-Rz gates, which of each can be further decomposed to two controlled-NOT and two R z gates, applied on a data qubit controlled by the index qubit for implementing x 1 and x 2 .In addition, the initial start value x 0 can be implemented by one R z gate (see Supplementary Information).
To demonstrate the proof-of-principle, we performed classical simulations and experiments of the quantum algorithm for calculating the Delta of an European call option on a IBM quantum device named ibmqx2.As an example, the underlying asset is given with µ = 0, σ = 0.02, r = 0.02, S 0 = 100, t = 1 and the time of maturity T = 10.The Fourier series approximation is performed by choosing L = 100 and P = 100.The standard error mitigation protocol available in qiskit [40] is applied to reduce experimental errors.We also performed classical simulation of the quantum algorithm with the standard noise model provided in qiskit.These results are compared with the theoretical values in Fig. 4. The comparison shows that the noise model provided in qiskit explains the experimental error reasonably well (see Methods for a brief comment on the remaining difference between two results).
Since Φ is smooth, the Fourier-series coefficients exponentially decays with respect to l.However, the choice of the approximation scheme for Φ comes with the consequence that its approximation Φ is periodic and continuous at the boundary with Φ(P ) = Φ(−P ) = 0.As P → ∞, the tails extends farther out and flattens to run parallel to the horizontal axis, thereby converging to Φ.This means that for small P depending on the dynamics of S t (and hence Sn ) the approximation may be inaccurate.Therefore, one should keep in mind the minimum/maximum values that Sn can attain and choose P accordingly.
FIG. 4. Evaluation of the Delta for European call option.(left) The Delta of an European call option is calculated for various strike prices.The underlying asset is defined with µ = 0, σ = 0.02, r = 0.02, S0 = 100, t = 1 and the time of maturity T = 10.The dotted black line is the true evaluation of the Delta.The green (orange) solid line is the real (imaginary) part of the theoretical Delta calculation obtained by a Fourier approximation with P = 100 and L = 100, which serves as the reference for the experimental validation.The red (brown) crosses are the real (imaginary) part of the Delta calculated with experiment on the IBM quantum computer with error mitigation applied.The dots are the simulation with noise model provided in qiskit with the same error mitigation applied.(right) This plot shows the characteristic function calculated in theory (×), by simulation with noise and error mitigation (dot) and the IBM quantum experiment with error mitigation (tri-down) for the example strike price of K = 110.

Correlated random walks
As another example, we consider a family of DSPs that strictly requires the ability to simulate path-dependent increments, for which the quantum advantage against classical methods manifests.We demonstrate the simulation of correlated random walks as one such example in the following section.
Correlated random walk (CRW) is a mathematical model that describes discrete random processes with correlations between successive steps [41].It has been a useful tool to study biological processes [42,43], and can also be used to approximate fractional Brownian motion [44,45], which has broad applications for example in mathematical finance [46][47][48] and data network [49][50][51][52].The correlation, often referred to as persistence, results in a local directional bias as the walk moves.More precisely, the CRW denoted by S n = n l=0 X l with n discrete random variables X l , l = 1, . . ., n and persistence parameters p l ∈ [0, 1] and q l ∈ [0, 1] has the following properties: (1) properties can be incorporated in quantum simulations easily by following the same procedure used in the previous example.Given p l and q l , the third property can be implemented as follows.First, all index qubits except the first one that encodes the probability distribution of X 1 are initialized in |0 .The first index qubit is prepared in (|0 + |1 )/ √ 2 in accordance with the second property.Then a controlled rotation gate is applied from each index ancilla qubit to an index qubit of the successive step.The controlled operation can be expressed as where θ p l = 2 cos −1 ( √ p l ) and θ q l = 2 cos −1 ( √ q l ).For example, after one application of the above controlled operation, the index qubits representing the probability distribution of the first two steps of the CRW is given as The above state shows that P = q 1 as required by the property (3) of the CRW.We demonstrate the proof-of-principle with an example designed as follows.The correlated random walk is given by increments x l1 = 1 and x l2 = −1 with an initial value x 0 = 0.The persistence parameters are p l = (1/2, 2/3, 5/6, 1) and q l = (1/2, 1/3, 1/6, 0) with l = 1, 2, 3, 4. Experiments were performed to calculate the characteristic function for v l = 2πl/P , l = −L, . . ., L with L = 100 and P = 100 on ibmqx2.The standard error mitigation protocol available in qiskit is applied to reduce experimental errors.We also performed classical FIG. 5. Characteristic functions of a correlated random walk.Characteristic functions of a random walk is calculated theoretically (×), by simulation (dot), and by experiment (tri-down).The simulation includes the standard noise model provided in qiskit.The experiment employs the standard error mitigation technique provided in qiskit.The parameters that define the correlated random walk is described in the main text.
simulation of the quantum algorithm with the standard noise model provided by qiskit.These results are compared with the theoretical values in Fig. 5.

III. DISCUSSION
We presented a quantum-classical hybrid framework for estimating an expectation value of any integrable function of a random variable in discrete stochastic process with independent increments.As the main ingredient of the framework, we developed a quantum algorithm for efficiently calculating point evaluations of the characteristic function of a random variable, which may also lead to other interesting applications since the probability distribution of a random variable can be completely defined by its characteristic function.More specifically, in the quantum part, the framework proposes a succinct representation of a classical DSP with independent increments as a quantum state in the form of equation (2).The joint probability and the value of each realization are encoded in an entangled state |Ψ f , and therefore all necessary information about a DSP is present.We also detailed the construction of a quantum circuit for preparing the quantum state |Ψ f using the number of circuit elements that only grows linearly with the total number of time steps.For a DSP of n total steps each consisting of k possibilities, there are k n paths.Such process can be encoded in a quantum state using n k-dimensional index qudits (or log 2 (k) qubits) and nk controlled gates.There is no need for sampling strategies as all realizations exist in quantum superposition.This fact is exploited in the steps afterwards.
Two different measurement schemes are proposed for the estimation of expectation values E [cos(S n )] and E [sin(S n )].The first scheme is to measure an expectation value of σ x and σ y directly on the data system of |Ψ f , resulting in a convergence error of O(1/ √ N ) for N repeated experiments for any N > 0. This shows that the quantum brute-force approach acquires the optimal importance distribution.Furthermore, the sampling convergence can be improved by utilizing the quantum amplitude estimation technique, which promises to reduce the approximation error by a factor of O(1/2 m ) using m additional qubits.However, the resource overhead of the latter approach is m ancilla qubits, O(poly(m)) Groverlike controlled operations, and O(m 2 ) one-and two-qubit gates for implementing QFT.Therefore, with the noisy intermediate-scale quantum (NISQ) devices, it may be desirable to use the former approach.
The advantage of this algorithm lies in two facts.First, it is a quantum-classical hybrid computation and thus is a viable candidate to be solved with near-term quantum devices.One can envisage multiple near-term devices running in parallel to calculate all 2L terms independently at the same time.Similar parallelization is also suitable for multiple partitions of a large quantum devices where the qubit connectivity is high within each partition but low among partitions.Second, given one type of process S n , we can pre-compute a number of evaluations of the characteristic functions ϕ Sn (±v i ) for With those evaluations at hand, it is possible to assemble the Fourier-series ondemand when given the Fourier-coefficients of a function f .This approach makes it possible to invest resources to maximize the precision of the characteristic function evaluation for regularly used DSPs.Furthermore, our framework promotes the idea for a co-design (specialpurpose) quantum computer [53][54][55][56], which focuses on optimizations and design decisions at the hardware level to particularly support the DSP simulations at hand.We underscore our findings with two interesting examples.First, we showed an application to finance, the calculation of the Delta of a European call option.The key idea was to model the stochastic behavior of an underlying asset as a Brownian motion, which is then approximated as a random walk by Donsker's invariance principle.Next, an application to DSPs with pathdependent increments is demonstrated by an example of correlated random walks.We performed and presented the results of proof-of-principle experiments for each example to demonstrate the validity and the feasibility of our method.
The framework can be extended to multi-variate func-tions f : R d → R. Given P 1 , . . ., P d ∈ R, the period of the respective argument of f and L 1 , . . ., L d ∈ N + , a multi-dimensional Fourier-series approximation xi can be applied like equation (11) and consequently equation (12).The primary harmonics to calculate are thus ϕ Sn (v) = E e iv•Sn , which can be achieved by increasing the dimension of the index systems and applying a new set of a unitary family V r (r = 1, . . ., d) to the same data system D.
The ability to simulate multi-variate stochastic dynamics with quadratic quantum speed-up without requiring any sampling strategies is highly beneficial for the study of discrete stochastic processes.The path-dependent and state-dependent DSP simulations are excellent fits for studying random walks with internal states [57], and discrete processes that converge to Ornstein-Uhlenbeck [58] and fractional Brownian motion [44].Due to the broad applicability of these mathematical models, the framework developed in this work presents tremendous opportunities for solving problems that arise from various disciplines, such as physics, biology, epidemiology, hydrology, engineering, and finance.
In future work, we plan to provide explicit quantum circuit design for simulating state-dependent DSPs.Interesting applications to accompany the future work include the analysis of stochastic epidemic models [59,60] and of hot streak [61].

METHODS
All experiments are performed using one of the publicly available IBM quantum devices consisting of five superconducting qubits, named as ibmqx2.In order to fully utilize the IBM quantum cloud platform, we used the IBM quantum information science kit (qiskit) framework [40].The versions-as defined by PyPi version numbers-used for this work were 0.20.0.
Superconducting quantum computing devices that are currently available via the cloud service, such as those used in this work, have limited coupling between qubits.Resolving coupling constraints as well as optimizations are done in qiskit with a preset of so-called pass managers.The optimization level ranges from 0 to 3. As we chose the ibmqx2 for our experiments the mapping of data register and index register follow the connectivity.For the family of devices that contains ibmq ourense, one swap operation must be used to exchange the physical qubit 3 and 4. For both layouts see Supp.Info. Figure 3a and 3b.During the compilation step of the experiments we first use a level 0 pass and then a subsequent level 3 pass to optimize and resolve connectivity constraints.To reduce experimental errors, we used the standard error mitigation functionality of qiskit.This requires an extra set of experiments to be executed prior to the main experiment.
In order to understand the source of experimental discrepancy, the experimental results were compared to simulation results that take a realistic noise model into account.During the execution of an experiment, the current device parameters were gathered and stored.Upon completion, a simulation was executed with the standard qiskit noise model, also applying error mitigation on the result.The remaining discrepancy between experimental and simulation results can be attributed to errors that are not included in the basic error model, such as various cross-talk effects, drift, and non-Markovian noise.The standard noise model is described in the supplementary information of Ref. [62] in detail.
The example for the Delta applied the strike price K from 10 to 220 in increments of 5 for the theoretical calculation, while the experiments were performed for K = 25, 55,85,105,110,115,120,125,130,160,190,220.The experiment for each K is executed for characteristic function evaluations at v l = 2πl/P with l = −L, . . ., L, L = 100 and P = 100, each with 8129 shots.For the correlated random walk experiment, we used 32786 shots for each of the evaluations of v l = 2πl/P with l = −L, . . ., L and P = L = 100.Note that this example did not use the Fourier approximation, and P has been chose to be the same as L so that v l goes through one period on each side, resulting in a symmetric picture.As the IBM quantum cloud platform allow for 8192 shots per execution, we created multiple identical experiments and manually added the results.
The measurement scheme of AE now will estimate a value a = Ψ 1 |Ψ 1 = j∈K n p 2 (j) sin 2 ( 1 2 n i=1 x i,ji ), by trigonometric identities.This is equivalent to .
Thus E [cos(S n )] = 1 − 2a, which was to be shown.The other result is given when setting V (x i,ji ) = R y (−x i,ji ) and having the initial state |ψ = R y (π/2).This results in a operator U (j) = R y (π/2 − n i=1 x i,ji ), and hence negates the angle and adds a pre-factor.Thus As a result, This concludes the proof: E [sin(S n )] = 1 − 2a .

IV. AMPLITUDE ESTIMATION
The technique of amplitude estimation (AE) uses phase estimation to find the amplitude of a state.It's use was originated by quantum search algorithms [1] in which the technique of Grover's search is extended to amplitude amplification to the "good" states.In short, say A is an algorithm/operator creating a final state A |0 = |Ψ f = |Ψ 0 + |Ψ 1 which is separated in two orthogonal subspaces.The operator applied is similar to the revert around the mean of Grover's search x i,ji |j |1 .
Therefore, the binary function is very simple; if we interpret x = j 1 • • • j n j n+1 as a binary number with the least significant bit being j n+1 , then we find the correct function to be f (x) = 1, x even 0, x odd .
The value to be under scrutiny is a = Ψ 1 |Ψ 1 and in particular sin 2 θ a = a.It is known that the eigenvalues of Q are λ ± = e ±i2θa .Then we define the operator Λ M (U ) : |j |y → |j (U j |y ) for 0 ≤ j < M .Given an M > 0, usually a power of two, the algorithm uses m = log 2 M qubits: where F m is the quantum Fourier transform.The implementation of Λ M (Q) is usually done by controlled operators with Π jl as given in equation ( 16 Proof.We will establish the proof by a simple inductive argument.To see that equation (S8) coincides with equation (3) of the main text, we look into the first level evolution: In summary, k orthogonal subspaces have been popularized by the application of A 1 and on each such subspace the data systems undergoes an evolution of V (x 1j ).Now, from earlier results, we know that exp −t 2 can be fourier-approximated as exp −t 2 ≈ a l cos(lπt/P ) with a l calculated in the previous section.Thus, the expectation value that we need to calculate becomes The expectation values in the second term of the above equation can be estimated by the standard procedure introduced in our work.On the other hand, the first term can be further approximated using the Fourier series, and the problem of finding the Fourier-coefficients depends on the function y.

IX. SIMULATION OF THE DELTA
We want to simulate with S t being a geometric Brownian motion as given in the main text.The simulation on a real quantum device, such as ibmqx2 (Suppl.Figure 3) puts the restriction on the maximum number of time steps one can simulate due to the number of qubits that are available.With this five-qubit device for example, one can study the random walk S4 = x 0 + X 1 + X 2 + X 3 + X 4 with X l i.i.d.given as in proposition 4 of the main text.The circuit to calculate the characteristic function ϕ S4 (v) is given by Suppl. Figure 1 and is therefore equal to the circuit to the numerical calculation example of proposition 3 of the main manuscript except the choices of x 0 and x lj (l = 1, 2, 3, 4, j = 1, 2).We calculate the characteristic function at L = 5 positions, so v 1 = π 10 , v 2 = π 5 , v 3 = 3π 10 , v 4 = 2π 5 , v 5 = π 2 .For each position, two expectation value measurements are needed, i.e.M = σ x and with M = σ y .We find that by setting together the coefficients and the characteristic function evaluations The corresponding quantum circuit is shown in Figure 1 and shows two important facts.The first is that the index registers are all initialized by a Hadamard gate to put each qubit in an equal superposition state.This corresponds to the probability of occurrence of each outcome.The second is, that controlled R z gates are easily implemented on a real quantum device, for example, a five-qubit quantum device ibmqx2 provided by IBM via the cloud, using two phase gates (these are usually virtual) and only two cnot gates each (refer to Figure 2).Looking at the connectivity of this particular quantum processor unit (QPU) as an example, refer to Figure 3, this also means that swap operations are not necessary.In total, there are 16 controlled-NOT gates, 4 single-qubit bit-flip (X) gates, 4 Hadamard gates and 16 Rz gates in use.

) and p 2 l
,j l = P [X l = x l,j l ].This structure allows to partition I into subsystems, so-called level-index (sub)systems: H I = n l=1 H I l , where each H I l = C k represents a qudit Hilbert space.Let |α(n) ∈ H I be the t e x i t s h a 1 _ b a s e 6 4 = " e o i 8 q S 2 1 3 b L g n Z 8 D u 8 e N A Y r / 4 Y b / 4 b u 7 A H B S d p M p l 5 L 2 8 6 Q c y Z N q 7 7 7 R T W 1 j c 2 t 4 r b p Z 3 d v f 2 D 8 u F R W 0 e J I r R F I h 6 p T o A 1 5 U z S l m G G 0 0 6 s K B Y B p 4 / B + C b z H y d U a R b J B z O N q S / w U L K Q E W y s 5 P c E N i O C e X o 3 6 9 f 6 5 Y p b d e d A q 8 T L S Q V y N P v l r 9 4 g I o m g 0 h C O t e 5 6 b m z 8 F C v D C K e z U i / R N M Z k j I e 0 a 6 n E g m o / n Y e e o T O r D 4 h l d 4 Q x P 0 g t 7 R x 2 K 0 g P K d Y / g D 9 P k D r J y S C w = = < / l a t e x i t > I n < l a t e x i t s h a 1 _ b a s e 6 4 = " G E p c Y B 6 N p 8 g p o f v Z O w 6 9 1 a x a 6 / s

FIG. 2 .
FIG. 2. Operator Tree.(a) The operator tree spanned by a binary example (k = 2) of V for three time steps with the edges named.Each edge represents a state transformation from a parent node to a child node, and each node represents a quantum state.Three states at the final leaves starting from an initial state |ψ at the root node are explicitly shown as an example.(b) The path j = (2, 1, 2) ∈ K 3 and the operator U (j) = V32V21V12 are shown as an example.