Note

Cette page a été générée à partir de tutorials/finance/06_basket_option_pricing.ipynb.

Exécuter en mode interactif dans le IBM Quantum lab.

Calcul du prix d’options panier

Introduction

Supposons une option panier avec prix d’exercice \(K\) et deux actifs sous-jacents dont les prix au comptant à l’échéance \(S_T^1\), \(S_T^2\) suivent des distributions aléatoires données. La fonction de gain correspondante est définie par :

\[\max\{S_T^1 + S_T^2 - K, 0\}\]

Dans ce qui suit, un algorithme quantique basé sur l’estimation d’amplitude est utilisé pour évaluer le gain attendu, c’est-à-dire le juste prix avant remise, pour l’option :

\[\mathbb{E}\left[ \max\{S_T^1 + S_T^2 - K, 0\} \right].\]

L’approximation de la fonction objectif ainsi qu’une introduction générale à la tarification des options et à l’analyse des risques sur les ordinateurs quantiques sont données dans les articles suivants :

• Quantum Risk Analysis. Woerner, Egger. 2018.

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

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 Aer, QuantumRegister, QuantumCircuit, execute, AncillaRegister, transpile
from qiskit.aqua.algorithms import IterativeAmplitudeEstimation
from qiskit.circuit.library import WeightedAdder, LogNormalDistribution, LinearAmplitudeFunction

Modèle d’Incertitude

Nous construisons un circuit afin de représenter une distribution aléatoire en log-normale multidimensionnelle dans un état quantique sur \(n\) qubits. Pour chaque dimension \(j = 1,\ldots,d\), la distribution est limitée à l’intervalle \([\text{low}_j, \text{high}_j]\) et discrétisée en utilisants \(2^{n_j}\) points, où \(n_j\) correspond au nombre de qubits utilisés pour représenter la dimension \(j\), c’est à dire : \(n_1+\ldots+n_d = n\). L’opérateur unitaire correspondant à ce circuit respecte la relation suivante :

\[\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},\]

où \(p_{i_1\ldots i_d}\) représentent les probabilités de chaque distribution tronquée et discrétisée et où \(i_j\) est mappé sur l’intervalle approprié en utilisant la carte affine suivante :

\[\{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].\]

Pour simplifier, nous supposons que les cours des deux actions sont indépendants et répartis de manière identique. Cette hypothèse simplifie la paramétrisation ci-dessous et peut être facilement étendue à des distributions multivariées plus complexes et également corrélées. La seule hypothèse importante pour l’implémentation suivante est que la grille de discrétisation des différentes dimensions a la même taille de pas.

# 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)))

# 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_T^1$ (\$)', size=15)
ax.set_ylabel('Spot Price $S_T^2$ (\$)', size=15)
ax.set_zlabel('Probability (\%)', size=15)
plt.show()

Fonction de Gain

La fonction de gain est nulle tant que la somme des prix au comptant à l’échéance \((S_T^1 + S_T^2)\) est inférieure au prix d’exercice \(K\) puis cette fonction augmente linéairement. L’implémentation utilise d’abord un opérateur de somme pondérée pour calculer la somme des prix au comptant dans un registre ancillaire, puis utilise un comparateur qui fait basculer un qubit ancillaire de \(\big|0\rangle\) à \(\big|1\rangle\) si \((S_T^1 + S_T^2) \geq K\). Ce registre ancillaire est utilisé pour contrôler la partie linéaire de la fonction de gain.

Une approximation de la partie linéaire proprement dite est faite comme indiqué ci-dessous. On exploite le fait que \(\sin^2(y + \pi/4) \approx y + 1/2\) pour \(|y|\) petit. Ainsi, pour un facteur d’échelle d’approximation \(c_\text{approx} \in [0, 1]\) et \(x \in [0, 1]\) on considère

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

pour de petites valeurs de \(c_\text{approx}\).

Nous pouvons facilement construire un opérateur qui agit ainsi

\[\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),\]

en utilisant des rotations-Y contrôlées.

Il est aussi intéressant de connaitre la probabilité de mesurer \(\big|1\rangle\) sur le dernier qubit, correspondant à \(\sin^2(a*x+b)\). En complément de l’approximation précédente, nous pouvons alors estimer la valeurs des intérêts. Plus \(c_\text{approx}\) est choisi petit, meilleure est l’approximation. Cependant, puisque nous estimons une propriété dépendant de \(c_\text{approx}\), le nombre de qubits d’évaluation \(m\) doit être choisi en conséquence.

Pour plus de détails sur l’approximation, nous nous référons à: Quantum Risk Analysis. Woerner, Egger. 2018.

