量子コンピュータと金融
~ポートフォリオのリスクを計算してみる~

導入
量子ビットと量子回路
1. 量子回路の測定
2. 量子振幅推定
量子回路でのVaR計算
量子回路でのVaR計算(改良版)
まとめ
参考資料

こんにちは、次世代システム研究室のT.I.です。

さて、世の中には様々なリスクに満ち溢れていますが、できることならリスクは回避したいですよね。(ハイリスク・ハイリターンという言葉もありますが、ローリスク・ハイリターンですませたいものです。)「卵は１つのカゴに盛るな」、「大事なものには予備が必要だ」といつも言われているように、リスクを回避する常套手段の１つが複数の選択肢を予め準備しておくことです。

昨今、世界中で量子コンピュータを巡っての研究開発が進められ、様々な分野への応用が期待されています。金融分野では、複数の金融資産の組み合わせ(ポートフォリオ)の最適化問題がまず挙げられる例ですが、今回はポートフォリオのリスクを量子コンピュータにより効果的に評価する方法を解説してみたいと思います。

TL;DR

量子コンピュータの金融分野への応用例としてポートフォリオのリスクを計算してみた。
損失発生率から合計損失額の集計と累積分布まで量子回路で計算し、量子振幅推定を利用し従来の手法よりも効率的に計算が可能となる。

導入

ポートフォリオとリスク分散

複数の資産の組み合わせをポートフォリオと呼びます。これにより単純に１つの金融資産に全資産を投資するよりもリスクを低減することができます。具体的に $\begin{align*} n\end{align*}$ 個の金融資産、それぞれの収益率を $\begin{align*} R_i\end{align*}$ とし、 $\begin{align*} w_i\end{align*}$ の割合で投資します。構成されたポートフォリオの収益は以下の通りです。

$\begin{align*} R_p = \sum_{i=1}^n w_i R_i \end{align*}$
ここで、割合の合計は $\begin{align*} \sum_{i=1}^n w_i = 1\end{align*}$ を満たすとします。

期待される収益はそれぞれの合算値
$\begin{align*} E(R_p) = \sum_{i=1}^n w_i E(R_i) \end{align*}$
となります。これだけでは単純に分けただけですが、重要なのはリスク(分散)を低減することです。収益の分散 $\begin{align*} \sigma_i = \sqrt{\mathrm{Var}(R_i)}\end{align*}$ と相関係数 $\begin{align*} \rho_{ij} = \mathrm{Corr}(R_i, R_j)/(\sigma_i \sigma_j) \end{align*}$ から、このポートフォリオの分散は以下となります。

$\begin{align*} \sigma_p^2 = \sum_{i=1}^n \sum_{j=1}^n w_i w_j \mathrm{Cov}(R_i, R_j) = \sum_{i=1}^n \sum_{j=1}^n w_i w_j \sigma_i \sigma_j \rho_{ij}\end{align*}$
ここで投資戦略として分散を最小化する組み合わせを決めてみます。

簡単のために２つの金融資産(A, B)に投資するケースを考え、各々の期待収益と分散、共分散を $\begin{align*} \mu_A, \mu_B, \sigma_A, \sigma_B, \rho_{AB}\end{align*}$ とし、A(B)の割合を $\begin{align*} w\end{align*}$ ( $\begin{align*} 1-w\end{align*}$ )とします。

$\begin{align*} \mu_p = w\mu_A + (1-w)\mu_B \end{align*}$
$\begin{align*} \sigma_p^2 = w^2\sigma_A^2 + (1-w)^2\sigma_B^2 + 2w(1-w)\sigma_A\sigma_B\rho_{AB} \end{align*}$
さて、 $\begin{align*} \sigma_p^2\end{align*}$ を最小化する $\begin{align*} w\end{align*}$ は、 $\begin{align*} \partial \sigma_p^2/\partial w = 0\end{align*}$ の条件から

$\begin{align*} \frac{\partial \sigma_p^2}{\partial w} = 2w \sigma_A^2 - 2(1-w)\sigma_B^2 + 2(1 - 2w) \sigma_A \sigma_B \rho_{AB} = 2w \left(\sigma_A^2 +\sigma_B^2 - 2 \sigma_A\sigma_B\rho_{AB}\right) + 2 \left(-\sigma_B^2 + \sigma_A \sigma_B \rho_{AB}\right) = 0 \end{align*}$ したがって

