# Derivative Pricing using Quantum Signal Processing

Nikitas Stamatopoulos  Goldman Sachs, New York, NY  William J. Zeng  Goldman
Sachs, New York, NY

###### Abstract

Pricing financial derivatives on quantum computers typically includes quantum
arithmetic components which contribute heavily to the quantum resources
required by the corresponding circuits. In this manuscript, we introduce a
method based on Quantum Signal Processing (QSP) to encode financial derivative
payoffs directly into quantum amplitudes, alleviating the quantum circuits
from the burden of costly quantum arithmetic. Compared to current state-of-
the-art approaches in the literature, we find that for derivative contracts of
practical interest, the application of QSP significantly reduces the required
resources across all metrics considered, most notably the total number of
T-gates by $\sim 16$x and the number of logical qubits by $\sim 4$x.
Additionally, we estimate that the logical clock rate needed for quantum
advantage is also reduced by a factor of $\sim 5$x. Overall, we find that
quantum advantage will require $4.7$k logical qubits, and quantum devices that
can execute $10^{9}$ T-gates at a rate of $45$MHz. While in this work we focus
specifically on the payoff component of the derivative pricing process where
the method we present is most readily applicable, similar techniques can be
employed to further reduce the resources in other applications, such as state
preparation.

##  1 Introduction

Financial derivatives can be priced on quantum computers using Quantum
Amplitude Estimation (QAE) [1] with a quadratic speedup compared to classical
Monte Carlo methods [2, 3, 4, 5, 6, 7]. However, the resources required by the
resulting quantum circuits have been estimated to be considerably larger than
what is expected to be possible with near term devices [8]. Possible reduction
of these resources has been proposed by replacing quantum arithmetic with
lookup tables [9], applying generative methods to prepare relevant probability
distributions [10], extracting multiple other useful quantities in derivative
pricing [11], and improving amplitude estimation [12, 13].

Quantum Signal Processing111This technique is also referred to as Quantum
Singular Value Transformation (QSVT) or Quantum Eigenvalue Transform (QET)
depending on the shape/dimensionality of the subsystem transformed. In this
manuscript we use QSP as an umbrella term encompassing all transformations of
this nature which rely on the same underlying process. (QSP) is a technique to
perform nearly arbitrary polynomial transformations to quantum (sub)systems
with low resource overhead [14]. Various quantum algorithms have been unified
under the QSP framework [15] and in multiple instances the QSP formulation has
led to a significant reduction in resource overheads [16, 17, 18, 19, 20]. In
particular, in Ref. [20] the authors show how to apply QSP for quantum state
preparation in cases where the target amplitudes are known functions of a
variable whose discretized values have been encoded as computational basis
states. Additionally, the authors demonstrated that this QSP approach
significantly outperforms the use of quantum (coherent) arithmetic for the
same task in terms of the total quantum resources required, with most
pronounced impact on the number of qubits.

Because the quantum oracles involved in derivative pricing inevitably include
components heavily relying on quantum arithmetic [8], in this manuscript we
show how QSP and the ideas from Ref. [20] can be used to reduce the quantum
resources needed for derivative pricing for certain very common types of
derivatives. We identify the specific circuit components where QSP can replace
quantum arithmetic routines, and provide implementations for all the unitaries
involved across the entire applicable domain. Moreover, we introduce a new
unitary to drive the QSP process, leading to further reduction in resources
compared to that considered in Ref. [20]. We examine the impact of this method
on concrete derivative applications of practical interest and demonstrate the
precise impact on quantum resources by explicitly constructing the
corresponding circuits.

The rest of the paper in structured as follows: In Sec. 2 we describe the core
components used for derivative pricing on a quantum computer and highlight the
bottlenecks we target in this work. In Sec. 3 we outline the QSP framework and
present the core ideas used to apply the framework to derivative pricing.
Then, in Sec. 4 we describe the implementation details of the proposed QSP
approach and in Sec. 5 we demonstrate how to apply the ideas developed to
derivative pricing instances of practical interest. In Sec. 6 we calculate the
quantum resources needed for quantum advantage in derivative pricing using the
method we develop, and contrast them with those of previous approaches in the
literature. Finally, in Sec. 7, we discuss the implications of this work,
highlight its limitations, and outline potential future research directions on
the topic.

##  2 Derivative Pricing with Quantum Computers

Financial derivatives are contracts whose value depends on the future
performance of one or more underlying assets. Derivative contracts are
typically defined by their _payoff_ , which determines the value of the
contract for a given future _path_ of the underlying assets. Classically,
these contracts are often priced by simulating a large number of possible such
paths according to an appropriate stochastic model, computing the payoff for
each simulated path, and calculating the expectation value across all future
payoffs. The _fair price_ of the derivative contract is then given by that
future expectation value, appropriately discounted to take into account the
difference between the time the contract is priced and its future value. For
an overview of derivative types, features and classical pricing methodologies,
we refer the reader to Ref. [21].

In a quantum setting, the quantum circuits which drive the pricing process
consist of two main components [8]: 1) loading the probability distribution
over paths of the (discretized) stochastic variables involved in the
calculation and 2) encoding the payoff of each path into the amplitude of an
ancilla qubit. The first component can be represented by an operator
$\mathcal{P}$ which, given stochastic paths $\omega\in\Omega$ where each path
occurs with probability $p(\omega)$, it prepares a probability-weighted
superposition of all paths

$\mathcal{P}:\ket{\vec{0}}\rightarrow\sum_{\omega}\sqrt{p(\omega)}\ket{\omega}.$
(1)

The second component is an operator $\mathcal{F}$ which computes the payoff
$f(\omega)$ of the derivative on each path $\ket{\omega}$, and encodes that
value in the amplitude of an ancilla qubit

$\mathcal{F}:\ket{\omega}\ket{0}\rightarrow\ket{\omega}\ket{f(\omega)}\left(\sqrt{f(\omega)}\ket{0}+\sqrt{1-f(\omega)}\ket{1}\right).$
(2)

The operator $\mathcal{A}=\mathcal{F}\mathcal{P}$

$\mathcal{A}:\ket{\vec{0}}\ket{0}\rightarrow\left(\sum_{\omega}\sqrt{p(\omega)}\sqrt{f(\omega)}\ket{\omega}\ket{f(\omega)}\ket{0}+\sum_{\omega}\sqrt{p(\omega)}\sqrt{1-f(\omega)}\ket{\omega}\ket{f(\omega)}\ket{1}\right),$
(3)

creates a state such that the probability of measuring $\ket{0}$ in the last
qubit is the price of the derivative, represented as the expected value of the
payoff, $\mathbb{E}[f]=\sum_{\omega\in\Omega}p(\omega)f(\omega)$.

The implementation of the $\mathcal{F}$ operator of Eq. (2) requires binary
quantum arithmetic to compute a digital representation of the payoff
$\ket{f(\omega)}$ from a path $\ket{\omega}$, before that value is encoded
into the amplitude of the last qubit. It is this component we will aim to
implement using QSP, such that the value $\sqrt{f(\omega)}$ is directly
encoded in the amplitude of an ancilla qubit without requiring the computation
of $\ket{f(\omega)}$, thus eliminating the costly quantum arithmetic involved.

##  3 Quantum Signal Processing

Given a quantum register $\ket{x}_{n}$ which stores an $n$-bit representation
of a scalar $x\in[0,1]$, a real function $g:[0,1]\rightarrow[0,1]$ and a
unitary which creates the state

$U\ket{x}_{n}\ket{0}_{n+1}=\ket{x}_{n}\left(g(x)\ket{\psi_{0}}_{n}\ket{0}+\sqrt{1-g(x)^{2}}\ket{\psi_{1}}_{n}\ket{1}\right),$
(4)

