# Quantum algorithm for
copula-based risk aggregation
using orthogonal series density estimation

1st Hitomi Mori  Graduate School of Engineering Science
Osaka University
Toyonaka, Osaka, Japan
<EMAIL_ADDRESS>2nd Koichi Miyamoto  Center for Quantum Information
and Quantum Biology
Osaka University
Toyonaka, Osaka, Japan
<EMAIL_ADDRESS>
###### Abstract

Quantum Monte Carlo integration (QMCI) provides a quadratic speed-up over its
classical counterpart, and its applications have been investigated in various
fields, including finance. This paper considers its application to risk
aggregation, one of the most important numerical tasks in financial risk
management. Risk aggregation combines several risk variables and quantifies
the total amount of risk, taking into account the correlation among them. For
this task, there exists a useful tool called copula, with which the joint
distribution can be generated from marginal distributions with a flexible
correlation structure. Classically, the copula-based method utilizes sampling
of risk variables. However, this procedure is not directly applicable to the
quantum setting, where sampled values are not stored as classical data, and
thus no efficient quantum algorithm is known. In this paper, we propose a
quantum algorithm for copula-based risk aggregation that is compatible with
QMCI. In our algorithm, we first estimate each marginal distribution as a
series of orthogonal functions, where the coefficients can be calculated with
QMCI. Then, by plugging the marginal distributions into the copula and
obtaining the joint distribution, we estimate risk measures using QMCI again.
With this algorithm, nearly quadratic quantum speed-up can be obtained for
sufficiently smooth marginal distributions.

###### Index Terms:

Quantum algorithm, Quantum Monte Carlo integration, Finance, Risk aggregation

##  I Introduction

Monte Carlo integration is a widely used method to perform numerical
integration using random sampling. It is known that the quantum algorithm for
Monte Carlo integration provides a quadratic speed-up over its classical
counterpart [1]. In finance, Monte Carlo integration is used in various tasks,
including risk analysis and derivative pricing. Studies have explored
potential applications of quantum Monte Carlo integration (QMCI) in these
domains [2].

Risk aggregation stands as a paramount objective within risk analysis and
plays a crucial role in the management of financial institutions [3]. Risk
aggregation is preceded by risk analysis, wherein substantial risks for a
company are identified and quantified. For example, insurance companies
typically have two major risk categories: underwriting risk and investment
risk. Underwriting risk reflects the uncertainty of the total amount of loss
that may occur from their insurance contracts, while investment risk arises
from the investment of insurance premiums. Furthermore, these categories are
broken down into more granular risk categories. Each risk is then quantified
using the loss distribution obtained through Monte Carlo simulation with an
appropriate stochastic model. The risk amount can be determined as the
quantile of the distribution. In risk aggregation, these simulated loss
distributions are combined to derive the total amount of risk, considering the
correlation among different risks111When risks have perfect correlation,
meaning the correlation coefficients are 1, the total risk amount is just the
sum of each risk amount. However, in general, the correlation coefficients are
within $[-1,1]$, and thus the total risk amount can be smaller than the simple
sum, benefiting from the so-called diversification effect.. Consequently, a
company can obtain the total risk amount representing a certain level (e.g.
1-in-100 year) of potential total loss. By comparing the total risk amount
with its capital, the company can assess the sufficiency of its capital.
Moreover, the company can calculate other related metrics222Such as return on
risk (the ratio of return to risk amount) and risk contribution (the amount
contributing to total risk amount). for each business unit to develop business
plans and determine capital allocation.

To incorporate the correlation structure into the modelling, a statistical
tool called copula [4] is widely used in finance. Copula is a function that
takes marginal distribution functions of multiple random variables and returns
the joint distribution function. This allows for the modeling of multiple risk
variables with a flexible correlation structure. When copula is employed
within Monte Carlo simulation, where computations are conducted based on
sample values, the correlation is reflected by rearranging the samples. In
this process, the initial step involves sorting the sample values of each risk
variable stored in memory. Then, the set of sample values is rearranged
according to the order specified by the copula. This rearranged set is
regarded as a set of samples from the joint distribution. However, this
approach is not directly applicable to QMCI, where sample values are not
stored in memory, and thus no efficient quantum algorithm is known.

In this paper, we propose a quantum algorithm for copula-based risk
aggregation that aligns with QMCI. In our algorithm, assuming access to the
oracles to generate quantum states encoding the possible values of the risk
variables, we first estimate the marginal distribution functions of each risk
variable using orthogonal series density estimation (OSDE) [5]. In OSDE,
distribution functions are approximated as a series of orthogonal functions
with coefficients represented as expected values of these functions. QMCI is
then employed to estimate these expected values. Subsequently, we derive the
joint distribution by plugging the estimated marginal distribution functions
into the copula. Finally, we estimate the risk measures from the joint
distribution using QMCI. Through this approach, we successfully obtain the
quantum speed-up. The number of queries to the state generation oracle of each
risk variable scales on the required accuracy $\epsilon$ as
$\widetilde{O}\left(\left(\frac{1}{\epsilon}\right)^{1+\frac{8}{2r-1}}\right)$,
where $r$ is a parameter characterizing the smoothness of the marginal
distribution functions. Thus, if the smoothness $r$ is high, the query
complexity approaches $\widetilde{O}(1/\epsilon)$, which represents a
quadratic improvement over the sample complexity of order
$\widetilde{O}(1/\epsilon^{2})$ in the classical algorithm.

The organization of this paper is as follows. In Section II, we introduce
notations and preliminary knowledge on risk aggregation, OSDE, and QMCI. We
present the quantum algorithm for risk aggregation in Section III along with
an analysis on the error and complexity of this algorithm, followed by the
conclusion in Section IV.

##  II Preliminaries

In this section, we introduce notations and preliminary knowledge used to
develop the quantum algorithm in Section III.

###  II-A Notations

Throughout this paper, we use the notations commonly used in statistics. We
use upper case letter (e.g. $X$) to denote random variable and use lower case
letter (e.g. $x$) to denote the value of the random variable. For a
$\mathbb{R}^{d}$-valued random variable $\bm{X}=(X_{1},\cdots,X_{d})$, its
probability density function (PDF) $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$
satisfies

$\mathrm{Pr}(\alpha_{1}\leq X_{1}\leq\beta_{1},\cdots,\alpha_{d}\leq
X_{d}\leq\beta_{d})\\\
=\int^{\alpha_{1}}_{\beta_{1}}dt_{1}\cdots\int^{\alpha_{d}}_{\beta_{d}}dt_{d}f(t_{1},\cdots,t_{d}),$
(1)

and the cumulative distribution function (CDF)
$F:\mathbb{R}^{d}\rightarrow[0,1]$ is

$\displaystyle F(x_{1},\cdots,x_{d})$ $\displaystyle=\mathrm{Pr}(X_{1}\leq
x_{1},\cdots,X_{d}\leq x_{d})$ (2)
$\displaystyle=\int_{-\infty}^{x_{1}}dt_{1}\cdots\int_{-\infty}^{x_{d}}dt_{d}f(t_{1},\cdots,t_{d}).$
(3)

$\mathbb{E}_{\bm{X}}[\cdot]$ denotes the expected value; for example,

$\displaystyle\mathbb{E}_{\bm{X}}[\phi(\bm{X})]=\int_{-\infty}^{\infty}dt_{1}\cdots\int_{-\infty}^{\infty}dt_{d}\phi(t_{1},\cdots,t_{d})f(t_{1},\cdots,t_{d}).$
(4)

Estimators are represented by the hat symbol, such as $\hat{a}$.

As a quantum computing platform, we hereafter consider a system of multiple
qubits. When we express a variable $x\in\mathbb{R}$ with $n$ qubits, for
$N=2^{n}$ and $j\in\\{0,1,\cdots,N-1\\}$, we represent the computational basis
as $\ket{j}$ and the value of $x$ on the $j$th grid point as $x^{(j)}$. We
assume that every quantum register has a sufficient number of qubits and thus
consider discretization errors to be negligible.

