††thanks: Current affiliation: BQP Advisors LLC††thanks: Current affiliation:
Quantonation

# Quantum Risk Analysis of Financial Derivatives

Nikitas Stamatopoulos  Goldman Sachs, New York, NY  B. David Clader  Goldman
Sachs, New York, NY  Stefan Woerner  IBM Quantum, IBM Research Europe – Zurich
William J. Zeng  Goldman Sachs, New York, NY

###### Abstract

We introduce two quantum algorithms to compute the Value at Risk (VaR) and
Conditional Value at Risk (CVaR) of financial derivatives using quantum
computers: the first by applying existing ideas from quantum risk analysis to
derivative pricing, and the second based on a novel approach using Quantum
Signal Processing (QSP). Previous work in the literature has shown that
quantum advantage is possible in the context of individual derivative pricing
and that advantage can be leveraged in a straightforward manner in the
estimation of the VaR and CVaR. The algorithms we introduce in this work aim
to provide an additional advantage by encoding the derivative price over
multiple market scenarios in superposition and computing the desired values by
applying appropriate transformations to the quantum system. We perform
complexity and error analysis of both algorithms, and show that while the two
algorithms have the same asymptotic scaling the QSP-based approach requires
significantly fewer quantum resources for the same target accuracy.
Additionally, by numerically simulating both quantum and classical VaR
algorithms, we demonstrate that the quantum algorithm can extract additional
advantage from a quantum computer compared to individual derivative pricing.
Specifically, we show that under certain conditions VaR estimation can lower
the latest published estimates of the logical clock rate required for quantum
advantage in derivative pricing by up to $\sim 30$x. In light of these
results, we are encouraged that our formulation of derivative pricing in the
QSP framework may be further leveraged for quantum advantage in other relevant
financial applications, and that quantum computers could be harnessed more
efficiently by considering problems in the financial sector at a higher level.

##  I Introduction

Ever since it was shown that quantum speedups were possible for Monte Carlo
methods [1], various methods have been proposed to apply the potential
benefits of quantum computing to the pricing of financial derivatives [2, 3,
4, 5]. While these quantum methods allow for a quadratic speedup over state of
the art classical approaches, the quantum resources required for practical
advantage are significant not just in terms of circuit depth and number of
qubits, but also in terms of the logical clock rate needed to match the
performance of modern classical computers [6]. On the other hand, the
application of quantum gradient methods utilizing the quantum oracles used for
pricing [7], showed that the logical clock rate required for advantage in the
task of computing the market risk of financial derivatives might be lower than
that of derivative pricing itself. As such, considering the quantum oracles
used in derivative pricing as components of higher-level algorithms which
compute other quantities of interest of financial derivatives could lower
certain quantum resources required for practical advantage.

Value at Risk (VaR) is a metric that estimates the maximum possible financial
depreciation of an asset (or collection of assets) within a specific time
frame at a certain confidence level [8, 9]. A related quantity, the
Conditional Value at Risk (CVaR) measures the expected depreciation of an
asset should losses exceed the Value at Risk, providing a risk metric for the
tail of the loss distribution. Financial firms make extensive use of the VaR
and CVaR metrics for their derivative holdings, aimed at assessing the
financial health of their balance sheet as well as to set trading limits in
order to minimize potential exposure to adverse market conditions. When
individual derivative contracts are priced using Monte Carlo simulation, the
computation of the VaR of portfolios of derivatives requires pricing each
contract multiple times over different possible market states and then
estimating the VaR from the resulting prices. For $N$ such market states, the
complexity of estimating the VaR thus scales as
$\mathcal{O}\left(N/\epsilon^{2}\right)$, where $\epsilon$ is the accuracy
with which each contract is priced. The same approach can be used to estimate
the VaR by pricing of each derivative on a quantum computer and the quadratic
speedup in the individual derivative pricing would carry through to the
overall complexity, improving it to $\mathcal{O}\left(N/\epsilon\right)$. In
Ref. [3, 10], a quantum method based on Quantum Amplitude Estimation (QAE) was
introduced to estimate the VaR in cases where potential future losses of a
financial asset can be modeled explicitly. While this method was shown to also
provide a quadratic advantage compared to classical VaR estimation using Monte
Carlo, it cannot be readily applied to cases where the financial assets in
question are derivative contracts which are individually priced using sampling
methods, and estimating the potential losses requires repeated pricing of the
contracts under different market conditions.

In this article, we introduce two quantum algorithms to estimate the VaR and
CVaR of derivative portfolios. The first is an extension of the QAE-based
method of Ref. [3, 10]. The second is a novel algorithm based on the Quantum
Signal Processing (QSP) framework [11, 12, 13, 14], which enables nearly
arbitrary polynomial transformations to quantum (sub)systems with minimal
overheads. Various quantum algorithms have been formulated in this framework
[15], and in most instances the QSP formulation has resulted in lower resource
requirements for the algorithms [16, 17, 18], including in quantum derivative
pricing [19]. Moreover, it has been shown that the quantum oracles employed in
derivative pricing naturally represent block-encodings of derivative prices
across different input parameters [20, 7]. Because the QSP framework relies on
efficiently constructing block-encodings of quantities of interest, we hence
examine how QSP can be harnessed for problems around derivative pricing.

Leveraging the QSP framework and derivative pricing oracles, we show how to
apply coherent transformations to the price of a financial derivative across
multiple inputs in superposition, and create tailored transformations that
allow us to estimate the VaR and CVaR of derivative portfolios. We perform
complexity and error analysis of the two quantum algorithms introduced and
compare their performance through explicit construction of the corresponding
quantum circuits and numerical simulations. Finally, we benchmark the QSP-
based algorithm against classical VaR estimation methods and determine that
quantum advantage is possible if suitable techniques to prepare and sample
from the required probability distributions are available.

In Sec. II we formally introduce the problem of estimating the VaR and CVaR of
derivatives and give an overview of the QAE-based VaR estimation method from
Ref. [3, 10]. In Sec. III we generalize this method to compute the VaR of a
portfolio of derivatives. In Sec. IV we give an overview of Quantum Signal
Processing, show how it can be used to perform transformations to
appropriately encoded derivative prices and introduce the QSP-based VaR and
CVaR estimation algorithm. In Sec. VI we compare the performance of the QAE
and QSP methods in terms of the quantum resources required for a target VaR
accuracy, and in Sec. VII we explore the possibility of quantum advantage in
estimating the VaR of derivative portfolios. We discuss our results and next
steps in Sec. VIII.

##  II Value at Risk

The VaR computation of derivative portfolios requires modeling the possible
future behavior of the market parameters underlying the derivatives. While the
pricing of derivatives relies on the value of underlying parameters today
(e.g. stock prices, interest rates, volatilities, etc.), which are available
(or can be estimated) from public markets, computing the VaR of derivatives
requires possible values of these parameters in the future. Two common
modeling methods for VaR computation consist of a) using Monte Carlo
simulation of the market parameters according to an appropriately defined
underlying stochastic process, and b) drawing historical samples of
(correlated) values for all relevant market parameters. Each such modeled
market state is also called a _scenario_ and can be defined as a collection of
_tweaks_ to today’s market parameters. To price a derivative under a scenario,
the tweaks are applied to the initial market state and the derivative is
priced as if the tweaked market state was today’s market. Examples of
scenarios defined this way could be

  1. 1.

Tweak asset X’s price up by 5%

  2. 2.

Tweak asset X’s price up by 5% and asset Y’s price down by 4%

  3. 3.

Tweak asset X’s price up by 5%, all points on its volatility surface up by 10%
and its correlation to asset Y by 20%

In both Monte Carlo and historical VaR methods, the value of the derivative
portfolio is computed under every scenario and from the resulting distribution
of gains/losses, maximum levels of potential loss can be derived at different
confidence levels.

More concretely, in order to compute the VaR of a derivative portfolio at a
confidence level $\alpha\in[0,1]$ (where $\alpha\geq 0.9$ is typically used),
we price the portfolio under a set of $N$ scenarios
$[s_{1},s_{2},\cdots,s_{N}]$ to get portfolio values
$\mathcal{V}=\\{V(s_{i})\\}$, and find the smallest portfolio value in the set
such that smaller values have a probability of occurrence greater than or
equal to $1-\alpha$

$V_{\alpha}=\textrm{inf }\\{x\;|\;\sum_{\\{V\in\mathcal{V}|V\leq
x\\}}\mathbb{P}[V]\geq 1-\alpha\\},$ (1)

where $\mathbb{P}[V]$ denotes the probability that the portfolio has value $V$
under the evaluated scenarios. The VaR of the portfolio at the $\alpha$
confidence level is then defined as $\textrm{VaR}_{\alpha}[V]\equiv
V_{0}-V_{\alpha}$, where $V_{0}$ denotes the current value of the portfolio.
This approach for computing the VaR of a portfolio can be used in both the
Monte Carlo and historical methods, the only difference being the way
probabilities are assigned to possible future market realizations.

The CVaR of a derivative portfolio measures the expected loss of the portfolio
beyond the VaR value and can be estimated using the same set of portfolio
values $\mathcal{V}$ calculated for the VaR, by computing the quantity

$C_{\alpha}=\frac{1}{1-\alpha}\sum_{\\{V\in\mathcal{V}|V\leq
V_{\alpha}\\}}\mathbb{P}[V]\cdot V,$ (2)

with $V_{\alpha}$ given by Eq. (1). The CVaR of the portfolio at the $\alpha$
confidence level is then given by $\textrm{CVaR}_{\alpha}[V]\equiv
V_{0}-C_{\alpha}$.

Classically, the $\textrm{VaR}_{\alpha}$ of a portfolio can be estimated by
pricing the portfolio under each of the $N$ scenarios, sorting the resulting
prices from smallest to largest, and picking the value at the $1-\alpha$
percentile as the VaR estimate. When the value of the derivative portfolio is
computed classically using Monte Carlo to accuracy $\epsilon_{p}$, the VaR
computation requires separate Monte Carlo pricings for each scenario, and
therefore the complexity of the algorithm scales as
$\mathcal{O}\left(N/\epsilon_{p}^{2}\right)$.

A quantum algorithm for computing the VaR of financial assets was first
presented in Ref. [3] and later refined and applied to the use case of
estimating capital requirements in order to manage counterparty credit risk in
Ref. [10]. The algorithm seeks to estimate the VaR of a financial asset at the
$\alpha$ confidence level, through direct access to a unitary modeling the
potential loss as a random variable $L$. This unitary $\mathcal{L}$ creates a
superposition

$\mathcal{L}:\ket{\vec{0}}\rightarrow\sum_{\ell}\sqrt{p_{\ell}}\ket{\ell}$ (3)

where the $\ket{\ell}$ register spans over possible values of $L$ encoded in
an appropriate binary representation, and $p_{\ell}$ the corresponding
probability of occurrence. Picking a VaR candidate $l$ and applying a unitary
$\mathcal{C}$ which performs the binary comparison

$\mathcal{C}:\ket{\ell}\ket{0}\rightarrow\ket{\ell}\ket{\ell<l},$ (4)

we get the state

$\sum_{\ell\leq
l}\sqrt{p_{\ell}}\ket{\ell}\ket{0}+\sum_{\ell>l}\sqrt{p_{\ell}}\ket{\ell}\ket{1}.$
(5)

The probability of measuring the last qubit in the $\ket{0}$ state gives us
the probability that the loss $L$ is smaller than $l$,
$\mathbb{P}[\ket{0}]=\mathbb{P}[L\leq l]$. Using QAE to estimate
$\mathbb{P}[\ket{0}]$ and applying a bisection search to find the smallest
$l_{\alpha}$ such that $\mathbb{P}[L\leq l_{\alpha}]\geq 1-\alpha$ gives us
the VaR of the asset at the $\alpha$ confidence level,
$l_{\alpha}=\textrm{VaR}_{\alpha}[L]$.