where $\ket{\psi_{0}}_{n}$ and $\ket{\psi_{1}}_{n}$ are normalized quantum
states, Quantum Signal Processing (QSP) can be used to apply polynomial
transformations to the amplitude $g(x)$.

In Ref. [14], the authors observe that given projectors
$\tilde{\Pi}\equiv\mathbb{I}^{\otimes 2n}\otimes\ket{0}\bra{0}$ and
$\Pi\equiv\mathbb{I}^{\otimes n}\otimes\left(\ket{0}\bra{0}\right)^{\otimes
n+1}$, $U$ is a block-encoding of the rank-1 matrix $A=\tilde{\Pi}U\Pi$

$\displaystyle U=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left[\kern
0.0pt\kern-2.5pt\kern-5.55557pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle\mbox{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\Pi$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\\\\\tilde{\Pi}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle A$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\cdot\\\$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\cdot$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\cdot$\hfil\kern 5.0pt\crcr}}}}\right]$}},$ (6)

with a single non-trivial singular value $g(x)$. They then show that
consecutive invocations of $U$ and $U^{\dagger}$, interleaved with projector-
controlled rotation operators $\Pi_{\phi}=e^{i\phi(2\Pi-\mathbb{I})}$ and
$\tilde{\Pi}_{\phi}=e^{i\phi(2\tilde{\Pi}-\mathbb{I})}$, can be used to apply
polynomial transformations to this singular value. Specifically, given phase
factors $\vec{\phi}=(\phi_{1},\phi_{2},\cdots,\phi_{d})$, define

$\displaystyle\mathcal{U}^{\vec{\phi}}=\begin{cases}\tilde{\Pi}_{\phi_{1}}U\displaystyle\prod_{k=1}^{(d-1)/2}\Pi_{\phi_{2k}}U^{\dagger}\tilde{\Pi}_{\phi_{2k+1}}U,&\text{for
odd $d$,}\\\
\displaystyle\prod_{k=1}^{d/2}\Pi_{\phi_{2k-1}}U^{\dagger}\tilde{\Pi}_{\phi_{2k}}U,&\text{for
even $d$}\end{cases}$ (7)

and

$\mathcal{U}_{C}^{\vec{\phi}}=\left(\mathcal{U}^{\vec{\phi}}\otimes\ket{0}\bra{0}+\mathcal{U}^{-\vec{\phi}}\otimes\ket{1}\bra{1}\right).$
(8)

Then, for any polynomial $P$ satisfying $\textrm{deg}(P)\leq d$, $|P(a)|\leq
1,\forall a\in[-1,1]$, and $P$ either even or odd, the authors of Ref. [14]
show that there exist phase factors $\vec{\phi}$ [22, 23, 24] such that the
unitary

$U^{\vec{\phi}}=(\mathbb{I}^{\otimes 2n+1}\otimes
H)\mathcal{U}_{C}^{\vec{\phi}}(\mathbb{I}^{\otimes 2n+1}\otimes H),$ (9)

applies the polynomial $P$ to the singular values of matrix $A$

$U^{\vec{\phi}}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left[\kern
0.0pt\kern-2.5pt\kern-5.55557pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle\mbox{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\big{(}\Pi\otimes\ket{0}\bra{0}\big{)}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\\\\\left(\Pi^{\prime}\otimes\ket{0}\bra{0}\right)$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle P(A)$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\cdot\\\$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\cdot$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\cdot$\hfil\kern 5.0pt\crcr}}}}\right]$}}$ (10)

where $\Pi^{\prime}=\tilde{\Pi}$ for odd $d$ and $\Pi^{\prime}=\Pi$ for even
$d$. Equivalently, $U^{\vec{\phi}}$ creates the state

$U^{\vec{\phi}}\ket{x}_{n}\ket{0}_{n+1}\ket{0}=\ket{x_{n}}\begin{cases}\left(P(g(x))\ket{\psi^{\prime}_{0}}_{n}\ket{0}_{2}+\sqrt{1-P(g(x))^{2}}\ket{\psi^{\prime}_{1}}_{n}\ket{0_{\perp}}_{2}\right),&\text{for
odd $d$,}\\\ &\\\
\left(P(g(x))\ket{0}_{n+2}+\sqrt{1-P(g(x))^{2}}\ket{0_{\perp}}_{n+2}\right),&\text{for
even $d$.}\end{cases}$ (11)

A circuit description of the operator $U^{\vec{\phi}}$ is shown in Fig. 1.

In Ref. [20], the authors consider the scenario where a target function $f(x)$
is computed by first constructing a unitary of the form of Eq. (4) and then
applying a polynomial transformation using QSP approximating the function
composition $\left(f\circ g^{-1}\right)(x)$. Specifically, the unitary in Eq.
(4) is taken to be

$U_{\textrm{sin}}\ket{x}_{n}\ket{0}=\ket{x}_{n}\left(\sin(x)\ket{0}+\cos(x)\ket{1}\right),$
(12)

and QSP transformations based on this unitary are used to prepare states whose
amplitudes are functions of $x$, without requiring quantum arithmetic. With
this unitary the polynomial transformation applied with QSP needs to
approximate $(f\circ\arcsin)(x)$. $U_{\textrm{sin}}$ can be constructed with
$n$ controlled-$R$y rotations, and the unitary $\mathcal{U}^{\vec{\phi}}$ of
Eq. (7) is slightly simpler to implement in this case because the
corresponding projectors satisfy $\Pi=\tilde{\Pi}\equiv\mathbb{I}^{\otimes
n}\otimes\ket{0}\bra{0}$. The choice of $U_{\textrm{sin}}$ additionally
simplifies the final state in Eq. (11), requiring $n$ fewer ancilla qubits,
leading to

$U^{\vec{\phi}}\ket{x}_{n}\ket{0}_{2}=\ket{x_{n}}\left(P(g(x))\ket{0}_{2}+\sqrt{1-P(g(x))^{2}}\ket{0_{\perp}}_{2}\right).$
(13)

In this manuscript we instead consider the unitary

$U_{\textrm{sqrt}}\ket{x}_{n}\ket{0}_{n+1}=\ket{x}_{n}\left(\sqrt{x}\ket{\psi_{0}}_{n}\ket{0}+\sqrt{1-x}\ket{\psi_{1}}_{n}\ket{1}\right),$
(14)

such that QSP needs to approximate the function $f(x^{2})$. The material
difference between these two approaches lies mainly in the amount of resources
required to construct $U_{\textrm{sin}}$ and $U_{\textrm{sqrt}}$, as well as
the polynomial degree needed to construct approximations of sufficient
accuracy to the target transformation in each case. Because the details of QSP
transformations based on $U_{\textrm{sin}}$ are analyzed in Ref. [20], for the
rest of this manuscript we focus on the implementation of the QSP
transformation driven by $U_{\textrm{sqrt}}$, and revisit the former when we
estimate the resources required to apply this technique in derivative pricing.

In this work we focus on target functions of the form
$\left(Ae^{x}-B\right)/C$ because this particular form is ubiquitous in
derivative pricing [21]. This motivation is further illustrated in practice in
Sec. 5 where we apply these methods to price concrete financial derivatives.
Namely, we are looking for a unitary which prepares the state

$\ket{x}_{n}\ket{0}_{n+q}\rightarrow\ket{x}_{n}\left(\sqrt{\frac{Ae^{x}-B}{C}}\ket{\psi^{\prime}_{0}}_{n}\ket{0}_{q}+\sqrt{1-\frac{Ae^{x}-B}{C}}\ket{\psi^{\prime}_{1}}_{n}\ket{0_{\perp}}_{q}\right),$
(15)