$\mathbbm{1}_{C}$ denotes the indicator function, which takes 1 if the
condition $C$ is satisfied and $0$ otherwise.

For $x\in\mathbb{R}$ and $\epsilon>0$, we say that $x^{\prime}\in\mathbb{R}$
is an $\epsilon$-approximation of $x$ if $|x-x^{\prime}|\leq\epsilon$.

###  II-B Risk aggregation

In financial mathematics, the term ”risk” indicates uncertainties regarding
potential profits or losses. Risks are quantified according to risk measures,
with two widely used types: Value at Risk (VaR) and Tail Value at Risk
(TVaR)333Also known as Conditional Value at Risk (CVaR).. VaR is a certain
quantile of the loss distribution, and it is most common in financial
regulations and disclosures. VaR of a random variable $X$ at a confidence
level $\alpha\in[0,1]$ is expressed as

$\displaystyle\mathrm{VaR}_{X}(\alpha):=\inf\\{x\in\mathbb{R}:F(x)\geq\alpha\\}=F^{-1}(\alpha).$
(5)

TVaR is the conditional expected value of the loss in the tail of the
distribution that exceeds the VaR. Because this can capture extreme loss
events included in the tail, it is preferred for natural catastrophe risks
such as hurricane risk. TVaR of $X$ at a confidence level $\alpha$ is
expressed as

$\displaystyle\mathrm{TVaR}_{X}(\alpha)$
$\displaystyle:=\mathbb{E}_{X}[X|X\geq\mathrm{VaR}_{X}(\alpha)]$
$\displaystyle=\frac{1}{1-\alpha}\int^{\infty}_{\mathrm{VaR}_{X}(\alpha)}tf(t)dt.$
(6)

Although risk measures for each risk variable $X_{1},\cdots,X_{d}$ are of
practical interest, those for the sum $S=X_{1}+\cdots+X_{d}$ are also crucial.
Risk aggregation estimates this total amount of risk by modeling the joint
distribution of $X_{1},\cdots,X_{d}$. For this purpose, copulas are often
used, because they enable the modeling of a flexible correlation structure. A
copula is a CDF for a $[0,1]^{d}$-valued random variable, where each variable
is uniformly distributed on $[0,1]$. By using a copula, we can handle the
individual univariate marginal distributions and their dependency separately,
thanks to Sklar’s theorem, which guarantees the consistency between the
copula-based joint distribution and each marginal distribution:

######  Theorem 1 (Sklar’s Theorem [4]).

For any joint distribution function $F$ for $\mathbb{R}^{d}$-valued random
variable $\bm{X}=(X_{1},\cdots,X_{d})$ with marginal distributions
$F_{1},...,F_{d}$, there exists a copula $C$ such that

$\displaystyle F(x_{1},...,x_{d})=C(F_{1}(x_{1}),...,F_{d}(x_{d}))$ (7)

for any $(x_{1},\cdots,x_{d})\in\mathbb{R}^{d}$. If $F_{1},...,F_{d}$ are
continuous, then $C$ is unique. Conversely, for any marginal distributions
$F_{1},...,F_{d}$ and $d$-variate copula $C$, the function $F$ defined as in
(7) is a CDF with marginal distributions $F_{1},...,F_{d}$.

Defining the density of the copula as

$\displaystyle c(u_{1},...,u_{d})=\frac{\partial}{\partial u_{1}...\partial
u_{d}}C(u_{1},...,u_{d}),$ (8)

the joint PDF can be expressed as

$\displaystyle
f(x_{1},...,x_{d})=c(F_{1}(x_{1}),...,F_{d}(x_{d}))\prod^{d}_{i=1}f_{i}(x_{i}).$
(9)

The popular choices of copula include Gaussian copula and T-copula, the
details of which can be found in appendix A.

In practice, we often rely on samples for risk calculation. In many cases,
explicit formulas for each marginal CDF $F_{i}$ may not be available, but we
can prepare samples of $X_{i}$444For example, if the random time evolution of
$X_{i}$ is modeled by a stochastic differential equation (SDE), which is
common for market risk variables, we generally do not have an explicit formula
for the CDF of $X_{i}$ at a future time $T$. However, we can generate sample
paths of $X_{i}$ from the present to $T$ using some discretization method of
the SDE, such as the Euler-Maruyama method [6].. Using samples, $F_{i}$ can be
estimated based on ordering [7], where we consider the $k$-th smallest sample
value as the $100\frac{k}{N}$-percentile point of $X_{i}$. Specifically,
denoting the $N$ sample values of $X_{i}$ sorted in the ascending order by
$x_{i}^{(1)},\cdots,x_{i}^{(N)}$, we define the percentile function
$\tilde{F}_{i}^{-1}:[0,1]\rightarrow\mathbb{R}$ of $X_{i}$ by

$\tilde{F}_{i}^{-1}(u):=\begin{cases}x_{i}^{(1)}&\mathrm{for}\ 0\leq
u\leq\frac{1}{N},\\\ x_{i}^{(2)}&\mathrm{for}\
\frac{1}{N}<u\leq\frac{2}{N},\\\ &\vdots\\\ x_{i}^{(N)}&\mathrm{for}\
\frac{N-1}{N}<u\leq 1\end{cases}.$ (10)

When we have $N$ samples for each $X_{i}$ from $F_{i}$ and
$(U_{1},\cdots,U_{d})\in[0,1]^{d}$ from $C$555The procedure for sampling from
a copula can be found in [3], we can generate samples from the joint
distribution in (7) as described in Algorithm 1. Given the sample values
$\bm{x}^{(1)},\cdots,\bm{x}^{(N)}$, we estimate $\mathrm{VaR}_{S}(\alpha)$ as
$\tilde{F}_{S}^{-1}(\alpha)$, which is defined similarly to (10), and
$\mathrm{TVaR}_{S}(\alpha)$ as

$\frac{\sum_{l=1}^{N}s^{(l)}\mathbbm{1}_{s^{(l)}\geq\mathrm{VaR}_{S}(\alpha)}}{\sum_{l=1}^{N}\mathbbm{1}_{s^{(l)}\geq\mathrm{VaR}_{S}(\alpha)}},$
(11)

where $s^{(l)}:=\sum_{i=1}^{d}\tilde{x}_{i}^{(l)}$. It is shown that samples
$\bm{x}^{(l)}$ generated by this method asymptotically obey $F$ in (7) with
the CDF converging at a rate $O(N^{-1/2})$, and so do the VaR and TVaR
estimations from the samples666See [7] for mathematically rigorous
statements..

Algorithm 1 Sampling from the copula-based joint distribution

1:

2:Samples from each marginal $F_{i}$ and copula $C$.

3:for $i=1,\cdots,d$ do

4: Get $N$ sample values of $X_{i}\sim F_{i}$

5: Sort samples in the ascending order $x_{i}^{(1)},\cdots,x_{i}^{(N)}$.

6:end for

7:for $l=1,\cdots,N$ do

8: Sample $(U_{1}^{(l)},...,U_{d}^{(l)})\sim C$.

9: Let $\bm{x}^{(l)}:=(\tilde{x}^{(l)}_{1},\cdots,\tilde{x}^{(l)}_{d})$, where
$\tilde{x}_{i}^{(l)}:=\tilde{F}_{i}^{-1}(u_{i}^{(l)})$.

10:end for

11:Output $\bm{x}^{(1)},\cdots,\bm{x}^{(N)}$.

Note that this procedure based on sampling and rearranging is not directly
applicable to the quantum algorithm, because QMCI does not use sample values
of $X_{i}$ but rather a superposition of quantum states that encode the
possible values of $X_{i}$, as discussed in Section II-D.

###  II-C Orthogonal series density estimation