This quantum method is not directly applicable to the calculation of
derivative VaR through the evaluation of scenario prices because in that case
we do not have direct access to the unitary of Eq. (3). In the following
section we show how to extended it in order to also be able to compute the VaR
of derivatives by incorporating one additional application of QAE.

##  III Value At Risk with Quantum Amplitude Estimation

Quantum methods for the VaR calculation of financial derivatives rely on the
existence of an oracle $\mathcal{A}$ that encodes the price of a derivative
$V$ into an amplitude

$\mathcal{A}:\ket{\vec{0}}\rightarrow\left(\sqrt{V}\ket{\psi_{0}}\ket{0}+\sqrt{1-V}\ket{\psi_{1}}\ket{1}\right),$
(6)

for arbitrary, normalized states $\ket{\psi_{0}}$, $\ket{\psi_{1}}$. Different
ways of constructing this oracle have been proposed in the literature [2, 4,
6]. In Appendix A, we describe in more detail the method proposed in [6] and
show how it can be extended to handle the case where we want Eq. (6) to encode
the value of a portfolio of derivatives instead of the price of a single
derivative.

In order to compute the VaR of financial derivatives in the way described in
Sec. II, we construct scenarios consisting of tweaks to (classical) market
data parameters which can be encoded in some appropriate way as a
computational basis state $\ket{s}$. For example, a scenario which
simultaneously tweaks $d$ underlying parameters (also called _risk factors_),
can be encoded as

$\ket{s}=\ket{t_{1}t_{2},\dots t_{d}},$ (7)

where $t_{i}$ is a binary encoding of the tweak amount to underlying $i$. For
such scenario encoding $\ket{s}$, we can generalize the form of operator
$\mathcal{A}$ of Eq. (6), such that it encodes the value of the portfolio
under scenario $s$

$\mathcal{A}:\ket{s}_{q}\ket{0}_{n}\ket{0}\rightarrow\ket{s}_{q}\left(\sqrt{V(s)}\ket{\psi_{0}^{s}}_{n}\ket{0}+\sqrt{1-V(s)}\ket{\psi_{1}^{s}}_{n}\ket{1}\right),$
(8)

for normalized quantum states $\ket{\psi_{0}^{s}}_{n}$ and
$\ket{\psi_{1}^{s}}_{n}$. Let $\mathcal{S}$ be an operator which loads $N$
scenarios $[s_{1},s_{2},\dots,s_{N}]$ encoded this way in superposition,
weighted by a probability denoting the likelihood of its occurrence

$\mathcal{S}:\ket{0}_{q}\rightarrow\sum_{i=0}^{N-1}\sqrt{p(s_{i})}\ket{s_{i}}_{q}.$
(9)

The operator $\mathcal{A}_{\mathcal{S}}\equiv\mathcal{A}\mathcal{S}$ then
encodes the portfolio value under all scenarios in superposition

$\mathcal{A}_{\mathcal{S}}\equiv\mathcal{A}\mathcal{S}:\ket{\vec{0}}\rightarrow\sum_{i=0}^{N-1}\sqrt{p(s_{i})}\ket{s_{i}}_{q}\left(\sqrt{V(s_{i})}\ket{\psi_{0}^{s_{i}}}_{n}\ket{0}+\sqrt{1-V(s_{i})}\ket{\psi_{1}^{s_{i}}}_{n}\ket{1}\right).$
(10)

Let us for now assume that each value $V(s_{i})$ can be exactly represented
using $m$ qubits. The circuitry of QAE (up to measurement) can then be applied
with the marked state identified as the last qubit being in the $\ket{0}$
state to get a binary representation of each $V(s_{i})$ into a quantum
register

$\sum_{i=0}^{N-1}\sqrt{p(s_{i})}\ket{s_{i}}_{q}\ket{V(s_{i})}_{m}\ket{\textrm{garbage}}$
(11)

where the state of the remaining registers $\ket{\textrm{garbage}}$ will
depend on the choice of QAE variant used, but is not involved in the remainder
of the process and can be ignored. Then the comparator of Eq. (4) can be
applied to the $\ket{V(s_{i})}$ register with a VaR candidate $\mu$ and an
ancilla qubit initially in the $\ket{0}$ state to obtain

$\sum_{\\{i\>|\>V(s_{i})\leq\mu\\}}\sqrt{p(s_{i})}\ket{s_{i}}_{q}\ket{V(s_{i})}_{m}\ket{0}+\sum_{\\{i\>|\>V(s_{i})>\mu\\}}\sqrt{p(s_{i})}\ket{s_{i}}_{q}\ket{V(s_{i})}_{m}\ket{1}.$
(12)

The probability of the last qubit being in the $\ket{0}$ state gives us the
probability that the portfolio value $V$ is smaller that $\mu$ across all
scenarios considered, $\mathbb{P}[\ket{0}]=\mathbb{P}[V\leq\mu]$. We can thus
employ consecutive rounds of QAE to find the smallest $\mu_{\alpha}$ such that
$\mathbb{P}[V\leq\mu_{\alpha}]\geq 1-\alpha$, which can be achieved with a
bisection search over values of $\mu$. Assuming the portfolio value today is
$V_{0}$, the VaR at the $\alpha$ confidence level will be given by
$\textrm{VaR}_{\alpha}[V]\equiv V_{0}-\mu_{\alpha}$. Note that if we use the
_canonical_ amplitude estimation [21] circuitry to generate Eq. (11), the
value encoded in the resulting register will not be $V(s_{i})$ but rather
$\theta_{i}=\arcsin(\sqrt{V(s_{i})})$ 111Canonical QAE gives an equal
superposition of two values $\theta_{1}=\arcsin(\sqrt{V})\in[0,\pi/2]$ and
$\theta_{2}=\pi-\theta_{1}\in[\pi/2,\pi]$. Quantum arithmetic can be used to
perform $\theta\rightarrow\pi-\theta$ if $\theta\in[\pi/2,\pi]$ which will
ensure the register represents a value $\theta\in[0,\pi/2]$ such that
$V=\sin^{2}(\theta)$.. While quantum arithmetic can be used to compute a
binary approximation to $V(s_{i})$ from $\theta_{i}$, since we are only
interested performing the comparison $V(s_{i})\leq\mu$ with
$V(s_{i})\in[0,1]$, we can adjust the comparator circuit to equivalently
perform $\theta_{i}\leq\arcsin(\sqrt{\mu})$. A diagram of the circuit required
to generate Eq. (12) using canonical QAE is shown in Fig. 1.

Figure 1: Quantum circuit required to prepare the state in Eq. (12). Given a
VaR candidate $\mu$, the probability that the portfolio value is smaller than
$\mu$ is encoded into the amplitude of a qubit being in the $\ket{0}$ state.
First, the $\mathcal{S}$ and $\mathcal{A}$ operators of Eq. (9) and Eq. (8)
respectively generate the superposition of all scenarios and encode the value
of the portfolio under each scenario into the amplitude of the bottom qubit
being in the $\ket{0}$ state. The $\oslash$ symbol indicates that an entire
register acts as control to a unitary. Then, repeated controlled applications
of an operator $\mathcal{Q}$ are applied with
$\mathcal{Q}=\mathcal{A}S_{0}\mathcal{A}^{\dagger}S_{1}$,
$S_{0}=\mathbb{I}-2\ket{\vec{0}}\bra{\vec{0}}$,
$S_{1}=\mathbb{I}-2\ket{1}\bra{1}$ where $S_{0}$ reflects about the bottom two
registers being all in the $\ket{0}$ state and $S_{1}$ reflects around the
bottom qubit in the $\ket{1}$ state. After an inverse Quantum Fourier
Transform, the $m$ qubits acting as controls to the $\mathcal{Q}$ applications
contain an $m$-qubit approximation to the portfolio value under each scenario
in superposition, as shown in Eq. (11). Finally, a comparator circuit
$C_{\mu}$ compares the binary value $V$ encoded in the $m$ qubits to $\mu$ and
flips the top qubit if $V>\mu$, generating Eq. (12). Estimating the
probability of the top qubit being in the $\ket{0}$ state using QAE and re-
running this circuit for different values of $\mu$ looking for the smallest
value which gives $\mathbb{P}[\ket{0}]\geq 1-\alpha$ will give us a VaR
estimate at the $\alpha$ confidence level.

In practice, we will not be able to represent all values $V(s_{i})$ exactly in
binary, and instead of Eq. (11), for each scenario value we will get a
superposition over all $M=2^{m}$ possible states representable using $m$
qubits, where each state will be weighed by a probability depending on the
value of $V(s_{i})$

$\sum_{i=0}^{N-1}\sqrt{p(s_{i})}\ket{s_{i}}_{q}\sum_{j=0}^{M-1}\sqrt{q_{j}(V(s_{i}))}\ket{\tilde{V}_{j}}_{m},$
(13)

with the probabilities given by [21]

$q_{j}(x)=\frac{\sin^{2}\left(M(\tilde{V}_{j}-x)\pi\right)}{M^{2}\sin^{2}\left((\tilde{V}_{j}-x)\pi\right)},$
(14)

which will peak around the values $\tilde{V}_{j}$ closest to each $V(s_{i})$.
Applying the comparator unitary to perform
$\ket{\tilde{V}_{j}}\ket{0}\rightarrow\ket{\tilde{V}_{j}}\ket{\tilde{V}_{j}\leq\mu}$
will give

$\sum_{i=0}^{N-1}\sqrt{p(s_{i})}\ket{s_{i}}_{q}\left(\sum_{\\{j\>|\>\tilde{V}_{j}\leq\mu\\}}\sqrt{q_{j}(V(s_{i}))}\ket{\tilde{V}_{j}}_{m}\ket{0}+\sum_{\\{j\>|\>\tilde{V}_{j}>\mu\\}}\sqrt{q_{j}(V(s_{i}))}\ket{\tilde{V}_{j}}_{m}\ket{1}\right),$
(15)

meaning that every scenario value $V(s_{i})$ will (incorrectly) contribute
something to the probability of $\ket{1}$ if $V(s_{i})\leq\mu$ and to
$\ket{0}$ if $V(s_{i})>\mu$. Increasing the number of qubits used for this
_inner_ application of QAE will limit the probability that is erroneously
“allocated” to the wrong flag state.

Note that the QAE method used to produce Eq. (13) needs to coherently
approximate $V(s_{i})$, meaning highly performant iterative QAE methods [22,
23, 24] cannot be employed for this task. However, the original QAE
formulation [21] as well as the QAE variant introduced in Ref. [25] can be
used for coherent estimation. While the former requires the Quantum Fourier
Transform (QFT) and overhead in terms of ancilla qubits, the latter can
produce the approximation coherently with comparable performance to the
iterative versions.

##  IV Value At Risk with Quantum Signal Processing

###  IV.1 Quantum Signal Processing

Quantum Signal Processing is a technique which performs polynomial
transformations to the singular values of a matrix $A$ that has been _block-
encoded_ into a unitary operator 222This technique was originally introduced
as Quantum Singular Value Transformation in Ref. [12], but here we use the
nomenclature of Ref. [15], which refers to the underlying method as Quantum
Signal Processing, encompassing a variety of applications.. Specifically,
given projectors $\Pi$ and $\tilde{\Pi}$, we say that a unitary $U$ acting on
$n$ qubits is a block-encoding of a matrix $A$ if $A=\tilde{\Pi}U\Pi$, such
that the projectors $\Pi$ and $\tilde{\Pi}$ determine the location of $A$ in
$U$