$\begin{align*} w^\ast = \frac{\sigma_B^2 - \sigma_A \sigma_B \rho_{AB}} {\sigma_A^2 + \sigma_B^2 - 2 \sigma_A\sigma_B\rho_{AB} } \end{align*}$
の重さで資産を組み合わせると、個別収益の多いAのみを運用するよりも期待収益は少なくなりますが、リスクである分散を減らした投資が可能となります。

A(B)の期待収益 8(2)、分散 8%(5%)、および相関係数 -0.5 の場合のポートフォリオの収益と標準偏差の関係

２つの資産の場合はこのように簡単に計算できますが、投資先が多い場合計算は複雑になります。量子コンピュータと金融分野という文脈では、このポートフォリオの最適化問題として紹介されることが多いですが、今回はポートフォリオのリスクの評価といった別の観点からの応用を紹介します。

VaR (Value at Risk)

金融資産のリスクを図る一つの指標が VaR (Value at Risk) です。これは保有している金融資産(ポートフォリオ)にある一定の確率で発生しうる損失額を指しています。
以下の図はある資産で想定される損益額の分布となっています。損益は０を中心に分布しており、可能性は低いですが100以上の損失額も生じ得ます。具体的に想定される損失の発生する確率とその金額がVaRとなります。例えば、1 – alpha = 0.95つまり、95%以内の水準でしたら、その損失はVaR(0.05) = -82.2となります。

Monte Carlo法

金融資産のリスク(VaR)の推定に利用される手法の１つにMonte Carlo法があります。これは乱数を用いて数値計算やシミュレーションする計算手法です。ちなみに、この名前はモナコ公国にあるカジノで有名な地名に由来しています。

Monte Carlo法の具体的な例として、円周率を計算してみます。半径１の円がぴったり収まる１辺が２の正方形を考え、その中でランダムに１点を選択します。その点 $\begin{align*} (x,y) \end{align*}$ が円の内部 $\begin{align*} x^2 + y^2 \leq 1 \end{align*}$ にある確率は円の面積と正方形の面積の比である

$\begin{align*} p = \frac{\pi}{4} \end{align*}$

となります。したがって、たくさんの点を抽出し、それが円の内部である確率を計算していくと円周率が推定できます。以下のグラフはサンプル数を徐々に増やして円周率を分析したものです。サンプル数を増やすほど精度が高くなっていきます。 Monte Carlo法の誤差はサンプル数(N)に対して、 $\begin{align*} N^{-1/2} \end{align*}$ で小さくなります。つまり、精度を10倍にしたいならば、100倍の計算量が必要となります。

また、別の例としてランダム・ウォークする金融資産価格の一定期間後の価格についての
Monte Carloシミュレーションを見てみます。一定期間後の価格はシミュレーションの度に様々な値となりますが、極端に変動するものは稀で、多くは元の価格付近に集まっています。その累積分布(CDF cumulative distribution function)は測定回数を増やすほど滑らかになり精度よくVaRを推定できます。ただし、分布のテールのように大きく価格が変動する場合のリスク推定にはより多くの測定が必要となります。

量子ビットと量子回路

まずは、量子ビット(キュービット)の復習です。古典的ビットとは０か１の一方の値しかとることができませんが、量子ビットは０でも１でもある”重ね合わせ”状態を取ることができます。

$\begin{align*} |\psi\rangle = \alpha |0 \rangle + \beta | 1 \rangle \end{align*}$
しかし、この重ね合わせ状態そのものは観測することができません。実際の量子ビットを読みだすと重ね合わせ状態は破壊され、結果として $\begin{align*} |\alpha|^2\end{align*}$ の確率で0、 $\begin{align*} |\beta|^2\end{align*}$ の確率で1が得られます(つまり、 $\begin{align*} |\alpha|^2 + |\beta|^2 = 1\end{align*}$ を満たします)。１度の観測ではどちらか一歩しか得られませんので、沢山の量子ビットを繰り返し測定することで統計的にこれらの確率を推定します。

量子回路の測定

まずは、簡単な回路を作成して実験してみます。

from qiskit import QuantumCircuit

circ = QuantumCircuit(1)
circ.ry(np.pi/4,0)
circ.measure_all()
circ.draw('mpl', scale=2.)

この回路は０の量子ビットをY軸を中心に $\begin{align*} \pi/4\end{align*}$ (45度)回転させるもので、数式では次のように表され

$\begin{align*} R_Y(\pi/4)|0\rangle = \cos\frac{\pi}{8} |0\rangle + \sin\frac{\pi}{8} | 1 \rangle \end{align*}$
０と１の重ね合わせ状態となっております。この量子ビットの状態を模式的に表現したのが下のBloch球です。