Orthogonal series density estimation [5] is a technique used to estimate the
PDF $f$ of a random variable $X\in\mathbb{R}$ using orthogonal series. As a
specific choice, we employ a series of Hermite functions as outlined in [8]:

$h_{k}(x)=\frac{1}{\sqrt{2^{k}k!\pi^{1/2}}}H_{k}(x)e^{-x^{2}/2},$ (12)

where

$H_{k}=(-1)^{k}e^{x^{2}}\frac{d^{k}}{dx^{k}}e^{-x^{2}}$ (13)

is the $k$-th Hermite polynomial for $k=0,1,\cdots$. The Hermite function
$h_{k}$ satisfies the following orthogonal relationship:

$\displaystyle\int_{-\infty}^{\infty}h_{k}(x)h_{l}(x)dx=\delta_{kl},$ (14)

where $k,l\in\mathbb{N}$ and $\delta_{kl}$ is the Kronecker’s delta. We
approximate $f$ by the Hermite series ${f}^{(K)}$ defined as

$\displaystyle{f}^{(K)}(x):=\sum^{K}_{k=0}a_{k}h_{k}(x)$ (15)

for $K\in\mathbb{N}$ with the coefficients given by

$a_{k}=\int_{-\infty}^{\infty}f(x)h_{k}(x)dx.$ (16)

Because $f$ is the PDF of the random variable $X$, $a_{k}$ can be expressed as
the expected value

$a_{k}=\mathbb{E}_{X}[h_{k}(X)].$ (17)

The error bound for the truncated Hermite series is given by the following
theorem, which is extracted from the discussion leading to Theorem 5 in [8].

######  Theorem 2 (Approximation by Hermite series).

Let $f:\mathbb{R}\rightarrow\mathbb{R}$ be in $L^{\infty}\cap L^{2}$ and of
bounded variation on any finite interval. Besides, for some $r\in\mathbb{N}$,
$\left(x-\frac{d}{dx}\right)^{r}f\in L^{2}$ holds. Then, for any $L>0$, there
exists $\gamma_{f,L}\in\mathbb{R}$ such that

$\sup_{x\in[-L,L]}\left|{f}^{(K)}(x)-f(x)\right|\leq\gamma_{f,L}K^{-\frac{r}{2}+\frac{1}{4}}$
(18)

holds for any $K\in\mathbb{N}$.

The Hermite function also has the following property: for any $k=0,1,\cdots$
and any $x\in\mathbb{R}$,

$|h_{k}(x)|\leq\pi^{-1/4}$ (19)

holds [9].

###  II-D Quantum Monte Carlo integration

We first introduce quantum amplitude estimation (QAE) [10], which serves as
the central subroutine in QMCI and is the source of the quadratic speed-up.
QAE is a technique to estimate the squared amplitude of a specific basis state
in a superposition, and the error bound of QAE is given by the following
theorem.

######  Theorem 3 (Quantum amplitude estimation [10]).

There is a quantum algorithm called amplitude estimation which takes as input
one copy of a quantum state $\ket{\psi}$, a unitary transformation
$U=2\ket{\psi}\bra{\psi}-I$, a unitary transformation $V=I-2P$ for some
projector $P$, and an integer $t$. The algorithm outputs $\tilde{a}$, an
estimate of $a=\bra{\psi}P\ket{\psi}$, such that

$\displaystyle\left|\tilde{a}-a\right|\leq
2\pi\frac{\sqrt{a(1-a)}}{t}+\frac{\pi^{2}}{t^{2}}$ (20)

with probability at least $8/\pi^{2}$, using $U$ and $V$ $t$ times each.

Using QAE, we can construct QMCI, the quantum algorithm to estimate
$\mathbb{E}_{\bm{Y}}[\phi(\bm{Y})]$ for a random variable
$\bm{Y}\in\mathbb{R}^{d}$ and a function
$\phi:\mathbb{R}^{d}\rightarrow[0,1]$. Because we assume that we have
sufficient qubits to represent real numbers with the negligible discretization
error, we can express the expected value in terms of summation:

$\mathbb{E}_{\bm{Y}}[\phi(\bm{Y})]=\frac{1}{\mathcal{N}_{\bm{Y}}^{2}}\sum_{\bm{j}\in\mathcal{G}_{\bm{Y}}}\phi(y_{1}^{(j_{1})},\cdots,y_{d}^{(j_{d})})f(y_{1}^{(j_{1})},\cdots,y_{d}^{(j_{d})}),$
(21)

where
$\mathcal{G}_{\bm{Y}}=\mathcal{G}_{Y_{1}}\times\cdots\times\mathcal{G}_{Y_{d}}$,
$\mathcal{G}_{Y_{i}}=\\{0,\cdots,N_{i}-1\\}$ (i.e.
$\sum_{\bm{j}\in\mathcal{G}_{\bm{Y}}}=\sum_{j_{1}=0}^{N_{1}-1}\cdots\sum_{j_{d}=0}^{N_{d}-1}$),
and $y_{i}^{(j_{i})}$ is the $j_{i}$-th discretization point in $y_{i}$ axis.
Here we define two types of oracles required to run QMCI.

######  Definition 1 (State preparation).

For a random variable $\bm{Y}=(Y_{1},\cdots,Y_{d})\in\mathbb{R}^{d}$ with the
PDF $f$, we have access to the unitary operation $\mathcal{A}_{\bm{Y}}$ that
acts as

$\mathcal{A}_{\bm{Y}}\ket{0}^{\otimes n_{1}}\cdots\ket{0}^{\otimes n_{d}}\\\
=\frac{1}{\mathcal{N}_{\bm{Y}}}\sum_{\bm{j}\in\mathcal{G}_{\bm{Y}}}\sqrt{f(\bm{y}^{(\bm{j})})}\ket{j_{1}}\cdots\ket{j_{d}},$
(22)

where $\bm{y}^{(\bm{j})}=(y_{1}^{(j_{1})},\cdots,y_{d}^{(j_{d})})$ and
$\mathcal{N}_{\bm{Y}}:=\sqrt{\sum_{\bm{j}\in\mathcal{G}_{\bm{Y}}}f(\bm{y}^{(\bm{j})})}$.

If we know the exact form of the PDF, then the state preparation can be
achieved using existing methods [11, 12, 13, 14, 15, 16]. Even if we do not
know the exact form, given a quantum circuit to simulate ${\bm{Y}}$, we can
perform the state preparation. For example, if ${\bm{Y}}$ represents a value
of a stochastic process at a specific time and the process is modeled by an
SDE, the state can be prepared as described in [17, 18].

######  Definition 2 (Controlled rotation).

For a function $\phi:\mathbb{R}^{d}\to[0,1]$, we have access to the unitary
operation $W_{\phi}$ that acts as

$W_{\phi}\ket{j_{1}}\cdots\ket{j_{d}}\ket{0}\\\
=\ket{j_{1}}\cdots\ket{j_{d}}\left(\sqrt{1-\phi(\bm{y}^{(\bm{j})})}\ket{0}+\sqrt{\phi(\bm{y}^{(\bm{j})})}\ket{1}\right).$
(23)

for any $\bm{y}^{(\bm{j})}\in\mathbb{R}^{d}$.

Methods to implement this for some specific cases can be found in [19, 20,
17]. If $\phi$ can be expressed as an elementary function, we can compute
$\theta=\arcsin\left(\sqrt{\phi(\bm{y})}\right)$ onto an ancillary register
using arithmetic circuits [21, 22] and then apply a rotation gate with the
angle controlled
$\ket{\theta}\ket{0}\rightarrow\ket{\theta}(\cos\theta\ket{0}+\sin\theta\ket{1})$.

Applying $W_{\phi}(\mathcal{A}_{\bm{Y}}\otimes I)$ to the initial state
$\ket{0}^{\otimes s}\ket{0}$ with $s=n_{1}+\cdots+n_{d}$, we can prepare the
state

