Note

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

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

Prix des options d’achat européennes

Introduction

Supposons une option d’achat européenne avec un prix d’exercice \(K\) et un actif sous-jacent dont le prix au comptant à l’échéance \(S_T\) suit une distribution aléatoire donnée. La fonction de gain correspondante est définie comme suit:

\[\max\{S_T - 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\{K - S_T, 0\} \right]\]

ainsi que le \(\Delta\) résultant, c’est-à-dire la différence entre le prix de l’option et le prix au comptant, défini comme :

\[\Delta = \mathbb{P}\left[S_T \geq K\right]\]

L’approximation de la fonction objectif et 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
%matplotlib inline
import numpy as np

from qiskit import Aer, QuantumCircuit
from qiskit.aqua.algorithms import IterativeAmplitudeEstimation
from qiskit.circuit.library import LogNormalDistribution, LinearAmplitudeFunction

Modèle d’incertitude

Nous construisons un circuit afin de représenter une distribution aléatoire en log-normale dans un état quantique. La distribution est limitée à l’intervalle \([\text{low}, \text{high}]\) et rendue discrète en utilisant \(2^n\) points, où \(n\) correspond au nombre de qubits utilisés. L’opérateur unitaire correspondant à ce circuit respecte la relation suivante :

\[\big|0\rangle_{n} \mapsto \big|\psi\rangle_{n} = \sum_{i=0}^{2^n-1} \sqrt{p_i}\big|i\rangle_{n},\]

où les \(p_i\) représentent les probabilités de la distribution tronquée et discrétisée et où les \(i\) sont mappés sur l’intervalle approprié en utilisant la carte affine suivante :

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

# number of qubits to represent the uncertainty
num_uncertainty_qubits = 3

# 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

# construct A operator for QAE for the payoff function by
# composing the uncertainty model and the objective
uncertainty_model = LogNormalDistribution(num_uncertainty_qubits, mu=mu, sigma=sigma**2, bounds=(low, high))

# plot probability distribution
x = uncertainty_model.values
y = uncertainty_model.probabilities
plt.bar(x, y, width=0.2)
plt.xticks(x, size=15, rotation=90)
plt.yticks(size=15)
plt.grid()
plt.xlabel('Spot Price at Maturity $S_T$ (\$)', size=15)
plt.ylabel('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\) est inférieure au prix d’exercice \(K\) puis cette fonction augmente linéairement. Cette implémentation repose sur un comparateur qui modifie un qubit ancillaire de \(\big|0\rangle\) à \(\big|1\rangle\) si \(S_T \geq K\) , puis cet ancillaire est utilisé pour contrôler la partie linéaire de la fonction de gain.

La partie linéaire proprement dite suit l’approximation ci-dessous. On exploite le fait que \(\sin^2(y + \pi/4) \approx y + 1/2\) lorsque \(|y|\) est 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 rotation 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.

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

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

# setup piecewise linear objective fcuntion
breakpoints = [low, strike_price]
slopes = [0, 1]
offsets = [0, 0]
f_min = 0
f_max = high - strike_price
european_call_objective = LinearAmplitudeFunction(
    num_uncertainty_qubits,
    slopes,
    offsets,
    domain=(low, high),
    image=(f_min, f_max),
    breakpoints=breakpoints,
    rescaling_factor=c_approx
)

# construct A operator for QAE for the payoff function by
# composing the uncertainty model and the objective
num_qubits = european_call_objective.num_qubits
european_call = QuantumCircuit(num_qubits)
european_call.append(uncertainty_model, range(num_uncertainty_qubits))
european_call.append(european_call_objective, range(num_qubits))

# draw the circuit
european_call.draw()

     ┌───────┐┌────┐
q_0: ┤0      ├┤0   ├
     │       ││    │
q_1: ┤1 P(X) ├┤1   ├
     │       ││    │
q_2: ┤2      ├┤2   ├
     └───────┘│    │
q_3: ─────────┤3 F ├
              │    │
q_4: ─────────┤4   ├
              │    │
q_5: ─────────┤5   ├
              │    │
q_6: ─────────┤6   ├
              └────┘

# plot exact payoff function (evaluated on the grid of the uncertainty model)
x = uncertainty_model.values
y = np.maximum(0, x - strike_price)
plt.plot(x, y, 'ro-')
plt.grid()
plt.title('Payoff Function', size=15)
plt.xlabel('Spot Price', size=15)
plt.ylabel('Payoff', size=15)
plt.xticks(x, size=15, rotation=90)
plt.yticks(size=15)
plt.show()

# evaluate exact expected value (normalized to the [0, 1] interval)
exact_value = np.dot(uncertainty_model.probabilities, y)
exact_delta = sum(uncertainty_model.probabilities[x >= strike_price])
print('exact expected value:\t%.4f' % exact_value)
print('exact delta value:   \t%.4f' % exact_delta)

exact expected value:   0.1623
exact delta value:      0.8098

Evaluation du 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=european_call,
                                  objective_qubits=[3],
                                  post_processing=european_call_objective.post_processing)

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