この回路の測定をシミュレーションしてみます。量子回路の計算の実行には主に２つのモードが存在しています。

状態ベクトル計算
- 量子回路の数理モデルを $\begin{align*} 2^n\end{align*}$ 次元のベクトルを計算します。実際には観測不可能な量ではありますが、正確な量子重ね合わせの状態を理解できます。
- 得られた状態ベクトルから、量子ビットの測定結果を誤差なく推定できます。
量子ビット測定
- 現実に実装されている量子コンピュータやその挙動をシミュレーションして、量子ビットを測定します。
- 多数の測定を繰り返すことで、各々のビットの測定確率を推定します。精度のよい推定には多くの測定が必要です。

以下はQASMバックエンド(量子ビット測定のシミュレーション)による計算の実行です。

from qiskit import Aer, execute

backend_sim = Aer.get_backend('qasm_simulator')
job_sim = execute(circ, backend_sim, shots=1)
result_sim = job_sim.result()
counts = result_sim.get_counts(circ)

測定の度に、０または１が得られますので、それを繰り返し確率を更新していく様子が以下の図です。

最初は不安定で誤差も大きいですが、徐々に厳密値である sin(π/8)^2 ~ 0.854, cos(π/8)^2 ~ 0.146に近づきます。先ほどMonte Carlo法で解説したように統計の揺らぎは測定回数の平方根の逆数に比例するので、誤差を半分にするには４倍の統計量が必要となります。

量子振幅推定

量子コンピュータのアルゴリズムにグローバーのアルゴリズムというものがあります。これは探索問題を解くアルゴリズムで要素の数Nの平方根の計算時間でスケールする効率の良いものです。そのアイデアを拡張したものが量子振幅推定です。

n+1個の量子ビットを準備して、ある操作Aを実行した後に最後のビットが１となる確率をaとすると以下のような重ね合わせ状態になります。

$\begin{align*} {A}|0\rangle_{n+1} = \sqrt{1-a}|\psi_0\rangle_n |0\rangle + \sqrt{a} |\psi_1\rangle_n | 1 \rangle \end{align*}$
この確率 a が、この問題で測定したい量です。そのために新たに m 個の量子ビットを追加して作られる量子振幅推定の回路は以下のような構造となっています。

量子振幅推定回路(図は論文 arXiv:1907.03044 より引用)

詳細は割愛しますが、QはAを含む状態を操作する演算で、回路の上側にm個追加された量子ビットで結果を読み込みます。測定用の量子ビット(m)を増やすほど有効桁数の多い数値が得られます。この回路をM回測定し推定した誤差は

$\begin{align*} |a - \tilde{a}| \leq \frac{2\sqrt{a(1-a)}\pi}{M} + \frac{\pi^2}{M^2} = \mathcal{O}\left(\frac{1}{M}\right) \end{align*}$
となります。これは、測定回数の平方根の逆数で誤差が減少するMonte Carlo計算よりも効率的です。特にVaRの推定では0.1%などの稀な確率に興味がありますから、更に近似して誤差は

$\begin{align*} \frac{1}{5M} + \frac{\pi^2}{M^2} \end{align*}$
となります。

Qiskitでは既に量子振幅推定は実装されており、次のように利用できます。１個の量子ビットをY軸に沿って回転させた先ほどの例と同じ回路ですが、今回の場合は量子ビットの状態を破壊する測定をせずに組み合わせる点に注意してください。

from qiskit.algorithms import IterativeAmplitudeEstimation, EstimationProblem
from qiskit.utils import QuantumInstance

# set target precision and confidence level
circ = QuantumCircuit(1)
circ.ry(np.pi/4, 0)
# circ.measure_all() Do not measure!

epsilon = 0.01
alpha = 0.05

qi = QuantumInstance(Aer.get_backend('qasm_simulator'), shots=100)
problem = EstimationProblem(state_preparation=circ, objective_qubits=[0])

ae = IterativeAmplitudeEstimation(epsilon, alpha=alpha, quantum_instance=qi)
result = ae.estimate(problem)

print(r'p_1 = {:.4f}'.format(result.estimation_processed))
print('Confidence interval: [{:.4f}, {:.4f}]'.format(
      result.confidence_interval_processed[0],
      result.confidence_interval_processed[1]))

# p_1 = 0.1460
# Confidence interva: [0.1421, 0.1500]

量子回路でのVaR計算