$\ket{\psi}:=\frac{1}{\mathcal{N}_{\bm{Y}}}\sum_{\bm{j}\in\mathcal{G}_{\bm{Y}}}\sqrt{f(\bm{y}^{(\bm{j})})}\ket{j_{1}}\cdots\ket{j_{d}}\\\
\left(\sqrt{1-\phi(\bm{y}^{(\bm{j})})}\ket{0}+\sqrt{\phi(\bm{y}^{(\bm{j})})}\ket{1}\right).$
(24)

Note that for $P=I^{\otimes s}\otimes\ket{1}\bra{1}$,

$\braket{\psi}{P}{\psi}=\frac{1}{\mathcal{N}_{\bm{Y}}^{2}}\sum_{\bm{j}\in\mathcal{G}_{\bm{Y}}}\phi(\bm{y}^{(\bm{j})})f(\bm{y}^{(\bm{j})})={\mathbb{E}}_{\bm{Y}}[\phi(\bm{Y})].$
(25)

Thus, we can estimate ${\mathbb{E}}_{\bm{Y}}[\phi(\bm{Y})]$ using QAE with

$\displaystyle U$ $\displaystyle=2\ket{\psi}\bra{\psi}-I$ (26)
$\displaystyle=W_{\phi}(\mathcal{A}_{\bm{Y}}\otimes I)\left(2\ket{0}^{\otimes
s}\bra{0}^{\otimes
s}\otimes\ket{0}\bra{0}-I\right)(W_{\phi}(\mathcal{A}_{\bm{Y}}\otimes
I))^{\dagger}.$ (27)

The complexity of this method, QMCI, is summarized as follows.

######  Theorem 4 (Quantum Monte Carlo integration; Theorem 2.3 in [1],
modified).

Let $\delta,\epsilon\in(0,1)$. Assuming we have access to the state
preparation oracle $\mathcal{A}_{\bm{Y}}$ for a random variable
$\bm{Y}\in\mathbb{R}^{d}$ and the controlled rotation oracle $W_{\phi}$ for a
function $\phi:\mathbb{R}^{d}\rightarrow[0,1]$, there exists a quantum
algorithm called $\textup{{QMCI}}(\mathcal{A}_{\bm{Y}},\phi,\epsilon,\delta)$
that, with probability at least $1-\delta$, outputs an
$\epsilon$-approximation of ${\mathbb{E}}_{\bm{Y}}[\phi(\bm{Y})]$, querying
$\mathcal{A}_{\bm{Y}}$ and $W_{\phi}$
$O\left(\frac{1}{\epsilon}\log\left(\frac{1}{\delta}\right)\right)$ times
each.

Note that the success probability is enhanced from $8/\pi^{2}$ in Theorem 3 to
$1-\delta$ by taking the median of outputs of multiple runs of QAE (Lemma 2.1
in [1]), which is based on Lemma 6.1 in [23].

##  III Quantum algorithm for risk aggregation

###  III-A Overview

In this section, we show our main contribution: the algorithm for risk
aggregation. First, we note that risk aggregation for a set of risk variables
$\bm{X}=(X_{1},\cdots,X_{d})$ can be formulated as the calculation of the
expected value in the following form:

$\displaystyle\mathbb{E}_{\bm{X}}\left[g\left(\bm{X}\right)\right]$
$\displaystyle=\frac{1}{\mathcal{N}_{\bm{X}}^{2}}\sum_{\bm{x}\in\mathcal{G}_{\bm{X}}}g(\bm{x}^{(\bm{j})})f(\bm{x}^{(\bm{j})})$
$\displaystyle=\frac{1}{\mathcal{N}_{\bm{X}}^{2}}\sum_{\bm{x}\in\mathcal{G}_{\bm{X}}}g(\bm{x}^{(\bm{j})})c(F_{1}(x_{1}^{(j_{1})}),...,F_{d}(x_{d}^{(j_{d})}))\prod^{d}_{i=1}f_{i}(x_{i}^{(j_{i})})$
$\displaystyle=\mathbb{E}^{\mathrm{ind}}_{\bm{X}}[g(\bm{X})c(F_{1}(X_{1}),...,F_{d}(X_{d}))].$
(28)

Here, $g:\mathbb{R}^{d}\rightarrow[0,1]$ is a function chosen according to the
type of the risk measure, as explained later.
$\mathbb{E}^{\mathrm{ind}}_{\bm{X}}[\cdot]$ is the expected value under the
assumption that $X_{1},\cdots,X_{d}$ are independent of each other, i.e. the
PDF is expressed as

$f^{\mathrm{ind}}(x_{1},\cdots,x_{d})=\prod^{d}_{i=1}f_{i}(x_{i}).$ (29)

When we have the state preparation oracles
$\mathcal{A}_{X_{1}},\cdots,\mathcal{A}_{X_{d}}$ for $X_{1},\cdots,X_{d}$, we
can generate the following quantum state

$\mathcal{A}_{\bm{X}}^{\rm ind}\ket{0}^{\otimes
d}=\frac{1}{\mathcal{N}_{\bm{X}}}\sum_{\bm{j}\in\mathcal{G}_{\bm{X}}}\sqrt{f^{\mathrm{ind}}(\bm{x}^{(\bm{j})})}\ket{j_{1}}\cdots\ket{j_{d}},$
(30)

where $\mathcal{A}_{\bm{X}}^{\rm
ind}:=\mathcal{A}_{X_{1}}\otimes\cdots\otimes\mathcal{A}_{X_{d}}$ and
$\mathcal{N}_{\bm{X}}=\prod_{i=1}^{d}\mathcal{N}_{X_{i}}$ with
$\mathcal{N}_{X_{i}}:=\sqrt{\sum_{j_{i}\in\mathcal{G}_{X_{i}}}f_{i}(x_{i}^{(j_{i})})}$.
This motivates us to estimate (28) by QMCI. Since we are given explicit
formulas for $c$ and $g$ in many cases, the remaining requirement for running
QMCI is the formulation of marginal CDFs $F_{1},\cdots,F_{d}$. To do this, we
use OSDE, which gives us explicit approximation formulas of the CDFs as
orthogonal series. In summary, our algorithm is composed of two steps:
estimating the marginal CDFs by OSDE and estimating the risk measures.

###  III-B Estimating the marginal distributions

Now, we explain the process to estimate the marginal CDFs. We assume the
availability of the state preparation oracle $A_{X_{i}}$ with $X_{i}$ and the
controlled rotation with $\phi=\bar{h}_{k}$ for $k=0,1,\cdots$, where

$\bar{h}_{k}(x)=\frac{1}{2}+\frac{\pi^{1/4}}{2}h_{k}(x),$ (31)

which satisfies $0\leq\bar{h}_{k}(x)\leq 1$. Then, we run QMCI to estimate
${\mathbb{E}}_{X_{i}}[\bar{h}_{k}(X_{i})]$ for $k=0,\cdots,K_{i}$ and get the
estimations $\hat{a}^{X_{i}}_{k}$ for

$a^{X_{i}}_{k}:=\mathbb{E}_{X_{i}}[h_{k}(X_{i})]=\frac{1}{\pi^{1/4}}\left(2{\mathbb{E}}_{X_{i}}[\bar{h}_{k}(X_{i})]-1\right),$
(32)

the exact coefficient in the Hermite expansion for $f$. Then, given the
estimated coefficients, we have a Hermite series approximating $f_{i}$ as

$\hat{f}_{i}(x)=\sum^{K_{i}}_{k=0}\hat{a}^{X_{i}}_{k}h_{k}(x).$ (33)

Furthermore, by integrating this, we get an approximation $\hat{F}_{i}$ for
$F_{i}$. In fact, since the accuracy of the Hermite series approximation is
guaranteed not on the entire real axis but in a finite interval, we set
$\hat{F}_{i}$ to 0 or 1 outside the interval. Namely, we define

