Autocallables with Integration Amplitude Loading
In this Notebook we will go through the implementation of the Integration Amplitude Loading Method for the autocallables based on https://arxiv.org/pdf/2402.05574.pdf and https://arxiv.org/pdf/2012.03819 using classiq platform 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 number of qubits, we can 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
1.9215788783046464
postprocessing((np.sin((compute_constant_rotation(bin_1_payoff) / 2))) ** 2)
1.9215788783046523
bin_2_payoff
2.7693490391599074
postprocessing((np.sin((compute_constant_rotation(bin_2_payoff) / 2))) ** 2)
2.7693490391599074
postprocessing((np.sin((compute_constant_rotation(0) / 2))) ** 2)
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))
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)
@qfunc
def add_minimum(x: Permutable[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),
)
@qfunc
def check_barrier_crossed(log_return: Const[QNum], barrier_crossed: Permutable[QBit]):
barrier_crossed ^= log_return < round_b_norm
@qfunc
def check_binary(
log_return: Const[QNum], round_K_bin_norm: CReal, binary_valid: Permutable[QBit]
):
binary_valid ^= log_return > round_K_bin_norm
@qfunc
def check_K_put(
log_return: Const[QNum[PRECISION + int_places, SIGNED, PRECISION]],
k_put_valid: Permutable[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)
@qfunc
def check_put_activate(
barriers_crossed: Const[QArray],
k_put_valid: Const[QBit],
put_alone_valid: Permutable[QBit],
):
put_alone_valid ^= (
(
(barriers_crossed[0] == 1)
| (barriers_crossed[1] == 1)
| (barriers_crossed[2] == 1)
)
& (k_put_valid == 1)
) == 1
@qfunc
def populate_put_alone_valid(
sum_log_return: Const[QNum],
barriers_crossed: Permutable[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: 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 sythesis
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)
Starting synthesis
Quantum program link: https://platform.classiq.io/circuit/2yrxvhXoIZN8IsaCrh1mIqYP0Py
print("Circuit width: ", qprog.data.width)
Circuit width: 24
Execution takes a lot of time
See the results below
# 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))
)
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