for normalized quantum states $\ket{\psi^{\prime}_{0}}_{n}$,
$\ket{\psi^{\prime}_{1}}_{n}$, and $\ket{0_{\perp}}_{q}$ denoting a normalized
state orthogonal to $\ket{0}_{q}$. We assume that constants $A,B,C$ are such
that $(Ae^{x}-B)/{C}\in[0,1]$.

Figure 1: Circuit diagram of the $U^{\vec{\phi}}$ operator of Eq. (11) based
on the unitary $U$ of Eq. (4). The choice of angles
$\vec{\phi}=(\phi_{1},\phi_{2},\cdots,\phi_{d})$ determines the exact
$d$-degree polynomial $P$ applied to the state, such that the probability of
measuring $\ket{00}$ in the bottom two qubits if $d$ is odd and
$\ket{0}_{n+2}$ in the bottom $n+2$ qubits if $d$ is even, will be given by
$|P(g(x))|^{2}$. We use the $\oslash$ symbol to indicate a unitary that is
controlled on every qubit of a quantum register being in the $\ket{0}$ state.

##  4 Implementation

In this section we show how to implement the unitary $U_{\textrm{sqrt}}$ of
Eq. (14) and the choice of polynomial $P$ needed in Eq. (11) in order to
construct the target state Eq. (15) under various scenarios. We assume that
the register $\ket{x}_{n}$ is an $n$-bit, binary, fixed-point representation
of a scalar $x$ with $p$ digits to the left of the binary point

$x=\underbrace{x_{n-1}\cdots x_{n-p}}_{p}.\underbrace{x_{n-p-1}\cdots
x_{0}}_{n-p}.$ (16)

Let us first consider the case where $x$ is positive with $p=0$, such that
$x<1$

$x=\sum_{i=0}^{N-1}\frac{x_{i}\cdot 2^{i}}{N},$ (17)

where $x_{i}\in{0,1}$ denotes the $i$-th bit of $\ket{x}_{n}$ and $N=2^{n}$.
The first step to implement the transformation of Eq. (15) is to encode the
value $\sqrt{x}$ to the amplitude of an ancilla qubit. We start with the state
$\ket{x}_{n}$ and add $n$ ancilla qubits in a uniform superposition

$\ket{x}_{n}\ket{+}_{n}=\sum_{j=0}^{N-1}\frac{1}{\sqrt{N}}\ket{x}_{n}\ket{j}_{n}.$
(18)

We then apply the unitary

$\mathcal{C}:\ket{a}\ket{b}\ket{0}\rightarrow\ket{a}\ket{b}\ket{a<b}$ (19)

to perform a binary comparison between the two qubit registers of Eq. (18).
This produces the state

$\sum_{j=0}^{N-1}\frac{1}{\sqrt{N}}\ket{x}_{n}\ket{j}_{n}\ket{N\cdot
x<j}=\ket{x}_{n}\otimes\left[\sum_{j=0}^{N\cdot
x-1}\frac{1}{\sqrt{N}}\ket{j}_{n}\ket{0}+\sum_{j=N\cdot
x}^{N-1}\frac{1}{\sqrt{N}}\ket{j}_{n}\ket{1}\right],$ (20)

with $x$ given by Eq. (17). By inspection, we observe that the probability of
measuring $\ket{0}$ in the comparator ancilla qubit is then $\frac{N\cdot
x}{N}=x$. The unitary $U_{\textrm{sqrt}}=\mathcal{C}(\mathbb{I}^{\otimes
n}\otimes H^{\otimes n}\otimes\mathbb{I})$ acting on the state
$\ket{x}_{n}\ket{0}_{n+1}$therefore creates the target state of Eq. (14)

$U_{\textrm{sqrt}}\ket{x}_{n}\ket{0}_{n+1}=\ket{x}_{n}\left(\sqrt{x}\ket{\psi_{0}}_{n}\ket{0}+\sqrt{1-x}\ket{\psi_{1}}_{n}\ket{1}\right),$
(21)

where $\ket{\psi_{0}}_{n}$ and $\ket{\psi_{1}}_{n}$ are normalized quantum
states which depend on the value of $x$. A polynomial $P$ which approximates
the function

$f(x)=\sqrt{\frac{Ae^{x^{2}}-B}{C}},$ (22)

in the domain $x\in[0,1]$ allows us to approximately prepare the target state
given by Eq. (15) using Eq. (11).

When $p>0$ such that $x\geq 1$, we can use the same unitary
$U_{\textrm{sqrt}}$ as above, ignoring the location of the binary point $p$,
which in this case generates the state

$U_{\textrm{sqrt}}\ket{x}_{n}\ket{0}_{n+1}=\ket{x}_{n}\left(\sqrt{\frac{x}{2^{p}}}\ket{\psi_{0}}_{n}\ket{0}+\sqrt{1-\frac{x}{2^{p}}}\ket{\psi_{1}}_{n}\ket{1}\right).$
(23)

Because we now have an extra factor of $2^{p}$, the polynomial $P$ in this
case needs to approximate the function

$f(x)=\sqrt{\frac{Ae^{x^{2}\cdot 2^{p}}-B}{C}},$ (24)

in order to generate the target state of Eq. (15).

In order to handle negative values of $x$, we first need to pick a value
$s\geq|x|\geq 0$ and then apply a binary addition circuit which performs
$\ket{x}_{n}\ket{0}_{m}\rightarrow\ket{x}_{n}\ket{x+s}_{m}$, where
$\ket{x+s}_{m}$ is again a fixed-point register with $m$ total digits and $p$
to the left of the binary point. The unitary $U_{\textrm{sqrt}}$ of Eq. (23)
is then applied with the second register as input to construct the state

$U_{\textrm{sqrt}}\ket{x+s}_{m}\ket{0}_{m+1}=\ket{x+s}_{m}\left(\sqrt{\frac{x+s}{2^{p}}}\ket{\psi_{0}}_{m}\ket{0}+\sqrt{1-\frac{x+s}{2^{p}}}\ket{\psi_{1}}_{m}\ket{1}\right),$
(25)

where the normalization factor $2^{p}$ arising from treating the fixed-point
register $\ket{x+s}_{m}$ as a binary string in the comparator circuit,
similarly handles the case $x\leq-1.0$. The polynomial $P$ which generates the
target state of Eq. (15) needs to approximate the function

$f(x)=\sqrt{\frac{Ae^{x^{2}\cdot 2^{p}}e^{-s}-B}{C}}.$ (26)

Note that for both positive and negative values of $x$, the transformation in
Eq. (25) guarantees that the input to the polynomial $P$ will only be the
domain $x\in[0,1]$. This way, while only polynomials with definite parity
(either even or odd) can be implemented using QSP, we can extend the function
we approximate with QSP to the domain $x\in[-1,1]$ using $f(-x)=f(x)$ or
$f(-x)=-f(x)$, even if $f(x)$ does not have definite parity. Additionally, the
unitary $U_{\textrm{sqrt}}$ of Eq. (25) and the function Eq. (26) allow us to
construct the target state for any input $x$, provided that $f(x)\in[0,1]$.

##  5 Applications

###  5.1 European Call Option

In this section we illustrate how the QSP method described in the previous
sections to prepare the state in Eq. (15) can be used in derivative pricing to
compute payoffs of the form $f(S)=\max(0,S-K)$, by first considering a
European call option. Given a strike price $K$, the holder of a European call
option contract receives the amount

$f(S_{T})=\max(0,S_{T}-K),$ (27)

at a future date $T$ (the _expiration date_), where $S_{T}$ is the market
price of an underlying asset at $T$. In order to determine the value of this
contract today, the future price of the underlying asset is modeled as a
stochastic variable with an appropriate probability distribution, and the
price of the option is calculated by computing the expectation value of its
_payoff_ in Eq. (27).

