Japanese
言語
English
Japanese
German
Korean
Portuguese, Brazilian
French
Table of Contents
Shortcuts

注釈

当ページは tutorials/finance/07_asian_barrier_spread_pricing.ipynb から生成されました。

IBM Quantum lab でインタラクティブに実行します。

アジアン・バリア・スプレッド・オプションの価格推定

はじめに

アジアン・バリア・スプレッドオプションは３つのオプション・タイプを組み合わせたもので、それゆえにQiskit Finance オプションの価格設定フレームワークがサポートする複数の機能を組み合わせています。

  • アジアン・オプション: ペイオフは、定められた期間内の原資産の平均価格によって決まります。

  • バリア・オプション: 定められた期間内で原資産価格があるしきい値（バリア）を超過すると、ペイオフはゼロになります。

  • ブル・スプレッド: ペイオフは（平均価格に応じて）ゼロから始まり、線形に増加し、一定に維持される部分線形関数に従います。

権利行使価格 \(K_1 < K_2\) と期間 \(t=1,2\) が，与えられた多変量分布 (確率過程によって生成されたものなど) に従うスポット価格 \((S_1, S_2)\) とバリア閾値 \(B>0\) に対応するとします。対応するペイオフ関数は次のように定義されます。

\[\begin{split}P(S_1, S_2) = \begin{cases} \min\left\{\max\left\{\frac{1}{2}(S_1 + S_2) - K_1, 0\right\}, K_2 - K_1\right\}, & \text{ if } S_1, S_2 \leq B \\ 0, & \text{otherwise.} \end{cases}\end{split}\]

以下では、振幅推定に基づく量子アルゴリズムを使用して、期待されるペイオフ、すなわちオプションの割引前の公正価格を見積もります。

\[\mathbb{E}\left[ P(S_1, S_2) \right].\]

目的関数の近似と量子コンピューターによる一般的なオプション価格設定とリスク分析は、以下の論文で紹介されています。

  • Quantum Risk Analysis. Woerner, Egger. 2018.

  • Option Pricing using Quantum Computers. Stamatopoulos et al. 2019.

[1]:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.interpolate import griddata
%matplotlib inline
import numpy as np

from qiskit import QuantumRegister, QuantumCircuit, Aer, execute, AncillaRegister, transpile
from qiskit.circuit.library import IntegerComparator, WeightedAdder, LogNormalDistribution, LinearAmplitudeFunction
from qiskit.aqua.algorithms import IterativeAmplitudeEstimation

不確実性モデル

多変量対数正規ランダム分布を \(n\) 量子ビットの量子状態にロードする回路ファクトリを構築します。 すべての次元 \(j = 1,\ldots,d\) について、分布は指定された区間 \([\text{low}_j, \text{high}_j]\) に切り捨てられ、 \(2^{n_j}\) グリッドポイントを使用して離散化されます。ここで、 \(n_j\) は次元 \(j\) を表すために使用される量子ビットの数を示します。つまり \(n_1+\ldots+n_d = n\) です。 回路ファクトリに対応するユニタリー演算子は、以下のように実装されます:

\[\big|0\rangle_{n} \mapsto \big|\psi\rangle_{n} = \sum_{i_1,\ldots,i_d} \sqrt{p_{i_1\ldots i_d}}\big|i_1\rangle_{n_1}\ldots\big|i_d\rangle_{n_d},\]

ここで \(p_{i_1\ldots i_d}\) は、切り捨て分布と離散化分布に対応する確率を表し、 \(i_j\) はアフィンマップを用いて右区間にマップされています。

\[\{0, \ldots, 2^{n_j}-1\} \ni i_j \mapsto \frac{\text{high}_j - \text{low}_j}{2^{n_j} - 1} * i_j + \text{low}_j \in [\text{low}_j, \text{high}_j].\]

簡単のため、ここでは両方の株価が独立しており、同一の分布を取るものと仮定します。この仮定は、以下のパラメトライゼーション簡略化し、より複雑なそして相関のある多変量分布に緩和することが容易に出来ます。現在の実装のための唯一の重要な仮定は、異なる次元の離散化グリッドが同じステップサイズを持つということです。

[2]:

# number of qubits per dimension to represent the uncertainty
num_uncertainty_qubits = 2

# parameters for considered random distribution
S = 2.0 # initial spot price
vol = 0.4 # volatility of 40%
r = 0.05 # annual interest rate of 4%
T = 40 / 365 # 40 days to maturity

# resulting parameters for log-normal distribution
mu = ((r - 0.5 * vol**2) * T + np.log(S))
sigma = vol * np.sqrt(T)
mean = np.exp(mu + sigma**2/2)
variance = (np.exp(sigma**2) - 1) * np.exp(2*mu + sigma**2)
stddev = np.sqrt(variance)

