Rainbow options with Direct Amplitude Loading

View on GitHub Experiment in the IDE In this Notebook we will go through the implementation of the Direct Amplitude Loading Method for the rainbow option presented in [1]

Introduction

In finance, a crucial aspect of asset pricing pertains to derivatives. Derivatives are contracts whose value is contingent upon another source, known as the underlying. The pricing of options, a specific derivative instrument, involves determining the fair market value (discounted payoff) of contracts affording their holders the right, though not the obligation, to buy (call) or sell (put) one or more underlying assets at a predefined strike price by a specified future expiration date (maturity date). This process relies on mathematical models, considering variables like current asset prices, time to expiration, volatility, and interest rates.

Data Definitions

The problem inputs are:

NUM_QUBITS: the number of qubits representing an underlying asset
NUM_ASSETS: the number of underlying assets
K: the strike price
S0: the arrays of underlying assets prices
dt: the number of days to the maturity date
COV: the covariance matrix that correlate the underlyings
MU_LOG_RET: the array containing the mean of the log return of each underlyings

import numpy as np
import scipy

NUM_QUBITS = 2
NUM_ASSETS = 2

K = 190
S0 = [193.97, 189.12]
dt = 250

COV = np.array([[0.000335, 0.000257], [0.000257, 0.000418]])
MU_LOG_RET = np.array([0.00050963, 0.00062552])

MU = MU_LOG_RET * dt
CHOLESKY = np.linalg.cholesky(COV) * np.sqrt(dt)
SCALING_FACTOR = 1 / CHOLESKY[0, 0]

from classiq import *

EPSILON_VALUE = 0.05
ALPHA_VALUE = 0.1
EPSILON = QConstant("EPSILON", float, EPSILON_VALUE)
ALPHA = QConstant("ALPHA", float, ALPHA_VALUE)

Gaussian State preparation

Encode the probability distribution of a discrete multivariate random variable $W$ taking values in ${w_{0}, . ., w_{N - 1}}$ describing the assets' prices at the maturity date. The number of discretized values, denoted as $N$ , depends on the precision of the state preparation module and is consequently connected to the number of qubits ( $n =$ ) according to the formula $N = 2^{n}$ .

\sum_{i = 0}^{N - 1} \sqrt{p (w_{i})} | w_{i} ⟩

def gaussian_discretization(num_qubits, mu=0, sigma=1, stds_around_mean_to_include=3):
    lower = mu - stds_around_mean_to_include * sigma
    upper = mu + stds_around_mean_to_include * sigma
    num_of_bins = 2**num_qubits
    sample_points = np.linspace(lower, upper, num_of_bins + 1)

    def single_gaussian(x: np.ndarray, _mu: float, _sigma: float) -> np.ndarray:
        cdf = scipy.stats.norm.cdf(x, loc=_mu, scale=_sigma)
        return cdf[1:] - cdf[0:-1]

    non_normalized_pmf = (single_gaussian(sample_points, mu, sigma),)
    real_probs = non_normalized_pmf / np.sum(non_normalized_pmf)
    return sample_points[:-1], real_probs[0].tolist()


grid_points, probabilities = gaussian_discretization(NUM_QUBITS)

STEP_X = grid_points[1] - grid_points[0]
MIN_X = grid_points[0]

SANITY CHECK

The process must be stopped if the strike price $K$ is greater than the maximum value reacheable by the assets during the simulation, to avoid meaningless results. The payoff is $0$ in this case, so there is no need to simulate.

from IPython.display import Markdown

if K >= max(S0 * np.exp(np.dot(CHOLESKY, [grid_points[-1]] * 2) + MU)):
    display(
        Markdown(
            "<font color='red'> K always greater than the maximum asset values. Stop the run, the payoff is 0</font>"
        )
    )

Maximum Computation

Precision utils

FRAC_PLACES = 2


def round_factor(a):
    precision_factor = 2**FRAC_PLACES
    return round(a * precision_factor) / precision_factor


def floor_factor(a):
    precision_factor = 2**FRAC_PLACES
    return np.floor(a * precision_factor) / precision_factor

Affine and maximum arithmetic definitions

Considering the time delta between the starting date ( $t_{0}$ ) and the maturity date ( $t$ ), we can express the return value $R_{i}$ for the $i$ -th asset as $R_{i} = μ_{i} + y_{i}$ . Where:

$μ_{i} = (t - t_{0}) {\tilde{μ}}_{i}$ , being ${\tilde{μ}}_{i}$ the expected daily log-return value. It can be estimated by considering the historical time series of log returns for the $i$ -th asset.

$y_{i}$ is obtained through the dot product between the matrix $L$ and the standard multivariate Gaussian sample:

y_{i} = Δ x \cdot \sum_{k} l_{i k} d_{k} + x_{m i n} \cdot \sum_{k} l_{i k}

$Δ x$ is the Gaussian discretization step, $x_{m i n}$ is the lower Gaussian truncation value and $d_{k} \in [0, 2^{m} - 1]$ is the sample taken from the $k$ -th standard Gaussian. $l_{i k}$ is the $i, k$ entry of the matrix $L$ , defined as $L = C \sqrt{(t - t_{0})}$ , where $C$ is the lower triangular matrix obtained by applying the Cholesky decomposition to the historical daily log-returns correlation matrix.

