Photo by Ariana Suárez on Unsplash
本稿は量子コンピューター Advent Calendar 2019 の記事です.
みなさん,資産運用してますか?
現金と違って,株などを持っていると常に価格変動のリスクに晒されます.
出来るだけリスクを抑えたい,でも収益も欲しい!このバランスは様々です.
どうやってそのリスクとリターンをコントロールするのか.
金融資産以外で考えてみるともっと身近なものに保険 という商品があります.例えば生命保険は保険料を払うことで将来の死亡や病気の(金銭的な)リスクをコントロールする商品と見ることもできます.
同じように,デリバティブ と呼ばれる価格変動リスクをコントロールするための金融商品が存在します.このような金融商品の価格はいくらが適正か?というのがデリバティブ価格決定問題 です.
デリバティブについての標準的な教科書としてはHull本 が挙げられます.
デリバティブの価格決定は,問題によっては非常に大量の計算量を必要とします.量子計算をうまく用いることで価格決定に必要な計算量を削減できるので嬉しい,という構造になっています.
そのため,本稿では量子計算を用いてデリバティブの価格を決定していこうと思います.紹介する量子アルゴリズムは に基づいています.
本稿では
デリバティブ
ヨーロピアンオプションの価格決定問題
モンテカルロ法
量子アルゴリズム
の順で紹介していきます.デリバティブの価格決定問題に詳しい方は4. 量子アルゴリズム
から読み始めることをおすすめします.
1. デリバティブ
デリバティブについてもう少し詳しく見ていきましょう.
デリバティブ(derivative)とは金融"派生"商品とも呼ばれます.
金融派生商品は「ある資産(例えば株)の価格を参照してその価値が決まる金融商品」という意味で現資産から派生しています.
これだと抽象的すぎてよくわからないので具体例をあげてみます.
「トヨタ株のヨーロピアンプットオプション」というデリバティブは予め決められた期日(満期という.例えば6ヶ月後の第二金曜日) に予め決められた価格(権利行使価格という.例えば8,000円) でトヨタ株を売る権利(プットオプションという.買う権利はコールオプション) です.
ここで,例えば6ヶ月後のトヨタ株が7,500円であれば500円もらえます.
これは,8,000円で売る権利を持っているので市場で7,500円でトヨタ株を買ってきて即8,000円で売れば500円儲かるという話です
一方,8,500円であれば0円です.
なぜかというと売る権利 なので行使しなくても良いからです.
オプションと呼ばれるデリバティブはこのような権利 の取引となっています.
オプションにも色々な種類があり,株価を参照する以外にも金利や債券の価格を参照したものや,権利行使のタイミングが複数存在するもの,払い戻しの条件がもっと複雑なものが存在します.
このような金融商品の価格を決定することは自明な問題ではなさそうです.
本稿ではデリバティブの中でも特にヨーロピアンオプションと呼ばれるオプションの価格決定問題を扱います.他のデリバティブとしては先物などがあり,またオプションの中でもアメリカンオプション等の興味深い問題もありますが,本稿では扱いません.
2. ヨーロピアンオプションの価格決定問題
オプションの中でも最も簡単なヨーロピアンオプションについてその価格決定を考えてみましょう.
ヨーロピアンオプションとは予め決められた期日(満期T T T )において権利行使可能なオプション の事を言います.
権利行使のタイミングが異なる他のオプションとしてはアメリカンオプションと呼ばれるものが存在し,アメリカンオプションは満期までの任意のタイミングにおいて権利行使可能なオプション のことを言います.
ヨーロピアンオプションとアメリカンオプションの中間の性質を持ったバミューダオプションは満期までの複数のあるタイミングにおいて権利行使可能なオプション のことを言います.
オプションの価格C C C は何に依存するかを考えてみると
権利行使価格K ( > 0 ) K(>0) K ( > 0 )
満期T ( > 0 ) T(>0) T ( > 0 )
金利R ( t ) , t ∈ [ 0 , T ] R(t), t \in [0, T] R ( t ) , t ∈ [ 0 , T ]
資産の(将来の)価格変動S ( t ) , t ∈ [ 0 , T ] S(t), t \in [0, T] S ( t ) , t ∈ [ 0 , T ]
が挙げられます.
ここで,K K K とT T T はパラメータです.
そうすると,金利と資産の(将来の)価格変動を考えていけば良さそうです.
なぜ金利に依存するかは直感的にはよくわからないと思うので,まずは金利について考えていきましょう.本稿では簡単のために金利が固定された世界R ( t ) = r [ 1 / y e a r ] R(t)=r[1/\mathrm{year}] R ( t ) = r [ 1 / y e a r ] を考えます.
オプションと金利
この世界では例えば100円を金利0.01で貸すと1年後に101円で返ってきます.
逆に考えると,一年後の101円の価値は現在の価値で100円となります.
同じように株の価値についても,1年後のメルカリ株の価値が5000円だった場合(!),現在の価値は5000 / 1.01 ≅ 4950 5000/1.01\cong4950 5 0 0 0 / 1 . 0 1 ≅ 4 9 5 0 円です.このように,将来時点での価値を現在時点での価値に換算する(割引きする)必要があります.
オプションの価格決定についても将来時点でのオプションの価値を割引いて,現在時点での価値を計算する必要があります .そのため,オプションの価格は金利に依存することになります.詳細な話は実際にオプション価格決定を行う段階ででてきます.
資産の価格変動モデル
資産の(将来の)価格変動について考えてみます.資産として株を考えてみると,株価は取引所での売買の結果によって決まります.資産の価格変動を完全な形で扱おうとすると,取引所で取引するすべての人々の資金,戦略等を考慮した上でシミュレーションを行う必要が出てきますが,東証一部だけで1日数兆円ほどの売買があります.これをすべて考慮することは実質上不可能なので,何かしらのモデル化を行う必要があります.
金融の問題では「資産の価格変動が確率過程に従う 」というモデル化が多く用いられています.価格変動は人間の意思決定(どの銘柄をどれくらい買う,売るなど)の結果なので,そもそも確率的ではない(それぞれの人間の意思決定は確率的ではないはず)ですが,全参加者の意思決定を観察することもできないので,確率的に変動するものとして扱うことにしています.
一期間二項モデル
資産の価格変動について確率過程でモデル化することにしたので,オプション価格を決定する要素は揃いました.そこで,簡単なモデルでオプションの価格決定を行ってみます.
ここでは,最も簡単なモデルとして一期間二項モデルを考えてみます.
これは時間がt ∈ { 0 , T } t \in \{0, T\} t ∈ { 0 , T } の2時点だけであり,t = T t=T t = T における資産価格S ( T ) S(T) S ( T ) がS ( T ) = S u , S d S(T)=S_u,S_d S ( T ) = S u , S d の2つしかないモデルです.
t = 0 t=0 t = 0 においてはS ( 0 ) = S 0 S(0)=S_0 S ( 0 ) = S 0 を取ります.
また,S u S_u S u へ遷移する確率をp p p ,S d S_d S d へ遷移する確率を1 − p 1-p 1 − p とします.
数値例として,S 0 = 100 , S u = 110 , S d = 95 , p = 1 / 2 S_0=100, S_u=110, S_d=95, p=1/2 S 0 = 1 0 0 , S u = 1 1 0 , S d = 9 5 , p = 1 / 2 とします.このとき,価格変動は下図のようになります.
さて,一期間二項モデルで価格変動をモデル化しました.
ここで求めたいのは満期がT T T で権利行使価格がK K K のオプションの価格でした.
ここでは一期間の金利がr = 0.01 r=0.01 r = 0 . 0 1 でK = 95 K=95 K = 9 5 のコールオプション(t = T t=T t = T において資産1単位を95 95 9 5 で買う権利)の価格を決定します.
オプションの価格決定の際に重要なのは「複数の資産がある場合に,現金の収入,支出(キャッシュフロー)が同一であればそれらの資産は同じ価値を持つ 」というものです.もし,収入,支出が同一にも関わらず,価値が異なるものが存在すると,元手が0でも必ず利益を得られることになります
これに基づくと,オプションと同一のキャッシュフローを実現する資産が存在すれば,その資産の価値はオプションの価値と同一になる ことがわかります.これを使ってt = 0 t=0 t = 0 でのオプションの価値を求めることでオプションの価格C C C が決定できます.
ということでまず,コールオプションのキャッシュフローを考えます.t = 0 t=0 t = 0 ではコールオプションを買うので支出がC C C となります.t = T t=T t = T では,S ( T ) = 110 S(T)=110 S ( T ) = 1 1 0 の場合に収入が110 − 95 = 15 110-95=15 1 1 0 − 9 5 = 1 5 で,S ( T ) = 90 S(T)=90 S ( T ) = 9 0 の場合には権利行使すると損するので権利行使しないので収支は0 0 0 です.
さて,ここでt = 0 t=0 t = 0 においてy y y 単位の現金を借り入れ,資産をx x x 単位買ったとします.このような,資産の組み合わせをポートフォリオ と呼びます.ポートフォリオ自体を一つの資産として扱うことが可能です.
では,このポートフォリオのキャッシュフローを考えます.t = 0 t=0 t = 0 では収支は− 100 x + y -100x+y − 1 0 0 x + y となります.t = T t=T t = T で,持っている資産をすべて売却し,借入金をすべて返済すると,S ( T ) = 110 S(T)=110 S ( T ) = 1 1 0 の場合に110 x − 1.01 y 110x-1.01y 1 1 0 x − 1 . 0 1 y で,S ( T ) = 90 S(T)=90 S ( T ) = 9 0 の場合に90 x − 1.01 y 90x-1.01y 9 0 x − 1 . 0 1 y となります.
各時点でのコールオプションとポートフォリオのキャッシュフローが同一になるのは
− 100 x + y = − C 110 x − 1.01 y = 15 90 x − 1.01 y = 0 \begin{aligned}
&-100x + y = -C\\
&110x-1.01y=15 \\
&90x-1.01y=0
\end{aligned} − 1 0 0 x + y = − C 1 1 0 x − 1 . 0 1 y = 1 5 9 0 x − 1 . 0 1 y = 0
の場合です.これを解くとx ≅ 0.75 , y ≅ 66.83 , C ≅ 8.17 x\cong0.75,y\cong66.83, C\cong8.17 x ≅ 0 . 7 5 , y ≅ 6 6 . 8 3 , C ≅ 8 . 1 7 となります.
これにより,ヨーロピアンコールオプションの価格を決定することができました.
一期間二項モデルの結果は以下のようにまとめることが可能です.
C u = max ( S 0 u − K , 0 ) C d = max ( S 0 d − K , 0 ) x = C u − C d S 0 ( u − d ) y = u C d − d C u ( 1 + r ) ( 1 − d ) C = 1 1 + r { ( 1 + r ) − d u − d C u + u − ( 1 + r ) u − d C d } \begin{aligned}
C_u&=\max(S_0u-K, 0) \\
C_d&=\max(S_0d-K, 0) \\
x&=\frac{C_u-C_d}{S_0(u-d)} \\
y&=\frac{uC_d-dC_u}{(1+r)(1-d)}\\
C&=\frac{1}{1+r}\left\{\frac{(1+r)-d}{u-d}C_u+\frac{u-(1+r)}{u-d}C_d\right\}
\end{aligned} C u C d x y C = max ( S 0 u − K , 0 ) = max ( S 0 d − K , 0 ) = S 0 ( u − d ) C u − C d = ( 1 + r ) ( 1 − d ) u C d − d C u = 1 + r 1 { u − d ( 1 + r ) − d C u + u − d u − ( 1 + r ) C d }
ただし,ここではS u = S 0 u , S d = S 0 d , u > d S_u=S_0u, S_d=S_0d,u > d S u = S 0 u , S d = S 0 d , u > d とおきました.
また,C u , C d C_u, C_d C u , C d はt = T t=T t = T でのオプションの価格を表しています.
ここで一つ気づくことは,オプションの価格が確率p p p よらないことです.一方で
p ∗ = ( 1 + r ) − d u − d q ∗ = u − ( 1 + r ) u − d \begin{aligned}
p^\ast = \frac{(1+r)-d}{u-d}\\
q^\ast = \frac{u-(1+r)}{u-d}
\end{aligned} p ∗ = u − d ( 1 + r ) − d q ∗ = u − d u − ( 1 + r )
とおくと,0 ≤ p ∗ , q ∗ ≤ 1 , p ∗ + q ∗ = 1 0\leq p^\ast,q^\ast\leq 1, p^\ast+q^\ast=1 0 ≤ p ∗ , q ∗ ≤ 1 , p ∗ + q ∗ = 1 であり,p ∗ , q ∗ p^\ast, q^\ast p ∗ , q ∗ は確率とみなすことができます.
このp ∗ , q ∗ p^\ast, q^\ast p ∗ , q ∗ はリスク中立確率 と呼ばれます.
これをつかうと
C = 1 1 + r ( p ∗ C u + q ∗ C d ) C=\frac{1}{1+r}\left(p^\ast C_u+q^\ast C_d\right) C = 1 + r 1 ( p ∗ C u + q ∗ C d )
と書くことができます.この結果を改めて解釈すると,オプションの価格は確率p ∗ , q ∗ p^\ast, q^\ast p ∗ , q ∗ のもとでのt = T t=T t = T におけるオプション価格C ( T ) = { C u , C d } C(T)=\{C_u, C_d\} C ( T ) = { C u , C d } の期待値を現在価値に割り引いたもの とみなすことができます.
一期間二項モデルについて,オプションと同じキャッシュフローを持つポートフォリオを作成することで,オプションの価格決定を行うことができました.この結果はm m m 期間n n n 項モデルにも拡張可能です.
連続型確率過程のもとでのオプションの価格決定
さて,ここまで考えてきた一期間二項モデルは現実とは程遠いモデルなので,もう少し現実の価格変動に近い確率過程を考えてみます.
デリバティブの価格決定問題では価格変動のモデルに幾何ブラウン運動 などがよく用いられます.
幾何ブラウン運動のサンプルパスは上図のようになります.
一期間二項モデルは時刻は2時点のみでT T T 時点でのとりうる資産価格は2つしかありません.
対して幾何ブラウン運動はt ∈ [ 0 , T ] t \in [0, T] t ∈ [ 0 , T ] の連続時間確率過程となっており,とりうる資産価格は連続量になっています.
ここでは連続時間確率過程を使って資産価格の変動をモデル化します.
連続時間確率過程においては,ヨーロピアンコールオプションの価格は
C = e − r T E Q [ f ( S ( T ) ) ∣ S ( 0 ) = S 0 ] C=e^{-rT}E_Q\left[f(S(T))|S(0)=S_0\right] C = e − r T E Q [ f ( S ( T ) ) ∣ S ( 0 ) = S 0 ]
と書くことができます.ここでE Q [ ⋅ ∣ S ( 0 ) = S 0 ] E_Q\left[\cdot|S(0)=S_0\right] E Q [ ⋅ ∣ S ( 0 ) = S 0 ] はリスク中立確率(測度)のもとでの条件付期待値を表し,f ( S ( T ) ) = max ( S ( T ) − K , 0 ) f(S(T))=\max(S(T)-K, 0) f ( S ( T ) ) = max ( S ( T ) − K , 0 ) は満期時点でのオプション価格を表します.e − r T e^{-rT} e − r T は割引のファクターで,一期間二項モデルの1 1 + r \frac{1}{1+r} 1 + r 1 に対応しています.
多少表記がややこしくなりましたが,やることは一緒で,リスク中立確率のもとで満期時点でのオプション価格の期待値を計算すれば良い ということです.
幾何ブラウン運動の場合には解析解が知られています(Black-Scholesモデル).
しかしながら,一般の連続型確率過程については数値計算 によりこの期待値を求める必要があります.
数値計算による方法としてはツリー法 や,モンテカルロ法が挙げられます.
3. モンテカルロ法
ここでは最もよく使われるモンテカルロ法を紹介します.
ここではt = T t=T t = T での資産価格のリスク中立確率測度下での確率分布が既知であるとします.
特に,幾何ブラウン運動の場合には
S ( T ) = S 0 exp [ ( r − σ 2 2 ) T + σ T W ] S(T) = S_0\exp\left[\left(r-\frac{\sigma^2}{2}\right)T+\sigma\sqrt{T}W\right] S ( T ) = S 0 exp [ ( r − 2 σ 2 ) T + σ T W ]
となることが知られています.ここでr r r は金利,σ \sigma σ は幾何ブラウン運動の標準偏差,W W W は標準正規分布に従う確率変数です.
モンテカルロ法による期待値計算では積分をする代わりに確率分布に従う乱数を発生させ,積分を標本平均で近似します .
ここではヨーロピアンコールオプションの価格決定をモンテカルロ法で行う手順を紹介します.
以下の2から4をN回繰り返す
標準正規分布に従う乱数を発生させてW i W_i W i とする
S i = S 0 exp [ ( r − σ 2 2 ) T + σ T W i ] S_i=S_0\exp\left[\left(r-\frac{\sigma^2}{2}\right)T+\sigma\sqrt{T}W_i\right] S i = S 0 exp [ ( r − 2 σ 2 ) T + σ T W i ] を計算する
f ( S i ) = max [ S i − K , 0 ] f(S_i)=\max\left[S_i-K, 0\right] f ( S i ) = max [ S i − K , 0 ] を計算する
N個のサンプルが得られるので,E Q [ f ( S ( T ) ) ∣ S ( 0 ) = S 0 ] ≅ 1 N ∑ i = 1 N f ( S i ) E_Q\left[f(S(T))|S(0)=S_0\right]\cong \frac{1}{N}\sum_{i=1}^Nf(S_i) E Q [ f ( S ( T ) ) ∣ S ( 0 ) = S 0 ] ≅ N 1 ∑ i = 1 N f ( S i ) として近似する
さて,自分の欲しい精度で価格決定を行いたい場合に必要なN N N を見積もりたいと思います.
ここでX = f ( S ( T ) ) X=f(S(T)) X = f ( S ( T ) ) とおきます.標本平均1 N ∑ i = 1 N X i \frac{1}{N}\sum_{i=1}^NX_i N 1 ∑ i = 1 N X i の分散を考えると
V [ 1 N ∑ i = 1 N X i ] = 1 N 2 ∑ i = 1 N V [ X i ] = σ 2 N \begin{aligned}
V\left[\frac{1}{N}\sum_{i=1}^NX_i\right]&=\frac{1}{N^2}\sum_{i=1}^NV[X_i]\\
&=\frac{\sigma^2}{N}
\end{aligned} V [ N 1 i = 1 ∑ N X i ] = N 2 1 i = 1 ∑ N V [ X i ] = N σ 2
したがって標本平均の標準偏差ϵ \epsilon ϵ はO ( 1 / N ) O(1/\sqrt{N}) O ( 1 / N ) で減少することがわかります.
逆に言うと欲しい精度ϵ \epsilon ϵ があった場合N = O ( 1 / ϵ 2 ) N=O(1/\epsilon^2) N = O ( 1 / ϵ 2 ) だけ繰り返せばいい という事がわかります.
4. 量子アルゴリズム
やっと本題に入ります .
ここでは量子アルゴリズムを用いてヨーロピアンオプションの価格決定を行います.
モンテカルロ法によるオプションの価格をϵ \epsilon ϵ の精度(標準偏差)で得たい場合にはN = O ( 1 / ϵ 2 ) N=O(1/\epsilon^2) N = O ( 1 / ϵ 2 ) のサンプルが必要になるのでした.
量子アルゴリズムをもちいることでこれがN = ( 1 / ϵ ) N=(1/\epsilon) N = ( 1 / ϵ ) に減らせます
これを見ていきましょう.
大まかに言うとオプション価格計算に必要なのは次の2ステップです
Distribution Loading: ∑ j = 0 2 n − 1 p j ∣ j > \sum_{j=0}^{2^n-1}\sqrt{p_j}\left| j \right> ∑ j = 0 2 n − 1 p j ∣ j ⟩ を作ります.ここでn n n は量子ビット数で,{ p i } , ( i = 0 , 1 , … , 2 n − 1 ) \{p_i\}, (i=0,1,\dots,2^n-1) { p i } , ( i = 0 , 1 , … , 2 n − 1 ) は確率密度関数p ( x ) p(x) p ( x ) を離散化したものです.
Quantum Amplitude Estimation : N = O ( 1 / ϵ ) N=O(1/\epsilon) N = O ( 1 / ϵ ) で期待値を計算します.
これらを詳しく説明したいと思います.
その前に,これらのステップに必要なモジュールを用意していきたいと思います.
Quantum Arithmetics
量子回路を用いた算術演算をQuantum Arithmeticsと呼んでいます.
これは既知の関数f f f に対して
∣ x > → ∣ f ( x ) > \left| x \right> \rightarrow \left| f(x) \right> ∣ x ⟩ → ∣ f ( x ) ⟩
を行う回路です.
重ね合わせられた状態に対しては
∑ i ∣ x i > → ∑ i ∣ f ( x i ) > \sum_i\left| x_i \right> \rightarrow \sum_i\left| f(x_i) \right> i ∑ ∣ x i ⟩ → i ∑ ∣ f ( x i ) ⟩
となります.
Quantum Phase Estimation
ユニタリ演算子U U U に対し,その固有値はe i θ e^{i\theta} e i θ という形で表されます.
このU U U と,固有値e i θ e^{i\theta} e i θ に対応する固有状態∣ ψ > \left| \psi \right> ∣ ψ ⟩ に対して
∣ ψ > → Q P E ∣ ψ > ∣ θ > \left| \psi \right> \stackrel{QPE}{\rightarrow} \left| \psi \right>\left| \theta \right> ∣ ψ ⟩ → Q P E ∣ ψ ⟩ ∣ θ ⟩
という形で固有値の位相を取り出すことができます.
詳細はNielsen Chuang 等を参照してください.
Distribution Loading
まず,仮定として任意のa ≤ b = 0 , … 2 n − 1 a\leq b=0,\dots2^{n-1} a ≤ b = 0 , … 2 n − 1 に対して,
∑ j = a b p j ≅ ∫ x a x b p ( x ) d x \sum_{j=a}^{b}p_j\cong\int_{x_a}^{x_b}p(x)dx j = a ∑ b p j ≅ ∫ x a x b p ( x ) d x
が古典で効率よく計算できるとします.
計算できる,というだけで2 n 2^n 2 n 個全てについて計算するわけではないことに注意します.例えば幾何ブラウン運動などではこれが満たされています.
ここで,m = 1 , … , n ; k = 0 , … , 2 m − 1 m=1,\dots,n; k=0,\dots,2^{m-1} m = 1 , … , n ; k = 0 , … , 2 m − 1 に対してp k ( m ) p_k^{(m)} p k ( m ) を
p k ( m ) = ∑ j = k 2 n − m ( k + 1 ) 2 n − m − 1 p j p_k^{(m)}=\sum_{j=k2^{n-m}}^{(k+1)2^{n-m}-1}p_j p k ( m ) = j = k 2 n − m ∑ ( k + 1 ) 2 n − m − 1 p j
と定義します.
ややこしいので,いくつかのm , k m,k m , k について見てみます.
m = 1 , k = 0 m=1, k=0 m = 1 , k = 0 でp 0 ( 1 ) = ∑ j = 0 2 n − 1 − 1 p j p_0^{(1)}=\sum_{j=0}^{2^{n-1}-1}p_j p 0 ( 1 ) = ∑ j = 0 2 n − 1 − 1 p j
m = 2 , k = 0 m=2, k=0 m = 2 , k = 0 でp 0 ( 2 ) = ∑ j = 0 2 n − 2 − 1 p j p_0^{(2)}=\sum_{j=0}^{2^{n-2}-1}p_j p 0 ( 2 ) = ∑ j = 0 2 n − 2 − 1 p j
m = 2 , k = 1 m=2, k=1 m = 2 , k = 1 でp 0 ( 2 ) = ∑ j = 2 n − 2 2 n − 1 − 1 p j p_0^{(2)}=\sum_{j=2^{n-2}}^{2^{n-1}-1}p_j p 0 ( 2 ) = ∑ j = 2 n − 2 2 n − 1 − 1 p j
m = 1 m=1 m = 1 では2 n − 1 2^{n-1} 2 n − 1 個の和ですが,m = 2 m=2 m = 2 ではk = 0 k=0 k = 0 とk = 1 k=1 k = 1 で2 n − 2 2^{n-2} 2 n − 2 個ずつの和に分けられます.
m = n m=n m = n ではp k ( n ) = p k p_k^{(n)}=p_k p k ( n ) = p k となります.
ここで,∣ ψ ( m ) > \left| \psi^{(m)} \right> ∣ ∣ ∣ ψ ( m ) ⟩ を
∣ ψ ( m ) > = ∑ k = 0 2 m − 1 p k ( m ) ∣ k > \left| \psi^{(m)} \right>=\sum_{k=0}^{2^m-1}\sqrt{p_k^{(m)}}\left| k \right> ∣ ∣ ∣ ∣ ψ ( m ) ⟩ = k = 0 ∑ 2 m − 1 p k ( m ) ∣ k ⟩
と定義します.
このとき∣ ψ ( m ) > \left| \psi^{(m)} \right> ∣ ∣ ∣ ψ ( m ) ⟩ から∣ ψ ( m + 1 ) > \left| \psi^{(m+1)} \right> ∣ ∣ ∣ ψ ( m + 1 ) ⟩ を計算できればm = n m=n m = n まで計算して
∣ ψ ( n ) > = ∑ k = 0 2 n − 1 p k ( n ) ∣ k > = ∑ k = 0 2 n − 1 p k ∣ k > \left| \psi^{(n)} \right> = \sum_{k=0}^{2^n-1}\sqrt{p_k^{(n)}}\left| k \right> = \sum_{k=0}^{2^n-1}\sqrt{p_k}\left| k \right> ∣ ∣ ∣ ∣ ψ ( n ) ⟩ = k = 0 ∑ 2 n − 1 p k ( n ) ∣ k ⟩ = k = 0 ∑ 2 n − 1 p k ∣ k ⟩
となり,所望の状態が得られます.
では,∣ ψ ( m ) > \left| \psi^{(m)} \right> ∣ ∣ ∣ ψ ( m ) ⟩ から∣ ψ ( m + 1 ) > \left| \psi^{(m+1)} \right> ∣ ∣ ∣ ψ ( m + 1 ) ⟩ を計算してみましょう.
まず,∑ j = a b p j \sum_{j=a}^{b}p_j ∑ j = a b p j が効率よく計算できるのでQuantum Arithmeticによって
∣ k > → ∣ k > ∣ cos − 1 p 2 k ( m + 1 ) p k ( m ) > \left| k \right> \rightarrow \left| k \right>\left| \cos^{-1}\sqrt{\frac{p_{2k}^{(m+1)}}{p_k^{(m)}}} \right> ∣ k ⟩ → ∣ k ⟩ ∣ ∣ ∣ ∣ ∣ ∣ ∣ cos − 1 p k ( m ) p 2 k ( m + 1 ) ⟩
が計算可能です.ここで,突然∣ cos − 1 p 2 k ( m + 1 ) p k ( m ) > \left| \cos^{-1}\sqrt{\frac{p_{2k}^{(m+1)}}{p_k^{(m)}}} \right> ∣ ∣ ∣ ∣ ∣ cos − 1 p k ( m ) p 2 k ( m + 1 ) ⟩ が現れましたが,これは補助量子ビットに関するものです.つまり,上の式は本当は
∣ k > ∣ 0 > → ∣ k > ∣ cos − 1 p 2 k ( m + 1 ) p k ( m ) > \left| k \right>\left| 0 \right>\rightarrow \left| k \right>\left| \cos^{-1}\sqrt{\frac{p_{2k}^{(m+1)}}{p_k^{(m)}}} \right> ∣ k ⟩ ∣ 0 ⟩ → ∣ k ⟩ ∣ ∣ ∣ ∣ ∣ ∣ ∣ cos − 1 p k ( m ) p 2 k ( m + 1 ) ⟩
なんですが,補助量子ビットの∣ 0 > \left| 0 \right> ∣ 0 ⟩ は省略して書く,という記法を採用しています.
さて,今得られた状態に対してcontrol rotationゲートを用いて
∣ k > ∣ p 2 k ( m + 1 ) p k ( m ) > ∣ 0 > → ∣ k > ∣ p 2 k ( m + 1 ) p k ( m ) > ( p 2 k ( m + 1 ) p k ( m ) ∣ 0 > + 1 − p 2 k ( m + 1 ) p k ( m ) ∣ 1 > ) \left| k \right>\left| \sqrt{\frac{p_{2k}^{(m+1)}}{p_k^{(m)}}}\right>\left| 0 \right>\rightarrow\left| k \right>\left| \sqrt{\frac{p_{2k}^{(m+1)}}{p_k^{(m)}}} \right>\left(\sqrt{\frac{p_{2k}^{(m+1)}}{p_k^{(m)}}}\left| 0 \right>+\sqrt{1-\frac{p_{2k}^{(m+1)}}{p_k^{(m)}}}\left| 1 \right>\right) ∣ k ⟩ ∣ ∣ ∣ ∣ ∣ ∣ ∣ p k ( m ) p 2 k ( m + 1 ) ⟩ ∣ 0 ⟩ → ∣ k ⟩ ∣ ∣ ∣ ∣ ∣ ∣ ∣ p k ( m ) p 2 k ( m + 1 ) ⟩ ⎝ ⎜ ⎛ p k ( m ) p 2 k ( m + 1 ) ∣ 0 ⟩ + 1 − p k ( m ) p 2 k ( m + 1 ) ∣ 1 ⟩ ⎠ ⎟ ⎞
という状態を用意できます.control rotationゲートはレジスタのビット列を参照して,k k k 番目のビットが1 1 1 であれば補助量子ビットにR Y ( 2 − k ) R_Y(2^{-k}) R Y ( 2 − k ) をかける,という演算です.今後∣ p 2 k ( m + 1 ) p k ( m ) > \left|\sqrt{\frac{p_{2k}^{(m+1)}}{p_k^{(m)}}}\right> ∣ ∣ ∣ ∣ ∣ p k ( m ) p 2 k ( m + 1 ) ⟩ は使わないのでこのレジスタを逆演算を用いて∣ 0 > \left| 0 \right> ∣ 0 ⟩ に戻してやります.ここまでの演算を∣ ψ ( m ) > \left| \psi^{(m)} \right> ∣ ∣ ∣ ψ ( m ) ⟩ に作用させると
∣ ψ ( m ) > = ∑ k = 0 2 m − 1 p k ( m ) ∣ k > → ∑ k = 0 2 m − 1 p k ( m ) ∣ k > ( p 2 k ( m + 1 ) p k ( m ) ∣ 0 > + 1 − p 2 k ( m + 1 ) p k ( m ) ∣ 1 > ) = ∑ k = 0 2 m − 1 ( p 2 k ( m ) ∣ k > ∣ 0 > + p k ( m ) − p 2 k ( m + 1 ) ∣ k > ∣ 1 > ) = ∣ ψ ( m + 1 ) > \begin{aligned}
\left| \psi^{(m)} \right> = \sum_{k=0}^{2^m-1}\sqrt{p_k^{(m)}}\left| k \right> &\rightarrow \sum_{k=0}^{2^m-1}\sqrt{p_k^{(m)}}\left| k \right>\left(\sqrt{\frac{p_{2k}^{(m+1)}}{p_k^{(m)}}}\left| 0 \right>+\sqrt{1-\frac{p_{2k}^{(m+1)}}{p_k^{(m)}}}\left| 1 \right>\right)\\
&=\sum_{k=0}^{2^m-1}\left(\sqrt{p_{2k}^{(m)}}\left|k\right>\left|0\right>+\sqrt{p_{k}^{(m)}-p_{2k}^{(m+1)}}\left|k\right>\left|1\right>\right) = \left|\psi^{(m+1)}\right>
\end{aligned} ∣ ∣ ∣ ∣ ψ ( m ) ⟩ = k = 0 ∑ 2 m − 1 p k ( m ) ∣ k ⟩ → k = 0 ∑ 2 m − 1 p k ( m ) ∣ k ⟩ ⎝ ⎜ ⎛ p k ( m ) p 2 k ( m + 1 ) ∣ 0 ⟩ + 1 − p k ( m ) p 2 k ( m + 1 ) ∣ 1 ⟩ ⎠ ⎟ ⎞ = k = 0 ∑ 2 m − 1 ( p 2 k ( m ) ∣ k ⟩ ∣ 0 ⟩ + p k ( m ) − p 2 k ( m + 1 ) ∣ k ⟩ ∣ 1 ⟩ ) = ∣ ∣ ∣ ∣ ψ ( m + 1 ) ⟩
となり,∣ ψ ( m + 1 ) > \left|\psi^{(m+1)}\right> ∣ ∣ ∣ ψ ( m + 1 ) ⟩ が得られます.
これで∣ ψ ( m ) > \left|\psi^{(m)}\right> ∣ ∣ ∣ ψ ( m ) ⟩ から∣ ψ ( m + 1 ) > \left|\psi^{(m+1)}\right> ∣ ∣ ∣ ψ ( m + 1 ) ⟩ を得られるので,順にm = n m=n m = n まで計算していけば,ほしかった状態∑ k = 0 2 n − 1 p k ∣ k > \sum_{k=0}^{2^n-1}\sqrt{p_k}\left| k \right> ∑ k = 0 2 n − 1 p k ∣ k ⟩ が得られます.
Quantum Amplitude Estimation
Distribution Loadingで∑ k = 0 2 n − 1 p k ∣ k > \sum_{k=0}^{2^n-1}\sqrt{p_k}\left| k \right> ∑ k = 0 2 n − 1 p k ∣ k ⟩ が得られたので,期待値計算を考えてみます.
まずはQuantum Amplitude Estimationを使わずに期待値計算をしてみます.
rotationを使って
∑ k = 0 2 n − 1 p k ∣ k > → ∑ k = 0 2 n − 1 p k ∣ k > ( 1 − v ( k ) ∣ 0 > + v ( k ) ∣ 1 > ) : = ∣ χ > \sum_{k=0}^{2^n-1}\sqrt{p_k}\left| k \right> \rightarrow \sum_{k=0}^{2^n-1}\sqrt{p_k}\left| k \right>\left(\sqrt{1-v(k)}\left| 0 \right>+\sqrt{v(k)}\left| 1 \right>\right):=\left| \chi \right> k = 0 ∑ 2 n − 1 p k ∣ k ⟩ → k = 0 ∑ 2 n − 1 p k ∣ k ⟩ ( 1 − v ( k ) ∣ 0 ⟩ + v ( k ) ∣ 1 ⟩ ) : = ∣ χ ⟩
を得ます.補助量子ビットの∣ 1 > \left| 1 \right> ∣ 1 ⟩ を測定することで
μ = < χ ∣ I ⊗ n ⊗ ∣ 1 > < 1 ∣ χ > = ∑ k = 0 2 n − 1 p k v ( k ) \mu=\left< \chi \right|I^{\otimes n}\otimes\left| 1 \right>\left< 1 \right.\left| \chi \right>=\sum_{k=0}^{2^n-1}p_kv(k) μ = ⟨ χ ∣ I ⊗ n ⊗ ∣ 1 ⟩ ⟨ 1 ∣ χ ⟩ = k = 0 ∑ 2 n − 1 p k v ( k )
となり,期待値を計算することができます.
このとき,必要な測定回数を考えてみます.
平均がμ \mu μ の二項分布ではN N N 回サンプルした場合の標本分散がϵ 2 = μ ( 1 − μ ) N \epsilon^2=\frac{\mu(1-\mu)}{N} ϵ 2 = N μ ( 1 − μ ) であるので,精度ϵ \epsilon ϵ が必要な場合にはN = O ( 1 / ϵ 2 ) N=O(1/\epsilon^2) N = O ( 1 / ϵ 2 ) が必要になります.
これだと古典でのモンテカルロ法と変わらない,という結果になります.
そこで,量子加速を得るためにQuantum Amplitude Estimation (QAE)を用いることにします.
QAEではある回転行列を構成します.
回転行列の固有値がその回転角θ \theta θ を用いてe ± i θ e^{\pm i\theta} e ± i θ となるので,Quantum Phase Estimationを用いてθ \theta θ を取り出し,μ = 1 2 ( 1 − cos θ / 2 ) \mu=\frac{1}{2}(1-\cos\theta/2) μ = 2 1 ( 1 − cos θ / 2 ) から期待値を計算することができます.
以下でQAEの手順を見ていきたいと思います.
まず,
V = I ⊗ n + 1 − 2 I ⊗ n ⊗ ∣ 1 > < 1 ∣ V=I^{\otimes n+1}-2I^{\otimes n}\otimes\left| 1 \right>\left< 1 \right| V = I ⊗ n + 1 − 2 I ⊗ n ⊗ ∣ 1 ⟩ ⟨ 1 ∣
というユニタリ演算子を用意します.
この演算子V V V の∣ χ > \left| \chi \right> ∣ χ ⟩ での期待値で考えると
< χ ∣ V ∣ χ > = 1 − 2 μ \left< \chi \right|V\left| \chi \right>=1-2\mu ⟨ χ ∣ V ∣ χ ⟩ = 1 − 2 μ
一方で,V V V を∣ χ > \left| \chi \right> ∣ χ ⟩ に作用させると,∣ χ > \left| \chi \right> ∣ χ ⟩ とそれに直交する∣ χ ⊥ > \left| \chi_\perp \right> ∣ χ ⊥ ⟩ の線形結合でかけます.これを
V ∣ χ > = cos ( θ / 2 ) ∣ χ > + e i ϕ sin ( θ / 2 ) ∣ χ ⊥ > V\left| \chi \right>=\cos(\theta/2)\left| \chi \right>+e^{i\phi}\sin(\theta/2)\left| \chi_\perp \right> V ∣ χ ⟩ = cos ( θ / 2 ) ∣ χ ⟩ + e i ϕ sin ( θ / 2 ) ∣ χ ⊥ ⟩
と表現すると,
1 − 2 μ = cos ( θ / 2 ) 1-2\mu=\cos(\theta/2) 1 − 2 μ = cos ( θ / 2 )
であるから,μ \mu μ を求める代わりにθ \theta θ を計算する問題に置き換えることができます.
また,V ∣ χ > V\left| \chi \right> V ∣ χ ⟩ は∣ χ > \left| \chi \right> ∣ χ ⟩ と∣ χ ⊥ > \left| \chi_\perp \right> ∣ χ ⊥ ⟩ で張られる二次元平面内で∣ χ > \left| \chi \right> ∣ χ ⟩ をθ / 2 \theta/2 θ / 2 だけ回転させたものとみなすことができます.
さらに
U = I ⊗ n + 1 − 2 ⊗ ∣ χ > < χ ∣ U=I^{\otimes n+1}-2\otimes \left| \chi \right>\left< \chi \right| U = I ⊗ n + 1 − 2 ⊗ ∣ χ ⟩ ⟨ χ ∣
S = I ⊗ n + 1 − 2 ⊗ V ∣ χ > < χ ∣ V S=I^{\otimes n+1}-2\otimes V\left| \chi \right>\left< \chi \right|V S = I ⊗ n + 1 − 2 ⊗ V ∣ χ ⟩ ⟨ χ ∣ V
というユニタリ演算子を定義します.− U -U − U は,∣ χ > \left| \chi \right> ∣ χ ⟩ と∣ χ ⊥ > \left| \chi_\perp \right> ∣ χ ⊥ ⟩ で張られる二次元平面内のあるベクトル∣ ψ > \left| \psi \right> ∣ ψ ⟩ を∣ χ > \left| \chi \right> ∣ χ ⟩ に関して対称なベクトルに写す操作に対応しており,
− S -S − S は∣ ψ > \left| \psi \right> ∣ ψ ⟩ をV ∣ χ > V\left| \chi \right> V ∣ χ ⟩ に関して対称なベクトルに写す操作に対応しています.
ここで,
Q = U S Q=US Q = U S
と定義すると,このQ Q Q はV ∣ χ > V\left| \chi \right> V ∣ χ ⟩ は∣ χ > \left| \chi \right> ∣ χ ⟩ と∣ χ ⊥ > \left| \chi_\perp \right> ∣ χ ⊥ ⟩ で張られる二次元平面内で∣ ψ > \left| \psi \right> ∣ ψ ⟩ をθ \theta θ だけ回転させる操作となります.
回転角θ \theta θ をもつ回転行列の固有値はe ± i θ e^{\pm i\theta} e ± i θ であるので,QPEを用いてこの回転行列の固有値を推定することができれば,期待値μ \mu μ を計算できます.
では計算のステップ数はどの程度になるでしょうか?
実はQAEを用いてU , V U,V U , V をO ( N log ( 1 / δ ) ) O(N\log(1/\delta)) O ( N log ( 1 / δ ) ) 回作用させたとき,μ \mu μ の推定値μ ^ \hat{\mu} μ ^ の推定誤差が
∣ μ ^ − μ ∣ ≤ C ( μ N + 1 N 2 ) |\hat{\mu}-\mu| \leq C\left(\frac{\sqrt{\mu}}{N}+\frac{1}{N^2}\right) ∣ μ ^ − μ ∣ ≤ C ( N μ + N 2 1 )
となる確率の下限が1 − δ 1-\delta 1 − δ であることが知られています.
δ \delta δ を固定してやったとき,ほしい精度ϵ \epsilon ϵ に対して
∣ μ ^ − μ ∣ ≤ ϵ = O ( 1 / N ) |\hat{\mu}-\mu|\leq\epsilon=O(1/N) ∣ μ ^ − μ ∣ ≤ ϵ = O ( 1 / N ) という評価ができるので
N = O ( 1 ϵ ) N=O\left(\frac{1}{\epsilon}\right) N = O ( ϵ 1 )
が得られることになります.
これでQAEを用いることでヨーロピアンオプションの価格計算に必要な計算ステップ数が
O ( 1 ϵ 2 ) → O ( 1 ϵ ) O\left(\frac{1}{\epsilon^2}\right)\rightarrow O\left(\frac{1}{\epsilon}\right) O ( ϵ 2 1 ) → O ( ϵ 1 )
となることがわかりました.
まとめと展望
本稿では金融派生商品であるデリバティブを紹介し,特にデリバティブの中でもヨーロピアンオプションの価格計算について紹介しました.古典の数値計算手法としてはモンテカルロ法を紹介しました.
さらに,量子振幅推定を用いることでヨーロピアンオプションの価格決定問題に対して2次加速が得られることを示しました.
Distribution loadingに関してはquantum arithmeticsを最適化する方法 や,qGANを用いる方法 が提案されています.
また,QAEをQPEなしで行う方法 も提案されています.
本稿で紹介したアルゴリズムは量子回路のサイズがかなり大きくなってしまうので,これらの提案を用いることで少しでも回路を減らしていこうという試みがなされています.
間違い等あれば是非コメントお待ちしております.