$\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]$}}.$ (17)

QSP consists of consecutive applications 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})}$, which for phase
factors $\vec{\phi}=(\phi_{1},\phi_{2},\cdots,\phi_{d})$ induce a polynomial
transformation on $A$ such that [15, 12]

$A=\sum_{k}\sigma_{k}\ket{w_{k}}\bra{v_{k}}\rightarrow
P(A)=\sum_{k}P(\sigma_{k})\ket{w_{k}}\bra{v_{k}},$ (18)

where $\sigma_{k}$ denote the singular values of $A$ while
$\\{\ket{w_{k}}\\}$, $\\{\ket{v_{k}}\\}$ are the left and right singular
vectors of $A$ respectively, and $P$ is a $d$-degree polynomial. When $A$ is
Hermitian, in which case the singular values of the matrix coincide with its
eigenvalues, Quantum Signal Processing effectively performs a polynomial
transformation to the eigenvalues of the matrix.

More formally, define the operators

$\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}$ (19)

$\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).$
(20)

If $U$ is an $n$-qubit unitary, the polynomial transformation achieved through
QSP is implemented with [15]

$U^{\vec{\phi}}=(\mathbb{I}^{\otimes n}\otimes
H)\mathcal{U}_{C}^{\vec{\phi}}(\mathbb{I}^{\otimes n}\otimes
H)=\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]$}}$ (21)

where $\Pi^{\prime}=\tilde{\Pi}$ for odd $d$ and $\Pi^{\prime}=\Pi$ for even
$d$.

In practice, $d$ invocations of $U$ perform a $d$-degree polynomial
transformation to $A$ and the choice of phase factors
$\vec{\phi}=(\phi_{1},\phi_{2},\cdots,\phi_{d})$ determine the exact
polynomial that is applied. Remarkably, the reverse is also true: for
$a\in[-1,1]$ and any real polynomial $P\in\mathbb{R}(a)$, there exists a
sequence of QSP phase factors $\vec{\phi}=(\phi_{1},\phi_{2},\cdots,\phi_{d})$
for which Eq. (21) holds, as long as $\textrm{deg}(P)\leq d$, $|P(a)|\leq
1,\forall a\in[-1,1]$ and $P$ either even or odd. Additionally, for any
function $f(a)$ satisfying these conditions we can first identify a polynomial
approximation to the function, determine the QSP phase factors $\vec{\phi}$
for that polynomial, and use QSP to approximately apply $f(a)$ to the singular
values of the block-encoded matrix. Various methods have been proposed to
generate polynomial approximations to generic functions, where the
approximation is either constructed analytically [15, 26, 27] or with
optimization-based numerical methods [28, 29].

###  IV.2 Algorithm for Value at Risk

Defining $\tilde{\Pi}\equiv I^{\otimes(n+q)}\otimes\ket{0}\bra{0}$ and
$\Pi\equiv I^{\otimes q}\otimes\left(\ket{0}\bra{0}\right)^{\otimes(n+1)}$,
observe that $\tilde{\Pi}\mathcal{A}\Pi$ with $\mathcal{A}$ given by Eq. (8),
is a rank-1 matrix with a single non-trivial singular value $\sqrt{V(s)}$
[12]. Therefore, we can use QSP with the operator $\mathcal{A}$ acting as the
block encoding unitary of Eq. (17) to apply polynomial transformations to the
values $\sqrt{V(s)}$ as discussed in the previous section. After the
superposition of scenarios is created with the unitary $\mathcal{S}$ of Eq.
(9), the unitary $U^{\vec{\phi}}$ of Eq. (21) based on
$\mathcal{U}^{\vec{\phi}}$ with $U=\mathcal{A}$ produces the state

$U^{\vec{\phi}}\mathcal{S}\ket{0}_{q}\ket{0}_{n+1}\ket{0}=\sum_{i=0}^{N-1}\sqrt{p(s_{i})}\ket{s_{i}}_{q}\begin{cases}\left(P(\sqrt{V(s_{i})})\ket{\psi_{0}}_{n}\ket{0}_{2}+\sqrt{1-P(\sqrt{V(s_{i})})^{2}}\ket{\psi_{1}}_{n}\ket{0_{\perp}}_{2}\right),&\text{for
odd $d$,}\\\ &\\\
\left(P(\sqrt{V(s_{i})})\ket{0}_{n+2}+\sqrt{1-P(\sqrt{V(s_{i})})^{2}}\ket{0_{\perp}}_{n+2}\right),&\text{for
even $d$,}\end{cases}$ (22)

for normalized quantum states $\ket{\psi_{0}}_{n}$, $\ket{\psi_{1}}_{n}$, and
$\ket{0_{\perp}}_{a}$ denoting a normalized state orthogonal to $\ket{0}_{a}$.
The polynomial $P(x)$ is determined by the choice of phase factors
$\vec{\phi}$.

If we can find phase factors such that $P(x)$ is an approximation to the
threshold function

$\theta_{\mu}(x)=\begin{cases}1,\;{x}\leq{\mu}\\\ 0,\;{x}>{\mu},\end{cases}$
(23)

Eq. (22) will give us

$\sum_{\\{i\>|\>\sqrt{V(s_{i})}\leq\mu\\}}\sqrt{p(s_{i})}\ket{G}+\sum_{\\{i\>|\>\sqrt{V(s_{i})}>\mu\\}}\sqrt{p(s_{i})}\ket{G_{\perp}},$
(24)

where $\ket{G}=\ket{0}_{2}$ for odd $d$ and $\ket{G}=\ket{0}_{n+2}$ for even
$d$, ignoring the (normalized) state of the remaining registers for clarity.
The probability of measuring $\ket{G}$ will then be equal to the probability
that the portfolio value across all scenarios is smaller than $\mu^{2}$

$\mathbb{P}(\ket{G})=\sum_{\\{s\>|\>V(s)\leq\mu^{2}\\}}p(s).$ (25)

Therefore, the smallest value of $\mu=\mu_{\alpha}$ for which the probability
in Eq. (25) satisfies $\mathbb{P}(\ket{G})\geq 1-\alpha$ gives us the value of
$V_{\alpha}=\mu_{\alpha}^{2}$ we are looking for in Eq. (1) in order to
estimate the VaR of the portfolio at the $\alpha$ confidence level. We can
find such $\mu_{\alpha}$ by performing bisection search over values of
$\mu\in[0,1]$ and applying QAE to estimate the probability in Eq. (25) in each
iteration. The bisection search ends when we find a value of $\mu$ for which
our estimate of $\mathbb{P}(\ket{G})$ is $\epsilon$-close to $1-\alpha$ for
some value of $\epsilon$ determined by the target accuracy of the QAE method
employed. By the union bound, if each round of QAE is performed at confidence
level $1-\alpha_{A}$, the confidence of our final VaR estimate after $k$
rounds will be $1-k\cdot\alpha_{A}$.

This approach can be optimized by observing that not every bisection step
needs to use the same accuracy $\epsilon$ or confidence level parameter
$\alpha_{A}$ for QAE. Because at every step we need to estimate whether the
value $\mathbb{P}(\ket{G})$ in Eq. (25) is $\epsilon$-close to $1-\alpha$, we
can employ QAE with accuracy $\epsilon_{k}>\epsilon$ in cases where
$\mathbb{P}(\ket{G})$ is significantly far from $1-\alpha$. Using iterative
amplitude estimation routines such as IQAE [22], we can adaptively increase
the accuracy of the estimate by taking progressively more samples until we can
determine whether $|\mathbb{P}(\ket{G})-(1-\alpha)|\leq\epsilon$ with some
confidence $1-\alpha_{A}$ and then decide on the next bisection search step
accordingly. This observation is illustrated in Fig. 2.

Figure 2: For each round of bisection search in Algorithm 1, we are trying to
determine whether for a given value $\mu$, the value $\mathbb{P}(\ket{G})$ of
Eq. (25) is $\epsilon$-close to $1-\alpha$. However, the accuracy with which
we estimate the value $\mathbb{P}(\ket{G})$ can be relaxed depending on how
far $\mathbb{P}(\ket{G})$ is from $1-\alpha$. For example, for the values
$\mathbb{P}_{\mu_{1}}(\ket{G})$ and $\mathbb{P}_{\mu_{2}}(\ket{G})$
corresponding to some thresholds $\mu_{1},\mu_{2}$ shown above, each
probability only needs to be estimated with accuracy
$\epsilon_{1},\epsilon_{2}>\epsilon$ respectively in order to proceed with the
bisection search.

The confidence level $1-\alpha_{A}$ of QAE can similarly be adjusted in each
round. If the target confidence level of the VaR estimation is denoted as CL,
by the union bound the confidence level $\alpha_{k}$ in the $k$-th QAE
invocation can be adjusted as long as $\textrm{CL}=1-\sum_{k}\alpha_{k}$ is
satisfied. Therefore, the total probability of failure $1-\textrm{CL}$ can be
distributed across each QAE round as desired. One possible choice is to
allocate higher probability of failure when the QAE target error
$\epsilon_{k}$ is small, which would decrease the constant factor of the
$\mathcal{O}\left(1/{\epsilon_{k}}\right)$ complexity of QAE when
$1/{\epsilon_{k}}$ is large. The pseudocode for this QSP-based algorithm is
given in Algorithm 1 and the circuit required to generate Eq. (25) is shown in
Fig. 3.

Figure 3: Circuit diagram depicting one iteration of the QSP-based VaR
estimation method described algorithm in Sec. IV.2 and outlined in Algorithm
1. Given a VaR candidate $\mu$, the probability that the portfolio value is
smaller than $\mu^{2}$ is encoded into the amplitude of a marked state. The
$\mathcal{S}$ operator of Eq. (9) is first applied to generate the
superposition of all scenarios considered. Then, a $d$-degree polynomial
approximation to the threshold function $\theta_{\mu}(x)$ of Eq. (23) along
with the corresponding phase factors
$\vec{\phi}^{\mu}=(\phi_{1}^{\mu},\phi_{2}^{\mu},\cdots,\phi_{d}^{\mu})$ are
computed, and the generated polynomial transformation is implemented using the
QSP unitary of Eq. (22), which consists of alternating applications of
controlled phase rotations and the block-encoding operator $\mathcal{A}$ of
Eq. (8). After the transformation, the probability of measuring
$\ket{G}=\ket{0}_{2}$ in the bottom two qubits if $d$ is odd or
$\ket{G}=\ket{0}_{n+2}$ in the bottom $n+2$ qubits if $d$ is even, gives us
the probability that the portfolio value across all scenarios is smaller than
$\mu^{2}$ as shown in Eq. (25). Estimating this probability using QAE and re-
running this circuit for different values of $\mu$ looking for the smallest
value which gives $\mathbb{P}[\ket{G}]\geq 1-\alpha$ will give us a VaR
estimate at the $\alpha$ confidence level.

Input: Amplitude Estimation accuracy $\epsilon_{A}$

Input: VaR confidence level $\alpha$

1 Set $\mu_{l}=0$ and $\mu_{h}=1$

2 Set $k=0$

3 while  _True_ :

4 Apply the $S$ operator of Eq. (9) to generate the superposition of all
scenarios used to compute the VaR of a derivative portfolio whose value today
is $V_{0}$ $\sum_{i=0}^{N-1}\sqrt{p(s_{i})}\ket{s_{i}}$

5 Pick $\mu=\mu_{l}+(\mu_{h}-\mu_{l})/2$

6 Generate a $d$-degree polynomial and corresponding phase factors
$\vec{\phi}^{\mu}=(\phi_{1}^{\mu},\phi_{2}^{\mu},\cdots,\phi_{d}^{\mu})$ such
that the polynomial transformation $P(A)$ in Eq. (21) is an approximation to
the threshold function $\theta_{\mu}(x)$ of Eq. (23).