$\displaystyle\hat{F}_{i}(x):=\begin{cases}0&;\ x<-L_{i}\\\
\sum^{K_{i}}_{k=0}\hat{a}^{X_{i}}_{k}\mathcal{H}_{k,L_{i}}(x)&;-L_{i}\leq
x<L_{i}\\\ 1&;\ x\geq L_{i}\end{cases}$ (34)

with sufficiently large $L_{i}>0$, where

$\mathcal{H}_{k,L_{i}}(x):=\int_{-L_{i}}^{x}h_{k}(t)dt.$ (35)

Now, we just assume that we can evaluate $\mathcal{H}_{k,L_{i}}$ with a
negligible error, postponing the discussion on this point to Appendix B.

The above method is summarized as Algorithm 2.

Algorithm 2 $\textup{{EstMarg}}(X_{i},K_{i},\epsilon,\delta)$: estimation of
$F_{i}$

1:

2:Access to the state preparation oracle $\mathcal{A}_{X_{i}}$.

3:Access to the rotation oracles $W_{\bar{h}_{k}}$ for $k=0,\cdots,K_{i}$.

4:$L_{i}>0$ and $\epsilon,\delta\in(0,1)$.

5:for $k=0,\dots,K_{i}$ do

6: Run
$\textup{{QMCI}}\left(\mathcal{A}_{X_{i}},\bar{h}_{k},\frac{\sqrt{\pi}\epsilon}{16L_{i}(K_{i}+1)},\frac{\delta}{K_{i}+1}\right)$
and let the output be the estimation
$\hat{\mathbb{E}}_{X_{i}}[\bar{h}_{k}(X_{i})]$.

7: Set
$\hat{a}^{X_{i}}_{k}:=\frac{1}{\pi^{1/4}}\left(2\hat{\mathbb{E}}_{X_{i}}[\bar{h}_{k}(X_{i})]-1\right)$.

8:end for

9:Define $\hat{F}_{i}$ as (34) and output it.

We give the following theorem on the accuracy and complexity of this method.

######  Theorem 5.

Let $\delta,\epsilon\in(0,1)$. Let $X_{i}$ be a real-valued random variable
following the distribution $f_{i}$. Assume the following:

  1. (i)

$f_{i}$ has the properties of Theorem 2 for $r=r_{i}$.

  2. (ii)

Access to the state preparation oracle $\mathcal{A}_{X_{i}}$.

  3. (iii)

Access to the controlled rotation oracle $W_{\phi}$ with $\phi=\tilde{h}_{k}$
for $k=0,\cdots,K_{i}$.

  4. (iv)

For some $L_{i}>0$, $F_{i}(-L_{i})\leq\frac{\epsilon}{2}$ and
$F_{i}(L_{i})\geq 1-\epsilon$ hold.

Then, with probability at least $1-\delta$, we get $\hat{F}_{i}$ such that

$\left|\hat{F}_{i}(x)-F_{i}(x)\right|\leq\epsilon$ (36)

holds for any $x\in\mathbb{R}$ by
$\textup{{EstMarg}}\left(X_{i},K_{i},\epsilon,\delta\right)$, which queries
$\mathcal{A}_{X_{i}}$ and $W_{\bar{h}_{0}},\cdots,W_{\bar{h}_{K_{i}}}$

$\displaystyle
O\left(\frac{L_{i}K_{i}^{2}}{\epsilon}\log\left(\frac{K_{i}}{\delta}\right)\right)$
(37)

times. Here we set

$K_{i}=\left\lceil\left(\frac{8L_{i}\gamma_{i}}{\epsilon}\right)^{4/(2r_{i}-1)}\right\rceil,$
(38)

where $\gamma_{i}$ is a real number such that

$\sup_{x\in[-L_{i},L_{i}]}\left|{f}^{(K)}_{i}(x)-f_{i}(x)\right|\leq\gamma_{i}K^{-\frac{r_{i}}{2}+\frac{1}{4}}$
(39)

holds for any $K\in\mathbb{N}$ 777Because of (i), such $\gamma_{i}$ exists, as
implied by Theorem 2..

###### Proof.

Because of the definition of $\hat{F}_{i}$ in (34) and (iv), it is immediately
seen that (36) holds for any $x\in(-\infty,-L_{i}]\cup[L_{i},\infty)$. Thus,
we hereafter focus on the case that $x\in(-L_{i},L_{i})$.

We start by evaluating $\left|\hat{f}_{i}(x)-f_{i}(x)\right|$. Decomposing it
as

$\left|\hat{f}_{i}(x)-f_{i}(x)\right|\leq\left|\hat{f}_{i}(x)-{f}^{(K_{i})}_{i}(x)\right|+\left|{f}^{(K_{i})}_{i}(x)-f_{i}(x)\right|,$
(40)

we bound the first and second terms separately. On the second term, for
$x\in(-L_{i},L_{i})$, we have

$\left|{f}^{(K_{i})}_{i}(x)-f_{i}(x)\right|\leq\gamma_{i}K_{i}^{-\frac{r_{i}}{2}+\frac{1}{4}}\leq\frac{\epsilon}{8L_{i}},$
(41)

where we use (38) and (39). To bound the first one, we temporarily assume that
$\textup{{QMCI}}\left(\mathcal{A}_{X_{i}},\bar{h}_{k},\frac{\sqrt{\pi}\epsilon}{16L_{i}(K_{i}+1)},\frac{\delta}{K_{i}+1}\right)$
for every $k=0,\cdots,K_{i}$ outputs the estimation
$\hat{\mathbb{E}}_{X_{i}}[\bar{h}_{k}(X_{i})]$ such that

$\left|\hat{\mathbb{E}}_{X_{i}}[\bar{h}_{k}(X_{i})]-{\mathbb{E}}_{X_{i}}[\bar{h}_{k}(X_{i})]\right|\leq\frac{\sqrt{\pi}\epsilon}{16L_{i}(K_{i}+1)}.$
(42)

This implies that

$\displaystyle\left|\hat{a}_{k}^{X_{i}}-{a}_{k}^{X_{i}}\right|$
$\displaystyle=\frac{2}{\pi^{1/4}}\left|\hat{\mathbb{E}}_{X_{i}}[\bar{h}_{k}(X_{i})]-{\mathbb{E}}_{X_{i}}[\bar{h}_{k}(X_{i})]\right|$
(43) $\displaystyle\leq\frac{\pi^{1/4}\epsilon}{8L_{i}(K_{i}+1)}.$ (44)

We then have

$\displaystyle\left|\hat{f}_{i}(x)-{f}^{(K_{i})}_{i}(x)\right|$
$\displaystyle=\left|\sum_{k=0}^{K_{i}}(\hat{a}_{k}^{X_{i}}-a_{k}^{X_{i}})h_{k}(x)\right|$
$\displaystyle\leq\sum_{k=0}^{K_{i}}\left|\hat{a}_{k}^{X_{i}}-a_{k}^{X_{i}}\right|\pi^{-1/4}$
$\displaystyle\leq\frac{\epsilon}{8L_{i}},$ (45)

where we use (19) at the first inequality, and (44) at the third inequality.
Combining (40), (41), and (45), we get

$\left|\hat{f}_{i}(x)-f_{i}(x)\right|\leq\frac{\epsilon}{4L_{i}}.$ (46)

Then, this yields

$\displaystyle\left|F_{i}(x)-\hat{F}_{i}(x)\right|$
$\displaystyle\leq\left|F_{i}(-L_{i})\right|+\left|\int_{-L_{i}}^{x}f_{i}(y)dy-\int_{-L_{i}}^{x}\hat{f}_{i}(y)dy\right|$
$\displaystyle\leq\frac{\epsilon}{2}+\int_{-L_{i}}^{x}\left|\hat{f}_{i}(y)-f_{i}(y)\right|dy$
$\displaystyle\leq\frac{\epsilon}{2}+\frac{\epsilon}{4L_{i}}\times(x+L_{i})$
$\displaystyle\leq\epsilon$ (47)