Étant donné que l’opérateur de somme pondérée (dans son implémentation actuelle) ne peut additionner que des entiers, nous devons mapper depuis les intervalles d’origine vers les intervalles représentables pour estimer le résultat, et inverser cette correspondance avant d’interpréter le résultat. La cartographie correspond essentiellement à la cartographie affine décrite dans le cadre du modèle d’incertitude ci-dessus.

# 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

# set the strike price (should be within the low and the high value of the uncertainty)
strike_price = 3.5

# map strike price from [low, high] to {0, ..., 2^n-1}
max_value = 2**n_s - 1
low_ = low[0]
high_ = high[0]
mapped_strike_price = (strike_price - dimension*low_) / (high_ - low_) * (2**num_uncertainty_qubits - 1)

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

# setup piecewise linear objective fcuntion
breakpoints = [0, mapped_strike_price]
slopes = [0, 1]
offsets = [0, 0]
f_min = 0
f_max = 2*(2**num_uncertainty_qubits - 1) - mapped_strike_price
basket_objective = LinearAmplitudeFunction(
    n_s,
    slopes,
    offsets,
    domain=(0, max_value),
    image=(f_min, f_max),
    rescaling_factor=c_approx,
    breakpoints=breakpoints
)

         ┌───────┐┌────────┐
state_0: ┤0      ├┤0       ├──────
         │       ││        │
state_1: ┤1      ├┤1       ├──────
         │  P(X) ││        │
state_2: ┤2      ├┤2       ├──────
         │       ││        │
state_3: ┤3      ├┤3       ├──────
         └───────┘│        │┌────┐
  obj_0: ─────────┤        ├┤3   ├
                  │        ││    │
  sum_0: ─────────┤4 adder ├┤0   ├
                  │        ││    │
  sum_1: ─────────┤5       ├┤1   ├
                  │        ││    │
  sum_2: ─────────┤6       ├┤2 F ├
                  │        ││    │
 work_0: ─────────┤7       ├┤4   ├
                  │        ││    │
 work_1: ─────────┤8       ├┤5   ├
                  │        ││    │
 work_2: ─────────┤9       ├┤6   ├
                  └────────┘└────┘
objective qubit index 4

# 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 = AncillaRegister(max(n_aux, basket_objective.num_ancillas), 'work')  # additional qubits

objective_index = u.num_qubits

basket_option = QuantumCircuit(qr_state, qr_obj, ar_sum, ar)
basket_option.append(u, qr_state)
basket_option.append(agg, qr_state[:] + ar_sum[:] + ar[:n_aux])
basket_option.append(basket_objective, ar_sum[:] + qr_obj[:] + ar[:basket_objective.num_ancillas])

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

# plot exact payoff function (evaluated on the grid of the uncertainty model)
x = np.linspace(sum(low), sum(high))
y = np.maximum(0, x - strike_price)
plt.plot(x, y, 'r-')
plt.grid()
plt.title('Payoff Function', size=15)
plt.xlabel('Sum of Spot Prices ($S_T^1 + S_T^2)$', size=15)
plt.ylabel('Payoff', size=15)
plt.xticks(size=15, rotation=90)
plt.yticks(size=15)
plt.show()

# evaluate exact expected value
sum_values = np.sum(u.values, axis=1)
exact_value = np.dot(u.probabilities[sum_values>= strike_price], sum_values[sum_values>= strike_price]-strike_price)
print('exact expected value:\t%.4f' % exact_value)

exact expected value:   0.4870

Evaluation du Gain Attendu

Nous commençons par vérifier le circuit quantique en le simulant et en analysant la probabilité de mesurer l’état \(|1\rangle\) sur le qubit objectif.

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

state qubits:  5
circuit width: 11
circuit depth: 529

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

# 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

# map value to original range
mapped_value = basket_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.3954
Mapped Operator value: 0.4969
Exact Expected Payoff: 0.4870

Puis nous utilisons une estimation d’amplitude pour déterminer le gain attendu.

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

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

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

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.4870
Estimated value:        0.4726
Confidence interval:    [0.4238, 0.5213]

import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright

Version Information

Qiskit Software	Version
Qiskit	None
Terra	0.16.0.dev0+28d8c6a
Aer	0.6.1
Ignis	0.5.0.dev0+470d8cc
Aqua	0.8.0.dev0+ce81016
IBM Q Provider	0.8.0
System information
Python	3.7.7 (default, May 6 2020, 04:59:01) [Clang 4.0.1 (tags/RELEASE_401/final)]
OS	Darwin
CPUs	2
Memory (Gb)	16.0
Fri Oct 16 11:48:53 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.