Skip to content

Rainbow options with bruteforce methodology

View on GitHub Experiment in the IDE

In this Notebook we will go through the implementation using QMOD for the rainbow option. This Notebook role is to verify result of different metodology on a smal scale problem, as it grows exponentially in the gate count.


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


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)

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)}\left|w_i\right\rangle\]
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]


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(, [grid_points[-1]] * 2) + MU)):
            "<font color='red'> K always greater than the maximum asset values. Stop the run, the payoff is 0</font>"

Maximum Computation

Precision utils


def round_factor(a):
    precision_factor = 2 ** (FRAC_PLACES)
    return np.floor(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 = \mu_i + y_i\). Where:

\(\mu_i= (t-t_0)\tilde{\mu}_i\), being \(\tilde{\mu}_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 \(\mathbf{L}\) and the standard multivariate Gaussian sample:

\[y_i = \Delta x \cdot \sum_kl_{ik}d_k + x_{min} \cdot \sum_k l_{ik}\]

\(\Delta x\) is the Gaussian discretization step, \(x_{min}\) 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_{ik}\) is the \(i,k\) entry of the matrix \(\mathbf{L}\), defined as \(\mathbf{L}=\mathbf{C}\sqrt{(t-t_0)}\), where \(\mathbf{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 import Output, QNum, get_expression_numeric_attributes, qfunc
from classiq.qmod.symbolic import max as qmax

b = np.log(S0[0]) + MU[0] + MIN_X * CHOLESKY[0].sum()

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]

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

c = round_factor(c)

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

Brute-Force Amplitude Loading Method

This type of amplitude loading has an exponential scale, is therefore used as a "sanity check" method for validating result from the direct method and integration method that are part of the paper [1]

from classiq import (
from classiq.qmod.symbolic import exp, sqrt

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

def get_payoff_expression_normalized(x: QNum, size, fraction_digits):
    x_max = 1 - 1 / (2**size)
    payoff_max = get_payoff_expression(x_max, size, fraction_digits)
    payoff = get_payoff_expression(x, size, fraction_digits)
    return payoff / payoff_max

def brute_force_payoff(max_reg: QNum, ind_reg: QBit):
    max_reg_fixed = QNum("max_reg_fixed", max_reg.size, False, max_reg.size)
    bind(max_reg, max_reg_fixed)
    ind_reg *= get_payoff_expression_normalized(
        max_reg_fixed, max_reg.size, max_reg.fraction_digits
    bind(max_reg_fixed, max_reg)
def rainbow_brute_force(
    x1: QNum,
    x2: QNum,
    ind_reg: QBit,
) -> None:
    inplace_prepare_state(probabilities, 0, x1)
    inplace_prepare_state(probabilities, 0, x2)

    max_out = QNum("max_out")
        lambda: affine_max(x1, x2, max_out),
        lambda: brute_force_payoff(max_out, ind_reg),

def main(
    x1: Output[QNum],
    x2: Output[QNum],
    ind_reg: Output[QBit],
) -> None:
    allocate(NUM_QUBITS, x1)
    allocate(NUM_QUBITS, x2)
    allocate(1, ind_reg)
    rainbow_brute_force(x1, x2, ind_reg)

constraints = Constraints(max_width=23)
qmod = create_model(main, constraints=constraints)
qprog = synthesize(qmod)

IQAE algorithm

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

def qmci_oracle(ind: QBit):

def cmain():
    iqae_res = iqae(epsilon=0.01, alpha=0.05)
    save({"iqae_res": iqae_res})
from classiq import CInt, QArray, QCallable, grover_operator, power

def grover_algorithm(
    k: CInt,
    oracle_operand: QCallable[QArray[QBit]],
    sp_operand: QCallable[QArray[QBit]],
    x: QArray[QBit],
    power(k, lambda: grover_operator(oracle_operand, sp_operand, x))
def get_main():
    def main(
        k: CInt,
        ind_reg: Output[QBit],
    ) -> None:
        full_reg = QArray("full_reg")
        allocate(2 * NUM_QUBITS + 1, full_reg)
            lambda x: qmci_oracle(x[x.len - 1]),
            lambda x: rainbow_brute_force(
                x[NUM_QUBITS : 2 * NUM_QUBITS],
                x[x.len - 1],
        state_reg = QArray("state_reg")
        bind(full_reg, [state_reg, ind_reg])

    return main
from classiq import execute, write_qmod
from classiq.execution import ExecutionPreferences

def synthesize_and_execute(post_process):
    constraints = Constraints(max_width=25)
    qmod = create_model(
    write_qmod(qmod, "rainbow_options_bruteforce_method")
    print("Starting synthesis")
    qprog = synthesize(qmod)
    print("Starting execution")
    res = execute(qprog).result()
    iqae_res = res[0].value
    parsed_res = post_process(res[0].value)

    return (qmod, qprog, iqae_res, parsed_res)

Post Process

We need to add to the post-processing function a term: \begin{equation} \begin{split} &\mathbb{E} \left[\max\left(e^{b \cdot z}, Ke^{-b'}\right) \right] e^{b'} - K \ = &\mathbb{E} \left[\max\left(e^{-a\hat{x}}, Ke^{-b'-ax_{max}}\right) \right]e^{b'+ ax_{max}} - K \end{split} \end{equation}

import sympy

def parse_result_bruteforce(iqae_res):
    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))

    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

qmod_brute_force, qprog_brute_force, iqae_res_brute_force, parsed_res_brute_force = (
result, conf_interval = parsed_res_brute_force
    f"raw iqae results: {iqae_res_brute_force.estimation} with confidence interval {iqae_res_brute_force.confidence_interval}"
print(f"option estimated value: {result} with confidence interval {conf_interval}")
Starting synthesis
Starting execution
raw iqae results: 0.5057390928129042 with confidence interval (0.5035438550940906, 0.5079343305317179)
option estimated value: 23.68809763877735 with confidence interval [22.76055184 24.61564344]
assert -5 <= result <= 45


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