7 Apply the generated polynomial transformation using QSP and the block-
encoding operator $\mathcal{A}$ of Eq. (8) such that the probability of the
last qubits being in the $\ket{G}$ state with $\ket{G}=\ket{0}_{2}$ for odd
$d$ and $\ket{G}=\ket{0}_{n+2}$ for even $d$ becomes
$\mathbb{P}(\ket{G})=\sum_{\\{s\>|\>V(s)\leq\mu^{2}\\}}p(s)$

8 Iteratively estimate the probability $\mathbb{P}(\ket{G})$ using QAE with
decreasing estimation error $\epsilon_{k}\geq\epsilon_{A}$ to determine bounds
satisfying $\mathbb{P}(\ket{G})\in[p_{l},p_{h}]$ with $p_{h}-p_{l}\leq
2\epsilon_{k}$ at confidence level $1-\alpha_{k}$, until we either determine
that $1-\alpha\notin[p_{l},p_{h}]$ or the target estimation error reaches
$\epsilon_{k}=\epsilon_{A}$.

9 Set $k=k+1$

10 if  _$p_{h} <1-\alpha$_ :

11 Set $\mu_{l}=\mu$

12 elif  _$p_{l} >1-\alpha$_ :

13 Set $\mu_{h}=\mu$

14 else:

15 Return $\textrm{VaR}_{\alpha}$ estimate $V_{0}-\mu^{2}$ with confidence
$1-\sum_{k}\alpha_{k}$

16

Algorithm 1 VaR Estimation using Quantum Signal Processing

###  IV.3 Conditional Value at Risk

Once we have estimated the VaR, the CVaR of derivatives can be also be
estimated using QSP by applying a modified threshold function which is linear
below the threshold and zero above

$\tilde{\theta}_{\mu}(x)=\begin{cases}x,\;{x}\leq{\mu}\\\
0,\;{x}>{\mu},\end{cases}$ (26)

at $\mu=\mu_{\alpha}$ where $\mu_{\alpha}$ is the threshold estimated in the
VaR calculation such that the probability in Eq. (25) equals $1-\alpha$. Using
QSP to apply the transformation of Eq. (22) with $P(x)$ now a polynomial
approximation to this function gives us

$\sum_{\\{i\>|\>\sqrt{V(s_{i})}\leq\mu_{\alpha}\\}}\sqrt{p(s_{i})}\sqrt{V(s_{i})}\ket{G}+\sum_{\\{i\>|\>\sqrt{V(s_{i})}\leq\mu_{\alpha}\\}}\sqrt{p(s_{i})}\sqrt{1-V(s_{i})}\ket{G_{\perp}}+\sum_{\\{i\>|\>\sqrt{V(s_{i})}>\mu\\}}\sqrt{p(s_{i})}\ket{G_{\perp}}.$
(27)

Applying QAE to estimate the probability of measuring $\ket{G}$ will give us
the expected value of the portfolio below the VaR level

$\mathbb{P}(\ket{G})=\sum_{\\{s\>|\>V(s)\leq\mu_{\alpha}^{2}\\}}p(s)V(s).$
(28)

Dividing this quantity by $\mathbb{P}[V\leq\mu_{\alpha}^{2}]$ will give us
$C_{\alpha}$ of Eq. (2) and the CVaR of the portfolio is then given by
$V_{0}-C_{\alpha}$, if the portfolio value today is $V_{0}$.

Because the threshold function of Eq. (26) is discontinuous at $x=\mu$, we
will need a high degree polynomial to generate a close approximation to the
function. We can however avoid the discontinuity at $x=\mu$ by inverting the
threshold function in Eq. (26) in the region below $\mu$ such that

$\hat{\theta}_{\mu}(x)=\begin{cases}\mu-x,\;{x}\leq{\mu}\\\
0,\;{x}>{\mu}.\end{cases}$ (29)

Applying this transformation instead at $\mu=\mu_{\alpha}$ gives us

$\mathbb{P}(\ket{G})=\sum_{\\{s\>|\>V(s)\leq\mu_{\alpha}^{2}\\}}p(s)\left(\mu_{\alpha}-V(s)\right),$
(30)

from which we compute

$\hat{C}_{\alpha}=\frac{1}{\mathbb{P}[V\leq\mu_{\alpha}^{2}]}\sum_{\\{s\>|\>V(s)\leq\mu_{\alpha}^{2}\\}}p(s)\left(\mu_{\alpha}-V(s)\right)=\left(\frac{1}{\mathbb{P}[V\leq\mu_{\alpha}^{2}]}\sum_{\\{s\>|\>V(s)\leq\mu_{\alpha}^{2}\\}}p(s)\mu_{\alpha}\right)-C_{\alpha}=\mu_{\alpha}-C_{\alpha}.$
(31)

The CVaR can then be calculated as
$\textrm{CVaR}_{\alpha}[V]=V_{0}-C_{\alpha}=V_{0}-(\mu_{\alpha}-\hat{C}_{\alpha})$.

###  IV.4 Threshold Function Approximation using QSP

In order to apply the threshold functions in Eq. (23) and Eq. (26) using QSP,
we first need to find polynomials which approximate them. The polynomial
transformations possible through the QSP framework must have definite parity
(even or odd), and because we aim to transform real positive amplitudes
representing derivative prices, we restrict our attention to the interval
$[0,1]$ due to symmetry. For the threshold function in in Eq. (23) used for
the VaR calculation, we follow the method in Ref. [29] and look for a real
polynomial $f(x)$ satisfying

$|f(x)-c|\leq\epsilon,\quad\forall
x\in[0,\mu-\Delta/2];\quad|f(x)|\leq\epsilon,\quad\forall
x\in[\mu+\Delta/2,1]$ (32)

where $c$ is chosen close to $1$ but preferably slightly smaller to avoid
overshooting, $\epsilon$ controls how far away the function is allowed to
deviate from the values $0$ and $1$, and $\Delta$ is a gap parameter which
controls the steepness of the jump from $1$ to $0$ at $\mu$ where
$\theta_{\mu}(x)$ is discontinuous. Smaller values of $\Delta$ will create a
better approximation to the function in the interval
$[\mu-\Delta/2,\mu+\Delta/2]$ and a smaller value of $\epsilon$ will
correspondingly allow for a better approximation in the interval
$[0,\mu-\Delta/2]\cup[\mu+\Delta/2,1]$.

While there are analytical methods of obtaining an approximation $f(x)$ to
this function with a polynomial of degree
$\mathcal{O}\left((1/\Delta)\log(1/\epsilon)\right)$ [30], we employ the
optimization-based method described in Ref. [29] 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}\\}$

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