なお、以下の解析は基本的には Qiskit で公開されている解説と
(https://qiskit.org/documentation/finance/tutorials/09_credit_risk_analysis.html)、
元の論文 Credit Risk Analysis using Quantum Computers [arXiv:1907.03044] (https://arxiv.org/abs/1907.03044)を参考に補足などを追加したものとなっております。

ポートフォリオと損失額

K 個の金融資産があり、それぞれの発生しうる損失額 $\begin{align*} (L_1, \cdots, L_K) \in \mathbf{R}_{>0}^K\end{align*}$ 、
ポートフォリオの損失額とその期待値は以下の通りです。
$\begin{align*} \mathcal{L} = \sum_{k=1}^K L_k, \quad E[\mathcal{L}] = \sum_{k=1}^K E[L_k] \end{align*}$
ここで問題となるVaRは数式では以下の要因なります。
$\begin{align*} \mathrm{Var}_\alpha[\mathcal{L}] = \mathrm{inf}_{x \geq 0}\ \{ x \ | P[\mathcal{L} \leq x] \geq \alpha \} \end{align*}$
$\begin{align*} X_k \in \{ 0, 1\}\end{align*}$ を $\begin{align*} p_k\end{align*}$ の確率で１を取る変数として、損失額を以下のように仮定します。

$\begin{align*} L_k = \lambda_k X_k \end{align*}$ すると合計損失額は
$\begin{align*} E[\mathcal{L}] = \sum_{k=1}^K \lambda_k p_k\end{align*}$
各々のパラメータを決定すればこの期待値の計算そのものは非常に簡単です。

Gaussian conditional independence (GCI) model

このモデルはまだ各々の資産のリスクが独立なので、より現実的なリスクが連動するモデルを考えます。元の損失が発生する確率 $\begin{align*} p_k \end{align*}$ を以下のように修正します。 $\begin{align*} p_k(z) = F\left( \frac{F^{-1}(p_k^0) - \sqrt{\rho_k}z}{\sqrt{1-\rho_k}} \right)\end{align*}$
ここで、zは正規分布に従う乱数で、Fは正規分布の累積分布関数、ρは分布の広がり方を調整するパラメータ、 $\begin{align*} p_k^0 \end{align*}$ はもともとの $\begin{align*} p_k \end{align*}$ の値です。次の図はこの $\begin{align*} p_k(z) \end{align*}$ の例です。

図から分かるようにρを大きくするほど確率の揺らぎが大きくなります。GCI modelでの損失額のMonte Carlo計算は以下の通りになります。

$\begin{align*} z \end{align*}$ を正規分布に従って生成
$\begin{align*} p_k(z) \end{align*}$ を計算
その確率で損失額 $\begin{align*} \lambda_k \end{align*}$ を計上し集計

これを繰り返すことになります。

GCIモデルの量子計算

まず、K個の量子ビットで、1が測定される確率 $\begin{align*} p_k\end{align*}$ として資産の損失の有無を表現します。
$\begin{align*} \sqrt{1-p_k}|0\rangle + \sqrt{p_k} | 1 \rangle \end{align*}$
これは量子ビットをy軸回りに特定の角度回転することで得られます。
$\begin{align*} \bigotimes_{k=1}^K R_Y(\theta_p^k)|0\rangle_k = \bigotimes_{k=1}^K \left(\sqrt{1-p_k}|0\rangle + \sqrt{p_k}|1\rangle\right) \end{align*}$
更に、z での確率の修正し、かつ、zが正規分布での重みつけした状態を作ります。
$\begin{align*} |\Psi\rangle = \sum_{i=0}^{2^{n_z}-1} \sqrt{p_z^i} |z_i\rangle \bigotimes_{k=1}^K \left( \sqrt{1-p_k(z_i)}|0\rangle + \sqrt{p_k(z_i)}|1\rangle \right) \end{align*}$
ここで、 $\begin{align*} z\end{align*}$ は $\begin{align*} n_z\end{align*}$ 個の量子ビットを新たに追加して表現します。

さて、qiskit-finance の library に、実装されている Gaussian Conditional Independence Model を使うと以下のように作成できます。

ここでは簡単のため２つの資産( $\begin{align*} p_0 = 0.15, p_1 = 0.25\end{align*}$ )で、それぞれの発生しうる損失額( $\begin{align*} L_1 = 1, L_2 = 2 \end{align*}$ )で、ｚのために２つの量子ビットを準備します。
(実はこれらの数値は今回紹介する量子回路でうまく計算できるように非常に巧みに調整されています。)

# set problem parameters
n_z = 2
z_max = 2
z_values = np.linspace(-z_max, z_max, 2**n_z)
p_zeros = [0.15, 0.25]
rhos = [0.1, 0.05]
lgd = [1, 2]
K = len(p_zeros)
alpha = 0.05

from qiskit_finance.circuit.library import GaussianConditionalIndependenceModel as GCI
u = GCI(n_z, z_max, p_zeros, rhos)

u.darw('mpl')

この回路で、q0, q1 は資産の状態で、q2, q3 は z の状態を表しています。最初のP(X)で初期の確率を与えて、次の２箇所のLinRotで、z の分布に従った状態の更新を行います。この回路を実行すると４つの量子ビットの状態の重ね合わせ状態が得られます。以下のグラフは、量子ビットの可能な組み合わせ(0000, 0001, 0010, …, 1111) の重ね合わせを可視化したものです。

from qiskit.visualization import plot_state_city
job = execute(u, backend=Aer.get_backend('statevector_simulator'))
plot_state_city(job.result().get_statevector(), figsize=(18, 8))

シミュレーションでは、このような量子ビットの状態ベクトルを見ることができますが、
実際の測定では、以下のような確率で量子ビットが観測されます。(これは状態ベクトルシミュレーションの結果から計算してあります。)

この最初の２つの量子ビットが１の際に、対応する資産が損失を発生させます。例えば、損失額が１となる確率は、zの量子ビット(3rd, 4th)について足し合わせ、つまり1000, 1001, 1010, 1011 が発生する確率の合計として計算できます。各々の損失に対応した確率は以下の通りになります。

量子回路でのVaR計算(改良版)

さて、実は上の回路だけで目的は達成しているのですが、量子ビットを観測して件数を数え上げて古典的に損失を集計する。うーん、せっかく量子回路を使っているのですから全部量子回路上で計算を済ませたいですよね。そのためにGCI回路を拡張したいと思います。

そのためには、まず以下の回路が必要となります。

個別の資産の損失発生確率と損失額を集計
損失額が閾値以下となる確率を集計

そして、先に紹介した量子振幅推定を利用して計算できれば、量子ビットを何度も測定し確率を推定するよりも少ない測定回数で精度の良い結果が得られます。

量子回路での合計損失額評価

K量子ビットで、それぞれの資産の状態を表しましたが、
これに新たに $\begin{align*} n_s\end{align*}$ 個の量子ビットを追加して、そこに合計損失額を書き込みます
(zの量子ビットと同様に $\begin{align*} n_s\end{align*}$ ビットまでの数値を表現できます)。
$\begin{align*} \mathcal{S} : |x_1, \cdots, x_K\rangle_K |0\rangle_{n_s} \rightarrow |x_1, \cdots, x_K\rangle_K |\lambda_1 x_1 + \cdots + \lambda_K x_K \rangle_{n_s} \end{align*}$
この計算は、qiskit の WeightedAdder で実装されています。

from qiskit.circuit.library import WeightedAdder
agg = WeightedAdder(n_z + K, [0]*n_z + lgd)

この際には、zに対応する量子ビットの重さをゼロとして加えないよう注意します。

実際にこの回路で合計損失額( $\begin{align*} L \in \{ 0, \cdots, 2^{n_s}-1 \} \end{align*}$ )が計算できているか、測定して試してみたいと思います。そのために新しい量子ビットをさらに追加し、その振幅に変換します。 $\begin{align*} |L\rangle_{n_s}|0\rangle \rightarrow |L\rangle_{n_s} \left( \sqrt{1 - L/(2^{n_s}-1)}|0\rangle + \sqrt{L/(2^{n_s}-1)}|1\rangle \right) \end{align*}$
最後に追加した量子ビットの振幅や0,1の観測される確率を測定すれば、合計損失額は計算できます。この計算には、LinearAmplitudeFunctionを使います。

from qiskit.circuit.library import LinearAmplitudeFunction

# define linear objective function
breakpoints = [0]
slopes = [1]
offsets = [0]
f_min = 0
f_max = sum(lgd)
c_approx = 0.25

objective = LinearAmplitudeFunction(
    agg.num_sum_qubits,
    slope=slopes,
    offset=offsets,
    # max value that can be reached by the qubit register (will not always be reached)
    domain=(0, 2**agg.num_sum_qubits-1),
    image=(f_min, f_max),
    rescaling_factor=c_approx,
    breakpoints=breakpoints
)

先程作成したGCI model の回路の合計損失額、合計損失振幅を計算する回路は全体で以下のようになります。

from qiskit import QuantumRegister
# define the registers for convenience and readability
qr_state = QuantumRegister(u.num_qubits, 'state')
qr_sum = QuantumRegister(agg.num_sum_qubits, 'sum')
qr_carry = QuantumRegister(agg.num_carry_qubits, 'carry')
qr_obj = QuantumRegister(1, 'objective')

# define the circuit
state_preparation = QuantumCircuit(qr_state, qr_obj, qr_sum, qr_carry, name='A')

# load the random variable
state_preparation.append(u.to_gate(), qr_state)

# aggregate
state_preparation.append(agg.to_gate(), qr_state[:] + qr_sum[:] + qr_carry[:])

# linear objective function
state_preparation.append(objective.to_gate(), qr_sum[:] + qr_obj[:])

# uncompute aggregation
state_preparation.append(agg.to_gate().inverse(), qr_state[:] + qr_sum[:] + qr_carry[:])

# draw the circuit
state_preparation.draw('mpl')

この回路の最初の部分、P(X)が GCI model の部分、次の adder が、合計損失の計算、F が合計損失振幅の計算、最後に adder の逆操作(uncompute)を行います。
(uncomputeとは、計算の補助に使った量子ビット(ancilla)を元に戻す操作にあたります)

この回路を実行して、objective bit が１となる確率を調べてみると

prob_obj = 0
for i, a in enumerate(job.result().get_statevector()):
b = ('{0:0%sb}' % (len(qr_state) + 1)).format(i)[-(len(qr_state) + 1):]
am = np.round(np.real(a), decimals=4)
if np.abs(am) > 1e-6 and b[0] == '1':
prob_obj += am**2

LinearAmplitudeFunctionは、実際には元の数値を規格化した結果を出力しているので、実際の結果と比較するためには補正が必要です。(https://qiskit.org/documentation/stubs/qiskit.circuit.library.LinearAmplitudeFunction.html)
以下のように確かに期待値が計算できていることが判ります。

print('Exact Expected Loss: %.4f' % expected_loss)
print('Mapped Operator value: %.4f' % objective.post_processing(prob_obj))

# Exact Expected Loss: 0.6409
# Mapped Operator value: 0.6640

さて、上記のやり方では、objective bit が１となる確率を測定する必要があるので、結局はMonte Carlo計算に帰着してしまいます。Qiskitに実装されている量子位相推定のアルゴリズムを使えば以下のように計算ができます。

# set target precision and confidence level
epsilon = 0.01
alpha = 0.05

qi = QuantumInstance(Aer.get_backend('qasm_simulator'), shots=100)
problem = EstimationProblem(state_preparation=state_preparation,
objective_qubits=[len(qr_state)],
post_processing=objective.post_processing)
# construct amplitude estimation
ae = IterativeAmplitudeEstimation(epsilon, alpha=alpha, quantum_instance=qi)
result = ae.estimate(problem)

# print results
conf_int = np.array(result.confidence_interval_processed)
print('Exact value: \t%.4f' % expected_loss)
print('Estimated value:\t%.4f' % result.estimation_processed)
print('Confidence interval: \t[%.4f, %.4f]' % tuple(conf_int))

# Exact value: 0.6409
# Estimated value: 0.7011
# Confidence interval: [0.6322, 0.7699]

量子回路での累積分布関数

さて、total loss がある閾値よりも小さい確率を計算したいと思います。
$\begin{align*} P[L \leq x] \end{align*}$
そのために以下のような演算を行います。
$\begin{align*} C: |L\rangle_n | 0 \rangle \rightarrow \begin{cases} |L\rangle_n | 1 \rangle \quad \text{if} \ L \leq x \\ |L\rangle_n | 0 \rangle \quad \text{if} \ L > x \\ \end{cases} \end{align*}$
この計算を様々な数値Lの重ね合わせ状態に $\begin{align*} \sum_{L=0}^x \sqrt{p_L} |L\rangle \end{align*}$ 作用させると $\begin{align*} \sum_{L=0}^{x} \sqrt{p_L} |L\rangle_{n_s}|1\rangle + \sum_{L=x+1}^{n_S-1} \sqrt{p_L} | L \rangle_{n_s}0\rangle \end{align*}$
(Qiskit tutorial で紹介されている式では第２項の読み取りビットが1になってますが、０のtypoのようです)
この回路はIntegerComparatorで計算できます。

from qiskit.circuit.library import IntegerComparator

x_eval = 2
comparator = IntegerComparator(agg.num_sum_qubits, x_eval + 1, geq=False)
comparator.draw('mpl')

閾値が２の場合の回路は上図のようになります。通常のプログラミングと異なり量子計算回路ではロジックが回路構造に反映されていますので、閾値に応じて回路の接続が変わります。入力の数値は state0, 1 の量子ビットで２進数で符号化され、compare の量子ビットが閾値以下だと 1 を出力します。

具体的に3( $\begin{align*} |11\rangle \end{align*}$ )入力された場合、計算の流れは以下の通りです。

補助量子ビット(a)が制御ノット(CNOT)で、制御ビット(state 0)が１なので反転して１
compare ビットがCCNOTで制御ビット(state 1, 補助ビット)が１なので反転し１
compare ビットがXゲートで反転し0が出力される
補助ビット(a)が最後にCNOTで初期化される

一方、compare ビットが反転するには、CCNOTの制御ビットに両方１が入力される必要があるため0, 1, 2 の入力では、Xゲートに0がそのまま入力し、最終的な出力は１であることがわかります。

def get_cdf_circuit(x_eval):
    # define the registers for convenience and readability
    qr_state = QuantumRegister(u.num_qubits, 'state')
    qr_sum = QuantumRegister(agg.num_sum_qubits, 'sum')
    qr_carry = QuantumRegister(agg.num_carry_qubits, 'carry')
    qr_obj = QuantumRegister(1, 'objective')
    qr_compare = QuantumRegister(1, 'compare')

    # define the circuit
    state_preparation = QuantumCircuit(qr_state, qr_obj, qr_sum, qr_carry, name='A')

    # load the random variable
    state_preparation.append(u, qr_state)

    # aggregate
    state_preparation.append(agg, qr_state[:] + qr_sum[:] + qr_carry[:])

    # comparator objective function
    comparator = IntegerComparator(agg.num_sum_qubits, x_eval + 1, geq=False)
    state_preparation.append(comparator, qr_sum[:] + qr_obj[:] + qr_carry[:])

    # uncompute aggregation
    state_preparation.append(agg.inverse(), qr_state[:] + qr_sum[:] + qr_carry[:])

    return state_preparation

state_preparation = get_cdf_circuit(x_eval=2)

最終的な回路は以上となります。先程 total loss を計算していた部分がこの比較器に変更となっています。今回、この objective bit が１となる確率(振幅)が求めたいCDFとなります。

# evaluate resulting statevector
var_prob = 0
for i, a in enumerate(job.result().get_statevector()):
    b = ('{0:0%sb}' % (len(qr_state) + 1)).format(i)[-(len(qr_state) + 1):]
    prob = np.abs(a)**2
    if prob > 1e-6 and b[0] == '1':
        var_prob += prob
print('Operator CDF(%s)' % x_eval + ' = %.4f' % var_prob)
# Operator CDF(2) = 0.9591

量子回路でのVaR評価

さて、最後にCDFを計算してVaRを求めます。今回の模型のパラメータで可能な損失額ごとにCDFを計算してみます。

def run_ae_for_cdf(x_eval, epsilon=0.01, alpha=0.05, simulator='qasm_simulator'):
    # construct amplitude estimation
    state_preparation = get_cdf_circuit(x_eval)
    qi = QuantumInstance(Aer.get_backend('qasm_simulator'), shots=100)
    problem = EstimationProblem(state_preparation=state_preparation,
                                objective_qubits=[len(qr_state)])
    ae_var = IterativeAmplitudeEstimation(epsilon, alpha=alpha, quantum_instance=qi)
    result_var = ae_var.estimate(problem)

    return result_var.estimation

losses = np.array([0,1,2,3])
cdfs = np.array([run_ae_for_cdf(x) for x in losses])

結果は以下の通りです。ほぼ95％(より厳密には95.86%)で、今回作成したポートフォリオの損失は2以下で収まることが判りました。

今回の場合は可能な損失額全ての累積分布関数を計算しました。実際の問題では、目標となるαの確率となる損失額を二分法で探索すれば効率良くVaRを評価できます。

まとめ

最後に今回の計算は特別な少数のパラメータしかない模型でしたが、より現実的な系について考えようと思います。約100万種類( $\begin{align*} K = 2^{20}\end{align*}$ 量子ビット)の金融資産を考え、そのVaRの計算に必要な量子ビット・計算のコストは次のようになります。まず、正規分布の乱数( $\begin{align*} z\end{align*}$ )として $\begin{align*} n_z=10\end{align*}$ とすれば1024点刻みのサンプル、合計損益を集計する量子ビットとして十分な領域を表現するため $\begin{align*} n_s=30\end{align*}$ 、量子振幅推定のために $\begin{align*} m=10\end{align*}$ の量子ビットが必要とします。それにより約100万銘柄からなるポートフォリオのVaR( $\begin{align*} \alpha=99.9\end{align*}$ )を0.06% ポイントの精度で計算することができます。計算の実行で必要なステップの深さとしては3億7千万のゲート演算が必要になりますが１つ実行が $\begin{align*} 10^{-4}\end{align*}$ 秒ですむなら１時間程度で計算できます。

さて、今後の量子コンピュータの発展の見込みとして、IBMは昨年に65量子ビットのチップを発表しており、2021年中には127量子ビット、更に2023年には1,121量子ビットまで到達することを計画しています。効率的に素因数分解を実行できるShorのアルゴリズムのような汎用的な量子コンピュータの登場までには、まだまだ遠いです(2048 bit RSA暗号の解読には、2千万量子ビットが必要な見込みです)。しかし、数十から数百個の誤り補正がないNISQ(Noisy Intermediate Scale Quantum)デバイスと今回紹介したようなアルゴリズムは、今後様々な分野で活用が期待されます。

次世システム研究室では、ビッグデータ解析プラットホームの設計・開発を行うアーキテクトとデータサイエンティストを募集しています。興味を持って頂ける方がいらっしゃいましたら、ぜひ募集職種一覧からご応募をお願いします。一緒に勉強しながら楽しく働きたい方のご応募をお待ちしております。

参考資料

Qiskit Finance (https://qiskit.org/documentation/finance/)
Quantum Amplitude Amplification (https://en.wikipedia.org/wiki/Amplitude_amplification)
S. Woerner and D. J. Egger, “Quantum Risk Analysis” [arXiv:1806.06893] (https://arxiv.org/abs/1806.06893)
D. J. Eggger, R. G. Guitierrez, J. C. Mestre, and S. Woerner, “Credit Risk Analysis using Quantum Computers” [arXiv:1907.03044] (https://arxiv.org/abs/1907.03044)
D. J. Egger, C. Gambella, J. Marecek, S. McFaddin, M. Mevissen, R. Raymond, A. Simonetto, S. Woerner, E. Yndurain, “Quantum Computing for Finance: State of the Art and Future Prospects”, [arXiv:2006.14510] (https://arxiv.org/abs/2006.14510)
IBM’s Roadmap For Scaling Quantum Technology, https://www.ibm.com/blogs/research/2020/09/ibm-quantum-roadmap/
C. Gidney and M. Ekera, “How to factor 2048 bit RSA integers in 8 hours using 20 million noisy qubits”, [arXiv:1905.09749] (https://arxiv.org/abs/1905.09749)

量子コンピュータと金融
~ポートフォリオのリスクを計算してみる~

TL;DR

導入

ポートフォリオとリスク分散

VaR (Value at Risk)

Monte Carlo法

量子ビットと量子回路

量子回路の測定

量子振幅推定

量子回路でのVaR計算

ポートフォリオと損失額

Gaussian conditional independence (GCI) model

GCIモデルの量子計算

量子回路でのVaR計算(改良版)

量子回路での合計損失額評価

量子回路での累積分布関数

量子回路でのVaR評価

まとめ

参考資料

関連記事

Zero-shot learningの紹介：見たことがない画像やニュースを予測してみました

Seaborn Objects
~ グラフィックの文法で強化された Python 可視化ライブラリの新形態 ~

FPGAに機械学習モデルを実装する –
その1：ランダムフォレストによるクラス分類

FX予測においてON-LSTMとGANはLSTMに勝てるのか？

量子コンピュータと金融~ポートフォリオのリスクを計算してみる~

TL;DR

導入

ポートフォリオとリスク分散

VaR (Value at Risk)

Monte Carlo法

量子ビットと量子回路

量子回路の測定

量子振幅推定

量子回路でのVaR計算

ポートフォリオと損失額

Gaussian conditional independence (GCI) model

GCIモデルの量子計算

量子回路でのVaR計算(改良版)

量子回路での合計損失額評価

量子回路での累積分布関数

量子回路でのVaR評価

まとめ

参考資料

関連記事

Zero-shot learningの紹介：見たことがない画像やニュースを予測してみました

Seaborn Objects ~ グラフィックの文法で強化された Python 可視化ライブラリの新形態 ~

FPGAに機械学習モデルを実装する – その1：ランダムフォレストによるクラス分類

FX予測においてON-LSTMとGANはLSTMに勝てるのか？

量子コンピュータと金融
~ポートフォリオのリスクを計算してみる~

Seaborn Objects
~ グラフィックの文法で強化された Python 可視化ライブラリの新形態 ~

FPGAに機械学習モデルを実装する –
その1：ランダムフォレストによるクラス分類