The re-parameterization method for derivative pricing described in Ref. [8],
implements the operator of Eq. (1) by loading a discretized normal
distribution and then performing an affine transformation to the registers
using quantum arithmetic to generate a superposition of possible asset log-
returns $r_{i}$ at time $T$, each weighed by the probability of occurrence

$\sum_{i=0}^{N-1}\sqrt{p_{i}}\ket{r_{i}}_{n},$ (28)

where $N=2^{n}$. The price of the underlying asset for each possible log-
return $r_{i}$ is given by $S_{i}=S_{0}e^{r_{i}}$, where $S_{0}$ is the price
of the underlying asset today. The payoff operator of Eq. (2) is implemented
by first calculating the price of the underlying asset from the log-return
into another register

$\sum_{i=0}^{N-1}\sqrt{p_{i}}\ket{r_{i}}_{n}\rightarrow\sum_{i=0}^{N-1}\sqrt{p_{i}}\ket{r_{i}}_{n}\ket{S_{i}=S_{0}e^{r_{i}}},$
(29)

computing the payoff of Eq. (27)

$\sum_{i=0}^{N-1}\sqrt{p_{i}}\ket{r_{i}}_{n}\ket{S_{i}}\rightarrow\sum_{i=0}^{N-1}\sqrt{p_{i}}\ket{r_{i}}_{n}\ket{S_{i}}\ket{f(S_{i})},$
(30)

and finally encoding the payoff into the amplitude of an ancilla qubit

$\sum_{i=0}^{N-1}\sqrt{p_{i}}\ket{r_{i}}_{n}\ket{S_{i}}\ket{f(S_{i})}\rightarrow\\\
\sum_{i=0}^{N-1}\sqrt{p_{i}}\ket{r_{i}}_{n}\ket{S_{i}}\ket{f(S_{i})}\left(\sqrt{\frac{\max(0,S_{i}-K)}{f_{\textrm{max}}}}\ket{0}+\sqrt{\frac{1-\max(0,S_{i}-K)}{f_{\textrm{max}}}}\ket{1}\right),$
(31)

where $f_{\textrm{max}}=e^{\textrm{max}\\{r_{i}\\}}-K$ is the maximum value of
the payoff possible by the discretization chosen in Eq. (28). The probability
of measuring $\ket{0}$ in the last qubit is

$\mathbb{P}[\ket{0}]=\sum_{i=0}^{N-1}p_{i}\frac{\max(0,S_{i}-K)}{f_{\textrm{max}}},$
(32)

which is the expectation value of the option’s payoff, normalized by
$f_{\textrm{max}}$. Thus, applying amplitude estimation to the last qubit
allows us to estimate the price of the option.

We now describe how QSP and the methods developed in the previous sections can
be used to perform the same computation. Let the $\ket{r_{i}}_{n}$ register of
Eq. (28) represent each log-return in fixed-point representation and let
$s=\ln(S_{0}/K)$. Apply a circuit which performs the comparison
$\ket{r_{i}}_{n}\ket{0}\rightarrow\ket{r_{i}}_{n}\ket{r_{i}\leq s}$ to Eq.
(28) to get

$\sum_{\\{i|r_{i}>s\\}}\sqrt{p_{i}}\ket{r_{i}}_{n}\ket{0}+\sum_{\\{i|r_{i}\leq
s\\}}\sqrt{p_{i}}\ket{r_{i}}_{n}\ket{1}.$ (33)

Add an $m$-qubit register and compute
$\ket{r_{i}}_{n}\ket{0}_{m}\rightarrow\ket{r_{i}}_{n}\ket{r_{i}+|s|}_{m}$,
where $m$ is large enough to hold the maximum value of $r_{i}+|s|$ in fixed-
point representation, to get

$\sum_{\\{i|r_{i}>s\\}}\sqrt{p_{i}}\ket{r_{i}}_{n}\ket{r_{i}+|s|}_{m}\ket{0}+\cdots,$
(34)

where we ignore terms with the last qubit being in the $\ket{1}$ state. Note
that $0\leq r_{i}+|s|<2^{p}$ if the register $\ket{r_{i}+|s|}_{m}$ has $p$
digits to the left of the binary point. The application of the
$U_{\textrm{sqrt}}$ operator of Eq. (25) to the $\ket{r_{i}+|s|}_{m}$ register
and a new qubit register $\ket{0}_{m+1}$ gives

$\sum_{\\{i|r_{i}>s\\}}\sqrt{p_{i}}\ket{r_{i}}_{n}\ket{r_{i}+|s|}_{m}\left(\sqrt{\frac{r_{i}+|s|}{2^{p}}}\ket{\psi_{0}}_{m}\ket{0}+\sqrt{1-\frac{r_{i}+|s|}{2^{p}}}\ket{\psi_{1}}_{m}\ket{1}\right)\ket{0},$
(35)

for normalized quantum states $\ket{\psi_{0}}$ and $\ket{\psi_{1}}$. Employing
the QSP unitary $U_{+}^{\vec{\phi}}$ of Eq. (11) with phases
$\vec{\phi}=(\phi_{1},\phi_{2},\cdots,\phi_{d})$ chosen such that the
polynomial $P$ approximates the function

$f(x)=\sqrt{\frac{S_{0}e^{x^{2}\cdot 2^{p}}e^{-|s|}-K}{f_{\textrm{max}}}}$
(36)

produces

$\sum_{\\{i|r_{i}>s\\}}\sqrt{p_{i}}\ket{r_{i}}_{n}\ket{r_{i}+|s|}_{m}\left(\sqrt{\frac{S_{0}e^{r_{i}}-K}{f_{\textrm{max}}}}\ket{\psi^{\prime}_{0}}_{m}\ket{G}\right),$
(37)

with $\ket{G}=\ket{0}_{3}$ for odd $d$, and

$\sum_{\\{i|r_{i}>s\\}}\sqrt{p_{i}}\ket{r_{i}}_{n}\ket{r_{i}+|s|}_{m}\left(\sqrt{\frac{S_{0}e^{r_{i}}-K}{f_{\textrm{max}}}}\ket{G}\right),$
(38)

with $\ket{G}=\ket{0}_{m+3}$ for even $d$, where we have ignored terms with
any qubit in the last register being in the $\ket{1}$ state. Applying
amplitude estimation to estimate the probability of measuring the $\ket{G}$
state gives us

$\mathbb{P}[\ket{G}]=\sum_{\\{i|r_{i}>s\\}}p_{i}\frac{S_{0}e^{r_{i}}-K}{f_{\textrm{max}}},$
(39)

which is the same as the value we estimate with the standard re-
parameterization method in Eq. (32). In this case however, we do not
explicitly compute the values $\ket{S_{i}=S_{0}e^{r_{i}}}$ or $\ket{f(S_{i})}$
using quantum arithmetic, but rather implicitly calculate them by applying
amplitude transformations using QSP. In Fig. 2, we see plot the function of
Eq. (36) with $S_{0}=K=100$, $p=3$, $s=0$ and
$f_{\textrm{max}}=S_{0}e^{2^{p}}-K$ along with polynomial approximations of
varying degrees constructed using the optimization-based method described in
Appendix A.

We note that Ref. [20] discussed the possibility of applying QSP to compute
common derivative payoffs such as that of Eq. (27). The authors pointed out
that polynomial approximations to functions of the form $f(x)=\sqrt{x}$ suffer
due to the fact that the function is not differentiable at the origin. By
including the exponential in the function we are approximating, not only do we
pass more computation to the QSP transformation, but we also address this
issue since Eq. (36) is differentiable in the entire domain of applicability.

