Note

Run interactively in jupyter notebook.

Pricing European Call Options¶

Introduction¶

Suppose a European call option with strike price $$K$$ and an underlying asset whose spot price at maturity $$S_T$$ follows a given random distribution. The corresponding payoff function is defined as:

$\max\{S_T - K, 0\}$

In the following, a quantum algorithm based on amplitude estimation is used to estimate the expected payoff, i.e., the fair price before discounting, for the option:

$\mathbb{E}\left[ \max\{S_T - K, 0\} \right]$

as well as the corresponding $$\Delta$$, i.e., the derivative of the option price with respect to the spot price, defined as:

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

The approximation of the objective function and a general introduction to option pricing and risk analysis on quantum computers are given in the following papers:

• Quantum Risk Analysis. Woerner, Egger. 2018.

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

[1]:

import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np

from qiskit import Aer, QuantumCircuit
from qiskit.utils import QuantumInstance
from qiskit.algorithms import IterativeAmplitudeEstimation, EstimationProblem
from qiskit.circuit.library import LinearAmplitudeFunction
from qiskit_finance.circuit.library import LogNormalDistribution


Uncertainty Model¶

We construct a circuit factory to load a log-normal random distribution into a quantum state. The distribution is truncated to a given interval $$[\text{low}, \text{high}]$$ and discretized using $$2^n$$ grid points, where $$n$$ denotes the number of qubits used. The unitary operator corresponding to the circuit factory implements the following:

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

where $$p_i$$ denote the probabilities corresponding to the truncated and discretized distribution and where $$i$$ is mapped to the right interval using the affine map:

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

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

[3]:

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


Payoff Function¶

The payoff function equals zero as long as the spot price at maturity $$S_T$$ is less than the strike price $$K$$ and then increases linearly. The implementation uses a comparator, that flips an ancilla qubit from $$\big|0\rangle$$ to $$\big|1\rangle$$ if $$S_T \geq K$$, and this ancilla is used to control the linear part of the payoff function.

The linear part itself is then approximated as follows. We exploit the fact that $$\sin^2(y + \pi/4) \approx y + 1/2$$ for small $$|y|$$. Thus, for a given approximation rescaling factor $$c_\text{approx} \in [0, 1]$$ and $$x \in [0, 1]$$ we consider

$\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}$$.

We can easily construct an operator that acts as

$\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),$

using controlled Y-rotations.

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.

For more details on the approximation, we refer to: Quantum Risk Analysis. Woerner, Egger. 2018.

[4]:

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

[4]:

     ┌───────┐┌────┐
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   ├
└────┘
[5]:

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

[6]:

# 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


Evaluate Expected Payoff¶

[7]:

european_call.draw()

[7]:

     ┌───────┐┌────┐
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   ├
└────┘
[8]:

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

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

[9]:

result = ae.estimate(problem)

[10]:

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

Exact value:            0.1623
Estimated value:        0.1722
Confidence interval:    [0.1668, 0.1776]


Instead of constructing these circuits manually, Qiskit’s finance module offers the EuropeanCallPricing circuit, which already implements this functionality as building block.

[11]:

from qiskit_finance.applications.estimation import EuropeanCallPricing

european_call_pricing = EuropeanCallPricing(num_state_qubits=num_uncertainty_qubits,
strike_price=strike_price,
rescaling_factor=c_approx,
bounds=(low, high),
uncertainty_model=uncertainty_model)

[12]:

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

qi = QuantumInstance(Aer.get_backend('qasm_simulator'), shots=100)
problem = european_call_pricing.to_estimation_problem()
# construct amplitude estimation
ae = IterativeAmplitudeEstimation(epsilon, alpha=alpha, quantum_instance=qi)
result = ae.estimate(problem)

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

Exact value:            0.1623
Estimated value:        0.1679
Confidence interval:    [0.1628, 0.1731]


Evaluate Delta¶

The Delta is a bit simpler to evaluate than the expected payoff. Similarly to the expected payoff, we use a comparator circuit and an ancilla qubit to identify the cases where $$S_T > K$$. However, since we are only interested in the probability of this condition being true, we can directly use this ancilla qubit as the objective qubit in amplitude estimation without any further approximation.

[13]:

from qiskit_finance.applications.estimation import EuropeanCallDelta

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

[14]:

european_call_delta._objective.decompose().draw()

[14]:

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

european_call_delta_circ = QuantumCircuit(european_call_delta._objective.num_qubits)
european_call_delta_circ.append(uncertainty_model, range(num_uncertainty_qubits))
european_call_delta_circ.append(european_call_delta._objective, range(european_call_delta._objective.num_qubits))

european_call_delta_circ.draw()

[15]:

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

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

qi = QuantumInstance(Aer.get_backend('qasm_simulator'), shots=100)
problem = european_call_delta.to_estimation_problem()

# construct amplitude estimation
ae_delta = IterativeAmplitudeEstimation(epsilon, alpha=alpha, quantum_instance=qi)

[17]:

result_delta = ae_delta.estimate(problem)

[18]:

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

Exact delta:            0.8098
Esimated value:         0.8101
Confidence interval:    [0.8054, 0.8148]

[19]:

import qiskit.tools.jupyter
%qiskit_version_table


Version Information

Qiskit SoftwareVersion
QiskitNone
Terra0.17.0
Aer0.8.0
IgnisNone
AquaNone
IBM Q ProviderNone
System information
Python3.8.8 (default, Feb 19 2021, 19:42:00) [GCC 9.3.0]
OSLinux
CPUs2
Memory (Gb)6.791343688964844
Fri Apr 02 21:33:18 2021 UTC

This code is a part of Qiskit

[ ]: