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> _.

Pricing Asian Barrier Spreads

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.

Suppose strike prices \(K_1 < K_2\) and time periods \(t=1,2\), with corresponding spot prices \((S_1, S_2)\) following a given multivariate distribution (e.g. generated by some stochastic process), and a barrier threshold \(B>0\). The corresponding payoff function is defined as

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

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

We construct a circuit factory to load a multivariate log-normal random distribution into a quantum state on \(n\) qubits. For every dimension \(j = 1,\ldots,d\), the distribution is truncated to a given interval \([\text{low}_j, \text{high}_j]\) and discretized using \(2^{n_j}\) grid points, where \(n_j\) denotes the number of qubits used to represent dimension \(j\), i.e., \(n_1+\ldots+n_d = n\). The unitary operator corresponding to the circuit factory implements the following:

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

where \(p_{i_1\ldots i_d}\) denote the probabilities corresponding to the truncated and discretized distribution and where \(i_j\) is mapped to the right interval using the affine map:

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

for small \(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.

Eventually, we are interested in the probability of measuring \(\big|1\rangle\) in the last qubit, which corresponds to \(\sin^2(a*x+b)\). Together with the approximation above, this allows to approximate the values of interest. The smaller we choose \(c_\text{approx}\), the better the approximation. However, since we are then estimating a property scaled by \(c_\text{approx}\), the number of evaluation qubits \(m\) needs to be adjusted accordingly.

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

Since the weighted sum operator (in its current implementation) can only sum up integers, we need to map from the original ranges to the representable range to estimate the result, and reverse this mapping before interpreting the result. The mapping essentially corresponds to the affine mapping described in the context of the uncertainty model above.

# 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

We first verify the quantum circuit by simulating it and analyzing the resulting probability to measure the \(|1\rangle\) state in the objective qubit.

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

Next we use amplitude estimation to estimate the expected payoff. Note that this can take a while since we are simulating a large number of qubits. The way we designed the operator (asian_barrier_spread) implies that the number of actual state qubits is significantly smaller, thus, helping to reduce the overall simulation time a bit.

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