and aims to find the coefficients $\\{c_{k}\\}$ which give the best
approximation to $\theta_{\mu}(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\max\left\\{\underset{x_{j}\in[0,\mu-\Delta/2]}{\max}|f(x_{j})-c|,\underset{x_{j}\in[\mu+\Delta/2,1]}{\max}|f(x_{j})|\right\\}$
$\displaystyle\textrm{s.t.}\quad
f(x_{j})=\sum_{k}c_{k}A_{jk},\quad|f(x_{j})|\leq c,\quad\textrm{for
}j=0,1,\ldots,M-1,$ (34)

where $A_{jk}$ is a matrix of coefficients defined as $A_{jk}=T_{2k}(x_{j})$
for $k=0,1,\ldots,d/2$.

We follow the same procedure to construct a polynomial approximation to the
inverted threshold function $\hat{\theta}_{\mu}$ of Eq. (29) used in the
estimation of CVaR. In this case, we solve the minimax optimization problem

$\displaystyle\underset{\\{c_{k}\\}}{\min}\quad\max\left\\{\underset{x_{j}\in[0,\mu-\Delta/2]}{\max}|f(x_{j})-(\mu-
x_{j})|,\underset{x_{j}\in[\mu+\Delta/2,1]}{\max}|f(x_{j})|\right\\},$ (35)

subject to the same constraints as Eq. (IV.4). The term $|f(x_{j})-(\mu-
x_{j})|$ now enforces that the polynomial approximation is close to $\mu-x$ in
the interval $[0,\mu-\Delta/2]$. Polynomial approximations generated using the
optimization-based method described in this section for the threshold
functions $\theta_{\mu}(x)$ and $\hat{\theta}_{\mu}(x)$ at $\mu=0.5$ are shown
in Fig. 4.

While the calculation of CVaR is just as (or more) important than VaR as a
business use case, in practice they are both formulated similarly in the QSP
framework. In fact, computing the CVaR using QSP requires first the
computation of VaR and the correct threshold to be identified. As such, in the
following sections we focus on the performance and estimation errors of the
VaR algorithm which is the core component in the estimation of both risk
metrics.

(a)

(b)

Figure 4: Polynomial approximations generated using the optimization-based
method described in Sec. IV.4. Using Quantum Signal Processing we can perform
the transformation $a\rightarrow P(a)$ for each singular value $a$ of a block-
encoded unitary, where increasing the degree of the polynomial $P(a)$ will
result in a better approximation to the target function. (a) Approximations to
the threshold function of Eq. (23) at $\mu=0.5$ used for the estimation of
VaR, generated by optimizing Eq. (IV.4) for different polynomial degrees $d$
and gap parameter $\Delta$, with $\epsilon=1-c=10^{-3}$. (b) Approximations to
the linear threshold function of Eq. (29) used for the estimation of CVaR,
generated using Eq. (35) for different polynomial degrees $d$,
$\Delta=10^{-3}$ and $\epsilon=1-c=10^{-3}$. Avoiding the discontinuity at
$a=\mu$ by using the function of Eq. (29) instead of Eq. (26) in the
calculation of CVaR results in significantly lower degree polynomials required
to approximate the function with high degree of accuracy.

##  V Complexity And Error Analysis

The VaR calculation methods discussed so far rely on estimating the
probability that a portfolio value $X$ will be less than some threshold $\mu$.
In the continuous setting, this can be written as

$P[X\leq\mu]=\int_{0}^{1}\theta_{\mu}(x)p(x)dx,$ (36)

where $p(x)$ denotes the probability of value $x$ occurring, $\theta_{\mu}(x)$
is the threshold function defined in Eq. (23) and the portfolio values have
been normalized to lie in $[0,1]$. The estimation of VaR is then the solution
to the inverse problem, where given a confidence level $\alpha$, we try to
determine the value $\mu_{\alpha}$ giving $P[X\leq\mu_{\alpha}]\geq 1-\alpha$,
which can be formulated as finding the root of

$\tilde{f}(\mu_{\alpha})=(1-\alpha)-P[X\leq\mu_{\alpha}].$ (37)

The error in the estimation of VaR will be determined by the accuracy with
which we can approximate Eq. (36) and the accuracy in finding the root of Eq.
(37), given the approximation to $P[X\leq\mu_{\alpha}]$. An error $\delta$ in
the estimate of $P[X\leq\mu_{\alpha}]$ will result in error $\epsilon$ in the
estimate of the VaR given by

$\epsilon=\delta\cdot\left.\frac{dF^{-1}(p)}{dp}\right\rvert_{p=1-\alpha}+\mathcal{O}\left(\delta^{2}\right),$
(38)

where $F^{-1}(p)$ is the inverse cumulative distribution function of the
distribution $p(x)$, assuming $F^{-1}(p)$ is analytic at $p=1-\alpha$. We now
examine the accuracy with which $P[X\leq\mu_{\alpha}]$ can be estimated by
various methods, which to first order gives the accuracy of the VaR estimate
(with a constant factor depending on the shape of the probability distribution
$p(x)$ at $\mu_{\alpha}$).

###  V.1 Classical and Semi-Classical Methods

Classical methods for the VaR estimation of financial derivatives use two
nested Monte Carlo invocations; an inner round to estimate the price of the
derivative $x$ (see Eq. (53)), and an outer round sampling scenarios to
estimate the integral of Eq. (36). When $N$ scenarios are used to approximate
the integral and the derivative price for each scenario is evaluated to
accuracy $\epsilon_{p}$, the approximation $P^{\prime}$ can be written as

$P^{\prime}[X\leq\mu]=\frac{1}{N}\sum_{i=1}^{N}\theta_{\mu}(x_{i}+\epsilon_{p}).$
(39)

Normalizing the prices $x_{i}$ to lie in the interval $[0,1]$, the error of
this approximation is

$|P^{\prime}-P|=\left|\frac{1}{N}\sum_{i=1}^{N}\theta_{\mu}(x_{i}+\epsilon_{p})-\int_{0}^{1}\theta_{\mu}(x)p(x)dx\right|\leq\epsilon_{\theta}+\epsilon_{S}$
(40) $\textrm{with
}\quad\epsilon_{\theta}\equiv\left|\frac{1}{N}\sum_{i=1}^{N}\theta_{\mu}(x_{i}+\epsilon_{p})-\frac{1}{N}\sum_{i=1}^{N}\theta_{\mu}(x_{i})\right|,\quad\epsilon_{S}\equiv\left|\frac{1}{N}\sum_{i=1}^{N}\theta_{\mu}(x_{i})-\int_{0}^{1}\theta_{\mu}(x)p(x)dx\right|,$
(41)

where $\epsilon_{\theta}$ is the error stemming from the evaluation of the
threshold function using an approximate derivative price, and $\epsilon_{S}$
the error from approximating the integral with Monte Carlo sampling. We can
derive the asymptotic dependence of $\epsilon_{\theta}$ on $\epsilon_{p}$ by
considering the expression for $\epsilon_{\theta}$ in Eq. (41) in the
continuous limit

$\displaystyle\epsilon_{\theta}$ $\displaystyle=$
$\displaystyle\left|\int_{0}^{1}\theta_{\mu}(x+\epsilon_{p})p(x)dx-\int_{0}^{1}\theta_{\mu}(x)p(x)dx\right|$
(42) $\displaystyle=$
$\displaystyle\left|\int_{0}^{\mu}p(x)dx+\int_{\mu}^{\mu+\epsilon_{p}}p(x)dx-\int_{0}^{\mu}p(x)dx\right|=\left|\int_{\mu}^{\mu+\epsilon_{p}}p(x)dx\right|\approx
p(\mu)\epsilon_{p}+O(\epsilon_{p}^{2}),$

where the last approximation follows by considering the left Riemann sum
approximation to the integral 333We can also consider the pricing error
contribution to the threshold function in Eq. (42) as
$\theta_{\mu}(x-\epsilon_{p})$ and the same asymptotic dependence can be
derived by considering the right Riemann sum approximation to the integral..
It then follows that $\epsilon_{\theta}=\mathcal{O}\left(\epsilon_{p}\right)$
with the constant factor depending on the shape of the probability
distribution $p(x)$ at $\mu$. Due to the two nested Monte Carlo evaluations,
the complexity of this method scales as
$\mathcal{O}\left(\epsilon_{S}^{-2}\epsilon_{p}^{-2}\right)$.

The _semi-classical_ approach proceeds similarly as in the purely classical
case, the only difference being that the derivative price for each scenario is
evaluated on a quantum computer using amplitude estimation. In this case, the
approximation error is the same, but the quadratic speedup of amplitude
estimation makes the approximation to $P[X\leq\mu]$ scale as
$\mathcal{O}\left(\epsilon_{S}^{-2}\epsilon_{p}^{-1}\right)$.

###  V.2 QAE VaR

With the QAE VaR method we create a superposition over scenarios
$\sum_{i}\sqrt{p(s_{i})}\ket{s_{i}}$ with the operator of Eq. (9), evaluate
the derivative price $x_{i}$ for each scenario $s_{i}$ in superposition using
amplitude estimation, and after applying a binary comparator at $\mu$, we read
out the approximate probability $P^{\prime}[X\leq\mu]$ using another round of
QAE (see Fig. 1). Contrary to the classical and semi-classical methods, in
this case we are not limited to classical sampling of scenarios if we have a
way to load the superposition of scenarios using a closed-form expression for
$p(s_{i})$. As such, the approximation $P^{\prime}[X\leq\mu]$ in this case can
be written as

$P^{\prime}[X\leq\mu]=\sum_{i=1}^{N}\theta_{\mu}(x_{i}+\epsilon_{p})p(x_{i})\delta
x+\epsilon_{{A}},$ (43)

where $\epsilon_{p}$ is the approximation error of using amplitude estimation
to approximate the derivative price $x_{i}$ and $\epsilon_{{A}}$ is the error
from the outer amplitude estimation. The the sum over $x_{i}$ can be a (multi-
dimensional) Riemann sum approximating an integral if $p(x)$ can directly
loaded into the superposition of scenarios. If the scenarios are instead
sampled from $p(x)$, the above expression becomes a Monte Carlo approximation
with $p(x_{i})\delta x=1/N$. In this case the error of the approximation can
be written as

$|P^{\prime}-P|=\left|\sum_{i=1}^{N}\theta_{\mu}(x_{i}+\epsilon_{p})p(x_{i})\delta
x-\int_{0}^{1}\theta_{\mu}(x)p(x)dx\right|+\epsilon_{{A}}\leq\epsilon_{\theta}+\epsilon_{S}+\epsilon_{{A}}$
(44)
$\textrm{where}\quad\epsilon_{\theta}\equiv\left|\sum_{i=1}^{N}\theta_{\mu}(x_{i}+\epsilon_{p})p(x_{i})\delta
x-\sum_{i=1}^{N}\theta_{\mu}(x_{i})p(x_{i})\delta
x\right|,\quad\epsilon_{S}\equiv\left|\sum_{i=1}^{N}\theta_{\mu}(x_{i})p(x_{i})\delta
x-\int_{0}^{1}\theta_{\mu}(x)p(x)dx\right|,$ (45)

and the asymptotic dependence of $\epsilon_{\theta}$ on $\epsilon_{p}$ follows
similarly to the classical case, with
$\epsilon_{\theta}=\mathcal{O}\left(\epsilon_{p}\right)$.

If the cost of loading the scenario superposition is $C(\epsilon_{S})$, the
total complexity of estimating $P[X\leq\mu]$ using QAE VaR is
$\mathcal{O}\left((C(\epsilon_{S})+\epsilon_{p}^{-1})\epsilon_{A}^{-1}\right)$.
When the individual scenarios in the superposition are selected by sampling an
appropriate probability distribution similarly to the classical and semi-
classical approach, a sampling error $\epsilon_{S}$ requires
$\epsilon_{S}=\mathcal{O}\left(1/\sqrt{N}\right)$ scenarios. As this
superposition can be loaded in $\mathcal{O}\left(\log N\right)$ depth444It is
important to highlight that while this superposition can generally be loaded
in $\mathcal{O}\left(\log N\right)$ depth, in practice, this complexity might
include very large constant factors, as discussed in [31]. [31], the total
complexity becomes
$\mathcal{O}\left((\log(1/\epsilon_{S}^{2})+\epsilon_{p}^{-1})\epsilon_{A}^{-1}\right)=\tilde{\mathcal{O}}(\epsilon_{p}^{-1}\epsilon_{A}^{-1})$,
where the tilde notation ignores poly-logarithmic terms.

###  V.3 QSP VaR

The QSP VaR algorithm follows similar steps as the QAE variant, but in this
case the derivative prices are computed with the accuracy allowed of the
$\mathcal{A}$ operator from Eq. (6) and do not incur an approximation error
from amplitude estimation. While there is still approximation error in the
implementation of $\mathcal{A}$, it can be made exponentially small by
increasing the number of qubits in the construction of the operator with only
logarithmically increasing the oracle complexity of the algorithm [6] and as
such we ignore it in this error analysis. On the other hand, because in this
case we use QSP to approximate the threshold function of Eq. (23), we do not
implement exactly the threshold function $\theta_{\mu}(x)$ but rather an
approximation, which we denote as $\theta_{\mu}^{\prime}(x)$. As in the case
of QAE VaR, the preparation of the scenarios is not limited to sampling if
there is a method to load the scenario superposition efficiently, so we use
the same notation as the previous section to write the approximation to
$P[X\leq\mu]$ in this case as

$P^{\prime}[X\leq\mu]=\sum_{i=1}^{N}\theta_{\mu}^{\prime}(x_{i})p(x_{i})\delta
x+\epsilon_{{A}}.$ (46)

The approximation error is then given by

$|P^{\prime}-P|=\left|\sum_{i=1}^{N}\theta_{\mu}^{\prime}(x_{i})p(x_{i})\delta
x-\int_{0}^{1}\theta_{\mu}(x)p(x)dx\right|+\epsilon_{{A}}\leq\epsilon_{\theta}+\epsilon_{S}+\epsilon_{{A}}$
(47)
$\textrm{where}\quad\epsilon_{\theta}\equiv\left|\sum_{i=1}^{N}\theta_{\mu}^{\prime}(x_{i})p(x_{i})\delta
x-\sum_{i=1}^{N}\theta_{\mu}(x_{i})p(x_{i})\delta
x\right|,\quad\epsilon_{S}\equiv\left|\sum_{i=1}^{N}\theta_{\mu}(x_{i})p(x_{i})\delta
x-\int_{0}^{1}\theta_{\mu}(x)p(x)dx\right|.$ (48)

Contrary to the classical, semi-classical and QAE cases, QSP allows us to
bound the error from the first term $\epsilon_{\theta}$. For a choice of the
gap parameter $\Delta$, picking a polynomial degree of
$d=\mathcal{O}\left(1/\Delta\right)$ will guarantee that
$\epsilon_{\theta}\leq\mathcal{O}\left(\Delta\right)\sim\mathcal{O}\left(1/d\right)$
[30, 29]. In Fig. 5 we show numerically how this approximation error scales
with the degree of the polynomial used to approximate
$\theta_{\mu}^{\prime}(x)$ for $\Delta=10^{-3}$ when the polynomial is
generated using the method described in Sec. IV.4. The complexity of
estimating $P[X\leq\mu]$ then scales as
$\mathcal{O}\left((C(\epsilon_{S})+\epsilon_{\theta}^{-1})\epsilon_{A}^{-1}\right)$,
where again $C(\epsilon_{S})$ denotes the cost of loading the superposition
over scenarios. In the case where the scenarios are selected by sampling a
distribution and loaded in superposition in $\mathcal{O}(\log(N))$ depth, we
get $\tilde{\mathcal{O}}\left(\epsilon_{\theta}^{-1}\epsilon_{A}^{-1}\right)$.

Figure 5: The approximation error $\epsilon_{\theta}$ in the estimation of
$P[X\leq\mu]$ due to the polynomial approximation used in the QSP VaR method
as shown in Eq. (48). We pick the gap parameter $\Delta=10^{-3}$ and
numerically compute the approximation error of $d$-degree polynomials
generated using the optimization-based method described in Sec. IV.4. This
error is expected to scale as $\mathcal{O}\left(1/d\right)$ and fitting a
function of the form $f(d)=a/d+b$ to the data in this case gives us $a\approx
0.31$.

From the analysis in this section, we see that for a target error $\epsilon$
in the estimation of the probability in Eq. (36), classical nested Monte Carlo
scales as $\mathcal{O}\left(\epsilon^{-4}\right)$ and the semi-classical
method scales as $\mathcal{O}\left(\epsilon^{-3}\right)$, where the advantage
stems from the quadratic advantage of QAE compared to classical sampling in
the inner (pricing) Monte Carlo step. The semi-classical method cannot provide
further advantage in the outer Monte Carlo sampling step, as both of these
methods are limited to sampling scenarios from an appropriate probability
distribution. When also restricted to scenario sampling, the QAE and QSP-based
quantum methods both scale as $\tilde{\mathcal{O}}(\epsilon^{-2})$. If the
relevant probability distribution of scenarios can be loaded efficiently in
superposition, we do not have to load individually sampled scenarios in
superposition, and the error $\epsilon_{S}$ in Eq. (45) and Eq. (48) can be
made exponentially small, which is not possible with the classical and semi-
classical approaches.

##  VI Comparison of the Quantum Methods for Estimating VaR

The QAE and QSP VaR quantum estimation methods described in the previous
sections are based on the same procedure and it is not obvious whether one
outperforms the other in terms of required quantum resources for a target
estimation accuracy. The only difference between the two methods is how the
probability that a portfolio value under a number of scenarios is smaller than
a VaR candidate $\mu$ is encoded in an amplitude so that we can perform QAE to
read it out. The QAE method performs an
analog$\rightarrow$digital$\rightarrow$analog transformation whereby the
portfolio values initially encoded as amplitudes, are digitized approximately
into a quantum register by the inner round of QAE before being transformed
appropriately back into amplitudes using a binary comparator circuit. On the
other hand, the QSP framework allows us to perform the necessary
transformation to the amplitudes without the need for an intermediate binary
representation which removes that source of error in the process. It does
however introduce an error from the approximate implementation of the
threshold function which encodes the value to be extracted with QAE.
Additionally, looking at the quantum circuits for each method in Fig. 1 and
Fig. 3, the QSP circuit is simpler in that a) it only requires one additional
qubit other than those required to implement the $\mathcal{S}$ and
$\mathcal{A}$ operators which load the scenario superposition and encode the
derivative price into an amplitude respectively, b) the $\mathcal{A}$ operator
is only controlled on the scenario qubits while it is controlled on one more
qubit in the QAE case by virtue of the controlled invocation $\mathcal{Q}$
operator, and c) while the QSP method contains single-qubit phase rotations,
the $\mathcal{Q}$ operators in the QAE method require reflections.

