Nota

Esta página foi gerada, a partir de tutorials/finance/06_basket_option_pricing.ipynb.

Execute interativamente no IBM Quantum lab.

Precificando Opções Basket

Introdução

Suponha uma opção basket com preço de exercício \(K\) e dois ativos subjacentes, cujos preços, à vista, no vencimento \(S_T^1\), \(S_T^2\) seguem uma dada distribuição aleatória. A função payoff correspondente é definida como:

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

A seguir, um algoritmo quântico baseado em estimativa de amplitude é usado para estimar o payoff esperado, ou seja, o preço justo antes de descontar para a opção:

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

A aproximação da função objetivo e uma introdução geral sobre precificação de opções e análise de risco em computadores quânticos são dadas nos seguintes artigos:

• Quantum Risk Analysis. Woerner, Egger. 2018.

• Precificação de Opções utilizando Computadores Quânticos. 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

Modelo de incerteza

Construímos uma fábrica de circuitos para carregar uma distribuição aleatória multivariada log-normal em um estado quântico de \(n\) qubits. Para cada dimensão \(j = 1,\ldots,d\), a distribuição é truncada para um determinado intervalo \([\text{low}_j, \text{high}_j]\) e discretizada usando \(2^{n_j}\) pontos da grade, onde \(n_j\) denota o número de qubits usados para representar a dimensão \(j\), ou seja, \(n_1+\ldots+n_d = n\). O operador unitário correspondente à fábrica de circuitos implementa o seguinte:

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

onde \(p_{i_1\ldots i_d}\) denotam as probabilidades correspondentes à distribuição truncada e discretizada e onde \(i_j\) é mapeado para o intervalo correto utilizando o mapeamento afim:

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

Para simplificar, assumimos que ambos os preços das ações são independentes e identicamente distribuídos. Esta suposição apenas simplifica a parametrização, abaixo, e pode ser facilmente relaxada para distribuições multivariadas mais complexas e também correlacionadas. A única suposição importante para a implementação atual é que a grade de discretização das diferentes dimensões tenham o mesmo tamanho de passo.

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

Função Payoff

A função payoff é igual a zero, enquanto a soma dos preços à vista no vencimento \((S_T^1 + S_T^2)\) é menor do que o preço de exercício \(K\), e, depois aumenta linearmente. A implementação usa inicialmente um operador de soma ponderada para computar a soma dos preços à vista no registrador auxiliar e, então, usa um comparador, que inverte o qubit auxiliar de \(\big|0\rangle\) para \(\big|1\rangle\) se \((S_T^1 + S_T^2) \geq K\). Este auxiliar é usado para controlar a parte linear da função payoff.

A própria parte linear é aproximada da seguinte maneira. Exploramos o fato de que \(\sin^2(y + \pi/4) \approx y + 1/2\) para \(|y|\) pequeno. Assim, para um dado fator de reescalonamento de aproximação \(c_\text{approx} \in [0, 1]\) e \(x \in [0, 1]\) consideramos

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

para \(c_\text{approx}\) pequeno.

Podemos facilmente construir um operador que atua como

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

usando rotações Y controladas.

Eventualmente, estamos interessados na probabilidade de medir \(\big|1\rangle\) no último qubit, o que corresponde a \(\sin^2(a*x+b)\). Juntamente com a aproximação, acima, isto permite aproximar os valores de interesse. Quanto menor o \(c_\text{approx}\) escolhido, melhor a aproximação. No entanto, uma vez que estamos estimando uma propriedade escalonada por \(c_\text{approx}\) em seguida, o número de qubits de avaliação \(m\) precisa ser ajustado de acordo.

Para mais detalhes sobre a aproximação, podemos consultar: Quantum Risk Analysis. Woerner, Egger. 2018.

Uma vez que o operador de soma ponderada (em sua implementação atual) pode apenas somar números inteiros, precisamos mapear dos intervalos originais até o intervalo representável, para estimar o resultado, e reverter este mapeamento, antes de interpretar o resultado. O mapeamento corresponde, essencialmente, ao mapeamento afim, descrito no contexto do modelo de incerteza acima.

# 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

Avaliando o Payoff Esperado

Verificamos, primeiramente, o circuito quântico, simulando-o e analisando a probabilidade resultante, para medir o estado \(|1\rangle\) no qubit objetivo.

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

Em seguida, utilizamos estimativa de amplitude para estimar o payoff esperado.

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