Option Pricing using Quantum Computers

TL;DR

This paper introduces a quantum computing approach for option pricing that leverages amplitude estimation for quadratic speedup, covering various option types and demonstrating implementation on real quantum hardware.

Contribution

It develops quantum circuits for option pricing using amplitude estimation and evaluates their performance on IBM's quantum hardware with error mitigation.

Findings

Quantum circuits successfully price different options.

Amplitude estimation provides quadratic speedup over classical methods.

Error mitigation improves results on noisy hardware.

Abstract

We present a methodology to price options and portfolios of options on a gate-based quantum computer using amplitude estimation, an algorithm which provides a quadratic speedup compared to classical Monte Carlo methods. The options that we cover include vanilla options, multi-asset options and path-dependent options such as barrier options. We put an emphasis on the implementation of the quantum circuits required to build the input states and operators needed by amplitude estimation to price the different option types. Additionally, we show simulation results to highlight how the circuits that we implement price the different option contracts. Finally, we examine the performance of option pricing circuits on quantum hardware using the IBM Q Tokyo quantum device. We employ a simple, yet effective, error mitigation scheme that allows us to significantly reduce the errors arising from noisy two-qubit gates.

Tables

Table 1. Table 1: Example of the different option types.

	Single-asset	Multi-asset
Path-independent	European put/call	Basket option
Path-dependent	Barrier & Asian options	Multi-asset
		barrier options

Table 2. Table 2: Single-qubit, CNOT, Toffoli gate counts and overall circuit depth required for the full amplitude estimation circuits for each instance in Fig.  8 , as a function of the number of sampling qubits m 𝑚 m . These figures assume all-to-all connectivity across qubits.

#	Single-qubit	CX	CCX	Depth
$m=3$	2,091	2,056	90	3,927
$m=5$	12,768	9,078	378	17,332
$m=7$	52,275	37,132	1,530	70,916
$m=9$	210,144	149,290	6,138	285,204

Equations

$${\hat{\mathbb{E}}}_{\mathbb{P}}[f(\mathbf{X})]=\frac{1}{M}\sum _{i=1}^{M}f({\mathbf{X}}_i)$$

$$\mathcal{A}{|0⟩}_{n+1}=\sqrt{1-a}{|{\psi }_0⟩}_n|0⟩+\sqrt{a}{|{\psi }_1⟩}_n|1⟩$$

$$\tilde{a}={\sin }^2(y\pi /M)\in [0,1].$$

$$|a-\tilde{a}|\le$$

$${\mathcal{Q}}^k\mathcal{A}{|0⟩}_n|0⟩$$

$$\sum _{i=0}^{2^n-1}\sqrt{1-f(S_i)}\sqrt{p_i}|S_i⟩|0⟩+\sum _{i=0}^{2^n-1}\sqrt{f(S_i)}\sqrt{p_i}|S_i⟩|1⟩.$$

$$a=\sum _{i=0}^{2^n-1}f(S_i)p_i=\mathbb{E}[f(S)],$$

$${|\psi ⟩}_n=\sum _{i=0}^{2^n-1}\sqrt{p_i}{|S_i⟩}_n.$$

$${|\psi (\theta )⟩}_n=\sum _{i=0}^{2^n-1}\sqrt{p_i(\theta )}{|i⟩}_n,$$

$${|i⟩}_n|0⟩\to {|i⟩}_n\left(\cos \left[f(i)\right]|0⟩+\sin \left[f(i)\right]|1⟩\right)$$

$$\sum _{i=0}^{2^n-1}\sqrt{p_i}{|i⟩}_n\left[\cos \left(c\tilde{f}(i)+\frac{\pi }{4}\right)|0⟩+\sin \left(c\tilde{f}(i)+\frac{\pi }{4}\right)|1⟩\right],$$

$$\tilde{f}(i)=2\frac{f(i)-f_{\text{min}}}{f_{\text{max}}-f_{\text{min}}}-1,$$

$$P_1=\sum _{i=0}^{2^n-1}{p}_i{\sin }^2\left(c\tilde{f}(i)+\frac{\pi }{4}\right),$$

$$P_1\approx \sum _{i=0}^{2^n-1}{p}_i\left(c\tilde{f}(i)+\frac{1}{2}\right)=c\frac{2\mathbb{E}[f(X)]-f_{\text{min}}}{f_{\text{max}}-f_{\text{min}}}-c+\frac{1}{2}.$$

$${\sin }^2\left(c\tilde{f}(i)+\frac{\pi }{4}\right)=c\tilde{f}(i)+\frac{1}{2}+\mathcal{O}(c^3{\tilde{f}}^3(i))$$

$${|0⟩}_n\stackrel{{\mathcal{P}}_X}{\to }{|\psi ⟩}_n=\sum _{i=0}^{2^n-1}\sqrt{p_i}{|i⟩}_n,$$

$$f_C(S_T)=\max (0,S_T-K).$$

$$f_P(S_T)=\max (0,K-S_T).$$

$${|\psi ⟩}_n|0⟩\to |{\varphi }_1⟩=\sum _{i<K}\sqrt{p_i}{|i⟩}_n|0⟩+\sum _{i\ge K}\sqrt{p_i}{|i⟩}_n|1⟩.$$