The total resources required for both methods can be summarized as

$C_{\textrm{VaR{}}}=C_{S}C_{\mu}C_{\textrm{QAE}}\cdot k,$ (49)

where $C_{S}$ is the cost of implementing the scenario superposition unitary
$\mathcal{S}$, $C_{\mu}$ denotes the resources required by either method to
generate a state where the probability of measuring a qubit in the $\ket{0}$
state gives the probability that the derivative portfolio value falls below a
value $\mu$, $C_{\textrm{QAE}}$ is the cost of QAE to read out that
probability and $k$ is the number of bisection search iterations required to
compute the VaR estimate to a desired accuracy. The cost component $C_{\mu}$
is the main difference between the two methods in terms of the overall
algorithmic cost so we focus on estimating the cost of this component in each
case. Specifically, we measure the cost by the number of oracle invocations to
the pricing operator $\mathcal{A}$ (and its inverse $\mathcal{A}^{\dagger}$)
which is the most expensive component in derivative pricing [6]. The QAE
method requires one invocation of $\mathcal{A}$ to set up the initial state
Eq. (10) and $m$ ancilla qubits to discretize the value of the derivative
portfolio under the scenarios. Each ancilla qubit $j\in[0,m-1]$ controls
$2^{j}$ invocations of the $\mathcal{Q}$ operator as shown in Fig. 1, and
$\mathcal{Q}$ is comprised of one instance of $\mathcal{A}$ and one of
$\mathcal{A}^{\dagger}$, giving a total of $N_{\mu}^{\textrm{QAE}}(m)=2^{m}+1$
oracle calls. The QSP method uses a $d$-degree polynomial to approximate the
threshold function $\theta_{\mu}(x)$ of Eq. (23) and the its implementation
requires $N_{\mu}^{\textrm{QSP}}(d)=d-1$ invocations of $\mathcal{A}$ as shown
in Eq. (21) and Fig. 3.

We estimate the values of $N_{\mu}^{\textrm{QAE}}$ and
$N_{\mu}^{\textrm{QSP}}$ numerically by sampling simulated scenario prices
from a normal distribution with mean $0.5$ and standard deviation $0.09$ and
look for the VaR at the $\alpha=99\%$ level. The mean and standard deviation
of the normal distribution were chosen so that the sampled values lie with
high probability in the interval $[0,1]$ and the $\alpha=99\%$ confidence
level is typical for VaR estimations in practice [32]. We simulate the QAE VaR
estimation algorithm as described in Sec. III, where for each sampled scenario
price we generate the superposition of binary strings returned by QAE as
possible approximations to the value, weighed by the corresponding
probabilities of Eq. (14). Starting with $\mu=0.5$, we calculate the aggregate
probability $P_{0}$ across all scenarios that the price is smaller than $\mu$.
Then, for a parameter $\epsilon_{A}$ we determine bounds
$[p_{l},p_{h}]=[P_{0}-\epsilon_{A},P_{0}+\epsilon_{A}]$ as the confidence
interval given by an amplitude estimation round to determine $P_{0}$ with
target accuracy $\epsilon_{A}$. We run a bisection search over $\mu$ until we
find that $1-\alpha=0.01\in[p_{l},p_{h}]$ and we return that value of $\mu$ as
the VaR estimate.

The QSP method is simulated similarly, the only difference being how we
calculate $P_{0}$ for a given value of $\mu$. In this case, we pick a value
for the polynomial degree $d$ and for each value of $\mu$ during bisection
search we generate a $d$-degree polynomial approximation to the threshold
function $\theta_{\mu}(x)$ using the method described in Sec. IV.4 and apply
the resulting function to the sampled scenario prices. Subsequently, we
generate the same confidence interval for the amplitude estimation round given
a target accuracy $\epsilon_{A}$ as in the QAE method, and impose the same
stopping condition. The equivalent classical estimate is computed by sorting
the sampled scenario prices and finding the value at the $1-\alpha$
percentile. Because the simulations of the quantum methods use sampled
scenarios, the classical estimate represents the lowest bound either quantum
method can achieve.

The results from the simulation for both QAE and QSP methods for different
values of oracle call invocations and the final amplitude estimation accuracy
$\epsilon_{A}$ are shown in Fig. 6. Each data point is comprised of $100$
independent samples, where for each sample we randomly generate $N=5k$
scenario prices from the normal probability distribution
$\mathcal{N}(\mu=0.5,\sigma^{2}=0.09^{2})$, simulate the QAE and QSP
algorithms and then compute the error of each sample as the absolute value of
the difference between the quantum VaR estimate and the classical VaR estimate
from the same generated scenario prices. In the figure we plot the average
error across all samples for each data point. One initial observation is that
the QSP method performs consistently better for the error range examined,
requiring on average 10x fewer oracle calls than QAE for the same target
error. In addition, we notice that as the number of oracle calls increase, the
probability $P^{\prime}$ we are estimating each time using amplitude
estimation is encoded with increasing accuracy, but clearly at some point this
increased encoding accuracy will fail to reduce the VaR estimation error if
the value of $\epsilon_{A}$ is not small enough to resolve the difference.
This is what we observe at the tail end of the plots in Fig. 6, where the VaR
error in the $\epsilon_{A}=5\times 10^{-4}$ plots starts to decrease more
slowly than the $\epsilon_{A}=10^{-4}$ plots. Until that point, we calculate
that for a fixed value of $\epsilon_{A}$, the VaR error decreases
approximately linearly as the number of oracle calls increases for both QAE
and QSP cases (the slope of a log-log fit is approximately $-0.9$ for all four
plots), as expected from the error analysis in Sec. V.2 and Sec. V.3
respectively where both errors scale as $\mathcal{O}(\epsilon^{-1}_{p})$ with
everything else kept constant.

Our choice of parameters $N$, $\mu$ and $\sigma$ was made to simplify the
numerical simulations, but we observe qualitatively similar results for
different choices of distribution parameters $\mu$ and $\sigma$ as well as the
number of scenario prices generated for each data point.

Figure 6: VaR estimation errors as a function of the number of oracle calls
for the QSP and QAE quantum methods. For each data point we generate $100$
samples, and for each sample we randomly generate $N=5k$ scenario prices from
the normal probability distribution $\mathcal{N}(\mu=0.5,\sigma^{2}=0.09^{2})$
and compute the VaR estimate by simulating the QAE and QSP algorithms as
described in Sec. VI. The error of each sample is calculated as the absolute
value of the difference between the quantum VaR estimate and the classical VaR
estimate from the same generated scenario prices. We simulate each algorithm
for two values of the final amplitude estimation accuracy $\epsilon_{A}$ and
plot the average error across all samples for that each data point. On
average, we observe that the QSP method requires 10x fewer oracle calls than
QAE for the same VaR error. While increasing number of oracle calls $N_{o}$
leads to the VaR error decreasing as $\mathcal{O}(N_{o}^{-1})$, notice that
this scaling starts to depart from linearity as the amplitude estimation
accuracy $\epsilon_{A}$ gets close to the accuracy $1/N_{o}$ with which the
amplitude is encoded. Eventually, the scaling plateaus as the increased
encoding accuracy becomes too fine to be resolved by an amplitude estimation
process of fixed accuracy.

##  VII Conditions for Quantum Advantage

On account of the numerical results in the previous section indicating that
the QSP-based VaR method outperforms the corresponding QAE variant, in this
section we turn our attention on how the VaR estimation with QSP fares against
classical methods and what conditions (if any) there are for potential quantum
advantage. As discussed in Sec. V, the estimation error of the QSP method is
lower-bounded by the classical sampling error of
$\mathcal{O}\left(\epsilon^{-2}\right)$ when the superposition of scenarios is
generated with samples from a probability distribution. Assuming no particular
structure in the representation of the sampled scenarios, the superposition of
$N$ scenarios can be loaded to accuracy $\epsilon$ using the QRAM-based state
preparation routine from Ref. [31] with T-depth
$\mathcal{O}\left(\log(N)+\log(1/\epsilon)\right)$, at a cost of
$\mathcal{O}\left(N\right)$ T-count. On the other hand, if the scenarios can
be prepared by discretizing a multivariate probability distribution of the
relevant risk factors and loaded efficiently in superposition, the lower bound
on the VaR estimation error can be made exponentially small by increasing the
number of qubits representing the distribution. In the case that the
probability distribution can be modeled explicitly with an analytic density
function, the scenario superposition could be prepared using the re-
parameterization method from Ref. [6]. Alternatively, if the scenario
probabilities have to be inferred from historical data, quantum
representations of financially relevant copulas [33, 34] could potentially be
employed if they can scale to scenarios of high dimensionality. In either
case, the resources required to generate the scenario superposition will be
highly dependent on the nature of the target probability distribution if
explicit loading is possible, and the number of risk factors that make up each
scenario as well as the required granularity for each risk factor (as shown in
Eq. (7)). As such, in this section we do not examine specific choices of
scenario preparation methods, and instead consider the cost of loading the
superposition to be a free variable which needs to be taken into account when
estimating the total resources for the algorithm.