# lowest and highest value considered for the spot price; in between, an equidistant discretization is considered.
low  = np.maximum(0, mean - 3*stddev)
high = mean + 3*stddev

# map to higher dimensional distribution
# for simplicity assuming dimensions are independent and identically distributed)
dimension = 2
num_qubits=[num_uncertainty_qubits]*dimension
low=low*np.ones(dimension)
high=high*np.ones(dimension)
mu=mu*np.ones(dimension)
cov=sigma**2*np.eye(dimension)

# construct circuit factory
u = LogNormalDistribution(num_qubits=num_qubits, mu=mu, sigma=cov, bounds=(list(zip(low, high))))

[3]:

# plot PDF of uncertainty model
x = [ v[0] for v in u.values ]
y = [ v[1] for v in u.values ]
z = u.probabilities
#z = map(float, z)
#z = list(map(float, z))
resolution = np.array([2**n for n in num_qubits])*1j
grid_x, grid_y = np.mgrid[min(x):max(x):resolution[0], min(y):max(y):resolution[1]]
grid_z = griddata((x, y), z, (grid_x, grid_y))
fig = plt.figure(figsize=(10, 8))
ax = fig.gca(projection='3d')
ax.plot_surface(grid_x, grid_y, grid_z, cmap=plt.cm.Spectral)
ax.set_xlabel('Spot Price $S_1$ (\$)', size=15)
ax.set_ylabel('Spot Price $S_2$ (\$)', size=15)
ax.set_zlabel('Probability (\%)', size=15)
plt.show()
../../_images/tutorials_finance_07_asian_barrier_spread_pricing_5_0.png

ペイオフ関数

簡単のため、平均の代わりにスポット価格の合計を考えます。結果を２で割るだけで平均にできます。

ペイオフ関数は、スポット価格合計 \((S_1 + S_2)\) が権利行使価格 \(K_1\) よりも小さい間ゼロであり、その後スポット価格合計が \(K_2\) になるまで線形に増加します。そしてペイオフはどちらかのスポット価格がバリアしきい値 \(B\) を超えるまで、一定 \(K_2 - K_1\) ですが、その後ペイオフは直ちにゼロに落ちます。実装は初め加重合計演算子を使って補助量子ビットにスポット価格合計を計算しますが、それからコンパレーターを使って、もし \((S_1 + S_2) \geq K_1\) ならば補助量子ビットを \(\big|0\rangle\) to \(\big|1\rangle\) にフリップします。そしてもう一つのコンパレーターと補助量子ビットを \((S_1 + S_2) \geq K_2\) のケースを記録するためにフリップします。これらの補助量子ビットはペイオフ関数の線形部分をコントロールするために使われます。

それに加え、別の補助量子ビットを時間ステップごとに追加します。そして追加のコンパレーターで \(S_1\) とそれぞれの \(S_2\) がバリアしきい値 \(B\) を超えていないかチェックします。ペイオフ関数は \(S_1, S_2 \leq B\) の場合のみ適用されます。

線形部分自体は次のように近似されます。 小さい \(|y|\) に対して \(\sin^2(y + \pi/4) \approx y + 1/2\) という事実を利用します。 したがって、与えられた近似スケーリング係数 \(c_\text{approx} \in [0, 1]\) および \(x \in [0, 1]\) に対して、次を考えます。

\[\sin^2( \pi/2 * c_\text{approx} * ( x - 1/2 ) + \pi/4) \approx \pi/2 * c_\text{approx} * ( x - 1/2 ) + 1/2\]

\(c_\text{approx}\) は十分小さいものとします。

次のように機能する演算子を簡単に構築できます

\[\big|x\rangle \big|0\rangle \mapsto \big|x\rangle \left( \cos(a*x+b) \big|0\rangle + \sin(a*x+b) \big|1\rangle \right),\]

制御Y-回転を使って。

最終的には、 \(\sin^2(a*x+b)\) に対応する最後の量子ビットで \(\big|1\rangle\) を測定する確率に関心があります。 上記の近似と合わせて、対象の値を近似できます。小さい \(c_\text{approx}\) を選択するほど、近似は良くなります。 ただし、 \(c_\text{approx}\) でスケーリングされたプロパティを推定しているため、それに応じて評価量子ビットの数 \(m\) を調整する必要があります。

近似の詳細については、以下を参照してください: Quantum Risk Analysis, Woerner, Egger. 2018.