$$f(i)=\{\begin{array}{ll}{a}_{<}\cdot i+b_{<}& i<K,\\ a_{\ge }\cdot i+b_{\ge }& i\ge K.\end{array}$$

$$c\tilde{f}(i)+\frac{\pi }{4}=\{\begin{array}{ll}{g}_0& i<K,\\ g_0+g(i)& i\ge K,\end{array}$$

$$\sum _{i<K}\sqrt{p_i}{|i⟩}_n|0⟩\left[\cos (g_0)|0⟩+\sin (g_0)|1⟩\right]+$$

$$\sum _{i\ge K}\sqrt{p_i}{|i⟩}_n|1⟩\left\{\cos [g_0+g(i)]|0⟩+\sin [g_0+g(i)]|1⟩\right\}.$$

$$P_1=\sum _{i<K}{p}_i{\sin }^2(g_0)+\sum _{i\ge K}{p}_i{\sin }^2[g_0+g(i)].$$

$$g(i)=\frac{2c(i-K)}{i_{\text{max}}-K},$$

$$g_0=\frac{\pi }{4}-c.$$

$$P_1$$

$$=\frac{1}{2}-c+\frac{2c}{i_{\text{max}}-K}\sum _{i\ge K}{p}_i(i-K),$$

$$f_s(S)=\max (0,S-K_1)-\max (0,S-K_2),$$

$$f_s(S)=\{\begin{array}{ll}{f}_0(S)& S<K_1,\\ f_0(S)+f_1(S)& K_1\le S<K_2,\\ f_0(S)+f_1(S)+f_2(S)& K_2\le S.\end{array}$$

$${|\psi ⟩}_n=\sum _{i_1,...,i_d}\sqrt{p_{i_1,...,i_d}}{|i_1⟩}_{n_1}\otimes ...\otimes {|i_d⟩}_{n_d},$$

Full text

1 Introduction

Options are financial derivative contracts that give the buyer the right, but not the obligation, to buy (call option) or sell (put option) an underlying asset at an agreed-upon price (strike) and timeframe (exercise window). In their simplest form, the strike price is a fixed value and the timeframe is a single point in time, but exotic variants may be defined on more than one underlying asset, the strike price can be a function of several market parameters and could allow for multiple exercise dates. As well as providing investors with a vehicle to profit by taking a view on the market or exploit arbitrage opportunities, options are core to various hedging strategies and as such, understanding their properties is a fundamental objective of financial engineering. For an overview of option types, features and uses, we refer the reader to Ref. [1].

Due to the stochastic nature of the parameters options are defined on, calculating their fair value can be an arduous task and while analytical models exist for the simplest types of options [2], the simplifying assumptions on the market dynamics required for the models to provide closed-form solutions often limit their applicability [3]. Hence, more often than not, numerical methods have to be employed for option pricing, with Monte Carlo being one of the most popular due to its flexibility and ability to generically handle stochastic parameters [4, 5]. However, despite their attractive features in option pricing, classical Monte Carlo methods generally require extensive computational resources to provide accurate option price estimates, particularly for complex options. Because of the widespread use of options in the finance industry, accelerating their convergence can have a significant impact in the operations of a financial institution.

By leveraging the laws of quantum mechanics a quantum computer [6] may provide novel ways to solve computationally intensive problems such as quantum chemistry [7, 8, 9, 10], solving linear systems of equations [11], and machine learning [12, 13, 14]. Quantitative finance, a field with many computationally hard problems, may benefit from quantum computing. Recently developed applications of gate-based quantum computing for use in finance [15] include portfolio optimization [16], the calculation of risk measures [17] and pricing derivatives [18, 19, 20]. Several of these applications are based on the Amplitude Estimation algorithm [21] which can estimate a parameter with a convergence rate of $1/M$, where $M$ is the number of quantum samples used. This represents a theoretical quadratic speed-up compared to classical Monte Carlo methods.

In this paper we extend the pricing methodology presented in [17, 18] and place a strong emphasis on the implementation of the algorithms in a gate-based quantum computer. We first classify options according to their features and show how to take the different features into account in a quantum computing setting. In Sec. 3, we review the quantum methodology to price options and discuss how to represent relevant probability distributions in a quantum computer. In Sec. 4, we show a framework to price vanilla options and portfolios of vanilla options, options with path-dependent dynamics and options on several underlying assets. In Sec. 5 we show results from evaluating our option circuits on quantum hardware, and describe the error mitigation scheme we employ to increase the accuracy of the estimated option prices. In particular, we employ the maximum likelihood estimation method introduced in [22] to perform amplitude estimation without phase estimation in option pricing using three qubits of a real quantum device.

2 Review of option types and their challenges

Option contracts are valid for a pre-determined period of time, and their value at the expiration date is called the payoff. The goal of option pricing is to estimate the option payoff at the expiration date in the future and then discount that value to determine its worth today. The discounted payoff is also called the fair value and indicates the amount of money one should pay to enter the option contract today, making it worthwile receiving the payoff value at the expiration date.

In practice, we price complex options numerically using Monte Carlo methods by following these steps:

Model the asset price of the option’s underlying(s) and any other sources of uncertainty as random variables $\mathbf{X}=\{X_{1},X_{2},\dots,X_{N}\}$ following a stochastic process. 2. 2.

Generate a large number $M$ of random price paths $\{\mathbf{X}_{1},\mathbf{X}_{2},\dots,\mathbf{X}_{M}\}$ for the underlying(s) from the probability distribution $\mathbb{P}$ implied by the stochastic process. 3. 3.

Calculate the option’s payoff $f(\mathbf{X}_{i})$ on each generated price path and compute an estimator for the expectation value of the payoff $\mathbb{E}_{\mathbb{P}}[f(\mathbf{X})]$ as an average across all paths

[TABLE]

By the Central Limit Theorem, the estimator $\hat{\mathbb{E}}_{\mathbb{P}}$ converges to the expectation value $\mathbb{E}_{\mathbb{P}}$ as the number of paths goes to infinity, with convergence $\mathcal{O}\left(M^{-1/2}\right)$ [23]. 4. 4.

Discount the calculated expectation value to get the option’s fair value.

The discounting process requires knowledge of interest rates at future dates which is itself an important question from a financial modelling perspective. However, for the types of options we consider, this process is not computationally challenging and can be performed classically after the payoff calculation. We therefore do not discount the expected payoff for simplicity.

We classify options according to two categories: path-independent vs path-dependent and options on a single asset or on multiple assets. Path-independent options have a payoff function that depends on an underlying asset at a single point in time. Therefore, the price of the asset up to the exercise date of the option is irrelevant for the option price. By contrast, the payoff of path-dependent options depends on the evolution of the price of the asset and its history up to the exercise date. Table 1 exemplifies this classification. Options that are path-independent and rely on a single asset are the easiest to price, and in most cases numerical calculation is straightforward and would likely not benefit by the use of a quantum computer. Path-independent options on multiple assets are only slightly harder to price since more than one asset is now involved and the probability distributions must account for correlations between the assets, but usually these can be priced quite efficiently on classical computers as well. Path-dependent options on the other hand are significantly harder to price than path-independent options since they require an often expensive payoff calculation at multiple time points on each path, therefore minimizing the number of paths required for this step would lead to a significant benefit in the pricing process. It is this last case where we envision the largest impact of quantum computing.

3 Quantum Methodology

Here we outline the building blocks needed to price options on a gate-based quantum computer. As discussed in the previous section, the critical components are 1) represent the probability distribution $\mathbb{P}$ describing the evolution of random variables $\mathbf{X}=\{X_{1},X_{2},\dots,X_{N}\}$ on the quantum computer, 2) construct the circuit which computes the payoff $f(\mathbf{X})$ and 3) calculate the expectation value of the payoff $\mathbb{E}_{\mathbb{P}}[f(\mathbf{X})]$. In Sec. 3.1 we show how to use Amplitude Estimation to calculate the expectation value of a function of random variables. In Sec. 3.2 we describe the process of loading the relevant probability distributions to a quantum register, and in Sec. 3.3 we construct the circuits to compute the payoff and set up Amplitude Estimation to estimate the expectation value of the payoff. We then have all the ingredients to price options on a quantum computer.

3.1 Amplitude Estimation

The advantage of pricing options on a quantum computer comes from the Amplitude Estimation (AE) algorithm [21] which provides a quadratic speed-up over classical Monte Carlo simulations [24, 25]. Suppose a unitary operator $\mathcal{A}$ acting on a register of $(n+1)$ qubits such that

[TABLE]

for some normalized states $\ket{\psi_{0}}_{n}$ and $\ket{\psi_{1}}_{n}$, where $a\in[0,1]$ is unknown. AE allows the efficient estimation of $a$, i.e., the probability of measuring $\ket{1}$ in the last qubit. This estimation is obtained with an operator $\mathcal{Q}=\mathcal{A}\mathcal{S}_{0}\mathcal{A}^{\dagger}\mathcal{S}_{\psi_{0}}$, where $\mathcal{S}_{0}=1-2\ket{0}\bra{0}$ and $\mathcal{S}_{\psi_{0}}=1-2\ket{\psi_{0}}\ket{0}\bra{\psi_{0}}\bra{0}$, which is a rotation of angle $2\theta_{a}$ in the two-dimensional space spanned by $\ket{\psi_{0}}_{n}\ket{0}$ and $\ket{\psi_{1}}_{n}\ket{1}$. From Eq. (1) we know that $a=\sin^{2}(\theta_{a})$. To obtain an approximation for $\theta_{a}$, AE applies Quantum Phase Estimation [26, 27] to approximate certain eigenvalues of $\mathcal{Q}$, which are classically mapped to an estimator for $a$. The Quantum Phase Estimation uses $m$ additional sampling qubits to represent result and $M=2^{m}$ applications of $\mathcal{Q}$, i.e., $M$ quantum samples. The $m$ qubits, initialized to an equal superposition state by Hadamard gates, are used to control different powers of $\mathcal{Q}$. After applying an inverse Quantum Fourier Transform (QFT), their state is measured resulting in an integer $y\in\{0,...,M-1\}$, which is classically mapped to the estimator for $a$, i.e.

[TABLE]

The full circuit for AE is shown in Fig. 1. The estimator $\tilde{a}$ satisfies

[TABLE]

with probability of at least $8/\pi^{2}$. This represents a quadratic speedup compared to the $\mathcal{O}\left(M^{-1/2}\right)$ convergence rate of classical Monte Carlo methods [28].

To reduce the required number of qubits and the resulting circuit depth, Suzuki et al. have shown that AE can be performed without requiring quantum phase estimation while still maintaining a quadratic speed-up [22]. To this extent, they exploit that

[TABLE]

and by measuring $\mathcal{Q}^{k}\mathcal{A}\ket{0}$ for $k=2^{0},...,2^{m-1}$ for a given $m$ and applying a maximum likelihood estimation, an approximation for $\theta_{a}$ (and hence $a$) can be recovered. If we define $M=2^{m}-1$, i.e. the total number of $\mathcal{Q}$-applications, and we consider $N$ shots for each experiment, it has been shown empirically that the resulting estimation error scales as $\mathcal{O}(1/(M\sqrt{N}))$, i.e., the algorithm achieves the quadratic speed-up in terms of $M$. We will use this approach to demontrate results from real quantum hardware in Sec. 5.