As a benchmark to examine the conditions for practical quantum advantage, we
consider the autocallable derivative product studied in Ref. [6]. We calculate
the resources required by the QSP-based VaR algorithm to estimate the VaR of
the product to the same accuracy as a classical Monte Carlo method, assuming
that the product is evaluated under scenarios such that the resulting prices
are distributed normally. We use the latest available resource estimates from
Ref. [19], whereby an autocallable derivative product can be priced to within
$\epsilon_{p}=2\times 10^{-3}$ in one second classically and the $\mathcal{A}$
operator required to price the contract on a quantum computer can be
constructed with a T-depth of $T_{A}=3900$.

The pricing of this product under $N$ scenarios is simulated classically by
randomly generating $N$ values from a normal distribution
$\mathcal{N}(\mu,\sigma^{2})$ representing $N$ scenario prices, and adding a
noise term sampled from $\mathcal{N}(0,\epsilon_{p}^{2})$. Then we compute
$V_{\alpha}$ from Eq. (1) by sorting the resulting simulated prices from
smallest to largest, and picking the value at the $1-\alpha$ percentile. The
error of the VaR estimate from the classical simulation can be calculated as
$|V_{\alpha}-\Phi^{-1}(1-\alpha)|$, where $\Phi(x)$ is the cumulative
distribution function of the distribution $\mathcal{N}(\mu,\sigma^{2})$. We
repeat this process $200$ times, and define the error of the classical
algorithm $\epsilon_{C}$ as the standard deviation of the resulting
distribution of errors. Since we assume the contract can be priced classically
in one second, for $N$ scenarios, we consider the runtime of the end-to-end
classical VaR estimation to be $N$ seconds.

To evaluate the performance of the quantum QSP VaR algorithm, we similarly
generate $N$ values from the normal distribution $\mathcal{N}(\mu,\sigma^{2})$
representing $N$ scenario prices and simulate the QSP algorithm as described
in Algorithm 1. The error of the quantum estimate is given by
$|\mu_{\alpha}^{2}-\Phi^{-1}(1-\alpha)|$, where $\mu_{\alpha}$ is the estimate
of $\sqrt{V_{\alpha}}$ returned by the QSP algorithm. Repeating the process
$200$ times, we define the error of the quantum algorithm $\epsilon_{Q}$ as
the value at the $68$th percentile of the resulting distribution of errors.
This way we can say that with probability $68\%$ the estimation error of the
algorithm does not exceed $\epsilon_{Q}$, which corresponds to the same
confidence in the estimate of the classical algorithm, defined as the error at
one standard deviation.

We then proceed to estimate the T-depth of the end-to-end process, as a proxy
for the total runtime. The T-depth of each iteration, which is the T-depth of
the circuit of Fig. 3 is given by

$T_{i}=T_{S}+dT_{A}+dT_{R},$ (50)

where $T_{S}$ is the T-depth of the scenario preparation unitary
$\mathcal{S}$, $T_{A}=3900$ is the T-depth of the $\mathcal{A}$ operator [19],
$d$ is the polynomial degree used for the approximation to the threshold
function and $T_{R}$ is the T-depth of implementing each controlled $R_{z}$
rotation. Using the method in [35], the $R_{z}$ rotation can be performed to
precision $\epsilon_{R}$ with a T-depth of approximately
$3\log_{2}(1/\epsilon_{R})$ and the decomposition in [36] allows us to perform
the controlled rotation with an $R_{z}$ depth of one, using one ancilla qubit.
We measure the cost of each round of QAE for target accuracy $\epsilon_{A}>0$
and confidence level $1-\alpha_{k},\alpha_{k}\in(0,1)$ by the bound derived in
[22]

$N_{\text{oracle}}^{\text{wc}}\leq\frac{1.4}{\epsilon_{A}}\log\left(\frac{2}{\alpha_{k}}\log_{2}\left(\frac{\pi}{4\epsilon_{A}}\right)\right),$
(51)

where $N_{\text{oracle}}^{\text{wc}}$ denotes the worst-case number of QAE
oracle calls. In our case, the QAE oracle call consists of one invocation of
the entire circuit in Fig. 3 and one to its inverse (where we ignore the cost
of reflections in the QAE oracle). Thus, the end-to-end T-depth of the QSP VaR
algorithm which goes through $k$ rounds of QAE is given by

$T_{\textrm{QSP}}=\frac{2.8k}{\epsilon_{A}}\log\left(\frac{2}{\alpha_{k}}\log_{2}\left(\frac{\pi}{4\epsilon_{A}}\right)\right)\left(T_{S}+dT_{A}+3d\log_{2}(1/\epsilon_{R})\right).$
(52)

Additionally, we choose $\epsilon_{R}$ such that the total error from all the
$R_{z}$ rotations is the same order of magnitude as the error from the
polynomial approximation to the threshold function, which is $\gtrsim 10^{-4}$
as shown in Fig. 5. In our numerical simulations the polynomial degrees we use
satisfy $d\leq 1000$, so we pick $\epsilon_{R}=10^{-7}$.

The remaining free parameters in the total cost of the QSP algorithm are the
scenario loading cost $T_{S}$, the polynomial degree $d$ used for QSP and the
target accuracy $\epsilon_{A}$ with which we estimate the encoded probability
in each round of QAE. Additionally, for a fixed number of scenarios $N$, the
VaR estimation error is only dependent on the choices of $d$ and
$\epsilon_{A}$. In order to benchmark the performance of the QSP algorithm,
for a choice of $N$, we numerically look for the optimal values of parameters
$(d,\epsilon_{A})$ such that the estimation error of the quantum algorithm
$\epsilon_{Q}$ is equal to that of the classical algorithm $\epsilon_{C}$ and
calculate the overall T-depth given by Eq. (52) with $T_{S}$ left as the only
remaining free cost parameter. Once we estimate the optimal values of
$(d,\epsilon_{A})$, we can additionally define the minimum rate at which that
a quantum processor would need to execute T-gates so that the physical runtime
of the QSP method would equal the physical runtime of a classical processor.
That _logical_ QPU clock rate is defined as $L\equiv T_{\textrm{QSP}}/N$,
where, as described above, we assume it takes $N$ seconds to classically price
$N$ scenarios.

We use $\mu=0.5$ and $\sigma=0.09$ as the parameters of the normal
distribution from which we draw scenario samples, and simulate the classical
VaR estimation method for $\alpha=99\%$. The QSP algorithm is simulated for
the same value of $\sigma$, but in this case we average the results over
$\mu\in[0.45,0.48,0.5,0.52,0.55]$. The reason for this is that the QSP method
can do better (or worse) depending on where the true VaR value lies, because
of the bisection search performed. Specifically, the estimation error of the
quantum method will depend on the distance between the true VaR value and the
values reachable by bisection search in the interval $[0,1]$. We therefore use
this averaging over values of $\mu$ to benchmark the average performance of
the algorithm. On the other hand, the estimation error has no dependence on
the true value of VaR classically, so we fix the value of $\mu$ in that case.
We chose $\sigma=0.09$ for the simulated scenario distribution such that the
sampled scenarios are in the interval $[0,1]$ with high probability. We have
also simulated the algorithms for $\sigma\in[0.05,0.07]$ and our results
remain qualitatively unchanged.

In Fig. 7 we show the logical clock rate required for quantum advantage as
estimated from the numerical simulations, as function of the scenario
preparation T-depth $T_{S}$, for different numbers of scenarios $N$. For
reference, we also show the logical clock rate (45Mhz) required for quantum
advantage in pricing a single derivative calculated in Ref. [19]. Because the
runtime of the classical and semi-classical methods scales linearly with $N$,
as $N$ increases, the QSP method starts to provide room for quantum advantage
as long as the scenario preparation can be prepared with a smaller T-depth
than a certain threshold, which can be derived from Eq. (52). For small
preparation T-depths, the value of $T_{S}$ is dominated by the resources
required for QSP, reducing the logical clock rate for advantage until the two
terms become comparable. For $N=50$k, we notice that there is an advantage
over single derivative pricing as long as the scenario superposition can be
prepared with a T-depth $T_{S}\lesssim 8\times 10^{7}$. In the case that the
superposition can be prepared with a T-depth $T_{S}\lesssim 3\times 10^{5}$,
the logical clock rate required for advantage is reduced by a factor of $\sim
30$x.

We note here that these estimates for quantum advantage are based on the
assumption that the superposition over scenarios is created from sampled
scenarios. In the case where an explicit probability distribution of scenarios
can be loaded efficiently as a discretized probability density function, the
error from the scenario sampling would decrease exponentially and allow the
possibility of further quantum advantage.

Figure 7: Minimum QPU logical clock rate required for the QSP VaR method to
match classical methods in total execution time as a function of the scenario
loading T-depth. For reference, the horizontal dashed line indicates the
minimum logical clock rate required for quantum advantage in the pricing of a
single derivative as estimated in Ref. [19]. Both classical and quantum
algorithms are simulated by sampling scenarios from a normal distribution as
described in Sec. VII, and for each choice of number of scenarios $N$, we
search for the parameters $(d,\epsilon_{A})$ (controlling the polynomial
approximation to the threshold function and the accuracy of QAE respectively)
such that the estimation error from the QSP algorithm matches that of the
classical algorithm. For $N=50$k, we see that there is a possibility for
reduction of the logical clock rate required for quantum advantage for up to a
factor $\sim 30$x compared to that required in the pricing of a single
derivative. For all choices of $N$ shown here, we numerically estimate that
the QSP method requires at least $d=600$ and $\epsilon_{A}=1.2\times 10^{-3}$
in order to match the corresponding estimation error of the classical method.

##  VIII Discussion

While the efficacy of quantum computers in derivative pricing has been studied
in the context of estimating the value of individual derivatives, in various
relevant contexts, the quantity of interest is not the individual derivative
price itself, but rather metrics which characterize the risk profile of
financial derivatives. Because calculating these metrics often requires
pricing derivatives repeatedly, the power of quantum computers could be
harnessed more effectively by considering the aggregate computation instead of
the application of quantum computing to the individual components. As such,
devising quantum algorithms which use the oracles constructed for the pricing
of derivatives as subcomponents, opens up a new avenue for quantum advantage
in the context of derivative pricing [7]. In this work, we consider the
problem of estimating the Value at Risk and the Conditional Value at Risk of
financial derivatives. We first extend the algorithm introduced in Ref. [3,
10] to estimate the VaR of financial derivatives and then introduce a quantum
algorithm based on the Quantum Signal Processing framework to estimate these
metrics. Even though we show that asymptotically these two algorithms have
similar performance, we show that the use of QSP in practice allows for a more
efficient algorithm, a conclusion which has also been observed in different
contexts [16, 17, 18]. This result also highlights the practical advantage of
QSP-based methods which not only allow the application of generic
transformations to appropriately encoded values of interest, but that they do
so with minimal overhead. While our formulation of derivative pricing in the
QSP framework was used in this case to calculate the VaR and CVaR risk
metrics, we expect that it will open up new research avenues in the larger
context of derivative pricing using quantum computers.

Additionally, we study the possibility of quantum advantage by using this
algorithm over equivalent classical methods. The performance of the algorithm
relies on the ability to prepare a superposition of scenarios efficiently and
different approaches to that end have been proposed [6, 33, 34]. During the
writeup of this manuscript, a state preparation method was introduced to load
a uniform superposition of subsets of computational basis states scaling as
$\mathcal{O}\left(\log_{2}M\right)$ for $M$ states in superposition, with low
overhead [37]. The superposition over scenarios encoded using Eq. (7) is
precisely a question of generating a uniform superposition over subsets of
computational basis states, and therefore determining the applicability of
this method to the VaR algorithms we describe here is a very promising
research direction. In this work, we have left the scenario preparation as an
open question and instead, through numerical simulations, we estimate upper
bounds on the resources required to prepare the necessary superposition in
order for quantum advantage to exist. We find that in certain settings, the
quantum algorithm can provide a reduction in the logical clock rate required
for quantum advantage estimated in Ref. [19] by $\sim 30$x. While this
estimate is based on pricing the derivative under scenarios sampled from a
relevant probability distribution, quantum computing allows the possibility of
using the full (discretized) multivariate probability distribution in the
calculation, an option not available to classical methods due to the curse of
dimensionality. An extension of this work would be to study which (relevant)
probability distributions can be loaded efficiently on a quantum computer and
the impact on the estimate of possible quantum advantage calculated in this
work.