for $x\in(-L_{i},L_{i})$.

To complete the proof, let us prove the statements on the success probability
and complexity. Since the probability that each of
$\textup{{QMCI}}\left(\mathcal{A}_{X_{i}},\bar{h}_{k},\frac{\sqrt{\pi}\epsilon}{16L_{i}(K_{i}+1)},\frac{\delta}{K_{i}+1}\right)$
for $k=0,\cdots,K_{i}$ outputs $\hat{\mathbb{E}}_{X_{i}}[\bar{h}_{k}(X_{i})]$
satisfying (42) is at least $1-\frac{\delta}{K_{i}+1}$, the probability that
all of them output such estimations is at least
$\left(1-\frac{\delta}{K_{i}+1}\right)^{K_{i}+1}\geq 1-\delta$. The number of
queries to $\mathcal{A}_{X_{i}}$ and $W_{\bar{h}_{k}}$ in
$\textup{{QMCI}}\left(\mathcal{A}_{X_{i}},\bar{h}_{k},\frac{\sqrt{\pi}\epsilon}{16L_{i}(K_{i}+1)},\frac{\delta}{K_{i}+1}\right)$
is

$O\left(\frac{L_{i}(K_{i}+1)}{\epsilon}\log\left(\frac{K_{i}+1}{\delta}\right)\right)$
(48)

as implied by Theorem 4. Summing this up for $k=0,\cdots,K_{i}$, we reach the
total query number estimation in (37). ∎

###  III-C Estimating the risk measure

Given the approximations of the CDFs $\hat{F}_{1},\cdots,\hat{F}_{d}$ by the
method above, we finally consider estimating the expected value
$\mathbb{E}_{\bm{X}}\left[g\left(\bm{X}\right)\right]$ in (28). We now assume
that we have access to the controlled rotation oracle $W_{\Phi}$ for the
function $\Phi:\mathbb{R}^{d}\rightarrow[0,1]$ expressed by

$\Phi(\bm{x})=\frac{1}{c_{\rm
max}}g(\bm{x})c(\hat{F}_{1}(x_{1}),...,\hat{F}_{d}(x_{d}))$ (49)

with $c_{\rm max}:=\max_{\bm{u}\in[0,1]^{d}}c(\bm{u})$ and any functions
$\hat{F}_{1},\cdots,\hat{F}_{d}$ written in the form of (34). Then, we present
Algorithm 3 to estimate
$\mathbb{E}_{\bm{X}}\left[g\left(\bm{X}\right)\right]$.

Algorithm 3 $\textup{{EstRM}}(c,g,\epsilon,\delta)$: estimation of
$\mathbb{E}_{F_{\bm{X}}}\left[g\left(\bm{X}\right)\right]$

1:

2:Access to the state preparation oracles
$\mathcal{A}_{X_{1}},\cdots,\mathcal{A}_{X_{d}}$.

3:Access to the rotation oracles $W_{\bar{h}_{k}}$ for $k=0,\cdots,K_{i}$ and
$W_{\Phi}$ for $\Phi$ in (49).

4:$L_{1},\cdots,L_{d}>0$ and $\epsilon,\delta\in(0,1)$.

5:for $i=1,\dots,d$ do

6: Set $K_{i}$ as
$K_{i}=\left\lceil\left(\frac{16dc_{\mathrm{max}}^{\prime}L_{i}\gamma_{i}}{\epsilon}\right)^{4/(2r_{i}-1)}\right\rceil$
(50) with $c_{\mathrm{max}}^{\prime}$ given in Theorem 6.

7: Run
$\textup{{EstMarg}}\left(X_{i},K_{i},\frac{\epsilon}{2dc_{\mathrm{max}}^{\prime}},\frac{\delta}{2d}\right)$
and let the output be $\hat{F}_{i}$.

8:end for

9:Run $\textup{{QMCI}}\left(\mathcal{A}^{\rm
ind}_{\bm{X}},\hat{\Phi},\frac{\epsilon}{2c_{\mathrm{max}}},\frac{\delta}{2}\right)$,
where  $\hat{\Phi}(\bm{x})=\frac{1}{c_{\rm
max}}g(\bm{x})c(\hat{F}_{1}(x_{1}),...,\hat{F}_{d}(x_{d})),$ (51) and let the
result be $\hat{\mathbb{E}}^{\mathrm{ind}}_{\bm{X}}(\hat{\Phi})$.

10:Output $c_{\rm max}\hat{\mathbb{E}}^{\mathrm{ind}}_{\bm{X}}(\hat{\Phi})$.

The following theorem is on the accuracy and complexity of this algorithm.

######  Theorem 6.

Let $\delta,\epsilon\in(0,1)$. Let $c:[0,1]^{d}\rightarrow\mathbb{R}$ be the
density of copula and suppose that there exists
$c^{\prime}_{\mathrm{max}}\in\mathbb{R}$ such that, for any $i=1,\cdots,d$ and
any $\bm{u}\in[0,1]^{d}$,

$\left|\frac{\partial}{\partial u_{i}}c(u_{1},\cdots,u_{d})\right|\leq
c^{\prime}_{\mathrm{max}}.$ (52)

Let $X_{1},\cdots,X_{d}$ be real-valued random variables and suppose that all
the assumptions in Theorem 7 are satisfied for every $X_{i}$ except that
$\epsilon$ is replaced with $\frac{\epsilon}{2dc_{\mathrm{max}}^{\prime}}$.
Suppose that we have access to the rotation oracle $W_{\Phi}$ for any function
in the form of (49) with $g:\mathbb{R}^{d}\rightarrow[0,1]$. Then, with
probability at least $1-\delta$, Algorithm 3 outputs an
$\epsilon$-approximation of
$\mathbb{E}_{\bm{X}}\left[g\left(\bm{X}\right)\right]$, querying
$W_{\bar{h}_{k}}$

$\displaystyle
O\left(\frac{d^{2}c^{\prime}_{\mathrm{max}}LK^{2}}{\epsilon}\log\left(\frac{dK}{\delta}\right)\right)$
(53)

times, $W_{\Phi}$

$O\left(\frac{c_{\mathrm{max}}}{\epsilon}\log\left(\frac{1}{\delta}\right)\right)$
(54)

times, and $\mathcal{A}_{X_{i}}$

$O\left(\frac{d^{2}c^{\prime}_{\mathrm{max}}LK^{2}}{\epsilon}\log\left(\frac{dK}{\delta}\right)+\frac{c_{\mathrm{max}}}{\epsilon}\log\left(\frac{1}{\delta}\right)\right)$
(55)

times, where $L:=\max_{i=1,\cdots,d}L_{i}$ and $K:=\max_{i=1,\cdots,d}K_{i}$.

###### Proof.

As implied by Theorem 7, $\hat{F}_{i}$ obtained in step 7 in Algorithm 3
satisfies

$\forall x\in\mathbb{R},\
\left|\hat{F}_{i}(x)-F_{i}(x)\right|\leq\frac{\epsilon}{2dc_{\mathrm{max}}^{\prime}}$
(56)

with probability at least $1-\frac{\delta}{2d}$. If this holds, for any
$\bm{x}\in\mathbb{R}^{d}$,

$\displaystyle\left|\hat{\Phi}(\bm{x})-\Phi(\bm{x})\right|$
$\displaystyle\leq\frac{g(\bm{x})}{c_{\rm
max}}\left|c(\hat{F}_{1}(x_{1}),...,\hat{F}_{d}(x_{d}))-c(F_{1}(x_{1}),...,F_{d}(x_{d}))\right|$
$\displaystyle\leq\frac{c_{\mathrm{max}}^{\prime}g(\bm{x})}{c_{\mathrm{max}}}\sum_{i=1}^{d}\left|\hat{F}_{i}(x)-F_{i}(x)\right|$
$\displaystyle\leq\frac{\epsilon}{2c_{\mathrm{max}}}$ (57)