conf_int = np.array(result['confidence_interval'])
print('Exact value:        \t%.4f' % exact_value)
print('Estimated value:    \t%.4f' % (result['estimation']))
print('Confidence interval:\t[%.4f, %.4f]' % tuple(conf_int))

Exact value:            0.1623
Estimated value:        0.1694
Confidence interval:    [0.1633, 0.1755]

Au lieu de construire ces circuits manuellement, le module finance de Qiskit propose le circuit EuropeanCallExpectedValue, qui implémente déjà cette fonctionnalité en tant que bloc de construction.

from qiskit.finance.applications import EuropeanCallExpectedValue

european_call_objective = EuropeanCallExpectedValue(num_uncertainty_qubits,
                                                    strike_price,
                                                    rescaling_factor=c_approx,
                                                    bounds=(low, high))

# append the uncertainty model to the front
european_call = european_call_objective.compose(uncertainty_model, front=True)

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

# construct amplitude estimation
ae = IterativeAmplitudeEstimation(epsilon=epsilon, alpha=alpha,
                                  state_preparation=european_call,
                                  objective_qubits=[3],
                                  post_processing=european_call_objective.post_processing)
result = ae.run(quantum_instance=Aer.get_backend('qasm_simulator'), shots=100)

conf_int = np.array(result['confidence_interval'])
print('Exact value:        \t%.4f' % exact_value)
print('Estimated value:    \t%.4f' % (result['estimation']))
print('Confidence interval:\t[%.4f, %.4f]' % tuple(conf_int))

Exact value:            0.1623
Estimated value:        0.1708
Confidence interval:    [0.1656, 0.1759]

Evaluation du Delta

Le Delta est un peu plus simple à évaluer que la fonction de gain. Tout comme pour la fonction de gain, nous utilisons un circuit comparateur et un qubit ancillaire pour identifier les cas où \(S_T>K\). Cependant, puisque nous ne sommes intéressés que par la probabilité que cette condition soit vraie, nous pouvons utiliser le qubit ancillaire directement comme qubit objectif pour obtenir l’estimation d’amplitude sans approximations supplémentaires.

from qiskit.finance.applications import EuropeanCallDelta

european_call_delta = EuropeanCallDelta(num_uncertainty_qubits, strike_price, bounds=(low, high))

european_call_delta.decompose().draw()

state_0: ───────■─────────────────────────────■──
                │                             │
state_1: ───────┼────■───────────────────■────┼──
         ┌───┐  │    │            ┌───┐  │    │
state_2: ┤ X ├──┼────┼─────────■──┤ X ├──┼────┼──
         ├───┤  │    │       ┌─┴─┐└───┘  │    │
state_3: ┤ X ├──┼────┼───────┤ X ├───────┼────┼──
         └───┘┌─┴─┐  │       └─┬─┘       │  ┌─┴─┐
 work_0: ─────┤ X ├──■─────────┼─────────■──┤ X ├
              └───┘┌─┴─┐┌───┐  │  ┌───┐┌─┴─┐└───┘
 work_1: ──────────┤ X ├┤ X ├──■──┤ X ├┤ X ├─────
                   └───┘└───┘     └───┘└───┘     

state_preparation = QuantumCircuit(european_call_delta.num_qubits)
state_preparation.append(uncertainty_model, range(uncertainty_model.num_qubits))
state_preparation.append(european_call_delta, range(european_call_delta.num_qubits))
state_preparation.draw()

     ┌───────┐┌──────┐
q_0: ┤0      ├┤0     ├
     │       ││      │
q_1: ┤1 P(X) ├┤1     ├
     │       ││      │
q_2: ┤2      ├┤2     ├
     └───────┘│  ECD │
q_3: ─────────┤3     ├
              │      │
q_4: ─────────┤4     ├
              │      │
q_5: ─────────┤5     ├
              └──────┘

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

# construct amplitude estimation
ae_delta = IterativeAmplitudeEstimation(epsilon=epsilon, alpha=alpha,
                                        state_preparation=state_preparation,
                                        objective_qubits=[num_uncertainty_qubits])

result_delta = ae_delta.run(quantum_instance=Aer.get_backend('qasm_simulator'), shots=100)

conf_int = np.array(result_delta['confidence_interval'])
print('Exact delta:    \t%.4f' % exact_delta)
print('Esimated value: \t%.4f' % result_delta['estimation'])
print('Confidence interval: \t[%.4f, %.4f]' % tuple(conf_int))

Exact delta:            0.8098
Esimated value:         0.8068
Confidence interval:    [0.8013, 0.8123]

import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright

Version Information

Qiskit Software	Version
Qiskit	None
Terra	0.17.0.dev0+4ada179
Aer	0.6.1
Ignis	0.5.0.dev0+470d8cc
Aqua	0.9.0.dev0+5a88b59
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
Tue Oct 20 10:57:12 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.