(現在実装中の) 加重合計演算子は整数のみを合計することができるため、結果を推定するために元の範囲から表現可能な範囲に写像し、結果を解釈する前に逆写像する必要があります。このマッピングは不確定性モデルのコンテキストに出てきたアフィン写像に一致します。

[4]:

# determine number of qubits required to represent total loss
weights = []
for n in num_qubits:
    for i in range(n):
        weights += [2**i]

# create aggregation circuit
agg = WeightedAdder(sum(num_qubits), weights)
n_s = agg.num_sum_qubits
n_aux = agg.num_qubits - n_s - agg.num_state_qubits  # number of additional qubits

[5]:

# set the strike price (should be within the low and the high value of the uncertainty)
strike_price_1 = 3
strike_price_2 = 4

# set the barrier threshold
barrier = 2.5

# map strike prices and barrier threshold from [low, high] to {0, ..., 2^n-1}
max_value = 2**n_s - 1
low_ = low[0]
high_ = high[0]

mapped_strike_price_1 = (strike_price_1 - dimension*low_) / (high_ - low_) * (2**num_uncertainty_qubits - 1)
mapped_strike_price_2 = (strike_price_2 - dimension*low_) / (high_ - low_) * (2**num_uncertainty_qubits - 1)
mapped_barrier = (barrier - low) / (high - low) * (2**num_uncertainty_qubits - 1)

[6]:

# condition and condition result
conditions = []
barrier_thresholds = [2]*dimension
n_aux_conditions = 0
for i in range(dimension):
    # target dimension of random distribution and corresponding condition (which is required to be True)
    comparator = IntegerComparator(num_qubits[i], mapped_barrier[i] + 1, geq=False)
    n_aux_conditions = max(n_aux_conditions, comparator.num_ancillas)
    conditions += [comparator]

[7]:

# set the approximation scaling for the payoff function
c_approx = 0.25

# setup piecewise linear objective fcuntion
breakpoints = [0, mapped_strike_price_1, mapped_strike_price_2]
slopes = [0, 1, 0]
offsets = [0, 0, mapped_strike_price_2 - mapped_strike_price_1]
f_min = 0
f_max = mapped_strike_price_2 - mapped_strike_price_1
objective = LinearAmplitudeFunction(
    n_s,
    slopes,
    offsets,
    domain=(0, max_value),
    image=(f_min, f_max),
    rescaling_factor=c_approx,
    breakpoints=breakpoints
)

[8]:

# define overall multivariate problem
qr_state = QuantumRegister(u.num_qubits, 'state')  # to load the probability distribution
qr_obj = QuantumRegister(1, 'obj')  # to encode the function values
ar_sum = AncillaRegister(n_s, 'sum')  # number of qubits used to encode the sum
ar_cond = AncillaRegister(len(conditions) + 1, 'conditions')
ar = AncillaRegister(max(n_aux, n_aux_conditions, objective.num_ancillas), 'work')  # additional qubits

objective_index = u.num_qubits

# define the circuit
asian_barrier_spread = QuantumCircuit(qr_state, qr_obj, ar_cond, ar_sum, ar)

# load the probability distribution
asian_barrier_spread.append(u, qr_state)

# apply the conditions
for i, cond in enumerate(conditions):
    state_qubits = qr_state[(num_uncertainty_qubits * i):(num_uncertainty_qubits * (i + 1))]
    asian_barrier_spread.append(cond, state_qubits + [ar_cond[i]] + ar[:cond.num_ancillas])

# aggregate the conditions on a single qubit
asian_barrier_spread.mcx(ar_cond[:-1], ar_cond[-1])

# apply the aggregation function controlled on the condition
asian_barrier_spread.append(agg.control(), [ar_cond[-1]] + qr_state[:] + ar_sum[:] + ar[:n_aux])

# apply the payoff function
asian_barrier_spread.append(objective, ar_sum[:] + qr_obj[:] + ar[:objective.num_ancillas])

# uncompute the aggregation
asian_barrier_spread.append(agg.inverse().control(), [ar_cond[-1]] + qr_state[:] + ar_sum[:] + ar[:n_aux])

# uncompute the conditions
asian_barrier_spread.mcx(ar_cond[:-1], ar_cond[-1])

for j, cond in enumerate(reversed(conditions)):
    i = len(conditions) - j - 1
    state_qubits = qr_state[(num_uncertainty_qubits * i):(num_uncertainty_qubits * (i + 1))]
    asian_barrier_spread.append(cond.inverse(), state_qubits + [ar_cond[i]] + ar[:cond.num_ancillas])