holds, where we use Taylor’s theorem with (52). We then have

$\displaystyle\left|{\mathbb{E}}^{\mathrm{ind}}_{\bm{X}}[\hat{\Phi}\left(\bm{X}\right)]-{\mathbb{E}}^{\mathrm{ind}}_{\bm{X}}\left[\Phi\left(\bm{X}\right)\right]\right|$
$\displaystyle=\left|\frac{1}{\mathcal{N}_{\bm{X}}^{2}}\sum_{\bm{j}\in\mathcal{G}_{\bm{X}}}\left(\hat{\Phi}(\bm{x}^{(\bm{j})})-\Phi(\bm{x}^{(\bm{j})})\right)f^{\mathrm{ind}}(\bm{x}^{(\bm{j})})\right|$
$\displaystyle\leq\frac{1}{\mathcal{N}_{\bm{X}}^{2}}\sum_{\bm{j}\in\mathcal{G}_{\bm{X}}}\left|\hat{\Phi}(\bm{x}^{(\bm{j})})-\Phi(\bm{x}^{(\bm{j})})\right|f^{\mathrm{ind}}(\bm{x}^{(\bm{j})})$
$\displaystyle\leq\frac{\epsilon}{2c_{\mathrm{max}}}.$ (58)

On the other hand,
$\textup{{QMCI}}\left(\mathcal{A}^{\mathrm{ind}}_{\bm{X}},\hat{\Phi},\frac{\epsilon}{2c_{\mathrm{max}}},\frac{\delta}{2}\right)$
outputs $\hat{\mathbb{E}}_{\bm{X}}[\hat{\Phi}(\bm{X})]$ such that

$\left|\hat{\mathbb{E}}^{\mathrm{ind}}_{\bm{X}}[\hat{\Phi}(\bm{X})]-{\mathbb{E}}^{\mathrm{ind}}_{\bm{X}}[\hat{\Phi}\left(\bm{X}\right)]\right|\leq\frac{\epsilon}{2c_{\mathrm{max}}}$
(59)

with probability at least $1-\delta$. Combining (58) and (59), we get

$\displaystyle\left|c_{\mathrm{max}}\hat{\mathbb{E}}^{\mathrm{ind}}_{\bm{X}}[\hat{\Phi}(\bm{X})]-\mathbb{E}_{\bm{X}}\left[g\left(\bm{X}\right)\right]\right|$
$\displaystyle\leq\
c_{\mathrm{max}}\left(\left|\hat{\mathbb{E}}^{\mathrm{ind}}_{\bm{X}}[\hat{\Phi}(\bm{X})]-{\mathbb{E}}^{\mathrm{ind}}_{\bm{X}}[\hat{\Phi}\left(\bm{X}\right)]\right|\right.$
$\displaystyle\qquad\qquad\left.+\left|{\mathbb{E}}^{\mathrm{ind}}_{\bm{X}}[\hat{\Phi}\left(\bm{X}\right)]-{\mathbb{E}}^{\mathrm{ind}}_{\bm{X}}[\Phi\left(\bm{X}\right)]\right|\right)$
$\displaystyle\leq\ \epsilon.$ (60)

This holds if every $\hat{F}_{i}$ satisfies (56) and (59) holds, whose
probability is at least

$\left(1-\frac{\delta}{2d}\right)^{d}\left(1-\frac{\delta}{2}\right)\geq
1-\delta.$ (61)

Lastly, let us evaluate the query complexity of the algorithm. In estimating
each $\hat{F}_{i}$, $\mathcal{A}_{X_{i}}$ and $\\{W_{h_{k}}\\}_{k}$ are
queried

$O\left(\frac{dc^{\prime}_{\mathrm{max}}L_{i}K_{i}^{2}}{\epsilon}\log\left(\frac{dK_{i}}{\delta}\right)\right)$
(62)

times, and, in estimating all of $\hat{F}_{1},\cdots,\hat{F}_{d}$, this is
multiplied by $d$. In $\textup{{QMCI}}\left(\mathcal{A}^{\rm
ind}_{\bm{X}},\hat{\Phi},\frac{\epsilon}{2c_{\mathrm{max}}},\frac{\delta}{2}\right)$,
$\mathcal{A}^{\rm ind}_{\bm{X}}$ and $W_{\hat{\Phi}}$ are called the number of
times of order (54). In $\mathcal{A}_{\bm{X}}^{\rm ind}$,
$\mathcal{A}_{X_{1}},\cdots,\mathcal{A}_{X_{d}}$ are called once each and $d$
times in total. Combining these observations, we reach the query number bounds
in (53), (54), and (55).

∎

Substituting
$K=\left\lceil(16dc_{\mathrm{max}}^{\prime}L\gamma/\epsilon)^{4/(2r-1)}\right\rceil$
for $\gamma:=\max_{i=1,\cdots,d}\gamma_{i}$ and $r:=\min_{i=1,\cdots,d}r_{i}$
into (55), we obtain the number of queries to $\mathcal{A}_{X_{i}}$ scaling on
accuracy $\epsilon$ as

$\displaystyle\tilde{O}\left((1/\epsilon)^{1+\frac{8}{2r-1}}\right).$ (63)

If the smoothness parameter $r$ is high, the query complexity asymptotically
approaches $\tilde{O}(1/\epsilon)$, and the algorithm provides nearly
quadratic quantum speed-up.

Finally, we discuss how to utilize the algorithm EstRM to calculate risk
measures. To calculate ${\rm VaR}_{S}(\alpha)$, we first construct a procedure
to compute $\mathrm{Pr}(S\geq l)$ for $S=X_{1}+\cdots+X_{d}$ and
$l\in\mathbb{R}$ based on Algorithm 3, with $g$ set to

$\displaystyle g_{\rm VaR}(\bm{x};l)=\mathbbm{1}_{s\geq l}.$ (64)

Then, for $\alpha\in(0,1)$, we search the smallest $l_{\alpha}$ satisfying

$\displaystyle\mathrm{Pr}(S\geq l_{\alpha})\leq 1-\alpha,$ (65)

which corresponds to VaR${}_{S}(\alpha)$, using some root-finding method. For
example, we can use binary search, which increases the total query complexity
only by a logarithmic factor. For TVaR, using the calculated VaR $l_{\alpha}$,
we run Algorithm 3 setting $g$ to

$\displaystyle g_{\rm
TVaR}(\bm{x};l_{\alpha})=\sum_{i=1}^{d}\frac{x_{i}}{x_{\rm
max}}\mathbbm{1}_{s\geq l_{\alpha}},$ (66)

where $x_{\rm max}:=\max_{i=1,...,d}x_{i}^{(N_{i})}$, and let the output
multiplied by $x_{\rm max}$ be the estimation of TVaR.

##  IV Conclusion

We developed a quantum algorithm to conduct risk aggregation, a fundamental
task in financial institutions. In particular, we considered a standard
approach that uses copula, a tool to generate a joint distribution from
marginal distributions, to express the structure of the correlations among the
risk variables. However, the common classical algorithm for this approach is
based on sampling and rearranging, which cannot be straightforwardly applied
in the quantum setting, and thus no efficient quantum algorithm was known. To
address this challenge, we developed a quantum method consisting of two steps.
First, we estimate each marginal CDF based on OSDE. As each coefficient in the
orthogonal series is expressed as the expected value of the orthogonal
function, we apply QMCI to estimate it. Then, we plug the estimated CDFs into
the copula to get the joint distribution, from which we calculate the risk
measure using QMCI once more. Using this method, we achieved a nearly
quadratic quantum speed-up of risk aggregation, compared to the classical
method for smooth marginal PDFs.