from functools import reduce

from classiq.qmod.symbolic import max as qmax

a = STEP_X / SCALING_FACTOR
b = np.log(S0[0]) + MU[0] + MIN_X * CHOLESKY[0].sum()
c = (
    SCALING_FACTOR
    * (
        np.log(S0[1])
        + MU[1]
        - (np.log(S0[0]) + MU[0])
        + MIN_X * sum(CHOLESKY[1] - CHOLESKY[0])
    )
    / STEP_X
)
c = round_factor(c)


def get_affine_formula(assets, i):
    return reduce(
        lambda x, y: x + y,
        [
            assets[j] * round_factor(SCALING_FACTOR * CHOLESKY[i, j])
            for j in range(NUM_ASSETS)
            if CHOLESKY[i, j]
        ],
    )


def calculate_max_reg_type():
    x1 = QNum("x1", NUM_QUBITS, False, 0)
    x2 = QNum("x2", NUM_QUBITS, False, 0)
    expr = qmax(get_affine_formula([x1, x2], 0), get_affine_formula([x1, x2], 1) + c)
    size_in_bits, sign, fraction_digits = get_expression_numeric_attributes(
        [x1, x2], expr
    )
    return size_in_bits, fraction_digits


MAX_NUM_QUBITS = calculate_max_reg_type()[0]
MAX_FRAC_PLACES = calculate_max_reg_type()[1]

@qfunc
def affine_max(x1: QNum, x2: QNum, res: Output[QNum]):
    res |= qmax(get_affine_formula([x1, x2], 0), get_affine_formula([x1, x2], 1) + c)

Direct Method