print(asian_barrier_spread.draw())
print('objective qubit index', objective_index)

              ┌───────┐┌──────┐             ┌───────────┐      ┌──────────────┐»
     state_0: ┤0      ├┤0     ├─────────────┤1          ├──────┤1             ├»
              │       ││      │             │           │      │              │»
     state_1: ┤1      ├┤1     ├─────────────┤2          ├──────┤2             ├»
              │  P(X) ││      │┌──────┐     │           │      │              │»
     state_2: ┤2      ├┤      ├┤0     ├─────┤3          ├──────┤3             ├»
              │       ││      ││      │     │           │      │              │»
     state_3: ┤3      ├┤      ├┤1     ├─────┤4          ├──────┤4             ├»
              └───────┘│      ││      │     │           │┌────┐│              │»
       obj_0: ─────────┤      ├┤      ├─────┤           ├┤3   ├┤              ├»
                       │      ││      │     │           ││    ││              │»
conditions_0: ─────────┤2     ├┤      ├──■──┤           ├┤    ├┤              ├»
                       │  cmp ││      │  │  │           ││    ││              │»
conditions_1: ─────────┤      ├┤2     ├──■──┤           ├┤    ├┤              ├»
                       │      ││  cmp │┌─┴─┐│   c_adder ││    ││   c_adder_dg │»
conditions_2: ─────────┤      ├┤      ├┤ X ├┤0          ├┤    ├┤0             ├»
                       │      ││      │└───┘│           ││    ││              │»
       sum_0: ─────────┤      ├┤      ├─────┤5          ├┤0   ├┤5             ├»
                       │      ││      │     │           ││    ││              │»
       sum_1: ─────────┤      ├┤      ├─────┤6          ├┤1 F ├┤6             ├»
                       │      ││      │     │           ││    ││              │»
       sum_2: ─────────┤      ├┤      ├─────┤7          ├┤2   ├┤7             ├»
                       │      ││      │     │           ││    ││              │»
      work_0: ─────────┤3     ├┤3     ├─────┤8          ├┤4   ├┤8             ├»
                       └──────┘└──────┘     │           ││    ││              │»
      work_1: ──────────────────────────────┤9          ├┤5   ├┤9             ├»
                                            │           ││    ││              │»
      work_2: ──────────────────────────────┤10         ├┤6   ├┤10            ├»
                                            └───────────┘│    │└──────────────┘»
      work_3: ───────────────────────────────────────────┤7   ├────────────────»
                                                         └────┘                »
«                              ┌─────────┐
«     state_0: ────────────────┤0        ├
«                              │         │
«     state_1: ────────────────┤1        ├
«                   ┌─────────┐│         │
«     state_2: ─────┤0        ├┤         ├
«                   │         ││         │
«     state_3: ─────┤1        ├┤         ├
«                   │         ││         │
«       obj_0: ─────┤         ├┤         ├
«                   │         ││         │
«conditions_0: ──■──┤         ├┤2        ├
«                │  │         ││  cmp_dg │
«conditions_1: ──■──┤2        ├┤         ├
«              ┌─┴─┐│  cmp_dg ││         │
«conditions_2: ┤ X ├┤         ├┤         ├
«              └───┘│         ││         │
«       sum_0: ─────┤         ├┤         ├
«                   │         ││         │
«       sum_1: ─────┤         ├┤         ├
«                   │         ││         │
«       sum_2: ─────┤         ├┤         ├
«                   │         ││         │
«      work_0: ─────┤3        ├┤3        ├
«                   └─────────┘└─────────┘
«      work_1: ───────────────────────────
«
«      work_2: ───────────────────────────
«
«      work_3: ───────────────────────────
«
objective qubit index 4

[9]:

# plot exact payoff function
plt.figure(figsize=(7,5))
x = np.linspace(sum(low), sum(high))
y = (x <= 5)*np.minimum(np.maximum(0, x - strike_price_1), strike_price_2 - strike_price_1)
plt.plot(x, y, 'r-')
plt.grid()
plt.title('Payoff Function (for $S_1 = S_2$)', size=15)
plt.xlabel('Sum of Spot Prices ($S_1 + S_2)$', size=15)
plt.ylabel('Payoff', size=15)
plt.xticks(size=15, rotation=90)
plt.yticks(size=15)
plt.show()
../../_images/tutorials_finance_07_asian_barrier_spread_pricing_12_0.png 
[10]:

# plot contour of payoff function with respect to both time steps, including barrier
plt.figure(figsize=(7,5))
z = np.zeros((17, 17))
x = np.linspace(low[0], high[0], 17)
y = np.linspace(low[1], high[1], 17)
for i, x_ in enumerate(x):
    for j, y_ in enumerate(y):
        z[i, j] = np.minimum(np.maximum(0, x_ + y_ - strike_price_1), strike_price_2 - strike_price_1)
        if x_ > barrier or y_ > barrier:
            z[i, j] = 0