For the option contracts we consider, the random variables involved represent the possible values $S_{i}$ the underlying asset can take, and the corresponding probabilities $p_{i}$ that those values will be realized. For an option with payoff $f$, the $\mathcal{A}$ operator will create the state

[TABLE]

Comparing Eq. (1) and Eq. (5), we can see that

[TABLE]

meaning AE allows us to compute the undiscounted price of an option given a way to represent the option’s payoff as a quantum circuit and create the state of Eq. (5). In the following sections, we describe the necessary components to achieve that.

3.2 Distribution loading

The first component of our option pricing model is the circuit that takes a probability distribution implied for possible asset prices in the future and loads it into a quantum register such that each basis state represents a possible value and its amplitude the corresponding probability. In other words, given an $n$-qubit register, asset prices $\{S_{i}\}$ for $i\in\{0,...,2^{n}-1\}$ and corresponding probabilities $\{p_{i}\}$, the distribution loading module creates the state:

[TABLE]

The analytical formulas used to price options in the Black-Scholes-Merton (BSM) model [2, 29] assume that the underlying stock prices at maturity follow a log-normal distribution with constant volatility. In [30], the authors show that log-concave probability distributions (such as the log-normal distribution of the BSM model) can be efficiently loaded in a gate-based quantum computer. The option types considered in this paper are modeled using the underlying BSM dynamics and thus loading the relevant probability distributions onto quantum registers does not require prohibitive complexity.

However, in option models of practical interest, the simplified assumptions in the BSM model fail to capture important market dynamics, limiting the model’s applicability in real-life scenarios. As such, the market-implied probability distribution of the underlying needs to be captured properly in order for valuation models to accurately estimate the intrinsic value of option contracts. To address these shortcomings, Artificial Neural Networks (ANN) are becoming increasingly more popular as a means to capture the real-life dynamics of the relevant market parameters, without the need to assume a simplified underlying model [31, 32]. It is thus important to be able to efficiently represent distributions of financial parameters on a quantum computer which might not have explicit analytical representations.

The loading of arbitrary states into quantum systems requires exponentially many gates [33], making it inefficient to model arbitrary distributions as quantum gates. Since the distributions of interest are often of a special form, the limitation may be overcome by using quantum Generative Adverserial Networks (qGAN). These networks allow us to load a distribution using a polynomial number of gates [19]. A qGAN can learn the random distribution $X$ underlying the observed data samples $\set{x^{0},\ldots,x^{k-1}}$ and load it directly into a quantum state. This generative model employs the interplay of a classical discriminator, a neural network [34], and a quantum generator (a parametrized quantum circuit). More specifically, the qGAN training consists of alternating optimization steps of the discriminator’s parameters $\phi$ and the generator’s parameters $\theta$. After the training, the output of the generator is a quantum state

[TABLE]

that represents the target distribution. The $n$-qubit state $\ket{i}_{n}=\ket{i_{n-1}...i_{0}}$ encodes the integer $i=2^{n-1}i_{n-1}+...+2i_{1}+i_{0}\in\{0,...,2^{n}-1\}$ with $i_{k}\in\{0,1\}$ and $k=0,...,n-1$. The probabilities $p_{i}(\theta)$ approximate the random distribution underlying the training data. We note that the outcomes of a random variable $X$ can be mapped to the integer set $\{0,...,2^{n}-1\}$ using an affine mapping. Furthermore, the approach can be easily extended to multivariate data, where we use a separate register of qubits for each dimension [19].

Another useful feature in the use of qGANs for loading probability distributions, is the fact that we can tailor the qGAN variational form to construct short-depth circuits for an acceptable degree of accuracy. This in turn allows us to evaluate the performance of option pricing quantum circuits in near-term quantum hardware where resources are still quite limited.

The use of ANNs to represent probability distributions inevitably imposes a cost associated with the training component, in both classical and quantum models. However, during common business practices such as overnight risk assessment of large portfolios which may consist of millions of option contracts, the same probability distributions can be used across a large number of pricing calls defined on the same underlyings. For example, it is quite common during risk assessment valuations to require pricing of several thousands of option contracts defined on the same underlying. As such, the effective training cost per pricing call can be efficiently amortized to represent only a small fraction of the total individual option pricing cost.

It is also noteworthy that the qGAN training performs better if the initial distribution is close to the target distribution [19]. Therefore, as new market data comes in which needs to be incorporated into the probability distribution (e.g. for overnight risk, a lot of the same options as the day before need to be priced, but there is one more day’s worth of market data) the previously trained qGAN can be used as the initial distribution, leading to faster convergence.

3.3 Computing the payoff

We obtain the expectation value of a linear function $f$ of a random variable $X$ with AE by creating the operator $\mathcal{A}$ such that $a=\mathbb{E}[f(X)]$, see Eq. (6). Once $\mathcal{A}$ is implemented we can prepare the state in Eq. (1), and the $\mathcal{Q}$ operator, which only requires $\mathcal{A}$ and generic quantum operations [26, 27]. In this section, we show how to create a close relative of $\mathcal{A}$ and how to combine it with AE to calculate the expectation value of $f$.

The payoff function for the option contracts of interest is piece-wise linear and as such we only need to consider linear functions $f:\{0,...,2^{n}-1\}\to[0,1]$ which we write $f(i)=f_{1}i+f_{0}$. We can efficiently create an operator that performs

[TABLE]

using controlled Y-rotations [17]. To implement the linear term of $f(i)$ each qubit $j$ (where $j\in\{0,\dotsc n-1\}$) in the $\ket{i}_{n}$ register acts as a control for a Y-rotation with angle $2^{j}f_{1}$ of the ancilla qubit. The constant term $f_{0}$ is implemented by a rotation of the ancilla qubit without any controls, see Fig. 2. The controlled Y-rotations can be implemented with CNOT and single-qubit gates [35].

We now describe how to obtain $\mathbb{E}[f(X)]$ for a linear function $f$ of a random variable $X$ which is mapped to integer values $i\in\{0,...,2^{n}-1\}$ that occur with probability $p_{i}$ respectively. To do this we use the procedure outlined immediately above to create the operator that maps $\sum_{i}\sqrt{p_{i}}\ket{i}_{n}\ket{0}$ to

[TABLE]

where $\tilde{f}(i)$ is a scaled version of $f(i)$ given by

[TABLE]

with $f_{\text{min}}=\min_{i}f(i)$ and $f_{\text{max}}=\max_{i}f(i)$, and $c\in[0,1]$ is an additional scaling parameter. The relation in Eq. (11) is chosen so that $\tilde{f}(i)\in[-1,1]$. Consequently, $\sin^{2}[c\tilde{f}(i)+\pi/4]$ is an anti-symmetric function around $\tilde{f}(i)=0$. With these definitions, the probability to find the ancilla qubit in state $\ket{1}$, namely

[TABLE]

is well approximated by

[TABLE]

To obtain this result we made use of the approximation

[TABLE]

which is valid for small values of $c\tilde{f}(i)$. With this first order approximation the convergence rate of AE is $\mathcal{O}(M^{-2/3})$ when $c$ is properly chosen which is already faster than classical Monte Carlo methods [17]. We can recover the $\mathcal{O}(M^{-1})$ convergence rate of AE by using higher orders implemented with quantum arithmetic. The resulting circuits, however, have more gates. This trade-off, discussed in Ref. [17], also gives a formula that specifies which value of $c$ to use to minimize the estimation error made when using AE. From Eq. (12) we can recover $\mathbb{E}[f(X)]$ since AE allows us to efficiently retrieve $P_{1}$ and because we know the values of $f_{\text{min}}$, $f_{\text{max}}$ and $c$.

4 Option pricing on a quantum computer

In this section we show how to price the different options shown in Tab. 1. We put an emphasis on the implementation of the quantum circuits that prepare the states needed by AE. We use the different building blocks reviewed in Sec. 3.

