#はじめに
某社でクオンツアナリストとして働いておりますhalと申します。
本記事は自分の勉強も兼ねて、金融実務者の方向けに典型的なポートフォリオ最適化を通じて、
"量子アニーリング"の紹介&デモンストレーションをしていきます。
pythonの導入されたPCがあればどなたでも無料で試せるので、興味がある方は実際に動かしてみてください。
本記事ではD-Wave社のLeapという無料の量子アニーリングサービスを使用してきます。
ことはじめ1で導入方法を紹介しているのでこちらもご覧ください
#ポートフォリオ最適化
では実際にポートフォリオ最適化を考えていきます。
典型的な効用関数最大化問題を考えます。最小化問題形式にしているので注意してください。
(文字の定義は一番下に記載しています)
\begin{align}
\mathrm{minimize}~ ~ &\mathcal{H}=-\sum_i r_i w_i +\gamma \sum_{i,j} \sigma_{i,j} w_i w_j=-\boldsymbol{ r}^\mathrm{T}\cdot \boldsymbol{w}+\gamma \boldsymbol{w}^\mathrm{T}\boldsymbol{C}\boldsymbol{w}\\
\mathrm{subject ~ to}~~&||\boldsymbol{w}||=1
\end{align}
第1項がポートフォリオのリターンを表し、第2項はポートフォリオのリスクに、リスク回避度$\gamma$を掛け合わせたものになっており、
リスクをある程度制限しながらリターンを最大化する最適化です。
このままではLeapで計算をさせることはできず、QUBO形式もしくは Ising形式という形に問題を変換する必要があります。
どちらの形式でも本質的には同じものですが、今回はQUBO形式に変換していきます。
#QUBO
QUBO(Quadratic Unconstrained Binary Optimization)形式とは$ q_i$を0 or 1 の2値変数として
\mathrm{minimize}~ ~ \mathcal{H}=-\sum_{i,j}Q_{i,j}q_iq_j=-\boldsymbol{q}^\mathrm{T}\boldsymbol{Q}\boldsymbol{q}\\
の形式で表示される問題です。Unconstrained なので制約条件はつけられませんがラグランジュの未定乗数法で対応できますね。
D-Waveで解ける問題はこのQUBO形式なので、問題を変換して行列$\boldsymbol{Q}$を求める必要があります。
変換する際に一番の問題点は変数がbinaryであり0,1でしか答えが得られないことです。
任意のウェイトを得るには工夫をする必要があるので、今回は単純化した例題を解いてみます。
#等ウェイトでの最適ポートフォリオ
N種類の資産を投資対象銘柄として、そのうちK種類の資産での最適ポートフォリオを求める。
ただしポートフォリオは等ウェイトとする。
##問題のQUBO形式変換
さて等ウェイトであればBinaryでの表示は比較的簡単です。
資産iに対応するビットを$q_i$とし、$q_i=1$であればポートフォリオに採用、$q_i=0$であれば不採用と定義します。
この時最適化の制約条件は
\sum_i^N q_i=K
となります。ラグランジュの未定乗数$\lambda$を用いて制約条件なしの最小化問題に変換すると
\mathrm{minimize}~ ~ \mathcal{H}=-\boldsymbol{ r}^\mathrm{T}\cdot \boldsymbol{q}+
\gamma \boldsymbol{q}^\mathrm{T}\boldsymbol{C}\boldsymbol{q}+
\lambda \left(K-\sum_i^N q_i\right)^2
となります。これを変換し行列$\boldsymbol{Q}$を求め、QUBO形式に変換してきましょう。
以下の分配法則を利用してそれぞれの項ごとに変換していきます。
\boldsymbol{q}^T \boldsymbol{A}\boldsymbol{q}+\boldsymbol{q}^T \boldsymbol{B}\boldsymbol{q}=\boldsymbol{q}^T (\boldsymbol{A+B})\boldsymbol{q}
###第1項(リターン)
$q_i=q_i\times q_i$である($q_i$は0か1なので)ことを利用して以下のように書き直せます
(2019/5/5 後の記事と整合させるため表記を少し変更しました)
\begin{aligned}
\boldsymbol{ r}^\mathrm{T}\cdot \boldsymbol{q}&=\boldsymbol{q}^\mathrm{T}\mathrm{diag}(\boldsymbol{r^\mathrm{T}})\boldsymbol{q}\\ \\
\mathrm{diag}(\boldsymbol{r}^\mathrm{T})&=\left( \begin{array}{cccc}
r_1 & 0 & \cdots & 0\\
0 & r_2 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & r_N
\end{array} \right)
\end{aligned}
###第2項（リスク)
これはそのままでいけますね.
###第3項(制約条件)
第1項と同様$q_i=q_i\times q_i$であることと、定数は無視できることを利用して以下のように書き直せます。
（最小化問題なので定数項は無視できる)
(2019/5/5 第1項と同様に表記を少し変更しました)
\begin{align}
\left(K-\sum_i^N q_i\right)^2 &=q_1^2+q_2^2+\cdots+q_N^2
+2q_1q_2+\cdots 2q_{N-1}q_N
-2Kq_1-2Kq_2-\cdots -2Kq_N
+K^2\\
&=(1-2K)q_1^2+(1-2K)q_2^2+\cdots+(1-2K)q_N^2
+2q_1q_2+\cdots 2q_{N-1}q_N\\
&=\boldsymbol{q}^\mathrm{T}\left[
\underbrace{\left( \begin{array}{cccc}
1 & 1 & \cdots & 1\\
1 & 1 & \cdots & 1 \\
\vdots & \vdots & \ddots & \vdots \\
1 & 1 & \cdots & 1
\end{array} \right)}_\boldsymbol{A}-2K\boldsymbol{E}
\right]\boldsymbol{q}\\
&=\boldsymbol{q}^\mathrm{T}(\boldsymbol{A}-2K\boldsymbol{E})\boldsymbol{q}
\end{align}
###QUBO化
さあこれで全部QUBO形式に変換できました。以上の式をまとめると最適化問題は
\mathcal{H}=\boldsymbol{q}^\mathrm{T} \boldsymbol{Q} \boldsymbol{q}\\
\boldsymbol{Q}=-\mathrm{diag}(\boldsymbol{r^\mathrm{T}})+\gamma\boldsymbol{C}+\lambda (\boldsymbol{A}-2K\boldsymbol{E})
となります。あとはこれをLeapで計算してもらうだけです。
(2019/5/6追記)ただし$\boldsymbol{Q}$は最終的に変数同士の"組み合わせ"で定義する必要があるので、得られた$\boldsymbol{Q}$の非対角成分を2倍し、上三角行列だけにしてあげる必要があります。
##実際に解いてみる
$N=5,K=3$すなわち5資産のなかから3資産選ぶ場合を考えます。
分析データは以下のように定義します。（適当です)
ただし、共分散行列は上三角部分だけにしています
（2019/5/2追記）正確には共分散行列の非体対角成分だけ2倍して上三角行列として表示しています
(2019/5/6 追記) 次回の記事と整合性を取るために共分散行列を通常の対角行列として表示します。問題としては同値です。
\begin{align}
\boldsymbol{r}&=(1\%,2\%,3\%,4\%,5\%)\\
\boldsymbol{C}&=\left( \begin{array}{ccccc}
0.01 & -0.0062 & 0.0238 & 0.0209 & 0.0096 \\
-0.0062 & 0.0225 & 0.01455 & 0.0264 & 0.0147 \\
0.0238 & 0.01455 & 0.04 & 0.0114 & 0.0049 \\
0.0209 & 0.0264 & 0.0114 & 0.0625 & 0.0057 \\
0.0096 & 0.0147 & 0.0049 & 0.0057 & 0.09
\end{array} \right)
\end{align}
また$\gamma=1.2,\lambda =1$ とします。
この問題をLeapで最適化を行うコードは以下のようになります。
(2019/5/6 追記) ここも共分散行列を対角行列として表示に、最後に非対角項を2倍し上三角行列に変換する処理を追加しました。問題としては同値です。
import numpy as np
from dwave.system.samplers import DWaveSampler
from dwave.system.composites import EmbeddingComposite
#リターンベクトル
r=np.array([0.01,0.02,0.03,0.04,0.05])
#共分散行列
C=np.array([[ 0.01 , -0.0062 , 0.0238 , 0.0209 , 0.0096 ],\
[-0.0062 , 0.0225 , 0.01455, 0.0264 , 0.0147 ],\
[ 0.0238 , 0.01455, 0.04 , 0.0114 , 0.0049 ],\
[ 0.0209 , 0.0264 , 0.0114 , 0.0625 , 0.0057 ],\
[ 0.0096 , 0.0147 , 0.0049 , 0.0057 , 0.09 ]])
#N=5資産中K=3資産の等ウェイトポートフォリオ
N=5
K=3
#ガンマとラムダを定義
gamma=1.2
lam=1
#QUBO matrix を計算
R=np.diag(r)
A=np.ones([N,N])-2*K*np.eye(N)
Q=-R+gamma*C+lam*A
Q=np.triu(Q)+np.triu(Q,k=1)#非対角成分を2倍して上三角行列に変換
Q=Q/np.mean(Q) #規格化しないと桁落ちにより正しい解が得られないことがあります
#問題を与える時は辞書にして渡す必要があるので変換
Qdict=dict()
for i in range(N):
for j in range(N):
if i<=j:
Qdict.update({("x"+str(i),"x"+str(j)):Q[i,j]})
#問題をapiに投げ、結果を取得
sampler = EmbeddingComposite(DWaveSampler())
response = sampler.sample_qubo(Qdict,num_reads=1000)
print(response)
実行結果は以下のようになります(num_oc.は毎回異なります)
energyは目的関数の値を示しているので、一番上のものが最適解となります。
x0 x1 x2 x3 x4 energy num_oc. chain_.
0 1 1 1 0 0 -14.420885 180 0.0
1 1 1 0 0 1 -14.410704 60 0.0
2 1 1 0 1 0 -14.358505 144 0.0
3 0 0 1 1 1 -14.324203 130 0.0
4 0 1 1 0 1 -14.322323 85 0.0
5 1 0 1 0 1 -14.314282 33 0.0
...
実行結果から最適解は第0、1、2銘柄の等ウェイトポートフォリオであることがわかりました。
(通常の最適化でも同様な結果が得られますので、確認してみてください)
#任意のウェイトでの最適ポートフォリオ
ことはじめ3に続く...
#変数定義
| 変数 | 内容 |
|---|---|
| $\boldsymbol{w}=\{w_i\}$ | 対象銘柄のWeight |
| $\boldsymbol{r=\{r_i\}}$ | 対象銘柄の期待リターン（ベクトル) |
| $\boldsymbol{C}=\{ \sigma_{i,j}\}$ | 対象銘柄の共分散行列 |
#リンク
金融のための量子コンピューターことはじめ1 | D-wave Leap導入編
金融のための量子コンピューターことはじめ2 | ポートフォリオ最適化 1（本記事)
金融のための量子コンピューターことはじめ3 | ポートフォリオ最適化 2