The direct exponential amplitude loading encodes in $\tilde{f}$ the following function: $\tilde{f} (x) = {\begin{cases} e^{- a \hat{x}}, & if \frac{x}{2^{P}} \geq \frac{\log (K) - b^{'}}{b} \\ K e^{- (b^{'} + a x_{m a x})}, & if \frac{x}{2^{P}} < \frac{\log (K) - b^{'}}{b} \end{cases}$

where $\hat{x}$ is the binary complement of $x$ ( $\hat{x} = x - x_{m a x}$ ) and $x_{m a x} = 2^{R} - 1$ , the maximum value that can be stored in $| x ⟩$ register. For loading $e^{- a \hat{x}}$ the $| r ⟩$ is initialized to all zeros and one controlled rotation for each qubit is performed. The rotations angles are: $θ_{i} = 2 \arccos (\sqrt{e^{- a 2^{i}}})$ . All the probabilities of getting a $| 0 ⟩^{\otimes R}$ in the $| r ⟩$ are then collected by a multi-controlled X (MCX) gate and stored in the $| 1 ⟩$ state of a target qubit.

from classiq.qmod.symbolic import acos, asin, exp, sqrt


@qfunc
def exponential_amplitude_loading(
    exp_rate: CReal, x: QArray[QBit], aux: QArray[QBit], res: QBit
) -> None:
    apply_to_all(X, x)
    repeat(
        x.len,
        lambda index: control(
            x[index],
            lambda: RY(2 * acos(1 / sqrt(exp(exp_rate * (2**index)))), aux[index]),
        ),
    )
    apply_to_all(X, x)

    aux_num = QNum("aux_num", aux.len, False, 0)
    bind(aux, aux_num)
    control(aux_num == 0, lambda: X(res))
    bind(aux_num, aux)

class EstimationVars(QStruct):
    x1: QNum[NUM_QUBITS, False, 0]
    x2: QNum[NUM_QUBITS, False, 0]
    aux: QNum[MAX_NUM_QUBITS, False, 0]
    ind: QBit


def get_payoff_expression(x, size, fraction_digits):
    payoff = sqrt(
        qmax(
            S0[0]
            * exp(
                STEP_X / SCALING_FACTOR * (2 ** (size - fraction_digits)) * x
                + (MU[0] + MIN_X * CHOLESKY[0].sum())
            ),
            K,
        )
    )
    return payoff


def get_strike_price_theta_direct(x: QNum):
    x_max = 1 - 1 / (2**x.size)
    payoff_max = get_payoff_expression(x_max, x.size, x.fraction_digits)

    return 2 * asin(np.sqrt(K) / payoff_max)


@qfunc
def direct_load_amplitudes(geq_reg: QNum, max_reg: QNum, aux_reg: QNum, ind_reg: QBit):
    exp_rate = (1 / (2**max_reg.fraction_digits)) * a
    control(
        geq_reg == 1,
        lambda: exponential_amplitude_loading(exp_rate, max_reg, aux_reg, ind_reg),
    )
    strike_price_theta = get_strike_price_theta_direct(max_reg)
    control(geq_reg == 0, lambda: RY(strike_price_theta, ind_reg))


@qfunc
def direct_payoff(max_reg: QNum, aux_reg: QNum, ind_reg: QBit):
    geq_reg = QBit("geq_reg")
    within_apply(
        lambda: asset_geq_strike_price(max_reg, geq_reg),
        lambda: direct_load_amplitudes(geq_reg, max_reg, aux_reg, ind_reg),
    )


@qfunc
def asset_geq_strike_price(
    x: QNum,
    res: Output[QBit],
) -> None:
    a = STEP_X / SCALING_FACTOR
    b = np.log(S0[0]) + MU[0] + MIN_X * CHOLESKY[0].sum()
    COMP_VALUE = (np.log(K) - b) / a
    res |= x > floor_factor(COMP_VALUE)


@qfunc
def rainbow_direct(qvars: EstimationVars) -> None:
    inplace_prepare_state(probabilities, 0, qvars.x1)
    inplace_prepare_state(probabilities, 0, qvars.x2)
    max_out = QNum("max_out")
    within_apply(
        lambda: affine_max(qvars.x1, qvars.x2, max_out),
        lambda: direct_payoff(max_out, qvars.aux, qvars.ind),
    )


@qfunc
def main(qvars: Output[EstimationVars]) -> None:
    allocate(qvars.size, qvars)
    rainbow_direct(qvars)


qmod = create_model(main, constraints=Constraints(max_width=24))
print("Starting synthesis")
qprog = synthesize(qmod)
show(qprog)

Starting synthesis

IQAE algorithm

from classiq.qmod.builtins.classical_execution_primitives import iqae, save


@qfunc
def qmci_oracle(qvars: EstimationVars):
    Z(qvars.ind)


@cfunc
def cmain():
    iqae_res = iqae(epsilon=EPSILON, alpha=ALPHA)
    save({"iqae_res": iqae_res})

@qfunc
def grover_algorithm(
    k: CInt,
    oracle_operand: QCallable[QArray[QBit]],
    sp_operand: QCallable[QArray[QBit]],
    x: QArray[QBit],
):
    sp_operand(x)
    power(k, lambda: grover_operator(oracle_operand, sp_operand, x))

@qfunc
def main(
    k: CInt,
    ind: Output[QBit],
) -> None:
    qvars = EstimationVars("qvars")
    allocate(qvars.size, qvars)
    grover_algorithm(k, qmci_oracle, rainbow_direct, qvars)
    state = QArray("state")
    bind(qvars, [state, ind])

qmod = create_model(
    main,
    constraints=Constraints(max_width=25),
    classical_execution_function=cmain,
    out_file="rainbow_options_direct_method",
)
print("Starting synthesis")
qprog = synthesize(qmod)
show(qprog)
print("Starting execution")
result = execute(qprog).result_value()

Post Process

We need to add to the post-processing function a term: $\begin{aligned} E [max (e^{b \cdot z}, K e^{- b^{'}})] e^{b^{'}} - K \\ = & E [max (e^{- a \hat{x}}, K e^{- b^{'} - a x_{m a x}})] e^{b^{'} + a x_{m a x}} - K \end{aligned}$

import sympy

payoff_expression = f"sqrt(max({S0[0]} * exp({STEP_X / SCALING_FACTOR * (2 ** (MAX_NUM_QUBITS - MAX_FRAC_PLACES))} * x + ({MU[0]+MIN_X*CHOLESKY[0].sum()})), {K}))"
payoff_func = sympy.lambdify(sympy.symbols("x"), payoff_expression)
payoff_max = payoff_func(1 - 1 / (2**MAX_NUM_QUBITS))


def parse_result_direct(iqae_res):
    option_value = iqae_res.estimation * (payoff_max**2) - K
    confidence_interval = np.array(iqae_res.confidence_interval) * (payoff_max**2) - K
    return (option_value, confidence_interval)

Run Method

parsed_result, conf_interval = parse_result_direct(result)
print(
    f"raw iqae results: {result.estimation} with confidence interval {result.confidence_interval}"
)
print(
    f"option estimated value: {parsed_result} with confidence interval {conf_interval}"
)

Starting synthesis
Starting execution
raw iqae results: 0.08014469753569903 with confidence interval (0.07186168028239096, 0.08842771478900711)
option estimated value: 24.26836709887192 with confidence interval [ 2.12356356 46.41317064]

Assertions

expected_payoff = 23.0238
ALPHA_ASSERTION = 1e-5
measured_confidence = conf_interval[1] - conf_interval[0]
confidence_scale_by_alpha = np.sqrt(
    np.log(ALPHA_VALUE / ALPHA_ASSERTION)
)  # based on e^2=(1/2N)*log(2T/alpha) from "Iterative Quantum Amplitude Estimation" since our alpha is low, we want to check within a bigger confidence interval
assert (
    np.abs(parsed_result - expected_payoff)
    <= 0.5 * measured_confidence * confidence_scale_by_alpha
), f"Payoff result is out of the {ALPHA_ASSERTION*100}% confidence interval: |{parsed_result} - {expected_payoff}| > {0.5*measured_confidence * confidence_scale_by_alpha}"

References

[1]: Francesca Cibrario et al., Quantum Amplitude Loading for Rainbow Options Pricing. Preprint