4.1 Path-independent options

The price of path-independent vanilla options (e.g. European call and put options) depends only on the distribution of the underlying asset price $S_{T}$ at the option maturity $T$ and the payoff function $f(S_{T})$ of the option. To encode the distribution of $S_{T}$ in a quantum state, we truncate it to the range $[S_{T,\text{min}},S_{T,\text{max}}]$ and discretize this interval to $\{0,...,2^{n}-1\}$ using $n$ qubits. In the quantum computer, the distribution loading operator $\mathcal{P}_{X}$ creates a state

[TABLE]

with $i\in\{0,...,2^{n}-1\}$ to represent $S_{T}$. This state, exemplified in Fig. 3, may be created using the methods discussed in Sec. 3.2.

We start by showing how to price vanilla call or put options, and then generalize our method to capture the payoff structure of portfolios containing more than one vanilla option.

4.1.1 Vanilla options

Vanilla call options are structured so that if the underlying asset price is larger than a fixed value $K$ (the strike price) at the expiration date, the contract pays the difference between the realized price and the strike. As such, the call option payoff $f_{C}(S_{T})$ can be written as

[TABLE]

Equivalently, the corresponding put option has a similar payoff but it pays if the asset price at expiry is smaller than the strike. That is, the put option payoff is

[TABLE]

The linear part of the payoff can be computed as a quantum circuit with the approach outlined in Sec. 3.3, but we need a way to express the $\max$ operation as a quantum circuit as well. We show how to achieve that for both call and put option types by implementing a comparison between the values encoded in the basis states of Eq. (14) and $K$.

A quantum comparator circuit sets an ancilla qubit $\ket{c}$, initially in state $\ket{0}$, to the state $\ket{1}$ if $i\geq K$ and $\ket{0}$ otherwise. The state $\ket{\psi}_{n}$ in the quantum computer therefore undergoes the transformation

[TABLE]

This operation can be implemented by a quantum comparator [36] based on CNOT and Toffoli gates. Since we know the value of the strike, we can implement a circuit tailored to the specific strike price. We use $n$ ancilla qubits $\ket{a_{1},...,a_{n}}$ and compute the two’s complement of the strike price $K$ in binary using $n$ bits, storing the digits in a (classical) array $t[n]$. For each qubit $\ket{i_{k}}$ in the $\ket{i}_{n}$ register, with $k\in\{0,...,n-1\}$, we compute the possible carry bit of the bitwise addition of $\ket{i_{k}}$ and $t[k]$ into $\ket{a_{k}}$. If $t[k]=0$, there is a carry qubit at position $k$ only if there is a carry at position $k-1$ and $\ket{i_{k}}=1$. If $t[k]=1$, there is a carry qubit at position $k$ if there is a carry at position $k-1$ or $\ket{i_{k}}=1$. After going through all $n$ qubits from least to most significant, $\ket{i}_{n}$ will be greater or equal to the strike price, only if there is a carry at the last (most significant) qubit. This procedure along with the necessary gate operations is illustrated in Fig. 4. An implementation for $K=1.9$ and a three-qubit register is shown in Fig. 6.

To represent the payoff function $f(i)$, we add to $\ket{\phi_{1}}$ a second ancilla qubit, which corresponds to the last qubit in Eq. (10). The payoff function of vanilla options is piece-wise linear

[TABLE]

We now focus on a European call option with payoff $f(i)=\text{max}(0,i-K)$, i.e., $a_{<}=b_{<}=0$, $a_{\geq}=1$, and $b_{\geq}=-K$. To prepare the operator that calculates the payoff in the form of Eq. (10) for use with AE, we set

[TABLE]

where $g(i)$ is a linear function of $i$, and $g_{0}$ is an angle whose value we will carefully select. With this setup, the payoff in Eq. (10) can be constructed, starting from the state $\ket{\phi_{1}}\ket{0}$, by first initializing the last ancilla qubit to the state $\cos(g_{0})\ket{0}+\sin(g_{0})\ket{1}$, and then performing a rotation of the last ancilla qubit controlled by the comparator qubit $\ket{c}$ and the qubits in $\ket{\psi}_{n}$. This rotation operation, implemented by the quantum circuit in Fig. 7, applies a rotation with an angle $g(i)$ only if $i\geq K$. The state of the $n+2$ qubits after this operation becomes

[TABLE]

The probability to find the second ancilla in state $\ket{1}$, efficiently measurable using AE, is

[TABLE]

Now, we must carefully choose the angle $g_{0}$ and the function $g(i)$ to recover the expected payoff $\mathbb{E}[f(X)]$ of the option from $P_{1}$ using the approximation in Eq. (12). To reproduce $f(i)=i-K$ for $i\geq K$ and simultaneously satisfy $c\tilde{f}(i)=g_{0}+g(i)-\pi/4\in[-c,c]$ (see Eq. (11)), we must set

[TABLE]

where $i_{\text{max}}=2^{n}-1$. This choice of $g(i)$ forces us to choose

[TABLE]

To see why, we substitute Eqs. (20) and (21) into Eq. (19) and use the approximation in Eq. (13), which leads to

[TABLE]

where we have used $\sum_{i}p_{i}=1$ in the last equality. Eq. (4.1.1) shows that by setting $g_{0}=\pi/4-c$, we could recover $\mathbb{E}[\max(0,i-K)]$ from $P_{1}$ up to a scaling factor and a constant, from which we can subsequently recover the expected payoff $\mathbb{E}[f(i)]$ of the option given the probability distribution of the underlying asset. We should note that the fair value of the option requires appropriately discounting the expected payoff of the option to today, but as the discounting can be performed after the expectation value has been calculated, we omit it from our discussion for simplicity. We demonstrate the performance of our approach by running amplitude estimation using Qiskit [37] on the overall circuit produced by the elements described in this section, and verifying the convergence to the analytically computed value or classical Monte Carlo estimate. An illustration of the convergence of a European call option with increasing evaluation qubits is shown in Fig. 8.

A straightforward extension of the analysis above yields a pricing model for a European put option, whose payoff $f(i)=\text{max}(0,K-i)$ is equivalent to Eq. (17) with $a_{>}=b_{>}=0$, $a_{\leq}=-1$, and $b_{\leq}=K$.

4.1.2 Portfolios of options

Various popular trading and hedging strategies rely on entering multiple option contracts at the same time instead of individual call or put options and as such, these strategies allow an investor to effectively construct a payoff that is more complex than that of vanilla options. For example, an investor who wants to profit from a volatile asset without picking a direction of where the volatility may drive the asset’s price, may choose to enter a straddle option strategy, by buying both a call and a put option on the asset with the same expiration date and strike. If the underlying asset moves sharply up to expiration date, the investor can make a profit regardless of whether it moves higher or lower in value. Alternatively, the investor may opt for a butterfly option strategy by entering four appropriately structured option contracts with different strikes simultaneously. Because these option strategies give rise to piecewise linear payoff functions, the methodology described in the previous section can be extended to calculate the fair values of these option portfolios.

In order to capture the structure of such option strategies, we can think of the individual options as defining one or more effective strike prices $K_{j}$ and add a linear function $f_{j}(S)=a_{j}S+b_{j}$ between each of these strikes. For example, to price an option strategy with the payoff function

[TABLE]

which corresponds to a call spread (the option holder has purchased a call with strike $K_{1}$ and sold a call with strike $K_{2}$), we use functions $f_{0}$, $f_{1}$, and $f_{2}$ such that

[TABLE]

To match Eq. (23) with Eq. (24), we set $f_{0}(S)=0$, $f_{1}(S)=S-K_{1}$ and $f_{2}(S)=-S+K_{2}$. In general, to price a portfolio of options with $m$ effective-strike prices $K_{1}$, …, $K_{m}$ and $m+1$ functions $f_{0}(S)$, …, $f_{m}(S)$, we need an ancilla qubit per strike to indicate if the underlying has reached the strike. This allows us to generalize the discussion from Sec. 4.1.1. We apply a multi-controlled Y-rotation with angle $g_{j}(i)$ if $i\geq K_{j}$ for each strike $K_{j}$ with $j\in\{1,...,m\}$. The rotation $g_{0}(i)$ is always applied, see the circuit in Fig. 9. The functions $g_{j}(i)$ are determined using the same procedure as in Sec. 4.1.1.

