QiskitのHHLアルゴリズムを使ってポートフォリオ最適化問題を実装してみる

金融分野への量子コンピュータの応用先の1つとして挙げられているポートフォリオ最適化問題をQiskitを使って実装してみます。

スポンサーリンク

はじめに

本記事はQuantum Native Dojo 7-3「HHLアルゴリズムを用いたポートフォリオ最適化」を元に執筆しています。詳細はこちらをご覧ください。

7-3. HHLアルゴリズムを用いたポートフォリオ最適化 — Quantum Native Dojo ドキュメント

ポートフォリオ最適化とは

資産投資をする際に株式など将来の金額が不確実である資産へ投資する場合、ポートフォリオを組んで複数の資産へ分散投資するのが一般的です。

個人投資家にも言えることですが、金融機関のように巨大な金額のお金を運用する機関投資家にとってリスクを最小限に抑えた上でリターンを得るための最適なポートフォリオを組むことが課題となっています。

どのように最適ポートフォリオを得るかは現代ポートフォリオ理論という理論的枠組が提唱されています。ロボアドバイザーのWealthNaviもこの理論に基づいて資産運用のモデルを組んでいます。

ポートフォリオ最適化の定式化

現代ポートフォリオ理論では投資家が望む期待値リターンを設定したときに、それを達成するポートフォリオの中で収益率の分散が、つまりリスクが最小となるものを求めます。定式化すると以下のように表すことができます。

$$\begin{aligned} \min &~\bm{w} \Sigma \bm{w} \\ s.t. &~\bm{R^T} \bm{w} =\mu \\ &~\bm{\Pi^T} \bm{w} = \zeta \end{aligned}$$

各記号は次のとおりです。

$\bm{w}$ポートフォリオの重み(各資産への割当)ベクトル
$\Sigma$資産同士の共分散行列
$\bm{R}$資産ごとの期待リターンベクトル
$\mu$トータル期待リターン
$\bm{\Pi}$資産ごとの価格ベクトル
$\zeta$投資家の投資予算

Quantum Native Dojoでは問題を簡単にするため投資家の予算は考慮せずに、ポートフォリオへの投資する割合が合計で1となるような制限にしています。つまり数式で表すと下のようになります。

$$\begin{aligned} \min &~\bm{w} \Sigma \bm{w} \\ s.t. &~\bm{R^T} \bm{w} =\mu \\ &~\bm{I} \bm{w} = 1 \end{aligned}$$

このような等式制約付きの最大最小問題はラグランジュの未定乗数法を利用することで解くことができます。ラグランジュ乗数$\eta, \theta$を用いると、次のような方程式で表すことができます。

$$\left( \begin{array}{ccc} 0 & 0 & \bm{R^T} \cr 0 & 0 & \bm{I^T} \cr \bm{R^T} & \bm{I^T} & \Sigma \end{array} \right) \left( \begin{array}{c} \eta \cr \theta \cr \bm{w^T} \end{array} \right) = \left( \begin{array}{c} \mu \cr \zeta \cr \bm{0} \end{array} \right) $$

よって求めたい最適$\bm{w}$がわかればいいので、表現行列$W$の逆行列が求まればよいということになります。そのような問題にはHHLアルゴリズムが有効に働きます。

$$ W = \left( \begin{array}{ccc} 0 & 0 & \bm{R^T} \\ 0 & 0 & \bm{I^T} \\ \bm{R^T} & \bm{I^T} & \Sigma \end{array} \right) $$

それではQiskitを使って実装してみましょう。

Qiskitで実装してみる

では実際にQiskitを使って実装してみましょう。題材はNative Dojoと変えずに2017年のGAFAの株価時系列データを使ってポートフォリオを最適化してみます。

時系列データを入手する

まずはGAFA4社の時系列データを入手します。どんな方法でもいいですが、一例としてpandas_datareaderを使ってデータを取得します。

# 銘柄選択
codes = ["GOOG", 'AAPL', 'FB', 'AMZN'] # GAFA

# 2017年の1年間のデータを使用
start = datetime.datetime(2017, 1, 1)
end = datetime.datetime(2017, 12, 31)

# Yahoo! Financeから日次の株価データを取得
data = web.DataReader(codes, 'yahoo', start, end)

df = data['Adj Close']

データの前処理

先程のラグランジュ方程式に必要な量を求めていきます。

# 期待値リターン
daily_return = df.pct_change()
expected_return = daily_return.dropna(how='all').mean() * 252 # 年率換算のため年間の営業日数252を掛ける
# 共分散行列
cov = daily_return.dropna(how='all').cov() * 252 # 年率換算のため

これらの量を使って表現行列$W$を定義します。

R = expected_return.values
Pi = np.ones(4)
S = cov.values

row1 = np.append(np.zeros(2), R).reshape(1,-1)
row2 = np.append(np.zeros(2), Pi).reshape(1,-1)
row3 = np.concatenate([R.reshape(-1,1), Pi.reshape(-1,1), S], axis=1)
W = np.concatenate([row1, row2, row3])