## Acknowledgment

This work is supported by MEXT Quantum Leap Flagship Program (MEXTQLEAP) Grant
No. JPMXS0120319794, and JST COI-NEXT program Grant No. JPMJPF2014. KM is also
supported by JSPS KAKENHI Grant no. JP22K11924.

##  Appendix A Copulas

For readers unfamiliar with copula, we introduce some widely used types:

  * •

Gaussian copula: for any $n\times n$ correlation matrix $R=(\rho_{ij})$, the
Gaussian copula is defined by

$C=\Phi_{n,R}(\Phi_{\mathrm{SN}}^{-1}(\cdot),\cdots,\Phi_{\mathrm{SN}}^{-1}(\cdot)),$
(67)

where $\Phi_{\mathrm{SN}}$ is the CDF for the standard normal distribution and
$\Phi_{n,R}$ is the CDF for the $n$-variate normal distribution with zero mean
vector and covariance matrix equal to $R$. This is one of the most widely used
copulas in practice because of its simplicity, but its tail dependence can be
too weak in some situations.

  * •

T-copula: with a parameter $\nu>0$ called the number of degrees of freedom
(DOF) and a positive-define matrix $\Sigma\in\mathbb{R}^{n\times n}$, a
T-copula is defined by

$C=t_{n,\nu,\Sigma}(t_{\nu}^{-1}(\cdot),\cdots,t_{\nu}^{-1}(\cdot)),$ (68)

where $t_{\nu}$ is the CDF for the Student’s $t$ distribution with the number
of DOF $\nu$, and $t_{n,\nu,\Sigma}$ is the CDF for the $n$-variate $t$
distribution with the number of DOF $\nu$ and the covariance matrix $\Sigma$.
This copula exhibits stronger tail dependence compared to Gaussian copulas and
is preferred in certain scenarios.

##  Appendix B Evaluation of $\mathcal{H}_{k,L}$

To evaluate $\mathcal{H}_{k,L}(x)$, we first consider how to evaluate

$\tilde{\mathcal{H}}_{l}(x):=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x}y^{l}e^{-y^{2}/2}dy$
(69)

for $l\in\mathbb{N}\cup\\{0\\}$. We can utilize a recurrence formula

$\tilde{\mathcal{H}}_{l}(x)=-\frac{1}{\sqrt{2\pi}}x^{l-1}e^{-x^{2}/2}+(l-1)\tilde{\mathcal{H}}_{l-2}(x),l\geq
2$ (70)

with

$\tilde{\mathcal{H}}_{1}(x)=-\frac{1}{\sqrt{2\pi}}e^{-x^{2}/2},\tilde{\mathcal{H}}_{0}(x)=\Phi_{\rm
SN}(x).$ (71)

Thus, we can compute $\tilde{\mathcal{H}}_{l}(x)$ by evaluating the terms in
the form of $x^{l}e^{-x^{2}/2}$ and $\Phi_{\rm SN}(x)$. We can evaluate the
former on a quantum computer using circuits for elementary functions [21]. On
the latter, there are some approximation formulas with elementary functions
(see, e.g., [24]), and thus it is also evaluated on a quantum computer.

## References

  * [1] Ashley Montanaro.  Quantum speedup of monte carlo methods.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, Vol. 471, No. 2181, p. 20150301, September 2015.
  * [2] Dylan Herman, Cody Googin, Xiaoyuan Liu, Yue Sun, Alexey Galda, Ilya Safro, Marco Pistoia, and Yuri Alexeev.  Quantum computing for finance.  Nature Reviews Physics, Vol. 5, No. 8, pp. 450–465, 2023.
  * [3] Paul Embrechts, Alexander McNeil, and Daniel Straumann.  Correlation and dependence in risk management: properties and pitfalls.  Risk management: value at risk and beyond, Vol. 1, pp. 176–223, 2002.
  * [4] Roger B Nelsen.  An Introduction to Copulas.  Springer, New York, NY, USA, second edition, 2006.
  * [5] Sam Efromovich.  Orthogonal series density estimation.  WIREs Computational Statistics, Vol. 2, No. 4, pp. 467–476, 2010.
  * [6] G. Maruyama.  Continuous markov processes and stochastic equations.  Rend. Circ. Mat. Palermo, Vol. 4, pp. 48–90, 1955.
  * [7] Georg Mainik.  Risk aggregation with empirical margins: Latin hypercubes, empirical copulas, and convergence of sum distributions.  Journal of Multivariate Analysis, Vol. 141, pp. 197–216, 2015.
  * [8] Wlodzimierz Greblicki and Miroslaw Pawlak.  Hermite series estimates of a probability density and its derivatives.  Journal of Multivariate Analysis, Vol. 15, No. 2, pp. 174–182, 1984.
  * [9] Gabor Szegö.  Orthogonal polynomials.  American Mathematical Soc., Providence, 1939.
  * [10] Gilles Brassard, Peter Høyer, Michele Mosca, and Alain Tapp.  Quantum amplitude amplification and estimation, 2002.
  * [11] Lov Grover and Terry Rudolph.  Creating superpositions that correspond to efficiently integrable probability distributions.  2002\.
  * [12] Yuval R. Sanders, Guang Hao Low, Artur Scherer, and Dominic W. Berry.  Black-box quantum state preparation without arithmetic.  Phys. Rev. Lett., Vol. 122, p. 020502, Jan 2019.
  * [13] Arthur G. Rattew and BÃ¡lint Koczor.  Preparing arbitrary continuous functions in quantum registers with logarithmic complexity, 2022.
  * [14] Adam Holmes and A. Y. Matsuura.  Efficient quantum circuits for accurate state preparation of smooth, differentiable functions, 2020.
  * [15] Sam McArdle, András Gilyén, and Mario Berta.  Quantum state preparation without coherent arithmetic.  2022\.
  * [16] Mudassir Moosa, Thomas W Watts, Yiyou Chen, Abhijat Sarma, and Peter L McMahon.  Linear-depth quantum circuits for loading fourier approximations of arbitrary functions.  Quantum Science and Technology, Vol. 9, No. 1, p. 015002, October 2023.
  * [17] Almudena Carrera Vazquez and Stefan Woerner.  Efficient state preparation for quantum amplitude estimation.  Physical Review Applied, Vol. 15, No. 3, March 2021.
  * [18] Kazuya Kaneko, Koichi Miyamoto, Naoyuki Takeda, and Kazuyoshi Yoshino.  Quantum pricing with a smile: implementation of local volatility model on quantum computer.  EPJ Quantum Technology, Vol. 9, No. 1, p. 7, 2022.
  * [19] Stefan Woerner and Daniel J. Egger.  Quantum risk analysis.  npj Quantum Information, Vol. 5, No. 1, feb 2019.
  * [20] Nikitas Stamatopoulos, Daniel J. Egger, Yue Sun, Christa Zoufal, Raban Iten, Ning Shen, and Stefan Woerner.  Option pricing using quantum computers.  Quantum, Vol. 4, p. 291, July 2020.
  * [21] Thomas Häner, Martin Roetteler, and Krysta M Svore.  Optimizing quantum circuits for arithmetic.  arXiv preprint arXiv:1805.12445, 2018.
  * [22] Edgard Muñoz Coreas and Himanshu Thapliyal.  Everything you always wanted to know about quantum circuits.  Wiley Encyclopedia of Electrical and Electronics Engineering, pp. 1–17, 2022.
  * [23] Mark R. Jerrum, Leslie G. Valiant, and Vijay V. Vazirani.  Random generation of combinatorial structures from a uniform distribution.  Theoretical Computer Science, Vol. 43, pp. 169–188, 1986.
  * [24] Milton Abramowitz and Irene A Stegun.  Handbook of mathematical functions with formulas, graphs, and mathematical tables.  Dover, New York, 1968.