4.2 Multi-asset and path-dependent options

In this section we show how to price options with path-dependent payoffs as well as options on more than one underlying asset. In these cases, the payoff function depends on a multivariate distribution of random variables $\{S_{j}\}$ with $j\in\{1,...,d\}$. The ${S_{j}}$’s may represent one or several assets at discrete moments in time or a basket of assets at the option maturity. In both cases, the probability distribution of the random variables $S_{j}$ are truncated to the interval $[S_{j,\text{min}},S_{j,\text{max}}]$ and discretized using $2^{n_{j}}$ points so that they can be represented by $d$ quantum registers where register $j$ has $n_{j}$ qubits. Thus, the multivariate distribution is represented by the probabilities $p_{i_{1},...,i_{d}}$ that the underlying has taken the values $i_{1}$, …, $i_{d}$ with $i_{j}\in\{0,...,2^{n_{j}}-1\}$. The quantum state that represents this probability distribution, a generalization of Eq. (14), is

[TABLE]

with $n=\sum_{j}n_{j}$. Various types of options, such as Asian options or basket options, require us to compute the sum of the random variables $S_{j}$. The addition of the values in two quantum registers $\ket{a,b}\to\ket{a,a+b}$ may be calculated in quantum computers with adder circuits based on CNOT and Toffoli gates [38, 39, 40]. To this end we add an extra qubit register with $n^{\prime}$ qubits to serve as an accumulator. By recursively applying adder circuits we perform the transformation $\ket{\psi}_{n}\ket{0}_{n^{\prime}}\to\ket{\phi}_{n+n^{\prime}}$ where $\ket{\phi}_{n+n^{\prime}}$ is given by

[TABLE]

Here circuit optimization may allow us to perform this computation in-place to minimize the number of qubit registers needed. Now, we use the methods discussed in the previous section to encode the option payoffs into the quantum circuit.

4.2.1 Basket Options

A European style basket option is an extension of the single asset European option discussed in Sec. 4.1, only now the payoff depends on a weighted sum of $d$ underlying assets. A call option on a basket has the payoff profile

[TABLE]

where $S_{\text{basket}}=\vec{w}\cdot\vec{S}$, for basket weights $\vec{w}=[w_{1},w_{2},\dotsc,w_{d}]$, $w_{i}\in[0,1]$, underlying asset prices at option maturity $\vec{S}=[S_{1},S_{2},\dotsc S_{d}]$ and strike $K$. In the BSM model, the underlying asset prices are described by a multivariate log-normal distribution with probability density function [41]

[TABLE]

where $\ln S=(\ln S_{1},\ln S_{2}\dotsc,\ln S_{d})^{T}$ and $\mu=(\mu_{1},\mu_{2}\dotsc\mu_{d})^{T}$, where each $\mu_{i}$ is the log-normal distribution parameter for each asset defined in the caption of Fig. 3. $\Sigma$ is the $d\times d$ positive-definite covariance matrix of the $d$ underlyings

[TABLE]

with $\sigma_{i}$ the volatility of the $i$th asset, and $-1\leq\rho_{ij}\leq 1$ the correlation between assets $i$ and $j$ and $T$ the time to maturity.

The quantum circuit for pricing a European style basket call option is analogous to the single asset case, with an additional unitary to compute the weighted sum of the uncertainty registers $\ket{i_{1}}_{n_{1}}\dotsc\ket{i_{d}}_{n_{d}}$ before applying the comparator and payoff circuits, controlled by the accumulator register $\ket{b}_{n^{\prime}}=\ket{i_{1}+...+i_{d}}_{n^{\prime}}$. A schematic of these components is shown in Fig. 10. The implementation details of the circuit that performs the weighted sum operator is discussed in Appendix A.

We use a basket option to compare the estimation accuracy between AE and classical Monte Carlo. From Eq. (2), we know that for $M$ applications of the $\mathcal{Q}$ operator (see Fig. 1), the possible values returned by AE will be of the form $\sin^{2}(y\pi/M)$ for $y\in\{0,...,M-1\}$ and the maximum distance between two consecutive values is

[TABLE]

This quantity determines how close $M$ operations of $\mathcal{Q}$ can get us to the amplitude which encodes the option price. Using $\sin^{2}(\pi/4+x)=x+1/2+\mathcal{O}(x^{3})$ for $x\ll 1$, we get $\Delta_{\text{max}}=\pi/M+\mathcal{O}(M^{-3})$ for $\pi/M\ll 1$. From Eq. (3) and Eq. (4.1.1), we can determine that with probability of at least $8/\pi^{2}$, our estimated option price using AE will be within

[TABLE]

of the exact option price, where $c$, $i_{\text{max}}$ and $K$ are the parameters used to encode the option payoff into our quantum circuit, discussed in Sec. 4.1.1. To compare this estimation error to Monte Carlo, we use the same number of samples to price an option classically, and determine the approximation error at the same $8/\pi^{2}\approx 81\%$ confidence interval by repeated simulations. The comparison for a basket option on three underlying assets shows that AE provides a quadratic speed-up, see Fig. 11. It is worth noting that for typical business cases, the number of paths required for acceptable pricing accuracy goes from tens of thousands to millions (depending on the option underlyings) [42, 43], so they are well within the range where amplitude estimation becomes more efficient than Monte Carlo, as shown in Fig. 11.

4.2.2 Asian Options

We now examine arithmetic average Asian options which are single-asset, path-dependent options whose payoff depends on the price of the underlying asset at multiple time points before the option’s expiration date. Specifically, the payoff of an Asian call option is given by

[TABLE]

where $K$ is the strike price, $\bar{S}$ is the arithmetic average of the asset’s value over a pre-defined number of points $d$ between [math] and the option maturity $T$

[TABLE]

The probability distribution of the asset price at time $t$ will again be log-normal with probability density function

[TABLE]

with $\mu_{t}=\left(r-0.5\sigma^{2}\right)\Delta t+\ln(S_{t-1})$ and $\Delta t=T/d$. We can then use the multivariate distribution in Eq. (28), with $\vec{S}$ now a $d$-dimensional vector of asset prices at time points $[t_{1}\dotsc t_{d}]$, instead of distinct underlying prices at maturity $T$. As we are not considering multiple underlying assets that could be correlated, the covariance matrix is diagonal $\Sigma=\Delta t[\text{diag}(\sigma^{2},...,\sigma^{2})]$. An illustration of the probability density function used for an asset defined on two time steps is shown in Fig. 12.

We now prepare the state $\ket{\psi}_{n}$, see Eq. (25), where each register represents the asset price at each time step up to maturity. Using the weighted sum operator of Appendix A with equal weights $1/d$, we then calculate the average value of the asset until maturity $T$, see Eq. (33), into a register $\ket{\bar{S}}$

[TABLE]

Finally, we use the same comparator and rotation circuits that we employed for the basket option illustrated in Fig. 10 to load the payoff of an arithmetic average Asian option into the payoff qubit $\ket{p}$.

4.2.3 Barrier Options

Barrier options are another class of popular option types whose payoff is similar to vanilla European options, but they become activated or extinguished if the underlying asset crosses a pre-determined level called the barrier. In their simplest form, there are two general categories of barrier options

$\bullet$** Knock-Out**

The option expires worthless if the underlying asset crosses a certain price level before the option’s maturity.

$\bullet$** Knock-In**

The option has no value unless the underlying asset crosses a certain price level before maturity.

If the required barrier event for the option to have value at maturity occurs, the payoff then depends only on the value of the underlying asset at maturity and not on the path of the asset until then. If we consider a Knock-In barrier option and label the barrier level $B$, we can write the option’s payoff as

[TABLE]

where $T$ is the time to maturity, $S_{t}$ the asset price at time $t$ with $0<t\leq T$ and $K$ the option strike.