Figure 2: The function $f(x)$ of Eq. (36) with $S_{0}=K=100$, $p=3$, $s=0$ and
$f_{\textrm{max}}=S_{0}e^{2^{p}}-K$ used to compute the payoff of a European
call option, and approximating polynomials $P(x)$ of degree $d\in[6,8,16]$,
which can be applied to a quantum system using the QSP unitary of Eq. (11).
The polynomial approximations are generated using the method described in
Appendix A.

###  5.2 Autocallable

For a single underlying asset, an _autocallable_ derivative contract consists
of

  1. 1.

A set of $k$ binary payoffs $\\{K_{i},t_{i},f_{i}\\}$, such that if the asset
returns at time $t_{i}$ exceed $K_{i}$, the holder of the contract receives a
return $f_{i}$ and the contract terminates. The asset’s return at time $t$ is
defined as $S_{t}/S_{0}$, where $S_{t}$ is the value of the asset at time $t$
and $S_{0}$ the value of the asset today.

  2. 2.

If none of the binary payoffs have been triggered, and the asset returns have
crossed a barrier value $B$ at any time before the contract’s expiration, the
holder receives a return given by $1-\max(0,K_{T}-f_{T})$, where
$K_{T}\in[0,1]$ is a pre-determined _strike_ return and $f_{T}$ is the asset’s
return at expiration date $T$.

Because an autocallable is typically defined in terms of an asset’s
performance, it also includes a (dollar) notional value $V$ on which the
payoff returns are evaluated to determine a dollar amount to be exchanged. For
example, if the contract’s notional value is $\$1$M, and a binary payoff with
$f_{i}=1.5$ is triggered, the holder of the contract receives $\$1.5$M in
total. We now describe how this contract can be priced by employing QSP at the
appropriate place.

Starting again from the re-parameterization method from [8], we first
initialize a superposition of probability-weighted possible paths for the
asset’s log-returns at each timestep

$\sum_{r}\sqrt{p_{r}}\ket{r}\quad\textrm{ with
}\quad\ket{r}=\bigotimes_{t=1}^{T}\ket{r_{t}},$ (40)

where the asset’s total return at timestep $t$ is given by $e^{r_{t}}$. We add
a $k$-qubit register $\ket{s}$ and on each qubit $s_{i}$ we calculate whether
the $k$-th binary payoff is triggered, by checking whether the condition
$r_{t_{i}}>\ln K_{t_{i}}$ is satisfied and no previous binary payoffs have
been triggered

$\ket{s}=\bigotimes_{i=1}^{k}\Ket{s_{i}=\left((r_{t_{i}}>\ln
K_{t_{i}})\cdot(\neg s_{i-1})\right)},$ (41)

where the $\neg$ symbol denotes logical negation and $s_{0}=0$. Each binary
variable $s_{i}$ will thus be $1$ if the asset return at time $t_{i}$ exceeds
the threshold amount $K_{i}$ and none of the previous payoffs have been
triggered, hence indicating whether the payoff at $t_{i}$ should be payed. To
determine whether the payoff from the second autocallable condition is
triggered, we must check if the barrier value $B$ has been crossed at any
timestep of the simulation $t$. This is achieved by first adding a
$(T+1)$-qubit register $\ket{b}$ and calculating $b_{t}=r_{t}\leq\ln B$ in
each qubit for $t\in[1,T]$, and then computing

$b_{T+1}=\left(\neg\prod_{t=1}^{T}b_{i}\right)\cdot\prod_{i=1}^{k}\left(\neg
s_{i}\right),$ (42)

such that $b_{T+1}=1$ if $r_{t}>B$ for at least one value of $t$, and every
qubit in $\ket{s}$ is zero. Therefore, $b_{T+1}=1$ indicates that the barrier
$B$ has been crossed at least once before the contract’s expiration, and no
binary payoffs have been triggered, meaning the second condition of the
autocallable payoff is satisfied. Note that all these logical operations above
can be computed using just CNOT and Toffoli gates.

Un-computing and dropping the first $T$ qubits of the $\ket{b}$ register, we
arrive at the state

$\sum_{r}\sqrt{p_{r}}\ket{r}\bigotimes_{i=1}^{k}\ket{s_{i}}\ket{b}.$ (43)

The $\mathcal{F}$ operator needed to encode the autocallable contract’s payoff
into the amplitude of an ancilla qubit as outlined in Eq. (2), should
construct the state (up to normalization)

$\sum_{r}\sqrt{p_{r}}\ket{r}\bigotimes_{i=1}^{k}\ket{s_{i}}\ket{b}\begin{cases}\left(\sqrt{f_{i}}\ket{0}+\sqrt{1-f_{i}}\ket{1}\right),&\text{if
$s_{i}=1$}\\\
\left(\sqrt{1-\max(0,K_{T}-f_{T})}\ket{0}+\sqrt{\max(0,K_{T}-f_{T})}\ket{1}\right),&\text{if
$b=1$}\\\ \ket{1},&\text{otherwise}\end{cases}$ (44)

where $f_{T}=e^{r_{T}}$ is the asset’s return at the final timestep of the
computation. The first payoff clause can be implemented using controlled Ry
rotations since the values $f_{i}$ are known in advance of the computation.
For the second clause, we first apply a comparator circuit with an ancilla
qubit to calculate $\ket{r_{T}<\ln K_{T}}$, and then controlled on this qubit,
we need to prepare the state

$\ket{r_{T}}_{n}\ket{0}_{n+q}\rightarrow\ket{r_{T}}_{n}\left(\sqrt{1-\left(K_{T}-e^{r_{T}}\right)}\ket{\psi_{0}}_{n}\ket{0}_{q}+\sqrt{\left(K_{T}-e^{r_{T}}\right)}\ket{\psi_{1}}_{n}\ket{0_{\perp}}_{q}\right),$
(45)

where $\ket{\psi_{0}}_{n}$, $\ket{\psi_{1}}_{n}$ are arbitrary normalized
states. Note that we have already imposed the condition
$1-\left(K_{T}-e^{r_{T}}\right)\in[0,1]$, so the state is appropriately
normalized. Since this transformation has the same form as Eq. (15), we can
approximately prepare it using the QSP method described in Sections 3 and 4,
similarly to the European call option described in Sec. 5.1. Assuming
$\ket{r_{T}}_{n}$ is an $n$-qubit fixed-point representation of $r_{T}$ with
$p$ digits to the left of the binary point and spans over both positive and
negative values, we employ the unitary of Eq. (25) with an appropriately
chosen value $s\geq\max(|r_{T}|)\geq 0$, to first generate the state

$U_{\textrm{sqrt}}\ket{r_{T}+s}_{m}\ket{0}_{m+1}=\ket{r_{T}+s}_{m}\left(\sqrt{\frac{r_{T}+s}{2^{p}}}\ket{\psi_{0}}_{m}\ket{0}+\sqrt{1-\frac{r_{T}+s}{2^{p}}}\ket{\psi_{1}}_{m}\ket{0_{\perp}}\right).$
(46)

Subsequently, the polynomial transformation $P$ we apply with QSP in order to
generate Eq. (45) approximates the function

$f(x)=\sqrt{1-\left(K_{T}-e^{x^{2}\cdot 2^{p}}e^{-s}\right)}.$ (47)

While we have described an autocallable contract defined only on a single
underlying asset, the same approach can be extended to multi-asset cases such
as _BestOf_ and _WorstOf_ variants. In those cases, we generate stochastic
paths for all underlying assets and apply a circuit to compare the relative
performance of all assets, identifying the best/worst performing asset. The
payoff is then implemented similarly to the single-asset case following Eq.
(44), based on the returns of the identified asset.

##  6 Resource Estimates for Quantum Advantage

We evaluate the efficiency of employing QSP in the evaluation of derivative
payoffs by estimating the resources needed to price an autocallable derivative
contract using the method described in Sec. 5.2. As a benchmark we use the
resource estimates from Ref. [9], based on the reparameterization method
introduced in Ref. [8] where coherent quantum arithmetic is used to compute
the derivative payoff.