In our comparison between the performance of quantum and classical methods, we
have assumed that classical methods for the VaR estimation of financial
derivatives require pricing derivative contracts across a fixed number of
scenarios $N$, such that the total complexity scales as
$\mathcal{O}\left(N/\epsilon^{2}\right)$. In Ref. [38], the authors argue that
for sufficiently smooth probability density functions, an adaptive sampling
technique combined with Multi-level Monte Carlo can reduce the complexity of
general VaR and CVaR estimation to $\tilde{\mathcal{O}}(\epsilon^{-2})$. One
core element in that approach is that the probability distribution behind the
scenarios used in the VaR estimation can be easily sampled, such that the
number of samples generated can be adjusted to fit the target accuracy of the
calculation. However, in practice, it is extremely difficult to generate such
probability distribution across (typically) thousands of market factors, and
most commonly a fixed number of scenarios is created based on historical data
that capture market features deemed important from a financial modeler’s
perspective. That said, Multi-level Monte Carlo can also be applied in a
quantum setting [39], and adaptive samples can similarly be incorporated in
the quantum process, by adjusting the accuracy with which the derivative price
is encoded as a quantum amplitude. We leave the detailed analysis of this
approach to a future study.

###### Acknowledgements.

We thank Bryce Fuller, Patrick Rall and Farrokh Labib for discussions
regarding Quantum Signal Processing and feedback on this manuscript, Kunal
Kishore for his technical and business insights regarding the VaR estimation
of financial derivatives, and Noelle Ibrahim for discussions on classical
adaptive methods for estimating VaR.

## References

  * 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).
  * 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).
  * 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).
  * Markowitz [1952] H. Markowitz, “ _Portfolio Selection_ ,” The Journal of Finance 7, 77 (1952).
  * Roy [1952] A. D. Roy, “ _Safety First and the Holding of Assets_ ,” Econometrica 20, 431 (1952).
  * Egger  _et al._ [2020] D. J. Egger, R. G. Gutierrez, J. C. Mestre, and S. Woerner, “ _Credit risk analysis using quantum computers_ ,” IEEE Transactions on Computers (2020).
  * Low and Chuang [2017a] G. H. Low and I. L. Chuang, “ _Optimal Hamiltonian Simulation by Quantum Signal Processing_ ,” Phys. Rev. Lett. 118, 010501 (2017a).
  * 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.
  * Low and Chuang [2019] G. H. Low and I. L. Chuang, “ _Hamiltonian simulation by qubitization_ ,” Quantum 3, 163 (2019).
  * Kikuchi  _et al._ [2023] Y. Kikuchi, C. Mc Keever, L. Coopmans, M. Lubasch, and M. Benedetti, “ _Realization of quantum signal processing on a noisy quantum computer_ ,” npj Quantum Information 9 (2023).
  * Martyn  _et al._ [2021a] J. M. Martyn, Z. M. Rossi, A. K. Tan, and I. L. Chuang, “ _Grand Unification of Quantum Algorithms_ ,” PRX Quantum 2 (2021a).
  * Martyn  _et al._ [2021b] J. M. Martyn, Y. Liu, Z. E. Chin, and I. L. Chuang, “ _Efficient Fully-Coherent Hamiltonian Simulation_ ,” arXiv preprint arXiv:2110.11327 (2021b), arXiv:2110.11327 [quant-ph] .
  * 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).
  * Stamatopoulos and Zeng [2023] N. Stamatopoulos and W. J. Zeng, “ _Derivative Pricing using Quantum Signal Processing_ ,” (2023), arXiv:2307.14310 [quant-ph] .
  * Gilyén  _et al._ [2019] A. Gilyén, S. Arunachalam, and N. Wiebe, “ _Optimizing quantum optimization algorithms via faster quantum gradient computation_ ,” Proceedings of the Thirtieth Annual ACM-SIAM Symposium on Discrete Algorithms , 1425–1444 (2019).
  * Brassard  _et al._ [2002] G. Brassard, P. Hoyer, M. Mosca, and A. Tapp, “ _Quantum Amplitude Amplification and Estimation_ ,” Contemporary Mathematics 305 (2002).
  * Grinko  _et al._ [2021] D. Grinko, J. Gacon, C. Zoufal, and S. Woerner, “ _Iterative quantum amplitude estimation_ ,” npj Quantum Information 7 (2021).
  * 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).
  * Giurgica-Tiron  _et al._ [2022] T. Giurgica-Tiron, I. Kerenidis, F. Labib, A. Prakash, and W. Zeng, “ _Low depth algorithms for quantum amplitude estimation_ ,” Quantum 6, 745 (2022).
  * Rall [2021] P. Rall, “ _Faster Coherent Quantum Algorithms for Phase, Energy, and Amplitude Estimation_ ,” Quantum 5, 566 (2021).
  * 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).
  * 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).
  * Low and Chuang [2017b] G. H. Low and I. L. Chuang, “ _Hamiltonian simulation by uniform spectral amplification_ ,” arXiv preprint arXiv:1707.05391 (2017b).
  * Clader  _et al._ [2022] B. D. Clader, A. M. Dalzell, N. Stamatopoulos, G. Salton, M. Berta, and W. J. Zeng, “ _Quantum Resources Required to Block-Encode a Matrix of Classical Data_ ,” IEEE Transactions on Quantum Engineering 3, 1–23 (2022).
  * O’Brien and Szerszeń [2017] J. O’Brien and P. J. Szerszeń, “ _An evaluation of bank measures for market risk before, during and after the financial crisis_ ,” Journal of Banking & Finance 80, 215 (2017).
  * Milek [2020] J. Milek, “ _Quantum Implementation of Risk Analysis-relevant Copulas_ ,” arXiv preprint arXiv:2002.07389 (2020).
  * Zhu  _et al._ [2022] D. Zhu, W. Shen, A. Giani, S. R. Majumder, B. Neculaes, and S. Johri, “ _Copula-based Risk Aggregation with Trapped Ion Quantum Computers_ ,” arXiv preprint arXiv:2206.11937 (2022).
  * 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).
  * Kim and Choi [2018] T. Kim and B. Choi, “ _Efficient decomposition methods for controlled-R n using a single ancillary qubit_ ,” Scientific Reports 8 (2018).
  * Shukla and Vedula [2024] A. Shukla and P. Vedula, “ _An efficient quantum algorithm for preparation of uniform quantum superposition states_ ,” Quantum Information Processing 23 (2024).
  * Giles and Haji-Ali [2019] M. B. Giles and A.-L. Haji-Ali, “ _Multilevel Nested Simulation for Efficient Risk Estimation_ ,” SIAM/ASA Journal on Uncertainty Quantification 7, 497 (2019).
  * An  _et al._ [2021] D. An, N. Linden, J.-P. Liu, A. Montanaro, C. Shao, and J. Wang, “ _Quantum-accelerated multilevel Monte Carlo methods for stochastic differential equations in mathematical finance_ ,” Quantum 5, 481 (2021).

##  Appendix A Derivative Pricing using Amplitude Estimation

The price of a derivative contract is calculated as the expectation value of
the derivative’s discounted payoff, evaluated over a suitable stochastic
process that is assumed to govern the dynamics of the underlying market
parameters. Given stochastic paths $\omega\in\Omega$, each of which occuring
with probability $p(\omega)$, the expectation value of the discounted payoff
$f(\omega)$ can be written as

$\mathbb{E}[f]=\sum_{\omega\in\Omega}p(\omega)f(\omega).$ (53)

Quantum Amplitude Estimation (QAE) [21, 1] can be used to estimate such
expectation values if the quantity of interest can be efficiently encoded as
the probability of a specific measurement outcome. QAE takes as input a
unitary operator $\mathcal{A}$ which produces the state

$\mathcal{A}:\ket{0}_{n+1}\rightarrow\left(\sqrt{V}\ket{\psi_{0}}_{n}\ket{0}+\sqrt{1-V}\ket{\psi_{1}}_{n}\ket{1}\right),$
(54)

where $\ket{\psi_{0}}$ and $\ket{\psi_{1}}$ are arbitrary, normalized quantum
states and estimates the derivative price $V$ to additive error $\epsilon$
using $\mathcal{O}\left(1/\epsilon\right)$ invocations of the operator
$\mathcal{Q}=\mathcal{A}S_{0}\mathcal{A}^{\dagger}S_{\psi_{1}}$, where
$S_{0}=\mathbb{I}-2|0\rangle_{n+1}\langle 0|_{n+1}$ and
$S_{\psi_{1}}=\mathbb{I}-2\ket{1}\ket{\psi_{1}}_{n}\bra{1}\bra{\psi_{1}}_{n}$.

A method to construct the $\mathcal{A}$ operator of Eq. (54) for typical
derivative contracts of practical interest and the corresponding required
quantum resources is shown in [6]. For a contract with $d$ stochastic market
parameters $X=[x_{1},x_{2},\cdots x_{d}]$, an operator $\mathcal{P}$
constructs a probability-weighted superposition of all possible combinations
of stochastic realizations of the market parameters,

$\mathcal{P}:\ket{\vec{0}}^{\otimes
d}\rightarrow\sum_{\omega}\sqrt{p(\omega)}\ket{\omega},$ (55)

with an appropriate discretization for each parameter chosen such that the
overall error from such discretization satisfies a desired accuracy.
Typically, each market parameter is modeled stochastically at pre-defined
points in time that are determined by the definition of the derivative
contract. As such, each $\ket{\omega}$ contains all the information about the
values of the market parameters $X$ for each _path_ of the stochastic process.
Then, an operator $\mathcal{F}$ computes the discounted payoff $f$ of the
derivative on each path $\ket{\omega}$, and encodes that value in the
amplitude of an ancilla qubit

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

Taking $\mathcal{A}=\mathcal{F}\mathcal{P}$

$\mathcal{F}\mathcal{P}:\ket{\vec{0}}^{\otimes
d}\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),$
(57)

gives us the desired operator of Eq. (54), noticing that the probability of
measuring $\ket{0}$ in the last qubit is the expectation value of the
discounted payoff, as defined in Eq. (53).

This approach can be generalized to compute the value $V$ of a derivative
portfolio consisting of $m$ derivatives. The $\mathcal{P}$ operator of Eq.
(55) in this case will construct the stochastic process of the $d$ market
parameters of the entire portfolio, and the $\mathcal{F}$ operator of Eq. (56)
will compute the sum of the payoffs of each derivative

$\mathcal{F}:\ket{\omega}\ket{\vec{0}}\ket{0}\rightarrow\ket{\omega}\ket{\sum_{i=1}^{m}f_{i}(\omega)}\left(\sqrt{\sum_{i=1}^{m}f_{i}(\omega)}\ket{0}+\sqrt{1-\sum_{i=1}^{m}f_{i}(\omega)}\ket{1}\right),$
(58)

where $f_{i}$ is the payoff of derivative $i\in[1,m]$. Similarly to the
single-derivative case, QAE can then be used to estimate the expectation value
$V$ of the portfolio.

$V=\sum_{\omega\in\Omega}p(\omega)\sum_{i=1}^{m}f_{i}(\omega).$ (59)