To construct a quantum circuit to price a Knock-In barrier option, we use the same method as for the Asian option where $T$ is divided into $d$ equidistant time intervals with $\Delta t=T/d$, and use registers $\ket{i_{1}}_{n_{1}}\ket{i_{2}}_{n_{2}}\dotsc\ket{i_{d}}_{n_{d}}$ to represent the discretized range of asset prices at time $t\in\{\Delta t,2\Delta t,\dotsc,d\cdot\Delta t=T\}$. The probability distribution of Eq. (34) is used again to create the state $\ket{\psi}_{n}$ in Eq. (25).

To capture the path dependence introduced by the barrier, we use an additional $d$-qubit register $\ket{b}_{d}$ to monitor if the barrier is crossed. Each qubit $\ket{b_{t}}$ in $\ket{b}_{d}$ is set to $\ket{1}$ if $\ket{i_{t}}_{n_{t}}\geq B$. An ancilla qubit $\ket{b_{|}}$ is set to $\ket{1}$ if the barrier has been crossed in at least one time step. This is done by computing the logical OR, see Fig. 5, of every qubit in $\ket{b}_{d}$ and storing the result in the ancilla

[TABLE]

This is computed with $X$ (NOT) and Toffoli gates and $d-2$ ancilla qubits. The ancilla qubit $\ket{b_{|}}$ is then used as a control for the payoff rotation into the payoff qubit, effectively knocking the option in. For Knock-Out barrier options, we can follow the same steps and apply an $X$ gate to the ancilla barrier qubit before using it as control, in this manner knocking the option out if the barrier level has been crossed. A circuit displaying all the components required to price a Knock-In barrier option is shown in Fig. 13. Results from amplitude estimation on a barrier option circuit using a quantum simulator are shown in Fig. 14.

Even though we have focused our attention on barrier options where the barrier event is the underlying asset crossing a barrier from below, we can use the same method to price barrier options where barrier events are defined as the asset crossing the value from above. This only requires changing the comparator circuits to compute $S_{t}\leq B$ in the barrier register $\ket{b}_{d}$.

For all path-dependent options, including the Asian and barrier options we have examined in this section, we note that the choice of time intervals on which we need to represent the probability distribution on quantum registers depends on the type of option in question and the type of underlying(s). For instance, when pricing barrier options, we only need to represent the probability distribution at the barrier dates and at maturity. However, as the choice of which time intervals need to be represented for option pricing is independent of whether we employ a quantum or a classical pricing model, a detailed analysis of this choice is beyond the scope of this work.

5 Quantum hardware results

In this section we examine the performance of European call option circuits evaluated on quantum hardware. Quantum circuits based on standard amplitude estimation are not promising candidates for near-term devices, given that they require extra qubits to control the accuracy of the calculation, and multi-controlled gate operations. Tab. 2 lists the number of single and multi-qubit gates required for the European call option examples in Fig. 8, all of which are far from the reach of current and near-term quantum hardware.

However, a quadratic speed-up is also possible by performing AE without quantum phase estimation (see Sec. 3.1). Even though removing phase estimation from amplitude estimation does not make the overall algorithm immediately compatible with noisy quantum computers, it does lead to significantly shorter circuits than the original implementation of AE, allowing us to examine how it performs on noisy quantum hardware.

We hence focus on the circuits required to perform AE without phase estimation and show results for a European call option, using three qubits, two of which represent the uncertainty and one encodes the payoff. We consider a log-normal random distribution with $S_{0}=2$, $\sigma=40\%$, $r=5\%$, and $T=40/365$, see Fig. 3, and truncate the distribution to the interval defined by three standard deviations around the mean. With two qubits encoding this distribution, the possible values are $[1.21,1.74,2.28,2.81]$, represented by $\ket{00},\ldots,\ket{11}$, with corresponding probabilities $0.1\%,55.4\%,42.5\%$, and $1.9\%$. We set the strike price to $K=1.74$.

To examine the behavior of the circuit for different input probability distributions, we run eight experiments that differ by the initial spot price $S_{0}$ and all other parameters are kept constant. The spot price is varied from $1.8$ to $2.5$ in increments of $0.1$. This way we can use the same circuit for all experiments and only vary the Y-rotation angles used to map the initial probability distribution onto the qubit register. This choice of input parameters allows us to evaluate our circuits for expected option prices in the range $[0.0754,0.7338]$.

Following the procedure detailed in Sec. 3.1, we construct the circuits for $\mathcal{A}\ket{0}_{3}$ and $\mathcal{QA}\ket{0}_{3}$, which correspond to $k=0,1$ (i.e. $m=1$) in Eq. (4), with $n=2$. We then perform repeated measurements of the circuits, and by combining the measured probabilities for all $k$ in a single likelihood function, we can perform a maximum likelihood estimation for $\theta_{a}$, hence obtaining an estimate for $a=\mathbb{E}[f(S)]$ (see Eq. (6)), i.e. the expected payoff. Each experiment is evaluated on the IBM Q Tokyo 20-qubit device with 8192 shots. We repeat each 8192-shot experiment 20 times and average over the 20 measured probabilities in order to limit the effect of any one-off issues with the device. The standard deviation of the measured probabilities across the 20 runs was always $<2\%$. The connectivity of IBM Q Tokyo allows to choose three fully connected qubits for the experiments, and thus, no swaps are required to implement any two-qubit gate in our circuits [37]. For all circuits described in the following sections, we used qubits 6, 10 and 11.

Note that even though we are only interested in the result of a single qubit, we always measure all three qubits to be able to apply readout error mitigation as discussed later in Sec. 5.2.

5.1 Algorithm and Operators

We now show how to construct the operator $\mathcal{A}$ that is required for AE. The log-normal distribution on two qubits can be loaded using a single CNOT gate and four single-qubit rotations [44]. To encode the payoff function, we also exploit the small number of qubits and apply a uniformly controlled Y-rotation instead of the generic construction using comparators introduced in Sec. 4. A uniformly controlled Y-rotation, i.e.

[TABLE]

implements a different rotation angle $\theta_{i}$, $i=0,...,2^{n}-1$ for each state of the $n$-control qubits. For $n=2$, this operation can be efficiently implemented using four CNOT gates and four single qubit Y-rotations [45, 46]. The resulting circuit implementing $\mathcal{A}$ is shown in Fig. 15. Note that in our setup, different initial distributions only lead to different angles of the first four Y-rotations and do not affect the rest of the circuit. Although we use a uniformly controlled rotation, the rotation angles are constructed in the same way described in Sec. 3.3. We use an approximation scaling of $c=0.25$ and the resulting angles are $[\theta_{0},\ldots,\theta_{3}]=[1.1781,1.1781,1.5708,1.9635]$, which shows the piecewise linear structure of the payoff function.

The total resulting circuit requires five CNOT gates and eight single-qubit Y-rotations, see Fig. 15. Since we use uniformly controlled rotations, we do not need any ancilla qubit. Note that if we want to evaluate the circuit for $\mathcal{A}$ alone, we can replace the last CNOT gate in Fig. 15 by classical post-processing of the measurement result: if $q_{1}$ is measured as $\ket{1}$, we flip $q_{2}$ and otherwise we do nothing. This further reduces the overall CNOT gate count to four.

In the remainder of this section, we focus on $\mathcal{Q}\mathcal{A}\ket{0}_{3}$, i.e., the outlined algorithm for $m=1$, to examine the reach of today’s quantum hardware in evaluating AE option pricing circuits which do not require phase estimation. We note that this evaluation is relevant to any quantum algorithm realizing AE without phase estimation and is independent of the approach described in [22] or any other particular implementation.

After optimzing the gate count, the resulting circuit for $\mathcal{QA}\ket{0}_{3}$ consists of $18$ CNOT gates and 33 single-qubit gates. The detailed circuit diagram and applied circuit optimization steps are provided in Appendix B.

5.2 Error mitigation and results

We run the circuits for $\mathcal{A}\ket{0}_{3}$ and $\mathcal{QA}\ket{0}_{3}$ on noisy quantum hardware. The results are affected by readout errors and errors that occur during the execution of the circuits.

