A Quantum algorithm for linear PDEs arising in Finance
AA QUANTUM ALGORITHM FOR LINEAR PDES ARISING IN FINANCE
FILIPE FONTANELA, ANTOINE JACQUIER, AND MUGAD OUMGARI
Abstract.
We propose a hybrid quantum-classical algorithm, originated from quantum chemistry, to priceEuropean and Asian options in the Black-Scholes model. Our approach is based on the equivalence betweenthe pricing partial differential equation and the Schr¨odinger equation in imaginary time. We devise astrategy to build a shallow quantum circuit approximation to this equation, only requiring few qubits. Thisconstitutes a promising candidate for the application of Quantum Computing techniques (with large numberof qubits affected by noise) in Quantitative Finance. Introduction
Pricing financial derivatives accurately and efficiently is one of the most difficult and exciting challengesin Mathematical Finance. While the existence of closed-form solutions are known in very few simple cases,the vast majority of models do not admit one, and computing derivatives rely on computational techniques,either using Monte Carlo simulation methods or via numerical methods for partial differential equations, suchas finite differences or finite elements. Unfortunately, these techniques can rapidly become computationallyintensive when the model becomes complicated, or when the dimension of the problem (in the case of Basketoptions for example) becomes large. While huge progress has been made over the past decades, thesetechniques have inherent limitations which cannot be overcome when run on classical computers.Recently, the emergence of small-scale quantum computers has caught the eyes of mathematicians andfinancial engineers [15, 9] as a a potential goose laying golden eggs. Indeed, while Quantum Computinghas been around since Benioff [1], Deutsch [7], Feynman [8] and Manin [12] suggested that a quantumcomputer could perform tasks out of reach for classical computers, it really started getting some tractionwhen Shor [17] unearthed a polynomial-time quantum algorithm to factor integers. Since then, huge effortshave been made to actually build quantum hardware; only recently though have researchers, in a partnershipbetween Google AI and the US NASA, managed to stabilise a quantum chip over a short period of time,effectively starting what has already been called the Quantum Revolution. Several companies are nowproposing online quantum capabilities, using a few qubits. We are still far from being able to use thosein production mode on a large scale, but the revolution is fast marching, and tools are urgently needed toembrace it. Quantitative Finance is traditionally quick to respond to such calls, and several attempts haverecently been made to use quantum techniques for option pricing, in particular [16, 18], focusing on theso-called Amplitude Estimation algorithm [4], which allows a quadratic speed-up for simulation methodscompared to classical Monte Carlo schemes.We propose a new route to investigate the use of Quantum techniques in Quantitative Finance, and developa hybrid quantum-classical algorithm to solve the Schr¨odinger equation, which by simple transformationsis equivalent to the PDE satisfied by the option price. The paper is organised as follows. In Section 2,we introduce the main tools and notations, clarifying the link between the Black-Scholes PDE [3] and itsSchr¨odinger counterpart along the imaginary axis, while providing a short reminder on Quantum mechanics.In Section 3, we borrow an idea developed in quantum chemistry [10], and propose an algorithm to solvethe aforementioned Schr¨odinger equation using a quantum computer. The methodology develops a hybridalgorithm, where part of the computations can be run on a quantum computer, while the remaining tasks aresolved on a classical machine. We describe in detail in Section 4 the actual implementation of the algorithm,
Date : December 6, 2019.2010
Mathematics Subject Classification.
Key words and phrases. quantum algorithms, option pricing, PDE, Schr¨odinger equation.The views and opinions expressed here are the authors’ and do not represent the opinions of their employers. They are notresponsible for any use that may be made of these contents. No part of this presentation is intended to influence investmentdecisions or promote any product or service. The authors would like to thank Alexei Kondratyev for stimulating discussions. a r X i v : . [ q -f i n . C P ] D ec FILIPE FONTANELA, ANTOINE JACQUIER, AND MUGAD OUMGARI with a particular emphasis on the quantum circuit. The numerical results are presented in Section 5, andwe show that the shallow quantum circuit developed here is able to price European and arithmetic Asianoptions in the Black-Scholes model with very good accuracy.2.
Black-Scholes and the Schr¨odinger equation
The Black-Scholes model is at the core of financial modelling and assumes that, under a given risk-neutralmeasure, the underlying stock price process ( S t ) t ≥ satisfies the stochastic differential equation(2.1) d S t S t = r d t + σ d W t , for t ≥ , where W is a standard Brownian motion on a given filtered probability space, r ∈ R is the instantaneousrisk-free rate, and σ > Pricing PDEs.
For a European Call option with payoff f ( S T ) = max( S T − K,
0) at maturity T , theFeynman-Kac formulation allows us to write the option price V ( t, s ) := E (cid:2) e − r ( T − t ) f ( S ( T ) | S t = s (cid:3) , for any t ∈ [0 , T ], obtained by no-arbitrage arguments, as the unique smooth solution to the PDE(2.2) (cid:18) ∂ t + σ s ∂ ss + rS∂ s − r (cid:19) V ( t, s ) = 0 , for all s > t ∈ [0 , T ), with terminal condition V ( T, s ) = f ( s ). The changes of variables S = e x and τ = σ ( T − t ) reduce this partial differential equation to the heat equation(2.3) ∂ τ u ( τ, x ) = 12 ∂ xx u ( τ, x ) , on (0 , σ T ] × R , with boundary condition u (0 , x ) = e − ax f (e x ), where b := σ (cid:0) − r σ − σ − r (cid:1) and a := − rσ . Another common financial derivative is the fixed-strike arithmetic Asian Call option, withterminal payoff(2.4) max (cid:32) T (cid:90) T S ( t )d t − K, (cid:33) , for some strike K >
0. Setting Y t := (cid:0) B t − K e − r ( T − t ) (cid:1) /S t with B t = q ( t ) S t + e − r ( T − t ) T (cid:90) t S u d u and q ( t ) = − e − r ( T − t ) rT , if r (cid:54) = 0 , − tT , if r = 0 , Brown [5] and Vecer [19] showed that the price at time t ∈ [0 , T ) of the option reads V ( t, S t ) = S t (cid:101) Q ( t, Y t ).With τ := σ ( T − t ) and Q ( τ, y ) = (cid:101) Q ( t, y ), the function Q uniquely solves the PDE(2.5) ∂ τ Q ( τ, y ) = ( q ( τ ) − y ) ∂ yy Q ( τ, y ) , for ( τ, y ) ∈ (0 , σ T ] × R with terminal condition Q (0 , y ) = max( y, y ↓−∞ Q ( τ, y ) = 0 and Q ( τ, q (0)) = q (0). One should note that (2.5) is very similar to the heatequation (2.3), albeit with an additional time-dependence on the right-hand side. QUANTUM ALGORITHM FOR LINEAR PDES ARISING IN FINANCE 3
Schr¨odinger’s formulation.
We now translate the two pricing PDEs (2.3) and (2.5) into a linearSchr¨odinger equation, classical in quantum mechanics. The Wick rotation ξ = − i τ , where i is the imaginaryunit, transforms (2.3) into − i ∂ ξ u ( ξ, x ) = 12 ∂ xx u ( ξ, x ) , which, in the quantum computing literature, is often written using Dirac’s notation (see Section 2.3) as(2.6) − i ∂∂ξ | ψ (cid:105) = (cid:98) H | ψ (cid:105) , where the wave function ψ (and its associated quantum state | ψ (cid:105) ) plays the role of the modified derivativeprice u ( · , · ), while the Hamiltonian operator is simply (cid:98) H = ∂ xx . For this special case of a time-independentHamiltonian, the solution to (2.6) is given explicitly by(2.7) | ψ ( ξ ) (cid:105) = exp (cid:16) i (cid:98) H ξ (cid:17) | ψ (0) (cid:105) , where exp (cid:16) i (cid:98) H ξ (cid:17) is the time evolution operator, and | ψ (0) (cid:105) a normalised initial state with (cid:104) ψ (0) | ψ (0) (cid:105) = 1.For Asian options, the Wick rotation ξ = − i τ now turns (2.5) into − i ∂ ξ Q ( τ, y ) = ( q ( τ ) − Y ) ∂ yy Q ( τ, y ) , which reads, in Dirac’s notations,(2.8) − i ∂∂ξ | ψ (cid:105) = (cid:98) H ( ξ ) | ψ (cid:105) . The main difference between (2.6) and (2.8) is the time-dependent Hamiltonian. Nonetheless, the solutioncan be written by computing the new evolution operator as(2.9) | ψ ( ξ ) (cid:105) = exp (cid:32) i (cid:90) ξ (cid:98) H ( χ )d χ (cid:33) | ψ (0) (cid:105) . Reminder on quantum notations.
In order to facilitate the integration of Quantum Mechanics toolsinto the realm of Quantitative Finance, we shall endeavour to combine notations from both fields in a clearand consistent manner. To do so, we recall some standard Quantum notations, originating from Dirac, whichwe will use throughout the paper. In a given (complex-valued) Hilbert space H , a vector v is representedvia the ket notation | v (cid:105) ∈ H . For u ∈ H , the bra (cid:104) u | belongs to the dual space H ∗ , so is a linear map from H to C , and we denote the action of u on v as the bracket (cid:104) u | v (cid:105) . In classical linear algebra notations, this canbe recast in the following form: suppose that H is of dimension n ∈ N , then u , v ∈ H can be represented asvectors in R n as u = u ... u n and v = v ... v n , or u = ( u , . . . , u n ) (cid:62) and v = ( v , . . . , v n ) (cid:62) . In that case, we can write | v (cid:105) = v ... v n , (cid:104) u | = ( u ∗ , . . . , u ∗ n ) , and (cid:104) u | v (cid:105) = n (cid:88) i =1 u ∗ i v i , where ∗ denotes complex conjugacy. A physical state in quantum mechanics is represented by a state vector,or in fact a ket. The general quantum state of a qubit, the basic unit of quantum information, is a linearsuperposition of its orthonormal basis. In particular, a single qubit state can be written as | ψ (cid:105) = α | (cid:105) + β | (cid:105) ,for some α, β ∈ C satisfying | α | + | β | = 1, with ( | (cid:105) , | (cid:105) ) the orthonormal basis. In classical linear algebranotations, this boils down to ψ = (cid:18) ψ ψ (cid:19) = α (cid:18) (cid:19) + β (cid:18) (cid:19) . FILIPE FONTANELA, ANTOINE JACQUIER, AND MUGAD OUMGARI
Generally speaking, an n-qubit quantum state corresponds to a vector in C n . For such a state | ψ (cid:105) , that is,a physical system which can be in n different, mutually exclusive classical states | (cid:105) , | (cid:105) , . . . , | n − (cid:105) , so that | ψ (cid:105) = α | (cid:105) + · · · + α n − | n − (cid:105) , for ( α , . . . , α n − ) ∈ C n , such that (cid:80) ni =1 | α i | = 1. The basis ( | (cid:105) , | (cid:105) , . . . , | n − (cid:105) ) forms an orthonormalbasis of an n -dimensional Hilbert space H n . Given H n and H m , we can define the tensor product H := H n ⊗ H m as the nm -dimensional Hilbert space spanned by {| i (cid:105) ⊗ | j (cid:105) : i = 0 , . . . , n − , j = 0 , . . . , m − } . For example,a 2-qubit system, corresponding to a Hilbert space of dimension 4 can be viewed as the tensor product of twoHilbert spaces, each of dimension 2, and its basis being spanned by {| (cid:105) ⊗ | (cid:105) , | (cid:105) ⊗ | (cid:105) , | (cid:105) ⊗ | (cid:105) , | (cid:105) ⊗ | (cid:105)} ,which can also be written as {| (cid:105) , | (cid:105) , | (cid:105) , | (cid:105)} .Quantum logic gates are reversible quantum circuits operating on quantum states, which can be repre-sented as unitary matrices. We shall here make use of several particular gates repeatedly, which can be easilyrepresented in a 1-qubit state as the following 2 × X : X-Pauli gate: (cid:18) (cid:19) Y : Y-Pauli gate: (cid:18) − ii 0 (cid:19) H : Hadamard gate: 1 √ (cid:18) − (cid:19) R y ( θ ) : Rotation gate with angle θ : (cid:18) cos (cid:0) θ (cid:1) − sin (cid:0) θ (cid:1) sin (cid:0) θ (cid:1) cos (cid:0) θ (cid:1) (cid:19) Quantum Imaginary Time Evolution
We now describe a strategy to solve the linear Schr¨odinger equation along its imaginary time axis usinga quantum computer. Our approach relies on the hybrid algorithm developed in [10], where part of thecomputation is performed on a quantum computer, and part on a classical machine. Consider first a time-independent Hamiltonian (cid:98) H with evolution operator (or propagator) exp(i (cid:98) H ξ ) evolving along real values ξ .In this case, the propagator can be implemented efficiently on a quantum computer by means of a Trotterdecomposition, since exp(i (cid:98) H ξ ) is a unitary matrix [14]. Along the imaginary axis however, the correspondingevolution operator exp( (cid:98) H τ ) is represented by a non-unitary matrix; It can be easily simulated on a classicalcomputer, however as the dimension of the wave function grows exponentially, it rapidly becomes infeasible;its Trotter decomposition using unitary gates is not straightforward, making its implementation on a quantumcomputer more challenging. We follow instead a recent idea [10] to solve an equivalent normalised imaginarytime evolution(3.1) | ψ ( τ ) (cid:105) = γ ( τ ) e (cid:98) H τ | ψ (0) (cid:105) , with γ ( τ ) := (cid:16) (cid:104) ψ (0) | e (cid:98) H τ | ψ (0) (cid:105) (cid:17) − / , indirectly. The parameter γ ( τ ) is a normalisation constant, and (cid:98) H a time-independent Hamiltonian. Wecall this approach indirect since | ψ ( τ ) (cid:105) is not computed by constructing the operator exp( (cid:98) H τ ) on a quantumcomputer. Instead, we approximate | ψ ( τ ) (cid:105) by a quantum circuit composed of a sequence of parameterisedgates such that | ψ ( τ ) (cid:105) ≈ | φ ( θ τ ) (cid:105) , where θ τ = ( θ τ , · · · , θ Nτ ) (cid:62) ∈ R N is a vector of time-dependent parameters.Therefore, knowing how to describe the evolution of θ τ makes it possible to reconstruct the imaginarySchr¨odinger time evolution | ψ ( τ ) (cid:105) in (3.1). We shall refer to the approximation | φ ( θ τ ) (cid:105) as the ansatz circuit.Assume that our initial quantum state is | ψ init (cid:105) , so that the ansatz at time τ is | φ ( θ τ ) (cid:105) = Φ( θ τ ) | ψ init (cid:105) ,where Φ( θ τ ) is sequence of unitary gates Φ( θ τ ) = S (cid:0) U N ( θ Nτ ) , . . . , U k ( θ kτ ) , . . . , U ( θ τ ) (cid:1) which we willspecify later. In this paper, each parameterised gate U k ( θ k ) is considered as a rotation or a controlledrotation gate. It is also possible to show that the evolution of θ τ , and thus the Schr¨odinger imaginary timedynamics, can be computed using McLachlan’s variational principle [11](3.2) δ (cid:13)(cid:13)(cid:13)(cid:16) ∂ τ + (cid:98) H (cid:17) | ψ ( τ ) (cid:105) (cid:13)(cid:13)(cid:13) = 0 , where (cid:107) v (cid:107) := (cid:104) v | v (cid:105) , and δ denotes the infinitesimal variation. This determines the values of θ τ minimisingthe distance (cid:107) | ψ ( τ ) (cid:105) − | φ ( θ τ ) (cid:105) (cid:107) . The variational principle (3.2) results in the system of first-order ODEs [10](3.3) A( τ ) ˙ θ τ = C τ , for each τ, QUANTUM ALGORITHM FOR LINEAR PDES ARISING IN FINANCE 5 where ˙ θ τ = ∂ τ θ τ , and the matrix A( τ ) and the vector C( τ ) readA( τ ) = (cid:18) (cid:60) (cid:18) ∂ (cid:104) φ ( τ ) | ∂θ i ∂ | φ ( τ ) (cid:105) ∂θ j (cid:19)(cid:19) i,j =1 ,...,N and C( τ ) = (cid:18) (cid:60) (cid:18) ∂ (cid:104) φ ( τ ) | ∂θ i (cid:98) H | φ ( τ ) (cid:105) (cid:19)(cid:19) i =1 ,...,N . One of the main advantages of this method is that both A and C can be measured efficiently using aquantum computer [2, 10]. In order to build the hybrid classical-quantum scheme, we rely on the followingassumptions:
Assumption 3.1. (i) Every unitary gate in the algorithm depends on a single parameter.(ii) The decomposition (cid:98) H = N (cid:88) i =1 λ i h i holds, where λ i are real coefficients and h i are tensor products of Pauli matrices.Assumption 3.1(i) is in fact purely for convenience, as unitary gates with multiple parameters can beeasily decomposed into single-parameter ones. It in particular allows us to write, for every k = 1 . . . , N , thederivative of the unitary gate U k as ∂ θ k U k ( θ k ) = N (cid:88) i =1 f k,i U k ( θ k ) σ k,i , where σ k,i is a single-qubit or a two-qubit unitary operator, while f k,i is a scalar parameter, which in turnimplies ∂ θ k | φ ( τ ) (cid:105) = N (cid:88) i =1 f k,i (cid:101) Φ k,i | ψ init (cid:105) , where (cid:101) Φ k,i = S (cid:0) U N ( θ N ) , . . . , U k ( θ k ) σ k,i , . . . , U ( θ ) (cid:1) . For example, if U k ( θ k ) is a single-qubit R y ( θ k ) rotation gate such that U k ( θ k ) = exp (cid:0) − i θ k σ Y (cid:1) , where σ Y is the Y-Pauli matrix, the derivative is simply ∂ θ k U k ( θ k ) = − i σ Y U k ( θ k ). Therefore, the state ∂ θ k | φ ( τ ) (cid:105) canbe prepared by adding the extra σ Y gate, together with a constant factor − i2 only, such that ∂ θ k | φ ( τ ) (cid:105) = − i2 (cid:101) Φ k,k ( θ ) | ψ init (cid:105) , where (cid:101) Φ k,k ( θ ) = S (U N ( θ N ) , . . . , U k ( θ k ) σ Y , · · · , U ( θ )). In general, the matrix A is thencomputed as [10](3.4) A ij = (cid:60) (cid:88) k,l =1 ,...,N f ∗ k,i f l,j (cid:104) ψ init | (cid:101) Φ † k,i (cid:101) Φ l,j | ψ init (cid:105) , for i, j = 1 , . . . , N, where the dagger † denotes the complex conjugate of the transpose. Moreover, an equivalent circuit designedto measure the matrix A using a quantum computer, directly, is discussed in [2, 10]. We will also present anexample of a quantum circuit to measure an element of A in Section 4.For the vector C, we use Assumption 3.1(ii); however, this decomposition usually scales polynomially withthe system size, since (cid:98) H is a sparse matrix obtained from the discretisation of a differential operator. Thecorresponding vector C then reads [10]C i = (cid:60) (cid:88) k,j =1 ,...,N f ∗ k,i λ j (cid:104) ψ init | (cid:101) Φ † k,i h j Φ | ψ init (cid:105) , for i = 1 , . . . , N. Once the matrix A and the vector C are obtained, the time evolution can be computed numerically using aclassical computer. We suggest here an Euler scheme, as in [10], and the evolution of θ τ +∆ τ is calculated as(3.5) θ τ +∆ τ = θ τ + ∆ τ ˙ θ τ = θ τ + ∆ τ A( τ ) − C( τ ) , for some small time step ∆ τ . Usually, the matrix A is not well-conditioned, and we solve, for each τ (3.6) arg min (cid:110)(cid:13)(cid:13)(cid:13) A( τ ) ˙ θ τ − C( τ ) (cid:13)(cid:13)(cid:13) : θ τ ∈ R N (cid:111) , instead of (3.3) directly, assuming some small cut-off ratio for the eigenvalues of A. This computation is straightforward in
Python using the package numpy.linalg.lstsq for example.
FILIPE FONTANELA, ANTOINE JACQUIER, AND MUGAD OUMGARI
For a time-dependent Hamiltonian, for Asian options, we are interested in determining the quantum state(3.7) | ψ ( τ ) (cid:105) = γ ( τ ) exp (cid:18)(cid:90) τ (cid:98) H ( χ )d χ (cid:19) | ψ (0) (cid:105) , with γ ( τ ) := (cid:18) (cid:104) ψ (0) | exp (cid:26) (cid:90) τ (cid:98) H ( χ )d χ (cid:27) | ψ (0) (cid:105) (cid:19) − instead of (3.1), with the new normalisation constant. In order to avoid the integral calculation in (3.7),we use, along a left-point approximation when discretising of the integral, effectively freezing the time-dependent Hamiltonian (cid:98) H ( τ ) on every time interval; thus the only difference with the time-independentframework above is the update of the Hamiltonian at every time step. The main drawback of the presentapproach is the requirement for an ansatz circuit. Ideally, the latter has to be complex enough to approximatethe underlying quantum state | ψ ( τ ) (cid:105) , but not unnecessarily deep in order to avoid the expensive computation(on a classical computer) of the linear problem (3.6), i.e. too many parameters θ . In the examples below, wedesign a quantum circuit strategy for simple financial derivatives, and obtain promising results for relativelyshallow quantum circuits. 4. Hybrid Quantum-Classical Algorithm
We now design the quantum imaginary time algorithm described in Section 3 to price financial derivatives.The main challenge is to design an ansatz circuit able to represent the price of the underlying option fromexpiry down to the valuation date. We here employ the sequence of gates described by the quantum circuitdepicted in Figure 1. .. | ⟩ ... l .. X .... l . R y ( θ ). | ⟩ .... H .... R y ( θ )... l .. R y ( θ ).... l .. R y ( θ ).. | ⟩ .... H .... R y ( θ )... l .. R y ( θ )..... R y ( θ ).. | ⟩ .... H .... R y ( θ )... l .. R y ( θ )..... R y ( θ )..... R y ( θ ). Figure 1.
Ansatz implemented for European and Asian Call options. The vector θ is composedby the angles of each rotation gate. The notations l , . . . , l indicate the different layers. Remark 4.1.
In Figure 1, a line represents a qubit, whereas a squared box denotes the action of a quantumgate on the qubit on this line. The circuit is to be read from left (inputs) to right (outputs). A controlledgate (for example R ( θ )) is shown with a filled black circle on the control qubit line, and the usual symbolfor the gate is written on the target qubit, with a line connecting the two.The quantum circuit has 4 qubits only, and thus our strategy is able to recover Call prices with a resolutionof 2 = 16 points. The ansatz in Figure 1 has a layer l composed by H and X gates (2.10). In the nextlayer l , a sequence of R y gates is applied to each qubit. Each such gate has a free parameter that needs tobe calculated given some desired initial condition, and these parameters vary in time, as described above.The layers l , l and l , composed by controlled R cy gates, and a new layer l , similar to l , are added inorder to give more flexibility to the approximation, and thus achieve different and more involved values of | φ ( θ τ ) (cid:105) . In the present analysis, we increase or reduce the complexity of the ansatz circuit by adding orremoving new unit cells, i.e. layers l , l , l and l . The optimal configuration is obtained when | φ ( θ (cid:105) isable to represent | ψ (0) (cid:105) with good accuracy, and with as few free parameters as possible. One should notethat the ansatz is composed of gates with rotations along the y axis only. The motivation for using R y gatesis due to its real-valued representation, which is desirable as | ψ ( τ ) (cid:105) in (3.1) is a real number. However, otherapproaches might also be valid, using for instance R x and R z gates, at the cost of computing the real partof | φ ( θ τ ) (cid:105) at the end of the simulations. QUANTUM ALGORITHM FOR LINEAR PDES ARISING IN FINANCE 7
Mathematically, a quantum circuit is described by means of a series of matrix operations; the ansatz inFigure 1 is composed of X , H , and R y single-qubits gates, and their mathematical representations are givenby the matrices in (2.10). The controlled rotation gate R cy is implemented using two qubits, and its matrixrepresentation reads R cy ( θ ) := (cid:0) θ (cid:1) − sin (cid:0) θ (cid:1) (cid:0) θ (cid:1) cos (cid:0) θ (cid:1) = (cid:18) I O O R y ( θ ) (cid:19) , where O denotes the 2 × l in Figure 1, which contains the gates R y ( θ ), R y ( θ ), R y ( θ ) and R y ( θ ), is computed as E l = R y ( θ ) ⊗ R y ( θ ) ⊗ R y ( θ ) ⊗ R y ( θ ) , where ⊗ represents the Kronecker product. The effect of R cy ( θ ) in l is calculated similarly as E l = I ⊗ I ⊗ R cy ( θ ) , with I the 2 × E l and E l is the standard matrix multiplication E l ,l = E l × E l . This process is then repeated until the whole quantum circuit is completed, and the resulting effect of thesix layers in Figure 1 reads E l ,...,l = (cid:89) i =1 E l i . With this circuit, the final output | ψ F (cid:105) of the quantum circuit due to the input | ψ init (cid:105) is then computed as | ψ F (cid:105) = E l ,...,l | ψ init (cid:105) . For example, the input | ψ init (cid:105) of the ansatz circuit in Figure 1 is a quantum state where each of the 4 qubitsis | (cid:105) . Its vectorial representation is thus | ψ init (cid:105) = | (cid:105) ⊗ | (cid:105) ⊗ | (cid:105) ⊗ | (cid:105) Finally, note that if θ = 0, the resulting state | φ ( θ ) (cid:105) becomes a step function, since each gate R y behavesas the identity operator I . As discussed above, the payoff of each financial product needs to be representedby the ansatz | φ ( θ ) (cid:105) before starting the actual simulation. In practice, this representation is implementedby computing the values of θ which best describes the initial conditions of the algorithm, that is(4.1) θ = arg min θ ∈ R N {(cid:107)| φ ( θ ) (cid:105) − | ψ (0) (cid:105)(cid:107)} . The quantum state | ψ (0) (cid:105) is calculated directly from the payoff f ( · ) of each financial derivative. In the caseof an European option, the state | ψ (0) (cid:105) is simply | ψ (0) (cid:105) ≡ γ (0)e − aX f (e X ), where γ (0) is the normalisationconstant guaranteeing (cid:104) ψ (0) | ψ (0) (cid:105) = 1, and a comes from the change of variables from the Black-ScholesPDE to the heat equation in (2.3). Some care should be taken when solving (4.1) since the cost function is notconvex, and algorithms such as Differential Evolution will avoid being trapped in a local minimum. Ideally,the design of the ansatz should be considered together with the optimisation problem (4.1). In practice, ifthe ansatz circuit is not able to represent the payoff | ψ (0) (cid:105) accurately, it becomes pointless to compute thetime-marching simulation | φ ( θ τ ) (cid:105) since the algorithm already starts from the wrong initial conditions. Inthis case, the quantum circuit in Figure 1 needs to be improved for the desired application. The completestrategy to define the depth of the ansatz in Figure 1 and its initial condition θ is summarised in Table 1. FILIPE FONTANELA, ANTOINE JACQUIER, AND MUGAD OUMGARI
Ansatz and initial condition algorithminput ← Payoff function initialisation:
Define the normalised initial condition | ψ (0) (cid:105) Define the maximum acceptable number of unit cells N maxcells Define the maximum acceptable error ε max = (cid:107) | φ ( θ ) (cid:105) − | ψ (0) (cid:105) (cid:107) Define an initial ansatz with N cells unit cells endwhile N cells < N maxcells Compute the optimal θ from (4.1) and the approximation error ε = (cid:107) | φ ( θ ) (cid:105) − | ψ (0) (cid:105) (cid:107) if ε > ε max Increase the number of unit cells N cells in the ansatz else output → Ansatz circuit and initial condition θ endend output → Inform that the algorithm has not converged and that a different ansatz is needed
Table 1.
Algorithm to determine the depth of the ansatz and the initial conditions θ . The main advantage of the hybrid scheme is that the matrix A and the vector C in (3.3) can each bemeasured on a quantum computer. For example, A , can be measured using the quantum circuit in Figure 2... | ⟩ .... H .... X . | ⟩ .... H .... Y ..... X . | ⟩ .... X .... R y ( θ ) .... R y ( θ ) .... R y ( θ ) ..... Y ..... H ..... R y ( θ ) .. | ⟩ .... H .... R y ( θ ) .... R y ( θ ) ..... R y ( θ ) .. | ⟩ .... H .... R y ( θ ) .... R y ( θ ) ..... R y ( θ ) ..... R y ( θ ) . Figure 2.
Quantum circuit able to measure A , in (3.3). The symbol on the top right of theFigure denotes measurement of the quantum state. In summary, we need to add an ancillary qubit to the quantum circuit of Figure 1 together with two H and X gates, and a controlled R cy gate before each respective gate R y ( θ ) and R y ( θ ). In the following, A , is obtained by measuring the expectation value on the added ancillary qubit, as in [2, 10]. The determinationof C is implemented on a quantum computer using a very similar approach. The main difference is thatthe underlying Hamiltonian (cid:98) H first needs to be decomposed according to Assumption 3.1(ii). Table 2 belowsummarises the structure of the hybrid quantum-classical algorithm, where the labels (CC) and (QC) indicatewhich part of the algorithm should be computed on a classical and on a quantum computer respectively. Remark 4.2.
The current methodology only checks the quantum approximation accuracy at the payofflevel, i.e. verifying (cid:107) | φ ( θ ) (cid:105) − | ψ (0)) (cid:105) (cid:107) = O ( ε ), for some tolerance level ε , assuming that the error doesnot propagate much along the grid, namely that (cid:107) | φ ( θ τ ) (cid:105) − | ψ ( τ )) (cid:105) (cid:107) remains of order ε for all τ . Thisassumption may fail in cases where the price | ψ ( τ ) (cid:105) varies significantly. One computational issue of thescheme is the resolution of the high-dimensional minimisation problem (4.1). Other imaginary-time evolution QUANTUM ALGORITHM FOR LINEAR PDES ARISING IN FINANCE 9
Hybrid quantum-classical algorithminput ← Ansatz circuit and initial condition θ initialisation: Define the initial condition for the ansatz, i.e. | φ ( θ ) (cid:105) :Define the Hamiltonian (cid:98) H for the problemDecompose (cid:98) H = (cid:80) i λ i h i endwhile τ < σ T : Compute the matrix A (QC) if (cid:98) H is time-dependent: Update the Hamiltonian (cid:98) H ( τ )Update the decomposition (cid:98) H ( τ ) = (cid:80) i λ i h i end Compute the vector C (QC)Update the new values of θ (CC) end output → Price of the underlying derivative (in the original coordinates)
Table 2.
Hybrid quantum-classical algorithm. algorithms [13] bypassing this optimisation might provide an efficient alternative for large-scale applications.We leave these two issues for further investigations.5.
Numerical Examples
We now employ the imaginary time evolution technique described above to compute the price of twodifferent financial derivatives. The goal is to show that the proposed ansatz circuit in Figure 1 is able toreconstruct the complete price evolution of the underlying financial instruments accurately. The results showthat good accuracy is obtained when the ansatz is composed by the layers in Figure 1. This configurationleads to a quantum circuit, fully described in Appendix A, composed of 25 R y gates, resulting in θ τ ∈ R foreach τ . Finally, we simulate the evolution of an European Call option in the Black-Scholes model using (2.3).The second example is that of a fixed-strike arithmetic Asian Call, which uses the same ansatz circuit butdifferent initial conditions.5.1. European Call Option.
In the Black-Scholes model (2.1), we consider σ = 20%, S = K = 100, T = 1, and for simplicity zero interest rates. We discretise the state space on an equidistant grid boundedby S min = 50 and S max = 150, or x min = ln( S min ) ≈ .
91 and x max = ln( S max ) ≈ .
01. Using only 4 qubits,the discretisation allows the representation of | ψ (cid:105) using 2 = 16 points, where the states | ψ F (cid:105) = | (cid:105) and | ψ F (cid:105) = | (cid:105) represent the solution respectively at x min and x max . The evolution of | φ ( θ τ ) (cid:105) from the expirydown to the initial date is computed using the approach described in Section 3. The Hamiltonian of thesystem can be written as (cid:98) H = ∂ xx (Section 2.2), and we discretise it using second-order finite differences as − b · · · x − x x · · · x − x x · · · · · · x − x x · · · − b , where ∆ x is the discretisation step in space, and the first and top rows of the matrix correspond to thebehaviour at the boundaries x min and x max , where b appears in the change of variables from the Black-Scholes PDE to the heat equation (2.3). Once the operator (cid:98) H is constructed, the evolution of θ τ is obtainedfrom the Euler scheme (3.5). We split the time domain [0 , T ] into 500 equidistant steps, and compute the matrix A and the vector C in (3.5) on this grid. Since the underlying linear system is ill-conditioned, we usethe least squares approximation (3.6) to compute ˙ θ in (3.3) assuming an eigenvalue cutoff ration of 10 − .Figure 3 plots | φ ( θ τ ) (cid:105) using the ansatz circuit in Figure 1 and the difference between this and the priceusing a classical algorithm. The very low errors (cid:107)| ψ ( τ ) (cid:105) − | φ ( θ τ ) (cid:105)(cid:107) show that the proposed ansatz is able toreconstruct the evolution of | ψ ( τ ) (cid:105) very accurately, and is a promising step for the use of quantum technologyin quantitative finance. Figure 3.
European Call option prices. The left plot depicts the corresponding results ob-tained from the quantum ansatz implementation. The right one shows the errors (cid:107)| ψ ( τ ) (cid:105) − | φ ( θ τ ) (cid:105)(cid:107) . Figure 4 compares the expected solutions to the results obtained from the simulation of the quantumalgorithm. The results are depicted in the original option price coordinate system V ( t, S ) of (2.2). Thescaling parameter γ ( τ ) in (3.1) is straightforward to compute for the present case since the values of V ( t, S )are known at the boundary S max . Therefore, it is possible to determine u ( τ, x max ) and, consequently, thescaling parameter between | φ ( θ τ ) (cid:105) and u ( τ, x ) in (3.1) is known. Figure 4 displays the option price at
60 80 100 120 140 S P a y o ff (a) Classical payoffQuantum approximation
60 80 100 120 140 S C a ll p r i c e (b) Classical simulationQuantum approximation
96 98 100 102 104 S P a y o ff
96 98 100 102 104 S C a ll p r i c e Figure 4.
Prices of European Calls using the hybrid quantum-classical algorithm. The left plotshows the price at maturity and the right one at inception. inception ( τ = σ T ). The hybrid quantum-classical algorithm is able to reconstruct the price function withgreat accuracy. We list the initial and final values of θ in Appendix B. QUANTUM ALGORITHM FOR LINEAR PDES ARISING IN FINANCE 11
Asian Call Option.
We use the same parameters as in Section 5.1, and the same ansatz circuit inFigure 1. The system is then discretised between y min = − . y max = 0 .
4, corresponding to S min ≈ . S max ≈ .
67. The underlying time-dependent Hamiltonian (cid:98) H ( τ ) is thus discretised as1∆ y · · · ( q ( τ ) − y ) − ( q ( τ ) − y ) q ( τ ) − y ) · · · · · · ( q ( τ ) − y ) − ( q ( τ ) − y ) q ( τ ) − y ) · · · , where { y min , y , · · · , y , y max } is the discretised space grid and ∆ y the discretisation step. The first andlast rows of the matrix are null since | ψ ( τ ) (cid:105) is assumed constant at the boundaries y min and y max . Figure 6displays the results obtained from the classical and the quantum simulations, in the original price coordinates.Again, since the solution Q ( τ, y ) in (2.5) is known at the boundary y max , the corresponding scaling parameterlinking the classical and the quantum solution is trivial to compute. Figure 6(a) depicts the payoff functions,similar to those of European Call options in Section 5.1; again, the respective errors due to the quantumapproximation of the initial conditions are negligible. The results in Figure 6(b) show the option price atinception ( τ = σ T ). The quantum approximation obtained from the ansatz accurately match the expectedvalues obtained by classical simulation. Finally, the initial and final values of θ are listed in Appendix B. Figure 5.
Asian Call option prices computed with the hybrid algorithm. Conclusions and Outlook
Quantum computers are now part of the daily headlines in technology-related and financial news, andpromise a new area where heavy computations can be carried out in the blink of an eye, relegating slowpricing and onerous calibration to the Middle Ages. While this promise is appealing, we are still far fromit, and actual quantum computers are still in their infancy, still deeply affected by decoherence. On thesoftware side, current quantum algorithms require a large number of qubits together with precise gateimplementations, which may limit their applications in the near future. We put forward here a hybridquantum-classical algorithm, originating in computational chemistry, to price financial derivatives. Thestrategy is based on the PDE representation of the pricing problem, and the link between the latter and itsSchr¨odinger counterpart from quantum mechanics. The results show that a shallow quantum circuit is ableto represent European and Asian Call option prices accurately, and therefore suggest that this approach isa promising candidate for the application of quantum computing in Finance.The main challenge of the present methodology is the requirement for an ansatz circuit and the corre-sponding solution of an optimisation problem. More work is needed in the future to design an efficient ansatzfor more complex financial products, or in the development of an ansatz-free approach. One of the mostpromising application of this technique is for basket options, based on several stocks, and hence multidi-mensional. Classical PDE methods suffer from the so-called curse of dimensionality, making such pricingproblem cumbersome, or at least computationally intensive. Since its mathematical formulation is similar
60 80 100 120 140 160 S P a y o ff (a) Classical payoffQuantum approximation
60 80 100 120 140 160 S C a ll p r i c e (b) Classical simulationQuantum approximation S P a y o ff S C a ll p r i c e Figure 6.
Arithmetic Asian Call option prices calculated using the hybrid quantum-classical al-gorithm. The left plot displays the expected solution and the quantum approximationfor the payoff function, while the right one shows the corresponding prices at inception. to the simulation of large-scale quantum mechanical systems governed by systems of Schr¨odinger equations,we believe that a modified version of our algorithm is applicable there, leveraging the power of quantumcomputing, and we leave this investigation for the next step.
QUANTUM ALGORITHM FOR LINEAR PDES ARISING IN FINANCE 13
Appendix A. Complete ansatz circuit .. | ⟩ .... X .... R y ( θ ) . | ⟩ .... H .... R y ( θ ) .... R y ( θ ) ..... R y ( θ ) . | ⟩ .... H .... R y ( θ ) .... R y ( θ ) ..... R y ( θ ) .... R y ( θ ) ..... R y ( θ ) . | ⟩ .... H .... R y ( θ ) .... R y ( θ ) ..... R y ( θ ) .... R y ( θ ) ..... R y ( θ ) .... R y ( θ ) ..... R y ( θ ) ..... R y ( θ ) .... R y ( θ ) ..... R y ( θ ) .... R y ( θ ) ..... R y ( θ ) ..... R y ( θ ) .... R y ( θ ) ..... R y ( θ ) ..... R y ( θ ) . Figure 7.
Ansatz circuit explained in Section 5. The circuit consists of 25 R y gates. Appendix B. Initial and terminal conditions of θ European Call Option Asian Call Option τ = 0 τ = σ T τ = 0 τ = σ Tθ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ Table 3.
Values of θ at initial and terminal times. QUANTUM ALGORITHM FOR LINEAR PDES ARISING IN FINANCE 15
References [1] P. Benioff. The computer as a physical system: a microscopic Quantum Mechanical Hamiltonian model of computers asrepresented by Turing machines.
Journal of Statistical Physics , (563), 1980.[2] S.C. Benjamin and Y. Li. Efficient variational quantum simulator incorporating active error minimization. Physical ReviewX , (2), 2017[3] F. Black and M. Scholes. The pricing of options and corporate liabilities. Journal of Political Economy , (3): 637-654,1973.[4] G. Brassard, P. Hoyer, M. Mosca and A. Tapp. Quantum amplitude amplification and estimation. AMS ContemporaryMathematics , : 53-74, 2002.[5] C. Brown, J.C. Handley, C.T. Lin and K.J. Palmer. Partial differential equations for Asian option prices. QuantitativeFinance , (3), 447-460, 2016.[6] S. Endo, I. Kurata and Y.O. Nakagawa. Calculation of the Green’s function on near-term quantum computers. arXiv:1909.12250, 2019.[7] D. Deutsch. Quantum theory, the Church-Turing principle and the universal quantum computer. Proceedings of the RoyalSociety A , (1818): 97-117, 1985.[8] R.P. Feynman. Simulating physics with computers. International Journal of Theoretical Physics , : 467-488, 1982.[9] A. Martin, B. Candelas, A. Rodriguez-Rozas, J.D. Martin-Guerrero, X. Chen, L. Lamata, R. Orus, E. Solano and M. Sanz.Towards pricing financial derivatives with an IBM quantum computer. arXiv: 1904.05803, 2019.[10] S. McArdle, T. Jones, S. Endo, Y. Li, S.C. Benjamin and X. Yuan. Variational ansatz based quantum simulation ofimaginary time evolution. Quantum Information , (1), 2019.[11] A. McLachlan. A variational solution of the time-dependent Schr¨odinger equation. Molecular Physics , (1), 39-44, 1964.[12] Y.I. Manin. Computable and Noncomputable. Sov.Radio : 13-15, 1980.[13] M. Motta, C. Sun, A. T.K. Tan, M.J. O’Rourke, E. Ye, A.J. Minnich, F. Brandao and G.K Chan. Determining eigenstatesand thermal states on a quantum computer using quantum imaginary time evolution.
Nature Physics , 2019.[14] M.A. Nielsen and I.L. Chuang. A variational solution of the time-dependent Schr¨odinger equation. Cambridge UniversityPress, 2011.[15] R. Orus, S. Mugel and E. Lizaso. Quantum computing for finance: Overview and prospects.
Reviews in Physics , , 2019.[16] P. Rebentrost, B. Gupt and T.R. Bromley. Quantum computational finance: Monte Carlo pricing of financial derivatives. Physical Review A , (2), 2018.[17] P.W. Shor. Algorithms for quantum computation: discrete logarithms and factoring. Proceedings 35th Annual Symposiumon Foundations of Computer Science , 1994.[18] N. Stamatopoulos, D.J. Egger, Y. Sun, C. Zoufal, R. Iten, N. Shen and S. Woerner. Option pricing using quantumcomputers. arXiv: 1905.02666, 2019.[19] J. Vecer. A new PDE approach for pricing arithmetic average Asian options.
Journal of Comp. Finance , : 105-113, 2001. Lloyds Banking Group plc, Commercial Banking, 10 Gresham Street, London, EC2V 7AE, UK
E-mail address : [email protected]
Department of Mathematics, Imperial College London, and Alan Turing Institute
E-mail address : [email protected] Lloyds Banking Group plc, Commercial Banking, 10 Gresham Street, London, EC2V 7AE, UK
E-mail address ::