Nota

Esta página foi gerada de tutorials/finance/07_asian_barrier_spread_pricing.ipynb.

Execute interativamente no IBM Quantum lab.

Precificando Spread asiático de barreira

Introdução

Um spread asiático de barreira é uma combinação de 3 tipos de opção diferentes, e como tal, combina múltiplos recursos possíveis, que o framework de precificação de opção Qiskit Finance suporta:

• Opção asiática: O payoff depende do preço médio sobre o horizonte de tempo considerado.

• Opção de Barreira: O payoff é zero, se um determinado limiar for ultrapassado, a qualquer momento, dentro do horizonte de tempo considerado.

• (Bull) Spread: O payoff segue uma função linear por partes (dependendo do preço médio) começando em zero, aumentando linearmente, permanecendo constante.

Suponhamos os preços de exercício \(K_1 < K_2\) e os períodos de tempo \(t=1,2\), com os preços à vista correspondentes \((S_1, S_2)\), seguindo uma dada distribuição multivariada (por exemplo, gerada por algum processo estocástico), e um limiar de barreira \(B>0\). A função payoff correspondente é definida como

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

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[ P(S_1, S_2) \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 QuantumRegister, QuantumCircuit, Aer, execute, AncillaRegister, transpile
from qiskit.circuit.library import IntegerComparator, WeightedAdder, LogNormalDistribution, LinearAmplitudeFunction
from qiskit.aqua.algorithms import IterativeAmplitudeEstimation

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_1$ (\$)', size=15)
ax.set_ylabel('Spot Price $S_2$ (\$)', size=15)
ax.set_zlabel('Probability (\%)', size=15)
plt.show()

Função Payoff

Para simplificar, consideramos a soma dos preços à vista, em vez de sua média. O resultado pode ser transformado na média dividindo-o por 2.

A função payoff é igual a zero, enquanto a soma dos preços à vista \((S_1 + S_2)\) é inferior ao preço de exercício \(K_1\) e, em seguida, aumenta, linearmente, até que a soma dos preços, à vista, chegue a \(K_2\). Em seguida, o payoff permanece constante em \(K_2 - K_1\), a menos que qualquer um dos dois preços, à vista, exceda o limiar de barreira \(B\). Neste caso, o payoff iria diretamente para zero. A implementação, inicialmente, utiliza um operador de soma ponderada para computar a soma dos preços, à vista, em um registrador auxiliar, e depois, usa um comparador, que inverte um qubit auxiliar de \(\big|0\rangle\) para \(\big|1\rangle\) se \((S_1 + S_2) \geq K_1\) e outro comparador/auxiliar para tratar o caso que \((S_1 + S_2) \geq K_2\). Estes auxiliares são utilizados para controlar a parte linear da função payoff.

Além disso, adicionamos outra variável auxiliar para cada etapa de tempo e utilizamos comparadores adicionais para verificar se \(S_1\), respectivamente, \(S_2\), excedem o limiar de barreira \(B\). A função payoff só é aplicada se \(S_1, S_2 \leq B\).

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

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 = 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

Em seguida, utilizamos estimativa de amplitude para estimar o payoff esperado. Note que isto pode levar um tempo, já que estamos simulando um grande número de qubits. A forma como projetamos o operador (spread asiático de barreira) implica que o número de qubits de estado reais é significativamente menor, assim, ajudando a reduzir um pouco o tempo de simulação geral.

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