Note

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

Exécuter en mode interactif dans le IBM Quantum lab <https://quantum-computing.ibm.com/jupyter/tutorial/finance/07_asian_barrier_spread_pricing.ipynb> _.

Evaluation des options asiatiques sur spread

Introduction

Une option asiatique sur spread est une combinaison de 3 types d’options différents, et en tant que tel, combine plusieurs fonctionnalités que le framework d’évaluation des prix des options Qiskit Finance supporte :

• Option asiatique: Le paiement dépend du prix moyen durant une période déterminée.

• Option à barrière: le paiement est égal à zéro si un certain seuil est dépassé à n” importe quel moment de la période de temps considérée.

• Option ( haussières ) sur spread: Le paiement suit une fonction linéaire par morceaux (selon le prix moyen) partant de zéro, augmentant linéairement, demeurant constants.

Supposons les prix d’exercice \(K_1 < K_2\) et les périodes \(t=1,2\), avec les cours de sous-jacents correspondants \((S_1, S_ 2)\) suivant une distribution multivariée donnée (générée par exemple par un processus stochastique), et un seuil de barrière \(B> 0\). La fonction de paiement correspondante est définie comme

\[\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{ si } S_1, S_2 \leq B \\ 0, & \text{sinon.} \end{cases}\end{split}\]

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

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

L’approximation de la fonction objectif et une introduction générale à l’évaluation des prix d’options et à l’analyse des risques sur les ordinateurs quantiques sont donnés 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 QuantumRegister, QuantumCircuit, Aer, execute, AncillaRegister, transpile
from qiskit.circuit.library import IntegerComparator, WeightedAdder, LogNormalDistribution, LinearAmplitudeFunction
from qiskit.aqua.algorithms import IterativeAmplitudeEstimation

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 est le relation suivant :

\[\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 relation 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 titres sont indépendants et distribués de manière identique. Cette hypothèse simplifie la paramétrisation ci-dessous et peut être facilement étendue à des distributions plus complexes multivariées 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 ait 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_1$ (\$)', size=15)
ax.set_ylabel('Spot Price $S_2$ (\$)', size=15)
ax.set_zlabel('Probability (\%)', size=15)
plt.show()

Fonction de Paiement

Pour simplifier, nous considérons la somme des cours des sous-jacents au lieu de leur moyenne. Le résultat peut être transformé en moyenne en le divisant par 2.

La fonction de paiement est égale à zéro tant que la somme des cours des sous-jacents :math:` (S_1 + S_ 2) ` est inférieure au prix d’exercice :math:` K_1 , puis augmente de façon linéaire jusqu’à ce que cette somme atteigne :math: K_2 . Alors, le paiement reste constant à :math: K_2-K_1 ` à moins que l’un des deux cours dépasse le seuil barrière :math:` B , et le paiement va immédiatement à zéro. L’implémentation utilise d’abord un opérateur de somme pondérée pour calculer la somme des prix spot dans un registre ancilla, puis utilise un comparateur, qui flip un qubit ancilla de :math: big | 0rangle ` à :math:` big | 1rangle ` si :math:` (S_1 + S_ 2) geq K_1 ` et un autre comparateur / ancilla pour capturer le cas ou :math:` (S_1 + S_ 2) geq K_2 `. Ces ancillas sont utilisées pour contrôler la partie linéaire de la fonction de paiement.

En outre, nous ajoutons une autre variable ancilla pour chaque pas de temps et utilisons des comparateurs supplémentaires pour vérifier si :math:` S_1 , respectivement :math: S_2 , dépasse le seuil barrière :math: B . La fonction de paiement n’est appliquée que si :math: S_1, S_2 leq B `.

Une approximation de la partie linéaire proprement dite est faite de la façon suivante. 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 Y-rotations contrôlées.

Il est aussi intéressant de connaitre la probabilité de mesurer \(\big|1\rangle\) sur le dernier qubit, qui correspond à \(\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 faire correspondre les intervalles d’origine avec les intervalles représentables pour estimer le résultat, puis inverser cette correspondance avant d’interpréter le résultat. Cela correspond essentiellement à la relation 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_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)

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

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

# 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

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

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

# 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

Evaluation du Paiement 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 = 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

job = execute(asian_barrier_spread, 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
    # 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

Nous utilisons ensuite l’estimation de l’amplitude pour obtenir une estimation du gain attendu. Notez que cela peut prendre un certain temps car nous simulons un grand nombre de qubits. La façon dont nous avons conçu l’opérateur (asian_barrier_spread) implique que le nombre de qubits utilisé est significativement plus petit, contribuant ainsi à réduire un peu le temps de simulation global.

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

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.8023
Estimated value:        0.8299
Confidence interval:    [0.8243, 0.8355]

import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright

Version Information

Qiskit Software	Version
Qiskit	None
Terra	0.17.0.dev0+8bac9c1
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 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.