Skip to main content

View on GitHub

Open this notebook in GitHub to run it yourself
This notebook covers the implementation of the rainbow option using Qmod. The role of the notebook is to verify the result of a different methodology on a small 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 that afford their holders the right, although 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 such as current asset prices, time to expiration, volatility, and interest rates.

​
Data Definitions

The problem inputs:
  • 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 asset prices
  • dt: the number of days to the maturity date
  • COV: the covariance matrix that correlates to the underlyings
  • MU_LOG_RET: the array containing the mean of the log return of each underlying
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 = 0.05
ALPHA = 0.1

​
Gaussian State Preparation

Encode the probability distribution of a discrete multivariate random variable WW taking values in {w0,..,wN−1}\{w_0, .., w_{N-1}\} describing the asset prices at the maturity date. The number of discretized values, denoted as NN, depends on the precision of the state preparation module and is consequently connected to the number of qubits (n=n=) according to the formula N=2nN=2^n: ∑i=0N−1p(wi)∣wi⟩\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]

​

​
Sanity Check

To avoid meaningless results, the process must stop if the strike price KK is greater than the maximum value reacheable by the assets during the simulation. In this case, the payoff is 00, 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 = 1


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 (t0t_0) and the maturity date (tt), express the return value RiR_i for the ii-th asset as Ri=μi+yiR_i = \mu_i + y_i where μi=(t−t0)μ~i\mu_i= (t-t_0)\tilde{\mu}_i, being μ~i\tilde{\mu}_i the expected daily log-return value. It can be estimated by considering the historical time series of log returns for the ii-th asset. yiy_i is obtained through the dot product between the matrix L\mathbf{L} and the standard multivariate Gaussian sample: yi=Δx⋅∑klikdk+xmin⋅∑kliky_i = \Delta x \cdot \sum_kl_{ik}d_k + x_{min} \cdot \sum_k l_{ik} Δx\Delta x is the Gaussian discretization step, xminx_{min} is the lower Gaussian truncation value, and dk∈[0,2m−1]d_k \in [0,2^m-1] is the sample taken from the kk-th standard Gaussian. likl_{ik} is the i,ki,k entry of the matrix L\mathbf{L}, defined as L=C(t−t0)\mathbf{L}=\mathbf{C}\sqrt{(t-t_0)}, where C\mathbf{C} is the lower triangular matrix obtained by applying the Cholesky decomposition to the historical daily log-return 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()


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 = (
    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 calculate_max_reg_type():
    x1 = QNum(size=NUM_QUBITS)
    x2 = QNum(size=NUM_QUBITS)
    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]


@qperm
def affine_max(x1: Const[QNum], x2: Const[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, and is therefore used as a “sanity check” method for validating the result from the direct method and integration method that are part of the paper [1]. 
def get_payoff_expression(x, size, fraction_digits):
    payoff = np.sqrt(
        max(
            S0[0]
            * np.exp(
                STEP_X / SCALING_FACTOR * (2 ** (size - fraction_digits)) * x
                + (MU[0] + MIN_X * CHOLESKY[0].sum())
            ),
            K,
        )
    )
    return payoff


def get_payoff_expression_normalized(x, 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


@qfunc
def brute_force_payoff(max_reg: Const[QNum], ind_reg: QBit):
    max_reg_fixed = QNum(
        size=max_reg.size, is_signed=False, fraction_digits=max_reg.size
    )
    bind(max_reg, max_reg_fixed)
    assign_amplitude_table(
        lookup_table(
            lambda n: get_payoff_expression_normalized(
                n, max_reg.size, max_reg.fraction_digits
            ),
            max_reg_fixed,
        ),
        max_reg_fixed,
        ind_reg,
    )
    bind(max_reg_fixed, max_reg)


class EstimationVars(QStruct):
    x1: QNum[NUM_QUBITS]
    x2: QNum[NUM_QUBITS]


@qfunc
def rainbow_brute_force(qvars: EstimationVars, ind: QBit) -> None:
    inplace_prepare_state(probabilities, 0, qvars.x1)
    inplace_prepare_state(probabilities, 0, qvars.x2)

    max_out = QNum()
    affine_max(qvars.x1, qvars.x2, max_out),
    brute_force_payoff(max_out, ind)


@qfunc
def main(qvars: Output[EstimationVars], ind: Output[QBit]) -> None:
    allocate(qvars)
    allocate(ind)
    rainbow_brute_force(qvars, ind)


MAX_WIDTH_1 = 23
qmod_1 = create_model(main)
qmod_1 = update_constraints(qmod_1, max_width=MAX_WIDTH_1)
qprog_1 = synthesize(qmod_1)
show(qprog_1)
Output:
/var/folders/h1/cd4f05zs1bv9k8xcn5y2x17m0000gn/T/ipykernel_11657/1054662268.py:31: ClassiqDeprecationWarning: The '*=' operator is deprecated and will no longer be supported starting on 2025-12-03 at the earliest.

Use the 'assign_amplitude_table' function instead
    ind_reg *= get_payoff_expression_normalized(
Output:

Quantum program link: https://platform.classiq.io/circuit/35601S9I1x9ngoD5kSqo49d60eV

​
Iterative Quantum Amplitude Estimation (IQAE) Algorithm

from classiq.applications.iqae.iqae import IQAE

MAX_WIDTH_2 = 25

iqae = IQAE(
    state_prep_op=rainbow_brute_force,
    problem_vars_size=NUM_QUBITS * NUM_ASSETS,
    constraints=Constraints(max_width=MAX_WIDTH_2),
    preferences=Preferences(optimization_level=1),
)


qmod_2 = iqae.get_model()

print("Starting synthesis")
qprog_2 = iqae.get_qprog()
show(qprog_2)
print("Starting execution")
result = iqae.run(EPSILON, ALPHA)
Output:
/var/folders/h1/cd4f05zs1bv9k8xcn5y2x17m0000gn/T/ipykernel_11657/1054662268.py:31: ClassiqDeprecationWarning: The '*=' operator is deprecated and will no longer be supported starting on 2025-12-03 at the earliest.

Use the 'assign_amplitude_table' function instead
    ind_reg *= get_payoff_expression_normalized(
Output:

Starting synthesis
  Quantum program link: https://platform.classiq.io/circuit/3560IQUU8u8M4IzsByiJOeDmb6N
  Starting execution

​
Post-process

Add a term to the post-processing function: E[max⁡(eb⋅z,Ke−b′)]eb′−K=E[max⁡(e−ax^,Ke−b′−axmax)]eb′+axmax−K\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

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_bruteforce(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_bruteforce(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}"
)
Output:
raw iqae results: 0.5048828125 with confidence interval [0.4732664268687129, 0.5364991981312871]
  option estimated value: 23.326296635630996 with confidence interval [ 9.96754032 36.68505295]
  


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