Toward pricing financial derivatives with an IBM quantum computer

PHYSICAL REVIEW RESEARCH 3, 013167 (2021) Toward pricing financial derivatives with an IBM quantum computer Ana Martin ,1,*Bruno Candelas ,1,*Ángel Rodríguez-Rozas,2José D. Martín-Guerrero ,3Xi Chen,1,4Lucas Lamata ,5 Román Orús ,6,7,8Enrique Solano,1,4,7,9and Mikel Sanz 1,7,9,† 1Department of Physical Chemistry, University of the Basque Country UPV/EHU, Apartado 644, 48080 Bilbao, Spain 2Risk Division, Banco Santander, Avenida de Cantabria S/N, 28660 Boadilla del Monte, Madrid, Spain 3IDAL, Electronic Engineering Department, University of Valencia, Avinguda de la Universitat s/n, 46100 Burjassot, Valencia, Spain 4International Center of Quantum Artificial Intelligence for Science and Technology (QuArtist) and Physics Department, Shanghai University, 200444 Shanghai, China 5Departamento de Física Atómica, Molecular y Nuclear, Universidad de Sevilla, 41080 Sevilla, Spain 6Donostia International Physics Center, Paseo Manuel de Lardizabal 4, 20018 San Sebastián, Spain 7IKERBASQUE, Basque Foundation for Science, Plaza Euskadi 5, 48009 Bilbao, Spain 8Multiverse Computing, Pio Baroja 37, 20008 San Sebastián, Spain 9IQM, Nymphenburgerstr. 86, 80636 Munich, Germany (Received 6 August 2019; revised 1 November 2019; accepted 18 December 2020; published 22 February 2021) Pricing interest-rate financial derivatives is a major problem in finance, in which it is crucial to accurately reproduce the time evolution of interest rates. Several stochastic dynamics have been proposed in the literature to model either the instantaneous interest rate or the instantaneous forward rate. A successful approach to model the latter is the celebrated Heath-Jarrow-Morton framework, in which its dynamics is entirely specified by volatility factors. In its multifactor version, this model considers several noisy components to capture at best the dynamics of several time-maturing forward rates. However, as no general analytical solution is available, there is a trade-off between the number of noisy factors considered and the computational time to perform a numerical simulation. Here, we employ the quantum principal component analysis to reduce the number of noisy factors required to accurately simulate the time evolution of several time-maturing forward rates. The principal components are experimentally estimated with the five-qubit IBMQX2 quantum computer for 2 ×2and3×3 cross-correlation matrices, which are based on historical data for two and three time-maturing forward rates. This paper is a step towards the design of a general quantum algorithm to fully simulate on quantum computers the Heath- Jarrow-Morton model for pricing interest-rate financial derivatives. It shows indeed that practical applications of quantum computers in finance will be achievable in the near future. DOI: 10.1103/PhysRevResearch.3.013167 I. INTRODUCTION In finance, derivatives are contracts whose value derives from the value of an underlying financial asset or a set of assets, like an index, bonds, currency rates, stocks, market indices, or interest rates. Typical financial derivatives contracts include forwards, futures, swaps (currency swaps or interest rate swaps), caps, floors, and swaptions, among many others. They are typically used either to manage (mitigate) risk exposure (hedging), or for pure speculation. In the case of pricing interest-rate financial derivatives under the riskneutral assumption [1,2], it is crucial to model accurately the time evolution of interest rates. Several stochastic dynamics *These authors contributed equally to this work. †Corresponding author: [email protected] Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. have been proposed in the literature to model either the instantaneous interest rate r(t) (also known as the instantaneous spot rate or, simply, as the short rate) or the instantaneous forward rate, which is the forward rate at a future, infinitesimal period (T,T+δt) forecasted at a previous time t, denoted by f(t,T)[1]. Simple dynamics based on one or two noisy (random) factors for modeling both the short rate and the forward rates have been proposed [1,3–5]. For short rates, one- and two-factor models became popular, such as the the Vasicek model, the Hull-White model, the Cox-Ingesroll-Ross (CIR) model and its CIR++ extension as one-factor models, and the Gaussian-Vasicek model and the Hull-White model, as two-factor models. Furthermore, their corresponding algorithms are straightforward to implement. However, these models suffer from the strong requirements which arise from the necessity to calibrate to market data and to capture, at the same time, correlation and covariance structures from the time evolution of different forward rates. A highly successful approach proposed to overcome these constraints is the celebrated Heath-Jarrow-Morton (HJM) framework [6–9], which directly models the time evolution of forward rates. Indeed, the HJM model is a general family of models from which 2643-1564/2021/3(1)/013167(12) 013167-1 Published by the American Physical Society ANA MARTIN et al. PHYSICAL REVIEW RESEARCH 3, 013167 (2021) most of the aforementioned models may be derived [1]. To get a deeper understanding of this model, we refer the reader to Appendix A. In the HJM model, the dynamics is entirely specified by its volatility factors. Although general, there is a trade-off between the number of noisy factors considered and the computational time when executing the algorithm. Therefore, the computational power limits the accuracy of the model. Here is where a quantum computer becomes a useful tool. A quantum computer has a significantly bigger computational capacity than a classical one; therefore, using it would allow us to increase the accuracy of the HJM model. Quantum computing (QC) has emerged in the past years as one of the most exciting applications of quantum technologies [10], which promises to revolutionize the computational power at our disposal. In QC, entanglement, probably the most characteristic signature of quantum physics, is employed as an extra resource to speed up the performance of the computation, since it allows us to parallelize the calculations. Multiple algorithms with provable quantum speedup with respect to their best classical counterparts have been proposed for prime factorization [11], searching in a list [12], solving systems of linear equations [13], and finding the largest eigenvalues and eigenvectors of a given matrix [14], among many others. A particularly relevant application is quantum simulation, in which a controllable quantum system simulates the dynamics of another quantum system of interest whose classical simulation would be highly inefficient. Examples of applications of quantum simulations can be found in spin systems [15,16], quantum chemistry [17–20], quantum field theories [21,22], fluid dynamics [23], and quantum artificial life [24–26]. Some applications of quantum technologies to finances have already been proposed [27–31], but only few experiments have been carried out so far [32–35]. State-of-the-art technology, however, only provides us with small noisy quantum chips, which limits the applicability of digital quantum simulations to toy models. In this article, we employ an efficient quantum principal component analysis (qPCA) algorithm to effectively reduce the number of noisy factors needed to accurately simulate the joint dynamics of several time-maturing forward rates, according to the multifactor HJM model. Although this is a general implementation that can be run in any superconductive-circuit-based quantum processor, we have implemented this algorithm in the five-qubit IBMQX2 superconducting quantum processor of IBM due to its easy accessibility. The volatility factors, ¯σi(τj)=λi(vi)j,(1) where λiand viare the eigenvalues and eigenvectors of the covariance matrix, respectively, are estimated from 2 ×2 and 3 ×3 cross-correlation matrices between different timematuring forward rates based on historical data. To illustrate our qPCA algorithm, we apply this technique to the covariance matrix appearing in Fig. 19.3 in Ref. [36], based on historical data for 1-, 3- and 6-month rates: σ3=⎛ ⎝ 0.000189 0.000097 0.000091 0.000097 0.000106 0.000101 0.000091 0.000101 0.000126⎞ ⎠.(2) This is, to our knowledge, both the first quantum computing experiment in financial option pricing and the largest implementation of the qPCA algorithm on a quantum platform. Although for small matrices the problem can be easily solved on a classical computer, this contribution represents a promising attempt towards the quantum computation of large-scale financial problems which are today prohibitively expensive. In the present noisy intermediate-scale quantum (NISQ) technology era [22], we extend the applications of quantum computers to the field of finance, paving the way for achieving useful quantum supremacy or advantage in the following years. Quantum principal component analysis Principal component analysis (PCA) is a mathematical technique which allows us to find the optimal low-rank approximation of a given matrix by computing its spectral decomposition in eigenvalues and eigenvectors. Indeed, this approximation discards the smallest eigenvalues of the matrix, keeping only the principal components of the spectral decomposition. This technique is of paramount importance for a variety of applications that go from dimensionality reduction problems to problems related to finding patterns in data of big dimension. Unfortunately, the computational cost is too high when the size of the matrix is elevated. It is in this context in which quantum algorithms and quantum computers may play a relevant role. Indeed, in Ref. [14] the authors provided an elegant quantum algorithm to perform PCA with an exponential speedup. The authors assumed that the matrix can be represented by a quantum state; i.e., it is a non-negative matrix with trace equal to 1, which covers a wide range of interesting cases, including the case under study in this paper of covariance matrices associated to volatilities. In this article, we explore the use of this qPCA technique to effectively reduce the dimensionality of the HJM model, by reducing the number of noisy factors without detriment to its accuracy. The qPCA algorithm proposed is specially suited for this financial problem since it is only applicable to density matrices, such as the one given by the correlation matrix used as input for the HJM model. Finally, succeeding in reducing the effective number of noisy factors turns out to be critical for the construction of any quantum Monte Carlo algorithm for this model, after reducing the potential number of quantum gates of the resulting quantum circuit, thus alleviating the limitation in terms of decoherence times of the quantum processor. We implement a slightly modified version of the aforementioned algorithm, which is better adapted to be run in a small and noisy quantum chip, typical in this NISQ technology era [37]. This allows us to reduce the number of noisy factors presented within the HJM model. This is the first step toward the construction of a general quantum computing algorithm to fully simulate the HJM model on the IBM quantum computer for pricing interest-rate financial derivatives. In the following section, we briefly describe the algorithm. II. QUANTUM CIRCUIT Let us consider a non-negative matrix σN∈RN×RNwith tr[σN]=1, which is the matrix whose principal components we want to compute. Let us assume that we can efficiently 013167-2 TOWARD PRICING FINANCIAL DERIVATIVES WITH AN … PHYSICAL REVIEW RESEARCH 3, 013167 (2021) FIG. 1. Quantum circuit implementation for n+log Nqubits. The first nqubits are dedicated to the binary codification of the maximum eigenvalue of the matrix σNand they are initialized in the site |0. The rest of the qubits, a total of log N, encode the estimation of the corresponding eigenvector and are initialized on a random state |b. The single qubit gate Hcorresponds to the Hadamard gate. The rest of the gates are controlled operations. The controlled U2k ρgate applies the matrix U=eitσN2ktimes on the last set of qubits. The controlled R† kgate applies the matrix (10 0e2πi/2k) to each target qubit. After performing all the operations, if the initial |bis appropriate, one gets the final state |b1b2···bn⊗|umax,where|b1b2···bnis the n-bit estimation of the eigenvalue λ(n) jand |umaxis the best estimation of the exact eigenvector of σN. generate the unitary eitσN. It has been proven in the literature that, under certain conditions such as sparsity of the matrix [38,39] or access to several copies of σN[14], this is possible. In our case the covariance matrix σNis not sparse, but we can access several copies of it codified in quantum states. Thus, the best way to efficiently generate the unitary operation eitσNwith accuracy in O(t2−1) steps is the one described by Lloyd et al in Ref. [14]. This matrix admits a spectral decomposition σN=N j=1λj|ujuj|, with 0 ⩽λj⩽1 and N j=1λj=1, and we assume that σNcan be very well approximated by a matrix ρr=r j=1λj|ujuj|with rank rN. Therefore, the goal of the algorithm is the determination of the rlargest eigenvalues of σNand their corresponding eigenvectors. If we want to determine the eigenvalues with an n-bit precision, we will need n+log Nqubits, as depicted in Fig. 1, which represents the gate decomposition of the algorithm. A priori, we do not know the eigenvectors of our algorithm. Hence, we cannot make use of quantum phase estimation to compute directly the corresponding eigenvalue. Consequently, we initialize our system in a random state |bwhose (unknown) decomposition in terms of the eigenbasis is given by |b=N j=1βj|uj. If we take a random vector, the probability that there exists a component βk=0 is zero. The quantum state after the quantum Fourier transform can be written as |b=N j=1βj|(n) j⊗|uj, so eigenvalues and eigenvectors are entangled. (n) j=0.b1b2···bnis the n-bit precision binary representation of the jth eigenvalue of ρr. However, if our assumption that σNis well approximated by the r-rank matrix ρris correct, then the highest eigenvalues should be around 1/r≈n k=1yk2−k. Calling |y(n)=|y1y2···ynthe vector of these components, it means that, by projecting the eigenvalue component |(n) jof the state |baround this component, one may obtain the eigenvector corresponding to the maximum eigenvalue, i.e., y(n)|⊗1|b≈|umax.Itis possible, especially when nis small, as may happen in the NISQ chips, that the n-bit approximation of the eigenvalue cannot be able to distinguish between two or more eigenvectors. In this case, the projection is not into the maximum eigenvalue, but into a K-dimensional subspace containing the indistinguishable components y(n)|⊗1|b=K j=1˜ βj|uj, where the ˜ βjare the normalized βjin the subspace. As we do not know aprioriwhether K>1 or not, we could start with a different random state |c=N j=1γj|uj, which leads to |c=N j=1γj|(n) j⊗|uj. After projecting into |y(n) the expected state is a different superposition L j=1˜γj|uj with high probability, which helps us to check whether we have actually identified the eigenvector corresponding to the maximum eigenvalue. Otherwise, we must increase the n-bit precision until a unique eigenvalue is identified. Let us assume now that the n-bit precision is sufficient to determine a unique eigenvector. Taking into account the constraints due to the small number of qubits and the noise of the chip and the operations, we can sequentially improve the result of the eigenvector. As described above, we start the protocol with a random quantum state |b0, to which the noisy algorithm is applied and the projection into the |y(n)subspace is performed. Let us call the result |b0, which is an approximation for the eigenvector. If we employ now this state as the initial state in the protocol, |b0=|b1, then one expects that the approximation for the eigenvector provided by |b1 improves the fidelity due to the cancellation of coherent errors associated to the βcomponents. Nonetheless, there is a limitation in this sequential improvement related to the decoherence of the qubits and the statistical error of the measurement. In any case, the result can be (slightly) further improved by performing measurements in different bases and averaging, since this cancels some systematic errors of the gates. III. RESULTS As described in the previous section, the protocol is divided into two parts. First, we estimate the eigenvector |umax corresponding to the largest eigenvalue λmax. We start with a random state, apply the circuit implementation shown in Fig. 2, project on the binary n-bit estimation for the largest eigenvalue |y(n), and use this state as the initial state of the process, which sequentially approaches the exact eigenvector. Afterwards, we use this eigenvector to get a more accurate 013167-3 ANA MARTIN et al. PHYSICAL REVIEW RESEARCH 3, 013167 (2021) FIG. 2. (a) Quantum circuit implementation for the 2 ×2ma- trix. The first two qubits encode the 2-bit estimation of the greatest eigenvalue of ρ2and are initialized on the state |0. The last qubit is dedicated to the estimation of the corresponding eigenvector. It is initialized on a random state |b. For the first iteration, we initialize it on the state |+by applying a Hadamard gate. The single-qubit gates represented by the letter Hrefer to the Hadamard gate. The controlled Ugates represent the unitary controlled operations called U2k ρ in Fig. 1. The last two-qubit gate is a controlled S†gate, controlling the first qubit and acting on the second. The final state of the system after running the circuit and taking measures is |11⊗(0.719|0+ 0.659|1). (b) Populations for each iteration. The graphic shows the experimental probabilities of finding the three qubits in each state and its corresponding errors for the four iterations of the algorithm. We have considered both statistical and experimental errors, assuming for the latter an error of 8% for each two-qubit gate. approximation for the eigenvalue λmax by means of quantum phase estimation. For the estimation of the eigenvector, we start with a random state |b0. Hence, the initial state of the system is |0⊗|0⊗|b0. After the first iteration and projecting on the computational basis the eigenvector, we obtain a first estimation, which we call |b1, and use it as the initial state of the system on the next iteration. This is |0⊗|0⊗|b1.We continue this process and iterate ktimes until |bk−1≈|bk. Once we reach that point, we can say that |bk≈|umax. Let us now estimate the eigenvalue λmax. Once the first part is finished and we have an accurate approximation for |umax, we can apply quantum phase estimation [10] to obtain λmax with n-bit precision. The precision is limited in this case by the size of the processor. Our aim is to apply the algorithm to the 3 ×3 matrix given in Eq. (2) in the five-qubit IBMQX2 quantum processor. First, we solve the 2 ×2 submatrix of σ3containing only two maturities, and afterwards, we solve the 4 ×4 expansion of the same matrix. Despite the small size of the problem, the volume of the quantum algorithms allowed in this processor is almost achieved, but we can still obtain relatively accurate results. We have run the algorithm in both the simulator provided by QISKIT [40] and the real IBM quantum processor, reaching accurate results in both cases. A. 2 ×2matrix First, we need to codify the covariance matrix in a quantum state, so we only need to normalize it with respect to its trace, ρ2=σ2 tr(σ2)=0.6407 0.3288 0.3288 0.3593,(3) whose spectral decomposition is given by λ1=0.8576,|u1=0.8347|0+0.5508|1,(4) λ2=0.1424,|u2=0.5508|0−0.8347|1.(5) Let us remark that λmax λ2, a usual characteristic of these correlation matrices, so we can apply the PCA technique to find the optimal low-rank approximation of ρ2. Let us now define the unitary Uρ2=e2πiρ2=0.6260 −0.3068i−0.7170i −0.7170i0.6260 +0.3068i. (6) For the first part of the protocol, we make use of three qubits, two for a 2-bit approximation of the eigenvalue, and a third one one to represent the eigenvector. We apply the first part of the protocol as described above, starting with a quantum state |b0= 1 √2(|0+|1) and projecting into the |y(n)=|11state. After the fourth iteration, each of them averaged over 8192 realizations, the outcome vector estimating the eigenvector stabilizes and we stop. With this final eigenvector, we also rotate the measurement basis in x,y, and an arbitrary direction r=(cos α, −eiβsin α;eiβsin α,eiγcos α) to compute the relative phase and to improve the accuracy of the solution provided. We have chosen the set of angles α= 1.00, β=0.80, and γ=0.16, but any other choice would be valid as long as all of the angles are different from zero. Our estimation for |umaxis consequently given by |umax=[(0.87 ±δ)−i(0.10 ±δ)]|0 +[(0.47 ±δ)+i(0.10 ±δ)]|1,(7) with δ=0.9 the error estimated from the two-qubit gates and measurement fidelity provided by IBM and the statistical error related to the number of repetitions. This δshould be understood as an upper bound of all possible errors, including not only the gate and qubit errors provided by IBMQ, but also other errors, which do not increase with the depth of the algorithm, such as T1and T2, measurement errors, crosstalk, non-Markovian errors, etc. As we had no direct access to the processor, it is not possible to distinguish among them, so we consider the worst possible case. In Appendix Bwe provide more information about how we have estimated the error δ. Let us remark that we have split the complex phase between both states using the global phase. The estimation for the coefficients after each iteration in the zbasis is provided in Table I. We can observe that the algorithm has already converged in the first iteration, and the variations are within the estimated error. We take the eigenvector produced in the last iteration and repeat the algorithm with this one as the initial state measuring in x,y, and rrandom directions to check possible relative phases and to try to remove systematic errors, which yields the states |bx= 013167-4 TOWARD PRICING FINANCIAL DERIVATIVES WITH AN … PHYSICAL REVIEW RESEARCH 3, 013167 (2021) TABLE I. Estimated coefficients of the eigenvector for consecutive iterations of the algorithm in modulus and measured in the z basis. Here, the state of the previous iteration is employed as the initial state in the following iteration until the values are stabilized. Measurements of the eigenvector are performed in the zbasis and repeated for 8192 realizations. Iteration cz 0cz 1 1 0.719 0.695 2 0.707 0.707 3 0.720 0.694 4 0.680 0.734 0.878|0+(0.421 +i0.230)|1,|by=0.878|0+(0.427 + i0.220)|1, and |br=0.985|0+0.175|1. Let us remark that the previous estimation of the eigenvector was performed by projecting into the subspace estimating the eigenvalue into the 2-bit string max =0.11. However, we can now apply quantum phase estimation to improve the estimation for the eigenvalue. We divide the problem into these two stages for two reasons. First, we do not know apriorithe value of the maximum eigenvalue, only the approximate rank, and hence a low n-bit approximation covers a larger range, as explained in the previous section. Additionally, we observe a lower error when the protocol is performed in this manner, probably due to the accumulation of two-qubit gates and the error in the projection for the eigenvalue estimation. However, we cannot be sure, since IBM does not provide the exact quantum circuit which they are performing in the processor. Let us now use three qubits for the eigenvalue estimation max =0.b1b2b3, keeping one qubit to encode the corresponding eigenvector. The depth of the circuit implementation grows and leads us to the decoherence of the system when we run it on the real quantum processor, as depicted in Fig. 3. However, the result provided by the QISKIT simulator, producing the quantum state |111⊗[0.8150|0+0.5794|1], is an almost ideal result for the 3-bit string estimation of the eigenvalue. Indeed, the predicted eigenvalue is =0.111 in binary representation and corresponds to the number λ= 0.875 and the fidelity between |umaxand the one obtained after performing the quantum phase estimation in the QISKIT simulator |uQPEis F=|uQPE|umax|2=0.965.(8) This shows that, with few improvements in the gates and chips or with a lower level programming in the chip, one could substantially improve the results. B. 4 ×4matrix In this case, the matrix σ3will be represented by the twoqubit quantum state ρ4=σ4 tr(σ4)=⎛ ⎜ ⎝ 0.4489 0.2304 0.2162 0 0.2304 0.2518 0.2399 0 0.2162 0.2399 0.2993 0 0000 ⎞ ⎟ ⎠.(9) Thus, the unitary generated, Uρ4=e2πiρ4, is given by Uρ4=⎛ ⎜ ⎝ 0.415 +0.048i−0.108 −0.566i−0.029 −0.702i0 −0.108 −0.566i0.744 −0.030i−0.285 −0.181i0 −0.029 −0.702i−0.285 −0.181i0.618 +0.099i0 0001 ⎞ ⎟ ⎠. The spectral decomposition of ρ4, for the sake of comparability, is given by λ1=0.000,|u1=(0.000,0.000,0.000,1.000), λ2=0.031,|u2=(−0.119,0.786,−0.607,0.000), λ3=0.169,|u3=(0.734,−0.342,−0.587,0.000), λ4=0.800,|u4=(0.669,0.516,0.536,0.000), where the vectors are expressed in the basis {|00,|01,|10,|11}. This problem is much more complicated than the previous one, since we do not implement Uρ4, but the controlled Uρ4. This matrix must be decomposed in terms of two-qubit gates, which dramatically increases the depth of the algorithm and, consequently, the decoherence and the errors. The quantum circuit implementation for this problem is shown in Fig. 4. Following the aforementioned protocol, we start with the state |b0=(|00+|01+|10+|11)/2 and provide the coefficients in the zbasis, for both the simulator and the real processor, in Table II. Afterwards, we measure in different bases in order to compute the relative phases, and take the average to cancel systematic errors. The estimation of the eigenvector is, therefore, |umax=(0.6287 +i0.3991)|00+(0.4010 +i0.0693i)|01 +(0.4807 −i0.1964)|10+(0.0305 +i0.0959)|11. The number of entangling gates performed for this algorithm is at least 18, so the total estimated error δin the coefficients, assuming the 8% error per gate observed in the previous section, is over 100%, which makes in principle the result meaningless. A possible way to address this issue is to make use of some error mitigation methods [41], similar to how it is done in Ref. [32]. By measuring in different bases, we have mitigated the systematic errors associated to measurements. However, the errors that occur during the quantum circuit can be mitigated by employing, for instance, Richardson’s extrapolation 013167-5 ANA MARTIN et al. PHYSICAL REVIEW RESEARCH 3, 013167 (2021) FIG. 3. (a) Quantum phase estimation circuit implementation for the 3-bit estimation of the greatest eigenvalue, (3). After estimating the eigenvector |umax, it is used to improve the estimation of the corresponding eigenvalue λby applying the quantum phase estimation algorithm. In this case we dedicate three qubits to the binary codification of λ. The fourth qubit is initialized on the estimated eigenvector |umax. The fifth two-qubit gate represents a controlled T†gate, controlling the first qubit and acting on the third one. The rest of the gates are the same that have been applied on the previous part of the algorithm. Finally, one takes measures and gets the final state of the system: |(3)⊗|umax=|111⊗[0.8150|0+0.5794|1]. (b) Populations of the 3-bit eigenvalue estimation. This chart shows the probabilities of each state in the QISKIT simulator (green bars) and in the real quantum processor (yellow bars) for the quantum phase estimation algorithm taking the previously obtained eigenvector given in Eq. (C1). The first qubit refers to the subspace of the eigenvector estimation. The next three qubits refers to the subspace of the binary estimation of the eigenvalue. The quantum circuit comprises at least six entangling gates, which leads the system to an almost total decoherence, as reflected in the homogeneous distribution of probabilities in the real chip. [42] to eliminate powers of the noise perturbation. Additionally, to improve the error cancellation during measurements, readout error mitigation [43] could be applied. IV. CONCLUSIONS We have proposed and implemented an efficient quantum algorithm to reduce the number of noisy factors present in the time evolution of forward rates according to the multifactor Heath-Jarrow-Morton model. Indeed, this model considers several noisy components to accurately describe the dynamics of several time-maturing forward rates, which can be gathered in a cross-correlation matrix. The eigenvectors corresponding to the largest eigenvalues of this matrix provide the principal components of the correlations. When the considered data set is large, this calculation turns out to be challenging. The principal components are experimentally estimated using a hybrid classical-quantum algorithm with the five-qubit IBMQX2 quantum computer for 2 ×2 and 3 ×3 cross-correlation matrices, which are based on historical data for two and three time-maturing forward rates. We have obtained a reasonable approximation for both the maximum eigenvalue and its corresponding eigenvector in the 2 ×2 case. For the 4 ×4 matrix, the depth of the algorithm is too high and the experimental errors in the quantum processor prevent us from extracting any useful information. For the 4×4 matrix, the depth of the algorithm is too large and the experimental errors in the quantum processor prevent us from extracting any useful information. The application of some error mitigation techniques could be a solution for this problem. Simultaneously, the simulation in QISKIT shows that it would be achievable in a better experimental set. This means that we have exhausted the computational power provided by the current quantum processor in terms of gate fidelities, connectivity, and number of qubits. Nonetheless, this paper is a step towards the design of a general quantum algorithm to fully simulate on quantum computers the HJM model for pricing interest-rate financial derivatives, and shows that practical applications of quantum computers in finance will be achievable even in the NISQ technology era. The codes we have used to run our experiments both on the real chip and on the QISKIT [40] built-in simulator are available on the following GitHub public repository [44]. These codes have been written according to the 0.7.1 version of QISKIT. ACKNOWLEDGMENTS The authors acknowledge the use of IBM QISKIT for this work. The views expressed are those of the authors and do not reflect the official policy or position of IBM. We also acknowledge funding from projects QMiCS (Grant No. 820505) and OpenSuperQ (Grant No. 820363) of the EU Flagship on Quantum Technologies; the FETOPEN project QuroMorphic, Spanish Government, Grant No. PGC2018-095113-B-I00 (MCIU/AEI/FEDER, UE); Basque Government, Grant No. IT986-16; Spanish Ramón y Cajal Grant No. RYC-2017-22482; Shanghai Municipal Science and Technology Commission (Grants No. 18010500400 and No. 18ZR1415500); and the Shanghai Program for Eastern Scholar. This work is supported by the US Department of Energy, Office of Science, Office of Advanced Scientific Computing Research (ASCR) quantum algorithm teams program, under field work Proposal No. ERKJ333. We also acknowledge support from NSFC (Grant No. 12075145), STCSM (Grants No. 2019SHZDZX01- ZX04, No. 18010500400, and No. 18ZR1415500), and from Spanish Government PID2019-104002GB-C21, PID2019- 104002GBC22 (MCIU/AEI/FEDER, UE). A.M., as the first author, has been responsible for the development of this work with the support of B.C. M.S. suggested the seminal ideas and the adaptation of the algorithm to a NISQ chip. A.R.-R. helped with the financial background and the comparison with the classical results. J.M.-G., X.C., L.L., R.O., and E.S. helped to improve the ideas and results shown in the paper. All authors have carefully proofread the manuscript. M.S. supervised the project throughout all stages. The authors declare that there are no competing interests. 013167-6 TOWARD PRICING FINANCIAL DERIVATIVES WITH AN … PHYSICAL REVIEW RESEARCH 3, 013167 (2021) FIG. 4. Quantum circuit implementation for the 4 ×4 matrix. The first qubit is the only one dedicated for the binary codification of the greatest eigenvalue λof the matrix ρ4. It is initialized on the state |0. The last two qubits encode the estimation of the corresponding eigenvector and are initialized on a random state |b. The final state of the system after measuring is |1⊗|umax. APPENDIX A: THE HEATH-JARROW-MORTON FRAMEWORK In Sec. I, we introduced the HJM model to model the time evolution of forward rates. Here we present a brief description of such a framework. Calculating the fair price of financial derivatives represents a challenge from both the mathematical modeling and the computational power perspectives. Traditionally, simplistic assumptions on the dynamics of financial markets in terms of mathematical pricing models, such as the Black-Scholes model, provided simple analytical tools for both pricing and hedging purposes. However, these simple models have limited features in terms of modeling and, in particular, do not properly account for all the relevant financial risk factors. To solve this issue, sophisticated models such as the HJM framework considered in this work gained increasing interest, requiring one, however, to resort to numerical techniques due to the fact that its solution (price) is no longer available analytically. Solving sophisticated pricing models that account for multiple risk factors using numerical techniques represents a formidable computational task. The most popular technique used by the global financial community is based on Monte Carlo, which, due to its low convergence rate, requires a huge amount of computational resources in order to meet TABLE II. Estimated coefficients of the eigenvector for consecutive iterations of the algorithm in modulus and measured in the z basis. Here, the state of the previous iteration is employed as the initial state in the following iteration until the values are stabilized. Measurements of the eigenvector are performed in the zbasis and repeated for 8192 realizations. Iteration cz 00 cz 01 cz 10 cz 11 Chip 1 0.542 0.503 0.466 0.487 2 0.531 0.498 0.493 0.477 3 0.543 0.493 0.494 0.468 4 0.502 0.492 0.523 0.482 Simulator 1 0.719 0.695 0.695 0.695 2 0.707 0.707 0.695 0.695 3 0.720 0.694 0.695 0.695 4 0.680 0.734 0.695 0.695 acceptable levels of accuracy. Generally speaking, designing an efficient Monte Carlo algorithm requires two subsequent stages: first, to reduce the computational complexity of the model when possible, resorting to dimensionality reduction techniques, and second, to accelerate the Monte Carlo convergence from the resulting reduced model. Quantum computing appears as a promising alternative to speed up both stages. A T-maturity zero-coupon bond (also known as a pure discount bond) is a contract that ensures its investor to accrue one unit of currency at time T(its maturity), whose price at a previous time tis denoted by P(t,T). From this definition, it is clear that at expiry of the contract we must have P(T,T)=1. This time-dependent curve represents a fundamental element in the theory of risk-neutral derivative pricing [1] and will extensively be used throughout this article. P(t,T)isalso known as the curve of discount factors, since it is employed to calculate the present value of future cash flows. The inverse of this amount is called the capitalization factor, providing the capitalization of a present quantity to a future time. In finance, the instantaneous interest rate rt(also known as the instantaneous spot rate, or simply the short rate) is the rate of return of a risk-free investment at time t(for example, a U.S. Treasury bond) (see Ref. [1]). This is also the interest rate applied when borrowing money from the money market and it is given as an annual percentage. We denote by B(t)the time-tvalue B(t) of the money market account, defined as B(t)=exp −t 0 r(s)ds.(A1) Of course, if today is the time t,thevalueofr(t) can be observed in the money market and, therefore, it constitutes a known value. However, for future times T>t,r(T) is uncertain and modeled through a stochastic process. In the risk-neutral framework, when using the money market account B(t) as numéraire, the link between the short rate and the zero coupon is indeed materialized by the risk-neutral pricing formula P(t,T)=EQBB(t) B(T)×1|Ft=EQBe(−T tr(s)ds)|Ft, (A2) where QBis the equivalent martingale measure associated to the numéraire B(t), and Ftdenotes the filtration of the information observed in the market until time t. Therefore, if 013167-7 ANA MARTIN et al. PHYSICAL REVIEW RESEARCH 3, 013167 (2021) TABLE III. Complete outcome of each iteration of the first part of the algorithm. The first two qubits correspond to the state |λ(n), which encodes the 2-bit estimation of the largest eigenvalue λmax. The last qubit corresponds to the estimation of the corresponding eigenvector |umax. Each experiment has been run a total of 8192 times. The table shows both the number of times the experiment ends up on each state (under the heading “Counts”) and the percentage it represents over the total number of runs (“%”). Simulator Chip State Counts % Counts % First iteration 000 2205 27 2474 30 001 428 5 793 10 010 181 2 808 10 011 558 7 1253 15 100 325 4 457 6 101 316 4 304 4 110 2692 33 1087 13 111 1487 18 1016 12 Second iteration 000 2190 27 2538 31 001 422 5 787 10 010 181 2 766 9 011 516 6 1257 15 100 297 4 443 5 101 264 3 277 3 110 2757 34 1063 13 111 1565 19 1061 13 Third iteration 000 2196 27 2401 29 001 418 5 746 9 010 169 2 779 10 011 510 6 1347 16 100 292 4 465 6 101 287 4 297 4 110 2755 34 1117 14 111 1565 19 1040 13 Fourth iteration 000 2179 27 2182 27 001 414 5 671 8 010 191 2 754 9 011 469 6 1100 13 100 317 4 605 7 101 306 4 400 5 110 2801 34 1145 14 111 1515 18 1335 16 today is the time t, then P(t,T) is deterministic and it should match the information observed in the market. However, at any future time tffrom today, t<tf<T,P(tf,T) represents a random variable whose value is model dependent. Models for short rate are typically classified depending upon the number of noisy factors that defines their dynamics. Popular one-factor short-rate models include the Vasicek, the Hull-White model, and the CIR model and its CIR++ extension, among others. They quickly became of lesser interest due to their limitation when pricing financial instruments whose payoffs involve the joint distribution of several of such rates at different maturities, mainly due to their incapability to exhibit the intrinsic decorrelation among them. Motivated from this observation, multifactor models appeared to enrich the correlation structure. As a result, several two-factor models were proposed, such as the Gaussian-Vasicek model and the Hull-White two-factor model. Despite the freedom when modeling the instantaneous short rate in the models mentioned above, some limitations may appear when attempting to calibrate a particular model to the current (observed) market curve of discount factors and to capture, at the same time, the correlation and covariance structure of forward rates. The first sound alternative to shortrate models was introduced by Heath, Jarrow, and Morton in 1992 [8], developing a general framework for modeling the instantaneous forward rates. In its multifactor version, determining the number of noisy factors needed becomes a trade-off between the ability, with increasing noisy factors, to better reproduce correlation and covariance structures while capturing market data, and the computational cost when performing a numerical simulation. The connection between the forward rates f(t,T) and the short rate r(t) is established through the bond price as f(t,T)=−∂ ∂Tlog P(t,T).(A3) When f(t,T)isknownforallT,wemusthave P(t,T)=e−T tf(t,s)ds.(A4) By differentiating Eq. (A2) with respect to Twe obtain −∂P(t,T) ∂T=EQBexp −T t r(s)dsr(T) Ft.(A5) By changing to the T-forward measure QTassociated to the bond price numéraire P(t,T)wehave −∂P(t,T) ∂T =EQTexp −T t r(s)dsr(T)P(t,T) exp −T tr(s)ds Ft =P(t,T)EQT[r(T)|Ft],(A6) and therefore, f(t,T)=EQT[r(T)|Ft].(A7) As such, in the HJM multifactor model, the evolution of a risk-neutral zero-coupon bond price satisfies the following equation: dP(t,T)=P(t,T)r(t)dt + N  i=1T t σi(t,s)dsdWi(t), (A8) where dWi,i=1,...,N, are the uncorrelated Brownian increments associated to the volatilities σi. Using the bond price dynamics (A8) and (A3), we have df(t,T)=α(t,T)dt + N  i=1 σi(t,T)dWi(t),(A9) where α(t,T)= N  i=1 σi(t,T)T t σi(t,s)ds.(A10) 013167-8 TOWARD PRICING FINANCIAL DERIVATIVES WITH AN … PHYSICAL REVIEW RESEARCH 3, 013167 (2021) TABLE IV. After the fourth iteration, we rotate the measurement basis in x,y, and an arbitrary direction r= (cosα, −eiβsin α;eiβsin α, eiγcos α), with α=1.00, β=0.80, and γ=0.16. Here we present the complete outcome for each rotation. The first two qubits correspond to the state |(n), which encodes the 2-bit estimation of the largest eigenvalue λmax. The last qubit corresponds to the estimation of the corresponding eigenvector |umax. Each experiment has been run a total of 8192 times. The table shows both the number of times the experiment ends up on each state (under the heading “Counts”) and the percentage it represents over the total number of runs (“%”). xrotation yrotation rrotation Simulator Chip Simulator Chip Simulator Chip State Counts % Counts % Counts % Counts % Counts % Counts % 000 2215 27 2642 32 2229 27 2503 31 519 6 752 9 001 340 4 599 7 436 5 570 7 2089 26 2341 29 010 633 8 2183 27 631 8 2072 25 236 3 810 10 011 77 1 297 4 145 2 360 4 542 7 1413 17 100 588 7 557 7 549 7 567 7 75 1 160 2 101 4 0 125 2 3 0 133 2 508 6 422 5 110 4249 52 1555 19 4126 50 1736 21 475 6 449 5 111 86 1 234 3 73 1 251 3 3748 46 1845 23 This unique choice for α(t,T) as a function of the volatility terms is what prevents arbitrage. As mentioned before, this is a general framework from which many short-rate models may be derived, upon the particular choices for the σterms. However, not every choice generates a Markovian dynamics. They must also be carefully selected in order to derive practical algorithms that are efficient in terms of computational times. One possibility that ensures Markovianity is to assume that the volatility factors only depend on the time to maturity, so σi=¯σi(T−t)=¯σi(τ). At this point, we can use time series data to calculate the functions ¯σi. For this purpose, we build the covariance matrix between the changes in the forward rates for different time maturities τj(typically for maturities at 1 month, 3 months, 6 monts, 1 year, 2 years, etc.). The result is a symmetric matrix whose diagonal terms are the variances of the rates, while the off-diagonal terms represent the covariances between each pair of rates. Considering all possible time-maturing forwards is computationally costly for the numerical simulations. Using PCA we can obtain the most relevant eigenvectors and their associated eigenvalues. As seen in the literature, most of the evolution of the curve can be explained by considering two or three of such factors. Typically, it is observed that whenever the entries of the first principal component are all similar, then the dominant movement of the curve will be a parallel shift. Also, the second component typically accounts for a twist in the curve. In general, if the eigenvalues are λiand the eigenvectors are vi, the volatility factors will be given by ¯σi(τj)=λi(vi)j.(A11) where λiand viare the eigenvalues and eigenvectors of the covariance matrix, respectively, and are estimated from 2 ×2 and 3 ×3 cross-correlation matrices between different timematuring forward rates based on historical data. APPENDIX B: ERROR δ As we mention in Sec. III, the error δis an upper bound of all possible errors. We cannot differentiate where the errors that we get are coming from since they are a mixture between the two-qubit gates and measurement errors and the inherent errors of the system, which does not increase with the depth of the circuit. Thus, we have to consider the worst possible case, where all of the possible errors are taking place. We first run the experiment for the 2 ×2 covariance matrix and compute the fidelity of our circuit. After running this first experiment successfully and assuming that all of the errors that we get are a consequence of the two-qubit gate error, we estimate the error per two-qubit gate. Thus, the total error δof the process might be described as δ=Fidelity number of two-qubit gates,(B1) and the error per gate, δtwo-qubit gate, is defined as δtwo-qubit gate =Fidelity number of two-qubit gates.(B2) To estimate the upper bound of the error that we expect to suffer when we run the experiment for the 4 ×4 covariance matrix, we simply multiply δtwo-qubit gate by the number of twoqubit gates of the corresponding quantum circuit. APPENDIX C: COMPLETE RESULTS 1. 2 ×2matrix For the first part of the circuit we use two qubits for a 2-bit approximation of the eigenvalue λmax and one to represent the corresponding eigenvector. We start this part by initializing the last qubit, which corresponds to the eigenvector estimation, on the state |+ = 1/√2(|0+|1) and project into the state |y(n)=|11. In Table III we show the complete outcome we get on each iteration. After the fourth iteration the outcome vector estimating the eigenvector stabilizes and we stop. In order to improve the accuracy of the solution provided and to compute the relative phase, we rotate the measurement basis in x,y, and an arbitrary direction r= (cos α, −eiβsin α;eiβsin α,eiγcos α), with α=1.00, β= 013167-9