The goal of the reparameterization method is to prepare a probability
distribution across $a$ stochastic underlying assets and $T$ timesteps under
the assumption that asset prices follow geometric Brownian motion. In this
setting, asset prices are described by a multivariate lognormal distribution,
or equivalently, the asset log-returns are distributed normally. The
reparameterization method prepares the distribution over log-returns by first
loading independent standard normal distributions $\mathcal{N}(0,1)$ in
$a\cdot T$ quantum registers, and then applying quantum arithmetic across the
quantum registers to generate the multivariate normal distribution
$\mathcal{N(\bm{\mu},\mathbf{\Sigma})}$ representing the distribution of log-
returns at each timestep. The parameters $(\bm{\mu}$, $\mathbf{\Sigma})$ are
typically determined by historical market values and characterize the
stochastic process of each asset. The price of each asset at each timestep is
then calculated by exponentiating the resulting quantum registers, and the
payoff is computed into another register based on the asset prices, before
encoding the result into the amplitude of an ancilla qubit.

In the QSP approach outlined in Sec. 5.2, we again employ the
reparameterization method to prepare the distribution
$\mathcal{N(\bm{\mu},\mathbf{\Sigma})}$ over the assets’ log-returns. We
subsequently evaluate the payoff conditions of Eq. (41) and Eq. (42) using
ancilla qubits, and encode the binary payoffs to the payoff ancilla qubit
controlled on the corresponding payoff condition qubits. Note that the above
does not require direct knowledge of asset prices, and can all be performed
using only the log-returns, thus we do not need to calculate exponentials of
values stored in binary representation in quantum registers. At the final
timestep, we use QSP to implement the second clause of Eq. (44) which encodes
the payoff directly in the amplitude of the payoff qubit, and as such we
eliminate all instances of quantum register exponentiation in the circuit. A
schematic representation of the quantum circuits implementing the standard
reparameterization method of Ref. [8] and the corresponding QSP approach
outlined in Sec. 5.2 in order to price an autocallable contract is shown In
Fig. 3.

(a)

(b)

Figure 3: Schematic representation of the quantum circuits used to price the
autocallable derivative contract using (a) the standard reparameterization
method of Ref. [8] and (b) the QSP approach described in Sec. 5.2. The two
approaches employ the same recipe to prepare the multivariate normal
distribution $\mathcal{N(\bm{\mu},\mathbf{\Sigma})}$ for each timestep
$t\in[1,T]$ across $d$ underlying assets. The application of QSP in the
evaluation of the payoff via Eq. (44) removes the necessity of directly
computing exponentials of quantum registers, a major contributor to the
quantum resources required to implement the circuit [8]. The quantum circuits
depict only the main quantum registers involved in the calculation and not
ancillary registers that may be required in the subcomponents.

The resource estimates in Ref. [8, 9] for quantum advantage in pricing an
autocallable contract are computed for an instance with $a=3$ assets and
$T=20$ timesteps, and the overall target approximation error in the estimate
is set to $\epsilon=2\times 10^{-3}$. The corresponding accuracy with which
the payoff must be calculated and encoded into the amplitude of an ancilla
qubit is $\epsilon_{p}=10^{-3}$, meaning that the unitary $\mathcal{F}$ of Eq.
(2) must be implemented such that the amplitude of the $\ket{0}$ is accurate
to within $\epsilon_{p}$. In the standard reparameterization circuit of Fig.
3a, this target accuracy dictates the resources required to compute the
exponentials, the quantum arithmetic in the computation of the payoff and the
final rotation into an ancilla qubit using controlled-$Ry$ gates. By virtue of
Eq. (45), the QSP process encompasses the computation of the exponential and
the payoff simultaneously, as well as the amplitude encoding once we have
access to the unitary of Eq. (46) which encodes the value $\sqrt{r_{T}}$ in
the amplitude of an ancilla qubit. The corresponding payoff resources of the
QSP circuit in Fig. 3b are determined by the accuracy of the transformation of
Eq. (45), which in turn depends on the degree $d$ of the polynomial $P(x)$
required to approximate the function $f(x)$ of Eq. (47) such that
$\max|f(x)-P(x)|\leq\epsilon_{p}=10^{-3}$.

We extend the Q# [25] implementation from Ref. [9] to construct the quantum
circuit of Fig. 3b and estimate the corresponding quantum resources when QSP
is used to compute the exponentials and derivative payoff, for both
$U_{\textrm{sin}}$ (Eq. (12)) and $U_{\textrm{sqrt}}$ (Eq. (14)) as the
underlying unitary. In the original implementation, the parameters of the
quantum circuit up to and including the block labeled
$\mathcal{N(\bm{\mu},\mathbf{\Sigma})}$ were chosen to satisfy the target
accuracy requirements and are thus kept the same. As such, the register
$\ket{r_{T}}$ uses the fixed-point representation of Eq. (16) with
$(n,p)=(15,5)$ such that $r_{T}\in[-R,R]$ with $R=2^{5}$. Choosing $s=R$ in
Eq. (46) makes sure the encoded amplitudes are positive, and the
$U_{\textrm{sqrt}}$ operator is constructed using the method described in Sec.
4 with a ripple-carry comparator circuit [26]. Using the polynomial
approximation construction described in Appendix A, we numerically determine
that a polynomial degree $d=20$ suffices to generate a polynomial $P(x)$ which
satisfies $\max|f(x)-P(x)|\leq 10^{-3}$ for all values of $K_{T}\in[0,1]$.
With these parameters, in Fig. 4, we plot the function $f(x)$ and the
generated polynomial approximation $P(x)$, as well as the corresponding
approximation error $|f(x)-P(x)|$ for $K_{T}=1$.

The resource estimation from Ref. [8] was carried out with the assumption that
piecewise polynomial approximation [27] is used to compute the exponentials
required in the circuit. In Ref. [9], the quantum circuits were explicitly
constructed, increasing the accuracy of the estimates, and the exponentials
were computed using a lookup table. The lookup table is parameterized by the
number of swap bits, offering tradeoffs between the circuit depth/width and
the number of qubits required. The quantum resources reported are that of the
$\mathcal{Q}$ operator driving amplitude estimation, where
$\mathcal{Q}=\mathcal{A}S_{0}\mathcal{A}^{\dagger}S_{1}$, for appropriate
reflection operators $S_{0}$ and $S_{1}$ [1], and the $\mathcal{A}$ operator
of Eq. (3). In Table 1, we show the $T$-depth, $T$-count and number of logical
qubits required to implement the $\mathcal{Q}$ operator with (a) the approach
from Ref. [9] where a lookup table with one swap qubit is used to compute the
exponentials, while arithmetic is used to calculate the payoff, and (b) the
approach developed in this work where the exponentials and payoff are computed
using QSP either with the unitary of Eq. (12) or the unitary of Eq. (14). We
notice that the QSP method based on either unitary significantly outperforms
the Arithmetic approach, but because the $U_{\textrm{sin}}$-based QSP includes
$d\cdot n$ total controlled-Ry rotations ($n$ controlled-Ry rotations for each
invocation of $U_{\textrm{sin}}$ and $d$ total invocations of
$U_{\textrm{sin}}$ and $U_{\textrm{sin}}^{\dagger}$ in the QSP circuit), the
total $T$-depth is significantly larger compared to using $U_{\textrm{sqrt}}$.
In the Arithmetic column we use the results from Ref. [9] with one swap qubit
for reference as it gives the most balanced resources for the metrics
considered, but note that the $U_{\textrm{sqrt}}$-based QSP method outperforms
the arithmetic method also in the case where five or ten swap qubits are used
for the lookup table across all three metrics.