plt.title('Payoff Function', size=15)
plt.contourf(x, y, z)
plt.colorbar()
plt.xlabel('Spot Price $S_1$', size=15)
plt.ylabel('Spot Price $S_2$', size=15)
plt.xticks(size=15)
plt.yticks(size=15)
plt.show()
../../_images/tutorials_finance_07_asian_barrier_spread_pricing_13_0.png 
[11]:

# evaluate exact expected value
sum_values = np.sum(u.values, axis=1)
payoff = np.minimum(np.maximum(sum_values - strike_price_1, 0), strike_price_2 - strike_price_1)
leq_barrier = [ np.max(v) <= barrier for v in u.values ]
exact_value = np.dot(u.probabilities[leq_barrier], payoff[leq_barrier])
print('exact expected value:\t%.4f' % exact_value)

exact expected value:   0.8023

期待ペイオフの評価

我々は、まず量子回路を検証します。シミュレートで目的量子ビットの \(|1\rangle\) 状態の観測確率を分析します。

[12]:

num_state_qubits = asian_barrier_spread.num_qubits - asian_barrier_spread.num_ancillas
print('state qubits: ', num_state_qubits)
transpiled = transpile(asian_barrier_spread, basis_gates=['u', 'cx'])
print('circuit width:', transpiled.width())
print('circuit depth:', transpiled.depth())

state qubits:  5
circuit width: 15
circuit depth: 8541

[13]:

job = execute(asian_barrier_spread, backend=Aer.get_backend('statevector_simulator'))

[14]:

# evaluate resulting statevector
value = 0
for i, a in enumerate(job.result().get_statevector()):
    b = ('{0:0%sb}' % num_state_qubits).format(i)[-num_state_qubits:]
    prob = np.abs(a)**2
    if prob > 1e-4 and b[0] == '1':
        value += prob
    # all other states should have zero probability due to ancilla qubits
    if i > 2**num_state_qubits:
        break

# map value to original range
mapped_value = objective.post_processing(value) / (2**num_uncertainty_qubits - 1) * (high_ - low_)
print('Exact Operator Value:  %.4f' % value)
print('Mapped Operator value: %.4f' % mapped_value)
print('Exact Expected Payoff: %.4f' % exact_value)

Exact Operator Value:  0.6303
Mapped Operator value: 0.8319
Exact Expected Payoff: 0.8023

次に振幅推定で期待ペイオフを推定します。多数の量子ビットをシミュレートするため、推定には時間がかかることに注意してください。我々が設計した演算子 (asian_barrier_spread) は実際の状態量子ビット数よりもかなり少なくなっていますが、これはシミュレーション時間を短縮するためです。

[15]:

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

# construct amplitude estimation
ae = IterativeAmplitudeEstimation(epsilon=epsilon, alpha=alpha,
                                  state_preparation=asian_barrier_spread,
                                  objective_qubits=[objective_index],
                                  post_processing=objective.post_processing)

[16]:

result = ae.run(quantum_instance=Aer.get_backend('qasm_simulator'), shots=100)

[17]:

conf_int = np.array(result['confidence_interval']) / (2**num_uncertainty_qubits - 1) * (high_ - low_)
print('Exact value:    \t%.4f' % exact_value)
print('Estimated value:\t%.4f' % (result['estimation'] / (2**num_uncertainty_qubits - 1) * (high_ - low_)))
print('Confidence interval: \t[%.4f, %.4f]' % tuple(conf_int))

Exact value:            0.8023
Estimated value:        0.8299
Confidence interval:    [0.8243, 0.8355]

[18]:

import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright

Version Information

Qiskit SoftwareVersion
QiskitNone
Terra0.17.0.dev0+8bac9c1
Aer0.6.1
Ignis0.5.0.dev0+470d8cc
Aqua0.8.0.dev0+ce81016
IBM Q Provider0.8.0
System information
Python3.7.7 (default, May 6 2020, 04:59:01) [Clang 4.0.1 (tags/RELEASE_401/final)]
OSDarwin
CPUs2
Memory (Gb)16.0
Fri Oct 16 13:15:57 2020 CEST

This code is a part of Qiskit

© Copyright IBM 2017, 2020.

This code is licensed under the Apache License, Version 2.0. You may
obtain a copy of this license in the LICENSE.txt file in the root directory
of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.

Any modifications or derivative works of this code must retain this
copyright notice, and modified files need to carry a notice indicating
that they have been altered from the originals.

[ ]: