Autocallables with Integration Amplitude Loading

View on GitHub

Open this notebook in GitHub to run it yourself

This notebook covers the implementation of the Integration Amplitude Loading Method for the autocallables based on [1] and [2] using the Classiq platform’s Qmod language.

Data Definitions

import numpy as np
import scipy

S = 18  # notional value (can be asset initial price)
dt = 1  # in year
NUM_MARKET_DAYS_IN_YEAR = 250

DAILY_SIGMA = 0.0150694  # daily std log return
DAILY_MU = 0.00050963  # daily mean log return

SIGMA = DAILY_SIGMA * np.sqrt(NUM_MARKET_DAYS_IN_YEAR)  # annual
MU = DAILY_MU * NUM_MARKET_DAYS_IN_YEAR  # annual

b = 0.7  # binary barrier (to check with returns)
K_comp = 1  # K put (to check with returns)
K = K_comp * S  # K put (to to use for the payoff)
r = 0.04  # annual risk free rate

K_bin_1 = 1.1  # K binaries (to check with returns)
K_bin_2 = 1.1  # K binaries (to check with returns)
K_bin_norm = [
    np.log(K_bin_1),
    np.log(K_bin_2),
]  # K binaries (to check with log returns)
bin_1_payoff = 2 * np.exp(-r * 1)  # payoff binaries already discounted
bin_2_payoff = 3 * np.exp(-r * 2)  # payoff binaries already discounted


NUM_QUBITS = 2
TIME_STEPS = 3  # 3 years


PRECISION = 2

Gaussian State Preparation

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, stds_around_mean_to_include=3
)

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

Rescalings

Compute

R_T^{max}

resulting from discretization:

R_T_MAX_PROP = 0

if MU > 0 and SIGMA > 0:
    R_T_MAX_PROP = TIME_STEPS * (MU * dt + SIGMA * np.sqrt(dt) * (grid_points[-1]))

R_T_MIN_PROP = 0

if MU > 0 and SIGMA > 0:
    R_T_MIN_PROP = TIME_STEPS * (MU * dt + SIGMA * np.sqrt(dt) * (grid_points[0]))

R_T_MAX_PROP = np.max([np.abs(R_T_MIN_PROP), np.abs(R_T_MAX_PROP)])

R_T_MAX = np.log(
    np.max(
        [
            (np.exp(TIME_STEPS * r) * bin_2_payoff) + K,
            (np.exp(TIME_STEPS * r) * bin_1_payoff) + K,
            K,
        ]
    )
    / S
)
R_T_MAX = max(R_T_MAX_PROP, R_T_MAX)

In two’s complement, given

N

as the number of qubits, represent from

-2^{N-1}

and

2^{N-1}-1

a = 1 / (2**PRECISION)

if R_T_MAX < 1:
    int_places = 1
else:
    int_places = np.ceil(np.log2(np.ceil(R_T_MAX))) + 1

num_norm_factor = np.exp(a * (2 ** (int_places + PRECISION))) - 1
den_norm_factor = np.exp(a * (2 ** (int_places + PRECISION - 1) + 1))


norm_factor = S * (num_norm_factor / den_norm_factor)

Compute Constant Rotations

def compute_constant_rotation(const_payoff):
    b1 = const_payoff * np.exp(r * TIME_STEPS)
    num1 = (b1 + K) * den_norm_factor
    den1 = S * num_norm_factor
    den2 = num_norm_factor
    num2 = 1
    angle = (num1 / den1) - (num2 / den2)
    return 2 * np.arcsin(np.sqrt(angle))

bin_1_payoff_rotation = compute_constant_rotation(bin_1_payoff)
bin_2_payoff_rotation = compute_constant_rotation(bin_2_payoff)
zero_rotation = compute_constant_rotation(0)

bit_T_normalize = (S * np.exp(-r * TIME_STEPS)) / den_norm_factor


def postprocessing(x):
    return (
        (x * norm_factor * np.exp(-r * TIME_STEPS))
        + bit_T_normalize
        - (K * np.exp(-r * TIME_STEPS))
    )

Verifications

bin_1_payoff

Output:

np.float64(1.9215788783046464)

postprocessing((np.sin((compute_constant_rotation(bin_1_payoff) / 2))) ** 2)

Output:

np.float64(1.9215788783046523)

bin_2_payoff

Output:

np.float64(2.7693490391599074)