The difference in T-depth between the two QSP approaches can be understood by
considering the Clifford+T decomposition of both $U_{\textrm{sin}}$ and
$U_{\textrm{sqrt}}$. The unitary $U_{\textrm{sqrt}}$ can be implemented with a
Toffoli-depth of $2\lfloor\log_{2}(n-1)\rfloor+5$ using a ripple-carry
comparator for an $n$-qubit register [26], which we can also regard as the
equivalent T-depth if we use an ancilla qubit for the T-gate decomposition of
each Toffoli [28]. Thus, $U_{\textrm{sqrt}}$ can be constructed with a T-depth
of $11$ for the $n=15$ case we consider here. On the other hand, the
$U_{\textrm{sin}}$ unitary requires Ry rotations, each of which can be
performed to precision $\epsilon_{R}$ with a circuit of T-depth
$3\log_{2}(1/\epsilon_{R})$ [29]. Since there are $d\cdot n$ total Ry
rotations in the QSP circuit, if we want the total error from all rotations to
be within our total target payoff error
$\epsilon_{t}\leq\epsilon_{p}=10^{-3}$, each Ry rotation needs to be performed
with precision $\epsilon_{R}=\epsilon_{t}/(d\cdot n)\sim 3\times 10^{-6}$.
Hence, each $U_{\textrm{sin}}$ unitary requires a T-depth of $15\cdot
3\log_{2}(1/\epsilon_{R})\sim 818$, about $74$ times larger than that of
$U_{\textrm{sqrt}}$.

Using the Iterative Quantum Amplitude Estimation method [13] to compute the
total resources needed for the end-to-end process to estimate the target
amplitude to within $\epsilon=10^{-3}$ at confidence level $1-\alpha=0.68$
[8], we find that the $U_{\textrm{sqrt}}$-QSP method requires $4.7$k logical
qubits, a total T-depth of $4.5\times 10^{7}$ and T-count of $2.4\times
10^{9}$. With the assumption that this derivative contract can be priced
classically in one second [8], a quantum processor would need to execute
T-gates at a rate of $45$MHz to match the classical performance if QSP with
$U_{\textrm{sqrt}}$ is used, compared to $207$MHz with the arithmetic method
and $430$MHz using QSP and $U_{\textrm{sin}}$.

(a)

(b)

Figure 4: (a) The function $f(x)$ of Eq. (47) with $p=5$, $s=2^{p}$ and $K_{T}=1$ used to compute the payoff of an autocallable contract, and a $d=20$-degree approximating Chebyshev polynomial $P(x)$ generated using the method described in Appendix A. The QSP transformation of Eq. (11) based on the unitary $U_{\textrm{sqrt}}$ of Eq. (46) can then be employed to approximately create the state in Eq. (45), encoding the autocallable payoff directly in a quantum amplitude. (b) The absolute error between $f(x)$ and the approximation $P(x)$ satisfies the target accuracy $\max|f(x)-P(x)|\leq 10^{-3}$. Resource |  Arithmetic [9] | QSP (this work) |  Improvement
---|---|---|---
|  |  $U_{\textrm{sin}}$ [Eq. (12)] |  $U_{\textrm{sqrt}}$ [Eq. (14)] |
$T$-depth |  $36$k |  $75$k |  $7.8$k |  4.7x
$T$-count |  $6.6$M |  $605$k |  $414$k |  $16$x
# Qubits |  $19.2$k |  $4.7$k |  $4.7$k |  $4.1$x
Table 1: Quantum resources required to implement the $\mathcal{Q}$ operator
driving Amplitude Estimation to price the autocallable contract from Sec. 6
with $a=3$ assets and $T=20$ timesteps. The Arithmetic column refers to the
results from Ref. [9] where the exponentials are computed using a lookup table
with one swap qubit, and the payoff using quantum arithmetic [8]. In the QSP
column, the corresponding resources are shown when QSP is employed to compute
both the exponentials and the payoff as described in Sec. 5.2, based on either
Eq. (12) or Eq. (14) as the underlying unitary. In the last column we
highlight the improvement of $U_{\textrm{sqrt}}$-QSP for each resource
compared to the Arithmetic column. For all cases, the resources are calculated
such that all encoded amplitudes are accurate to within $\epsilon_{p}\leq
10^{-3}$.

##  7 Discussion

It has been previously highlighted that quantum arithmetic is a major
bottleneck in pricing financial derivatives on quantum computers [8, 9], and
the QSP-based method we have introduced to construct derivative payoffs
significantly reduces the quantum resources required for quantum advantage
compared to previous methods, in both circuit width and depth. This approach
seems more natural in the following sense: When an application requires the
encoding of a function of a variable into a quantum amplitude, it should be
more efficient to directly transform the amplitude, rather first applying a
digital transformation through arithmetic on computational basis states and
subsequently encoding the result into the amplitude. While QSP restricts the
domain and range of transformations to the interval $[0,1]$, we observe that
in cases like the one studied in this manuscript it is still advantageous
compared to utilizing additional quantum resources to increase the
range/precision of digital operations. We are therefore hopeful that similar
techniques can be explored in other use cases, potentially leading to
additional improvements to the method. While the topic we are addressing in
this manuscript is quite different from that of Ref. [20], the unitary we
introduced in Eq. (14) can also be employed for the state preparation method
introduced in that work, potentially further reducing the resources in that
case as well. Additionally, the availability of open source implementations
for QSP phase factor calculations [30, 31] enables the exploration of QSP for
specific applications more accessible, without requiring deep expertise in the
topic.

One limitation of the QSP method in its current form is that it can only
provide approximations to functions of a single variable. While for the
applications considered in this work this is enough to provide an advantage,
the ability to handle multivariate functions would greatly increase its
applicability. For example, multivariable QSP (M-QSP) could potentially reduce
the quantum resources in derivative pricing further by replacing the quantum
arithmetic in the circuit component labeled
$\mathcal{N(\bm{\mu},\mathbf{\Sigma})}$ in Fig. 3. Moreover, it would enable
the computation of more complicated derivative payoffs which are non-trivial
functions of multiple stochastic variables. Some properties and aspects of
M-QSP have been discussed in Ref. [32], and additional research on this topic
would have potentially far-reaching practical consequences for quantum
applications. We believe that this is a very worthwhile topic for further
research.

##  8 Acknowledgments

We thank Rajiv Krishnakumar, Dave Clader and Farrokh Labib for useful
discussions on derivative pricing and QSP, and Lin Lin for pointing us to
QSPPACK to numerically generate polynomial approximations and phase factors
used in QSP.