np.set_printoptions(linewidth=200)

右辺のベクトルを定義します。ここでハイパーパラメータであるトータル期待リターン$\mu$は10%に設定します。

mu = 0.1 # ポートフォリオのリターン
xi = 1.0 
# 右辺ベクトル
mu_xi_0 = np.append(np.array([mu, xi]), np.zeros_like(R))

ここで今用意したベクトルは6次元のベクトルですが、量子ビットで扱うために行列の次元を$2^n$に拡張する必要があります。QiskitではHHLモジュールの.matrix_resize()メソッドを利用することで簡単に拡張することができます。

matrix, vector, truncate_powerdim, truncate_hermitian = HHL.matrix_resize(W, mu_xi_0)

HHLアルゴリズムの準備

それではHHLアルゴリズムを組み立てていきます。といってもQiskitにはHHLアルゴリズムのもモジュールが用意されているのでこちらを利用します。

処理プロセスである固有値探索と制御回転の条件を設定します。今回固有値探索に用いる量子位相推定には補助量子ビットを7量子ビット利用します。

def create_eigs(matrix, num_ancillae, num_time_slices, negative_evals):
    ne_qfts = [None, None]
    if negative_evals:
        num_ancillae += 1
        # The QFT and IQFT circuits for handling negative eigenvalues
        ne_qfts = [QFT(num_ancillae - 1), QFT(num_ancillae - 1).inverse()]

    """
    Specifically, this class is based on PhaseEstimationCircuit with no measurements 
    and has additional handling of negative eigenvalues
    """
    return EigsQPE(MatrixOperator(matrix=matrix), # The Hamiltonian Operator object
                   QFT(num_ancillae).inverse(),
                   num_time_slices=num_time_slices,
                   num_ancillae=num_ancillae,
                   expansion_mode='trotter', 
                   expansion_order=2, # The suzuki expansion order, has a minimum value of 1
                   evo_time=None,  # This is t, can set to: np.pi*3/4
                   negative_evals=negative_evals, # Set True to indicate negative eigenvalues need to be handled
                   ne_qfts=ne_qfts)

# 扱う問題の次元
orig_size = len(vector)

# 固有値探索の初期化
eigs = create_eigs(matrix, 7, 50, True)
num_q, num_a = eigs.get_register_sizes()

# 制御回転の初期化
c = 0.5*(2 * np.pi * (1. / 2**(num_a) ))
reci = LookupRotation(scale=c, negative_evals=eigs._negative_evals, evo_time=eigs._evo_time)

最後に初期状態、つまり右辺ベクトルをセットしてHHLアルゴリズムのインスタンスを生成します。

# 初期状態の設定
init_state = Custom(num_q, state_vector=vector)

algo = HHL(
    matrix,
    vector, 
    truncate_powerdim, 
    truncate_hermitian, 
    eigs, # The eigenvalue estimation instance
    init_state, # The initial quantum state preparation
    reci, # The eigenvalue reciprocal and controlled rotation instance
    num_q, # Number of qubits required for the matrix Operator instance
    num_a, # Number of ancillary qubits for Eigenvalues instance
    orig_size # The original dimension of the problem (if truncate_powerdim)
    )

HHLアルゴリズムの実行

それでは組み立てたアルゴリズムを実行してみましょう。比較のためNumPyを用いて計算した厳密解と比較してみます。

## 結果
x_HHL = algo.run(QuantumInstance(Aer.get_backend('statevector_simulator')))
## 厳密解
# x_exact = np.linalg.lstsq(matrix, vector, rcond=0)[0]
x_exact = NumPyLSsolver(matrix, vector).run()

結果から最適ポートフォリオの重みベクトルを取り出した結果がこちらのグラフです。青が厳密解、赤がQiskitで計算した結果です。プラスは買い、マイナスは空売りを表していて概ねよく結果を再現できていると言えるでしょう。

最適ポートフォリオ

ためしに解の信頼度を求めてみると97%の信頼度と求まりました。この信頼度は補助ビットの数を変化させることで上下しますが、私のPC環境では補助ビット8以上では計算が終わりませんでした。。

# 解の信頼度を求める
def fidelity(hhl, ref):
    solution_hhl_normed = hhl / np.linalg.norm(hhl)
    solution_ref_normed = ref / np.linalg.norm(ref)
    fidelity = state_fidelity(solution_hhl_normed, solution_ref_normed)
    print("Fidelity:\t\t %f" % fidelity)

fidelity(x_HHL.solution, x_exact.solution) # => Fidelity: 0.968245

ちなみに今回組み立てた回路の情報はこちらです。使用した総量子ビットは14、量子ゲートは1252個、そのうちCNOTゲートは404個と、この程度の問題でも実機で動かすにはまだまだ大きな回路となってしまいました。

circuit_width:	 14
circuit_depth:	 1252
CNOT gates:	 404

コメント

タイトルとURLをコピーしました