postprocessing((np.sin((compute_constant_rotation(bin_2_payoff) / 2))) ** 2)

Output:

np.float64(2.7693490391599074)

postprocessing((np.sin((compute_constant_rotation(0) / 2))) ** 2)

Output:

np.float64(1.7763568394002505e-15)

Compute constants for comparisons

b_norm = np.log(b)
K_norm = np.log(K_comp)

Classical Payoff

from itertools import product

import pandas as pd

simulation = pd.DataFrame.from_dict(
    {
        "quantum_samples": list(range(len(grid_points))),
        "classical_samples": grid_points,
        "probabilities": probabilities,
    }
)

binary_combinations = list(product(range(2**NUM_QUBITS), repeat=TIME_STEPS))
sim = pd.DataFrame(binary_combinations)
new_col = ["time_" + str(i) for i in range(TIME_STEPS)]
sim.columns = new_col

def round_down(a):
    precision_factor = 2 ** (PRECISION)
    return np.floor(a * precision_factor) / precision_factor

def round_factor(a):
    # precision_factor = 2 ** (PRECISION)
    # return np.round(a * precision_factor) / precision_factor
    return floor_factor(a)


def floor_factor(a):
    precision_factor = 2 ** (PRECISION)

    if a >= 0:
        return np.floor(a * precision_factor) / precision_factor
    else:
        return np.ceil(a * precision_factor) / precision_factor

for i in range(TIME_STEPS):
    sim = sim.merge(simulation, left_on="time_" + str(i), right_on="quantum_samples")
    sim = sim.drop("quantum_samples", axis=1)
    new_col = new_col + ["c_" + str(i), "q_prob_" + str(i)]
    sim.columns = new_col
    sim["log_ret_val_" + str(i)] = (
        round_factor(MU * dt + np.sqrt(dt) * SIGMA * MIN_X)
        + round_factor(SIGMA * np.sqrt(dt) * STEP_X) * sim["time_" + str(i)]
    )
    new_col = new_col + ["log_ret_val_" + str(i)]
    sim.columns = new_col
    if i != 0:
        sim["log_ret_val_" + str(i)] = (
            sim["log_ret_val_" + str(i)] + sim["log_ret_val_" + str(i - 1)]
        )
    sim["ret_val_" + str(i)] = np.exp(sim["log_ret_val_" + str(i)])
    new_col = new_col + ["ret_val_" + str(i)]

sim["prob"] = 1
for i in range(TIME_STEPS):
    sim["prob"] = sim["prob"] * sim["q_prob_" + str(i)]

for i in range(TIME_STEPS):
    sim["b_crossed_" + str(i)] = sim["log_ret_val_" + str(i)] < round_factor(b_norm)

sim["bin_1_activate"] = sim["log_ret_val_0"] > round_factor(K_bin_norm[0])
sim["bin_2_activate"] = sim["log_ret_val_1"] > round_factor(K_bin_norm[1])

sim["K_put"] = sim["log_ret_val_" + str(TIME_STEPS - 1)] < round_factor(K_norm)

sim["payoff"] = 0.0

sim.loc[sim["bin_1_activate"], "payoff"] = bin_1_payoff
sim.loc[(~sim["bin_1_activate"]) & (sim["bin_2_activate"]), "payoff"] = bin_2_payoff

barrier_crossed_once = sim["b_crossed_0"]
for i in range(1, TIME_STEPS):
    barrier_crossed_once = barrier_crossed_once | sim["b_crossed_" + str(i)]

put_condition = (
    (~((sim["bin_1_activate"]) | (sim["bin_2_activate"])))
    & (barrier_crossed_once)
    & (sim["K_put"])
)

sim.loc[put_condition, "payoff"] = (
    S * (sim[put_condition]["ret_val_2"] 

- K_comp) * np.exp(-r * TIME_STEPS)
)


expected_payoff = sum(sim["prob"] * sim["payoff"])


print("expected payoff classical: " + str(expected_payoff))

Output:

expected payoff classical: -3.5032073946209352

Integration Method Circuit Synthesis

from classiq import *
from classiq.qmod.symbolic import sqrt


@qfunc
def integrator(exp_rate: CReal, x: Const[QNum], ref: QNum, res: QBit) -> None:
    prepare_exponential_state(-exp_rate, ref)
    res ^= x >= ref