## References

  * Brassard  _et al._ [2002] G. Brassard, P. Hoyer, M. Mosca, and A. Tapp, “ _Quantum Amplitude Amplification and Estimation_ ,” Contemporary Mathematics 305 (2002).
  * Montanaro [2015] A. Montanaro, “ _Quantum speedup of Monte Carlo methods_ ,” Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 471 (2015).
  * Rebentrost  _et al._ [2018] P. Rebentrost, B. Gupt, and T. R. Bromley, “ _Quantum computational finance: Monte Carlo pricing of financial derivatives_ ,” Phys. Rev. A 98, 022321 (2018).
  * Woerner and Egger [2019] S. Woerner and D. J. Egger, “ _Quantum risk analysis_ ,” npj Quantum Information 5 (2019).
  * Stamatopoulos  _et al._ [2020] N. Stamatopoulos, D. J. Egger, Y. Sun, C. Zoufal, R. Iten, N. Shen, and S. Woerner, “ _Option Pricing using Quantum Computers_ ,” Quantum 4, 291 (2020).
  * Doriguello  _et al._ [2022] J. a. F. Doriguello, A. Luongo, J. Bao, P. Rebentrost, and M. Santha, “ _Quantum Algorithm for Stochastic Optimal Stopping Problems with Applications in Finance_ ,” in 17th Conference on the Theory of Quantum Computation, Communication and Cryptography (TQC 2022), Leibniz International Proceedings in Informatics (LIPIcs), Vol. 232 (Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl, Germany, 2022) pp. 2:1–2:24.
  * Herbert [2022] S. Herbert, “ _Quantum Monte Carlo Integration: The Full Advantage in Minimal Circuit Depth_ ,” Quantum 6, 823 (2022).
  * Chakrabarti  _et al._ [2021] S. Chakrabarti, R. Krishnakumar, G. Mazzola, N. Stamatopoulos, S. Woerner, and W. J. Zeng, “ _A Threshold for Quantum Advantage in Derivative Pricing_ ,” Quantum 5, 463 (2021).
  * [9] “ _Using Q# to estimate resources needed for quantum advantage in derivative pricing_ ,” Accessed: 2023-06-21.
  * Zoufal  _et al._ [2019] C. Zoufal, A. Lucchi, and S. Woerner, “ _Quantum Generative Adversarial Networks for learning and loading random distributions_ ,” npj Quantum Information 5 (2019).
  * Stamatopoulos  _et al._ [2022] N. Stamatopoulos, G. Mazzola, S. Woerner, and W. J. Zeng, “ _Towards Quantum Advantage in Financial Market Risk using Quantum Gradient Algorithms_ ,” Quantum 6, 770 (2022).
  * Suzuki  _et al._ [2020] Y. Suzuki, S. Uno, R. Raymond, T. Tanaka, T. Onodera, and N. Yamamoto, “ _Amplitude estimation without phase estimation_ ,” Quantum Information Processing 19, 75 (2020).
  * Grinko  _et al._ [2021] D. Grinko, J. Gacon, C. Zoufal, and S. Woerner, “ _Iterative quantum amplitude estimation_ ,” npj Quantum Information 7 (2021).
  * Gilyén  _et al._ [2019] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, “ _Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics_ ,” in Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing (2019) pp. 193–204.
  * Martyn  _et al._ [2021] J. M. Martyn, Z. M. Rossi, A. K. Tan, and I. L. Chuang, “ _Grand Unification of Quantum Algorithms_ ,” PRX Quantum 2 (2021).
  * Low and Chuang [2017] G. H. Low and I. L. Chuang, “ _Optimal Hamiltonian Simulation by Quantum Signal Processing_ ,” Phys. Rev. Lett. 118, 010501 (2017).
  * Martyn  _et al._ [2023] J. M. Martyn, Y. Liu, Z. E. Chin, and I. L. Chuang, “ _Efficient fully-coherent quantum signal processing algorithms for real-time dynamics simulation_ ,” The Journal of Chemical Physics 158, 024106 (2023).
  * Lin and Tong [2020] L. Lin and Y. Tong, “ _Optimal polynomial based quantum eigenstate filtering with application to solving quantum linear systems_ ,” Quantum 4, 361 (2020).
  * Rall and Fuller [2023] P. Rall and B. Fuller, “ _Amplitude Estimation from Quantum Signal Processing_ ,” Quantum 7, 937 (2023).
  * McArdle  _et al._ [2022] S. McArdle, A. Gilyén, and M. Berta, “ _Quantum state preparation without coherent arithmetic_ ,” arXiv preprint arXiv:2210.14892 (2022).
  * Hull [2006] J. Hull,  _Options, futures, and other derivatives_, 6th ed. (Pearson Prentice Hall, Upper Saddle River, NJ [u.a.], 2006).
  * Haah [2019] J. Haah, “ _Product Decomposition of Periodic Functions in Quantum Signal Processing_ ,” Quantum 3, 190 (2019).
  * Chao  _et al._ [2020] R. Chao, D. Ding, A. Gilyen, C. Huang, and M. Szegedy, “ _Finding Angles for Quantum Signal Processing with Machine Precision_ ,” arXiv preprint arXiv:2003.02831 (2020), arXiv:2003.02831 [quant-ph] .
  * Dong  _et al._ [2021] Y. Dong, X. Meng, K. B. Whaley, and L. Lin, “ _Efficient phase-factor evaluation in quantum signal processing_ ,” Physical Review A 103, 042419 (2021).
  * Microsoft [2020] Microsoft,  _Q# Language Specification_ (2020).
  * Draper  _et al._ [2006] T. G. Draper, S. A. Kutin, E. M. Rains, and K. M. Svore, “ _A logarithmic-depth quantum carry-lookahead adder_ ,” Quantum Information and Computation 6, 351 (2006).
  * Häner  _et al._ [2018] T. Häner, M. Roetteler, and K. M. Svore, “ _Optimizing quantum circuits for arithmetic_ ,” arXiv preprint arXiv:1805.12445 (2018).
  * Selinger [2013] P. Selinger, “ _Quantum circuits of T-depth one_ ,” Physical Review A 87 (2013).
  * Ross and Selinger [2016] N. J. Ross and P. Selinger, “ _Optimal Ancilla-Free Clifford+T Approximation of z-Rotations_ ,” Quantum Info. Comput. 16, 901 (2016).
  * [30] “ _QSPPACK_ ,” Accessed: 2023-06-21.
  * [31] “ _pyqsp_ ,” Accessed: 2023-06-21.
  * Rossi and Chuang [2022] Z. M. Rossi and I. L. Chuang, “ _Multivariable quantum signal processing (M-QSP): prophecies of the two-headed oracle_ ,” Quantum 6, 811 (2022).
  * Dong  _et al._ [2022] Y. Dong, L. Lin, and Y. Tong, “ _Ground-State Preparation and Energy Estimation on Early Fault-Tolerant Quantum Computers via Quantum Eigenvalue Transformation of Unitary Matrices_ ,” PRX Quantum 3 (2022).

##  Appendix A Construction of Polynomial Approximations for QSP

In order to approximately apply a singular value transformation given by a
function $f(x)$ using QSP, we first need to find a polynomial which
approximates it. The polynomial transformations possible through the QSP
framework must have definite parity (even or odd), but because we only aim to
transform real positive amplitudes created by the unitary of Eq. (25), we
restrict our attention to the interval $[0,1]$ due to symmetry.

We employ the optimization-based method described in Ref. [33] which generates
near-optimal approximations without relying on any analytic computation. This
method constructs a linear combination of even Chebyshev polynomials with some
unknown coefficients $\\{c_{k}\\}$

$P(x)=\sum_{k=0}^{d/2}c_{k}T_{2k}(x),$ (48)

and aims to find the coefficients $\\{c_{k}\\}$ which give the best
approximation to $f(x)$. This process is formulated as a discrete optimization
problem by discretizing the interval $[0,1]$ using $M$ grid points generated
by the roots of Chebyshev polynomials
$\left\\{x_{j}=-\cos\frac{j\pi}{M-1}\right\\}_{j=0}^{M-1}$ and solving the
following minimax optimization problem for the coefficients $\\{c_{k}\\}$

$\displaystyle\underset{\\{c_{k}\\}}{\min}\quad\left\\{\underset{x_{j}\in[0,1]}{\max}|f(x_{j})-P(x_{j})|\right\\}$
$\displaystyle\textrm{s.t.}\quad
P(x_{j})=\sum_{k}c_{k}A_{jk},\quad|P(x_{j})|\leq c,\quad\textrm{for
}j=0,1,\ldots,M-1,$ (49)

where $c$ is chosen close to $1$ but preferably slightly smaller to avoid
overshooting, and $A_{jk}$ is a matrix of coefficients defined as
$A_{jk}=T_{2k}(x_{j})$ for $k=0,1,\ldots,d/2$.