To mitigate readout errors we run a calibration sequence in which we individually prepare and measure all eight basis states [47, 37]. The result is a $8\times 8$ readout-matrix $\mathcal{R}$ that holds the probability of measuring each basis state as function of the basis state in which the system was prepared. We use $\mathcal{R}$ to correct all subsequent measurements. The error we measure on $P_{1}$ for $\mathcal{A}\ket{0}_{3}$ was reduced from $\sim 6\%$ to $\sim 4\%$ using readout error mitigation.

Errors occuring during the quantum circuit can be mitigated using Richardson extrapolation [48]. First, the quantum circuit is run using a rescaled Hamiltonian to amplify the effect of the noise. Second, a Richardson extrapolation is used to extract the result of the quantum circuit at the zero noise limit. In hardware, error mitigation is done by stretching the duration of the gates. For each stretch factor the qubit gates need to be recalibrated [8]. Here, we use a simplified error mitigation protocol that circumvents the need to recalibrate the gates but still allows us to increase the accuracy of the quantum hardware [49]. Since the single-qubit and CNOT gates have an average randomized benchmarking fidelity of $99.7\%$ and $97.8$%, respectively, we restrict our error mitigation to the CNOT gates. Furthermore, because the optimized circuit for $\mathcal{A}\ket{0}_{3}$ contains only 4 CNOT gates, we employ the error mitigation protocol only when evaluating $\mathcal{QA}\ket{0}_{3}$ which consists of 18 CNOT gates.

We run the circuit for $\mathcal{QA}\ket{0}_{3}$ three times. In each run we replace the CNOT gates of the original circuit by one, three and five CNOT gates for a total of 18, 54, and 90 CNOT gates, respectively. Since a pair of perfect CNOT gates simplifies to the identity these extra gates allow us to amplify the error of the CNOT gate without having to stretch the gate duration, thus, avoiding the need to recalibrate the gate parameters. As the number of CNOT gates is increased the probability of measuring $\ket{1}$ tends towards 0.5 for all initial spot prices, see Fig. 16(b). After applying a second-order Richardson extrapolation, i.e quadratic extrapolation, we recover the same behavior as the option price obtained from classical simulations, see Fig. 16(c). Our simple error mitigation scheme dramatically increased the accuracy of the calculated option price: it reduced the error, averaged over the initial spot price, from $62\%$ to $21\%$.

6 Conclusion

We have presented a methodology and the quantum circuits to price options and option portfolios on a gate-based quantum computer. We showed how to account for some of the more complex features present in exotic options such as path-dependency with barriers and averages. The results that we show are available in the finance module in Qiskit [37]. Future work may involve calculating the price derivatives [50] with a quantum computer. Pricing options relies on AE. This quantum algorithm allows a quadratic speed-up compared to traditional Monte Carlo simulations, but will most likely require a universal fault-tolerant quantum computer [51]. However, research to improve the algorithms is ongoing [52, 53, 54]. Here we have used a new algorithm [22] that retains the AE speed-up but that uses less gates to measure the price of an option. Furthermore, we employed a simple error mitigation scheme that allowed us to greatly reduce the errors from the noisy quantum hardware. However, larger quantum hardware capable of running deeper quantum circuits with more qubits than the currently available quantum computers is needed to price the typical portfolios seen in the financial industry. Future work could focus on reducing the number of quantum registers in our implementation by performing some of the computation in-place.

7 Acknowledgments

The authors want to thank Abhinav Kandala for the very constructive discussions on error mitigation and real quantum hardware experiments. C.Z. and R.I. acknowledge the support of the National Centre of Competence in Research Quantum Science and Technology (QSIT).

Opinions and estimates constitute our judgment as of the date of this Material, are for informational purposes only and are subject to change without notice. This Material is not the product of J.P. Morgan’s Research Department and therefore, has not been prepared in accordance with legal requirements to promote the independence of research, including but not limited to, the prohibition on the dealing ahead of the dissemination of investment research. This Material is not intended as research, a recommendation, advice, offer or solicitation for the purchase or sale of any financial product or service, or to be used in any way for evaluating the merits of participating in any transaction. It is not a research report and is not intended as such. Past performance is not indicative of future results. Please consult your own advisors regarding legal, tax, accounting or any other aspects including suitability implications for your particular circumstances. J.P. Morgan disclaims any responsibility or liability whatsoever for the quality, accuracy or completeness of the information herein, and for any reliance on, or use of this material in any way. Important disclosures at: www.jpmorgan.com/disclosures

IBM, IBM Q, Qiskit are trademarks of International Business Machines Corporation, registered in many jurisdictions worldwide. Other product or service names may be trademarks or service marks of IBM or other companies.

Appendix A Circuit implementation of weighted sum operator

A.1 Weighted sum of single qubits

In this appendix, we demonstrate an implementation of the weighted sum operator on a quantum circuit. The weighted sum operator $\mathcal{S}$ computes the arithmetic sum of the values of $n$ qubits $\ket{a}_{n}=\ket{a_{1}\dotsc a_{n}}$ weighted by $n$ classically defined non-negative integer weights $\mathbf{\omega}=\left(\omega_{1},\omega_{2},\dotsc,\omega_{n}\right)$, and stores the result into another $m$-qubit register $\ket{s}_{m}=\ket{s_{1}\dotsm s_{m}}$ initialized to $\ket{0}_{m}$. In other words,

[TABLE]

where

[TABLE]

The choice of $m$ ensures that the sum register $\ket{s}_{m}$ is large enough to hold the largest possible weighted sum, i.e. the sum of all weights. Alternatively, we can write the weights in the form of a binary matrix $\Omega=\left(\Omega_{i,j}\right)\in\left\{0,1\right\}^{n\times n^{*}}$, where the $i$-th row in $\Omega$ is the binary representation of weight $\omega_{i}$ and $n^{*}=\max_{i=1}^{d}n_{i}$. We use the convention that less significant digits have smaller indices, so $\ket{s_{1}}$ and $\Omega_{i,1}$ are the least significant digits of the respective binary numbers. Using this binary matrix representation, $\mathcal{S}$ is to add the $i$-th qubit $\ket{a_{i}}$ of the state register to the $j$-th qubit $\ket{s_{j}}$ of the sum register if and only if $\Omega_{i,j}=1$. Depending on the values of the weights, an additional quantum register may be necessary to temporarily store the carries during addition operations. We use $\ket{c_{j}}$ to denote the ancilla qubit used to store the carry from adding a digit to $\ket{s_{j}}$. These ancilla qubits are initialized to $\ket{0}$ and will be reset to their initial states at the end of the computation.

Based on the above setup, we build quantum circuits for the weighted sum operator using three elementary gates: X (NOT), CNOT, and the Toffoli gate (CCNOT). These three gates suffice to build any Boolean function [38]. Starting from the first column in $\Omega$, for each column $j$, we find all elements with $\Omega_{i,j}=1$ and add the corresponding state qubit $\ket{a_{i}}$ to $\ket{s_{j}}$. The addition of two qubits involves three operations detailed in Fig. 17: (a) computation of the carry using a Toffoli gate ($\mathcal{M}$), (b) computation of the current digit using a CNOT ($\mathcal{D}$), (c) reset of the carry computation using two X gates and one Toffoli gate ($\widetilde{\mathcal{M}}$). When adding $\ket{a_{i}}$ to the $j$-th qubit of the sum register, the computation starts by applying $\mathcal{M}$ and then $\mathcal{D}$ to $\ket{a_{i}}$, $\ket{s_{j}}$ and $\ket{c_{j}}$, which adds $\ket{a_{i}}$ to $\ket{s_{j}}$ and stores the carry into $\ket{c_{j}}$. Then, using the same two operations, it adds the carry $\ket{c_{j}}$ to the next sum qubit $\ket{s_{j+1}}$ with carry recorded in $\ket{c_{j+1}}$. The process is iterated until all carries are handled. Finally, it resets the carry qubits by applying $\widetilde{\mathcal{M}}$ in reverse order of the carry computation. We reset the carry qubits in order to reuse them in later computations if necessary.

In general, we need $\max(k-2,0)$ carry qubits to compute the addition of $\ket{a_{i}}$ on $\ket{s_{j}}$, where $k\geq 1$ is the smallest integer satisfying

[TABLE]