def affine_py(x: QNum):
    return MU * dt + SIGMA * sqrt(dt) * (x * STEP_X + MIN_X)

@qperm
def add_minimum(x: QArray):
    X(x[x.size - 1])

round_K_bin_norm = []
round_b_norm = round_factor(b_norm)
round_K_bin_norm.append(round_factor(K_bin_norm[0]))
round_K_bin_norm.append(round_factor(K_bin_norm[1]))
round_k_norm = round_factor(K_norm)

@qfunc
def integration_load_amplitudes(y: Const[QNum], aux_reg: QNum, ind_reg: QBit):
    exp_rate = 1 / (2**PRECISION)
    integrator(exp_rate, y, aux_reg, ind_reg)


@qfunc
def integration_payoff(log_return: Const[QNum], aux_reg: QNum, ind_reg: QBit):
    log_return_unsigned = QNum(size=log_return.size)

    within_apply(
        lambda: [
            bind(log_return, log_return_unsigned),
            add_minimum(log_return_unsigned),
        ],
        lambda: integration_load_amplitudes(log_return_unsigned, aux_reg, ind_reg),
    )


@qperm
def check_barrier_crossed(log_return: Const[QNum], barrier_crossed: QBit):
    barrier_crossed ^= log_return < round_b_norm


@qperm
def check_binary(log_return: Const[QNum], round_K_bin_norm: CReal, binary_valid: QBit):
    binary_valid ^= log_return > round_K_bin_norm


@qperm
def check_K_put(
    log_return: Const[QNum[PRECISION + int_places, SIGNED, PRECISION]],
    k_put_valid: QBit,
):
    k_put_valid ^= log_return < round_k_norm


@qfunc
def binary_payoff(binary_payoff: CReal, target: QBit):
    RY(binary_payoff, target)


@qfunc
def zero_payoff(target: QBit):
    RY(zero_rotation, target)


@qperm
def check_put_activate(
    barriers_crossed: Const[QArray],
    k_put_valid: Const[QBit],
    put_alone_valid: QBit,
):
    put_alone_valid ^= (
        (
            (barriers_crossed[0] == 1)
            | (barriers_crossed[1] == 1)
            | (barriers_crossed[2] == 1)
        )
        & (k_put_valid == 1)
    ) == 1


@qperm
def populate_put_alone_valid(
    sum_log_return: Const[QNum],
    barriers_crossed: QArray,
    put_alone_valid: Output[QBit],
):
    k_put_valid = QBit()

    allocate(put_alone_valid)

    within_apply(
        lambda: [
            allocate(k_put_valid),
            check_K_put(sum_log_return, k_put_valid),
            check_barrier_crossed(
                sum_log_return, barriers_crossed[2]
            ),  # magari qui si può condizionare
        ],
        lambda: check_put_activate(barriers_crossed, k_put_valid, put_alone_valid),
    )


@qfunc
def final_payoff(
    check_all_zeros: Const[QNum],
    sum_log_return: Const[QNum],
    aux_reg: QNum,
    ind_reg: QBit,
):
    control(check_all_zeros == 0, lambda: zero_payoff(ind_reg))
    control(
        check_all_zeros == 1,
        lambda: integration_payoff(sum_log_return, aux_reg, ind_reg),
    )


@qfunc
def autocallable_integration(
    x: QArray[QNum, TIME_STEPS],
    aux_reg: QNum,
    ind_reg: QBit,
    sum_log_return: QNum,
    first_bin_valid: QBit,
    second_bin_valid: QBit,
    barriers_crossed: QArray[QBit],
) -> None:
    repeat(x.len, lambda i: inplace_prepare_state(probabilities, 0, x[i]))

    sum_log_return ^= affine_py(x[0])
    check_binary(sum_log_return, round_K_bin_norm[0], first_bin_valid)
    control(first_bin_valid, lambda: binary_payoff(bin_1_payoff_rotation, ind_reg))
    # check_barrier_crossed(sum_log_return, barriers_crossed[0])
    check_barrier_crossed(sum_log_return, barriers_crossed[0])

    sum_log_return += affine_py(x[1])
    check_binary(sum_log_return, round_K_bin_norm[1], second_bin_valid)
    control(
        ((first_bin_valid == 0) & (second_bin_valid == 1)) == 1,
        lambda: binary_payoff(bin_2_payoff_rotation, ind_reg),
    )  # check order
    # check_barrier_crossed(sum_log_return, barriers_crossed[1]) #magari qui si può condizionare
    check_barrier_crossed(sum_log_return, barriers_crossed[1])

    sum_log_return += affine_py(x[2])

    put_alone_valid = QBit()
    within_apply(
        lambda: populate_put_alone_valid(
            sum_log_return,
            barriers_crossed,
            put_alone_valid,
        ),
        lambda: final_payoff(
            [put_alone_valid, first_bin_valid, second_bin_valid],
            sum_log_return,
            aux_reg,
            ind_reg,
        ),
    )

IQAE Functions and QStruct

from classiq.applications.iqae.iqae import IQAE


class OracleVars(QStruct):
    x: QArray[QNum[NUM_QUBITS, False, 0], TIME_STEPS]
    aux_reg: QNum[int_places + PRECISION]
    sum_log_return: QNum[int_places + PRECISION, SIGNED, PRECISION]
    first_bin_valid: QBit
    second_bin_valid: QBit
    barriers_crossed: QArray[QBit, TIME_STEPS]


@qfunc
def iqae_state_preparation(state: OracleVars, ind: QBit):
    autocallable_integration(
        state.x,
        state.aux_reg,
        ind,
        state.sum_log_return,
        state.first_bin_valid,
        state.second_bin_valid,
        state.barriers_crossed,
    )

Base Simulator Synthesis

iqae = IQAE(
    state_prep_op=iqae_state_preparation,
    problem_vars_size=NUM_QUBITS * TIME_STEPS
    + 2 * (int_places + PRECISION)
    + 2
    + TIME_STEPS,
    constraints=Constraints(optimization_parameter="width"),
    preferences=Preferences(
        optimization_level=1, machine_precision=PRECISION, timeout_seconds=2000
    ),
)

qmod = iqae.get_model()

print("Starting synthesis")
qprog = iqae.get_qprog()
show(qprog)

Output:

Starting synthesis
  Quantum program link: https://platform.classiq.io/circuit/3564N6oLHen6uvAwSGk0stVkwht
  

print("Circuit width: ", qprog.data.width)

Output:


Circuit width:  24

Execution takes a lot of time. Examine the results:

# takes a lot of time

# EPSILON = 0.001
# ALPHA = 0.002
# res_iqae= iqae.run(EPSILON, ALPHA, execution_preferences=ExecutionPreferences(shots=100000))

# print("the expected result from classical computation is: " + str(expected_payoff))
# print("the result from IQAE is: " + str(postprocessing(res_iqae.estimation)))
# print("the confidence interval of the quantum estimation is: " + str(postprocessing(res_iqae.confidence_interval[0]))+"," +str(postprocessing(res_iqae.confidence_interval[1])) )
# print(f"Synthesis and execution time: {time() - start_time} seconds")

print("the expected result from classical computation is: " + str(expected_payoff))
print("the result from IQAE is: " + str(-3.5079217726515903))
print(
    "the confidence interval of the quantum estimation is: "
    + str(postprocessing(0.119158741))
    + ","
    + str(postprocessing(0.1197666))
)

Output:

the expected result from classical computation is: -3.5032073946209352
  the result from IQAE is: -3.5079217726515903
  the confidence interval of the quantum estimation is: -3.535334467336666,-3.4805134318653383
  

References

[1] Francesca Cibrario et al. (2024). Quantum Amplitude Loading for Rainbow Options Pricing. Preprint. [2] Shouvanik Chakrabarti et al. (2021). A Threshold for Quantum Advantage in Derivative Pricing, Quantum 5, 463.

Algorithms

Applications

Tutorials

Autocallables with Integration Amplitude Loading

View on GitHub

Data Definitions

Gaussian State Preparation

Rescalings

Compute Constant Rotations

Verifications

Classical Payoff

Integration Method Circuit Synthesis

IQAE Functions and QStruct

Base Simulator Synthesis

References

View on GitHub

​Data Definitions

​Gaussian State Preparation

​Rescalings

​Compute Constant Rotations

​Verifications

​Classical Payoff

​Integration Method Circuit Synthesis

​IQAE Functions and QStruct

​Base Simulator Synthesis

​References

Data Definitions

Gaussian State Preparation

Rescalings

Compute Constant Rotations

Verifications

Classical Payoff

Integration Method Circuit Synthesis

IQAE Functions and QStruct

Base Simulator Synthesis

References