where $\rho^{s}_{j,j+k-1}$ is the density matrix corresponding to $\ket{s_{j}\dotsm s_{j+k-1}}$. In the $k=1$ case, i.e. $\ket{s_{j}}=0$, the computation is reduced to “copying” $\ket{a_{i}}$ to $\ket{s_{j}}$ using the bit addition operator $\mathcal{D}$, and no carries would be produced. For $k\geq 2$, Eq. (41) guarantees no carries from $\ket{s_{j+k-1}}$ and beyond. Therefore we can directly compute the carry from $\ket{s_{j+k-2}}$ into $\ket{s_{j+k-1}}$ without worrying about additional carries. This eliminates the need for an ancilla qubit $\ket{c_{j+k-2}}$, and hence the number of carry qubits needed is $k-2$. To further reduce the number of ancilla qubits, we can use any sum qubit $\ket{s_{j}}=\ket{0}$ during the computation. In our case, since we are processing $\Omega$ column by column, all sum qubits more significant than $\ket{s_{j+k-1}}$ would be $\ket{0}$. In other words, we have the last $m-(j+k-1)$ sum qubits usable as carry qubits in the computation described above.

As the weights are known at the time of building the circuit, the possible states that $\ket{s}_{m}$ can have before each addition of a state qubit $\ket{a_{i}}$ are also computable. Since we are adding $\ket{a_{i}}$ to $\ket{s}_{m}$ starting from the least significant bit, $k$ equals the bit length of the maximum possible sum on $\ket{s_{j}\dotsc s_{m}}$ after adding $\ket{a_{i}}$ to $\ket{s_{j}}$. In other words,

[TABLE]

Therefore, the number of carry operations and additional ancilla qubits required for each addition of $\ket{a_{i}}$ can be determined. The term in the $\left\lfloor\cdot\right\rfloor$ in Eq. (42) is upper-bounded by

[TABLE]

where $n_{\mathrm{max}}=\max_{j=1}^{m}\sum_{i=1}{n}\Omega_{i,j}$ is the maximum number of $1$’s in a column of $\Omega$. It immediately follows that the number of non-trivial carry operations (i.e. carry operations that requires $\widetilde{\mathcal{M}}$) required to add $\ket{a_{i}}$ to $\ket{s_{j}\dotsc s_{m}}$ is upper-bounded by

[TABLE]

and the number of ancilla qubits required for the entire implementation of $\mathcal{S}$ is at most the upper bound for $k-2$, since we may use some sum qubits as carries. In other words, the number of ancilla qubits required for $\mathcal{S}$ grows at most logarithmically with the number of state qubits $n$.

A.2 Sum of multi-qubit integers

The weighted sum operator $\mathcal{S}$ can be used to calculate the sum of $d$ multi-qubit positive integers on a quantum register. To do that we first prepare the input register in the state

[TABLE]

where $\ket{a^{(i)}_{1}\dotsc a^{(i)}_{n_{i}}},i\in[1,d]$ is the binary representation of the $i$-th integer to sum with $n_{i}$ qubits, least significant figure first. Then we set the weights as

[TABLE]

or equivalently,

[TABLE]

where $I_{n_{i}\times n^{*}}=\left(I_{n_{i}},0_{n_{i}\times(n^{*}-n_{i})}\right),i\in[1,d]$ and $I_{n_{i}}$ is the $n_{i}$-dimensional identity matrix. Now if we build a weighted sum operator based on the weights in Eq. (46) and apply it on the input state qubits in Eq. (45), we would have the sum of the $d$ integers in $\ket{s}_{m}$.

Fig. 18 shows an example circuit computing the sum of two $3$-digit binary numbers represented on a $6$-qubit quantum register $\ket{a}_{3}\ket{b}_{3}$, and storing the result into a $4$-qubit register $\ket{s}_{4}$. The circuit is implemented by a weighted sum operator $\mathcal{S}$ with weights $\omega=(1,2,4,1,2,4)$. The addition of each qubit onto the sum qubits requires one carry gate ($\mathcal{M}$) followed by one addition gate ($\mathcal{D}$), except for the first bit $\ket{a_{1}}$ which does not have any carries before its addition. This results in a total of $6$ CNOT ($\mathcal{D}$) gates and $5$ Toffoli ($\mathcal{M}$) gates. The $11$ gates are grouped in three groups, as is shown in Fig. 18 by dashed boxes. Each group computes the sum of the bits $\ket{a_{j}}$ and $\ket{b_{j}}$ into $\ket{s_{j}}$ and the carry into $\ket{s_{j+1}}$. Note that separate carry qubits are not required, therefore no carry reset operators $\widetilde{\mathcal{M}}$ are used. In fact, using the above construction for $\mathcal{S}$, no extra carry qubits will be required for the addition of any two binary numbers. In general, $\mathcal{S}$ requires at most $\left\lfloor\log_{2}d\right\rfloor$ ancilla qubits for carrying operations, which directly comes from Eq. (44).

A.3 Weighted sum of multi-qubit integers

In addition to summing up $d$ integers equally, a weight $w_{i}$ may also be added to each integer $a^{(i)}$. In that case, the weight matrix would be

[TABLE]

In the case where $w_{i}$ are not integers, we can rescale the values represented on the quantum register by a common factor to make all weights integers. For example, if we are adding two numbers with weights $0.2$ and $0.8$, we could use integer weights of $w_{1}=1$ and $w_{2}=4$ instead, and reinterpret the resulting sum in postprocessing by dividing it by $5$.

Appendix B Optimized Circuit for $\mathcal{QA}\ket{0}_{3}$

In the following, we describe the circuit used for $\mathcal{QA}\ket{0}_{3}$ requiring only 18 CNOT gates. We have that $\mathcal{Q}=-\mathcal{A}\mathcal{S}_{0}\mathcal{A}^{\dagger}\mathcal{S}_{\psi_{0}}$, where $\mathcal{S}_{0}=1-2\ket{0}\bra{0}$ and $\mathcal{S}_{\psi_{0}}=1-2\ket{\psi_{0}}\ket{0}\bra{\psi_{0}}\bra{0}$ perform reflections on $\ket{0}_{3}$ and $\ket{\psi_{0}}_{2}\ket{0}$, respectively. $S_{\psi_{0}}$ can be implemented up to a global phase using a single-qubit $Z$-gate on the last qubit, which is sufficient to differentiate between $\ket{\psi_{0}}\ket{0}$ and $\ket{\psi_{1}}\ket{1}$. $\mathcal{S}_{0}$ is a bit more difficult and we use circuit synthesis for diagonal unitary matrices to achieve an efficient decomposition into gates [55]. This construction lead to 16 CNOT gates for $\mathcal{Q}$ and 21 for $\mathcal{QA}$, which was still a bit too much to run on real hardware.

To further reduce the CNOT count, we look at the full circuit $\mathcal{QA}\ket{0}_{3}$ and we applied the following optimization steps. The last part in $\mathcal{Q}$ is the application of $\mathcal{A}$. As mentioned in Sec. 5, we can drop the very last CNOT gate and apply it in a classical postprocessing. Furthermore, in $\mathcal{QA}\ket{0}_{3}$, we have $S_{\psi_{0}}$ between $\mathcal{A}$ and $\mathcal{A}^{\dagger}$, i.e. $\mathcal{A}^{\dagger}S_{\psi_{0}}\mathcal{A}$, where the uniformly controlled $Y$-rotations in $\mathcal{A}$ ($\mathcal{A}^{\dagger}$) are right before (after) $S_{\psi_{0}}$. On the other hand, the $Z$-gate that implements $S_{\psi_{0}}$ can be decomposed into an $X$-rotation and a $Y$-rotation. The $Y$-rotation can subsequently be absorbed into one of the uniformly controlled $Y$-rotations in $\mathcal{A}$ or $\mathcal{A}^{\dagger}$, changing the angles accordingly. Since the remaining $X$-rotation commutes with the two neighboring CNOT gates from $\mathcal{A}$ and $\mathcal{A}^{\dagger}$, we can move the $X$-rotation so that the two CNOT gates cancel each other. This reduces the CNOT gate count for $\mathcal{QA}\ket{0}_{3}$ to 18 and the resulting circuit is reported in Fig. 19.