diff --git a/research/rainbow_options/RainbowOptions-DirectMethod.ipynb b/research/rainbow_options/RainbowOptions-DirectMethod.ipynb
new file mode 100644
index 00000000..baaa0eac
--- /dev/null
+++ b/research/rainbow_options/RainbowOptions-DirectMethod.ipynb
@@ -0,0 +1,629 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "1197af31",
+   "metadata": {},
+   "source": [
+    "# Rainbow options with Direct Amplitude Loading\n",
+    "In this Notebook we will go through the implementation of the Direct Amplitude Loading Method for the rainbow option presented in [[1]](#QALROP)"
+   ]
+  },
+  {
+   "metadata": {},
+   "cell_type": "markdown",
+   "source": [
+    "### Introduction\n",
+    "\n",
+    "In finance, a crucial aspect of asset pricing pertains to derivatives.\n",
+    "Derivatives are contracts whose value is contingent upon another source, known as the underlying.\n",
+    "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).\n",
+    "This process relies on mathematical models, considering variables like current asset prices, time to expiration, volatility, and interest rates."
+   ],
+   "id": "f63599e2e3afbd49"
+  },
+  {
+   "cell_type": "markdown",
+   "id": "95faeb3c",
+   "metadata": {},
+   "source": [
+    "## Data Definitions\n",
+    "\n",
+    "The problem inputs are:\n",
+    "- `NUM_QUBITS`: the number of qubits representing an underlying asset\n",
+    "- `NUM_ASSETS`: the number of underlying assets\n",
+    "- `K`: the strike price\n",
+    "- `S0`: the arrays of underlying assets prices\n",
+    "- `dt`: the number of days to the maturity date\n",
+    "- `COV`: the covariance matrix that correlate the underlyings\n",
+    "- `MU_LOG_RET`: the array containing the mean of the log return of each underlyings"
+   ]
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "import numpy as np\n",
+    "import scipy\n",
+    "\n",
+    "NUM_QUBITS = 2\n",
+    "NUM_ASSETS = 2\n",
+    "\n",
+    "K = 190\n",
+    "S0 = [193.97, 189.12]\n",
+    "dt = 250\n",
+    "\n",
+    "COV = np.array([[0.000335, 0.000257], [0.000257, 0.000418]])\n",
+    "MU_LOG_RET = np.array([0.00050963, 0.00062552])"
+   ],
+   "id": "31301483d80050bc",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "code",
+   "id": "d2d99529",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "MU = MU_LOG_RET * dt\n",
+    "CHOLESKY = np.linalg.cholesky(COV) * np.sqrt(dt)\n",
+    "SCALING_FACTOR = 1 / CHOLESKY[0, 0]"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4be6b880",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-04-01T09:40:28.838996Z",
+     "start_time": "2024-04-01T09:40:28.226590Z"
+    }
+   },
+   "source": [
+    "## Gaussian State preparation\n",
+    "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$. \n",
+    "\n",
+    "$$\n",
+    "\\sum_{i=0}^{N-1} \\sqrt{p(w_i)}\\left|w_i\\right\\rangle \n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "172e7a00",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "def gaussian_discretization(num_qubits, mu=0, sigma=1, stds_around_mean_to_include=3):\n",
+    "    lower = mu - stds_around_mean_to_include * sigma\n",
+    "    upper = mu + stds_around_mean_to_include * sigma\n",
+    "    num_of_bins = 2**num_qubits\n",
+    "    sample_points = np.linspace(lower, upper, num_of_bins + 1)\n",
+    "\n",
+    "    def single_gaussian(x: np.ndarray, _mu: float, _sigma: float) -> np.ndarray:\n",
+    "        cdf = scipy.stats.norm.cdf(x, loc=_mu, scale=_sigma)\n",
+    "        return cdf[1:] - cdf[0:-1]\n",
+    "\n",
+    "    non_normalized_pmf = (single_gaussian(sample_points, mu, sigma),)\n",
+    "    real_probs = non_normalized_pmf / np.sum(non_normalized_pmf)\n",
+    "    return sample_points[:-1], real_probs[0].tolist()\n",
+    "\n",
+    "\n",
+    "grid_points, probabilities = gaussian_discretization(NUM_QUBITS)\n",
+    "\n",
+    "STEP_X = grid_points[1] - grid_points[0]\n",
+    "MIN_X = grid_points[0]"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "markdown",
+   "id": "151bb5cb",
+   "metadata": {},
+   "source": [
+    "### SANITY CHECK"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "bbb4f025",
+   "metadata": {},
+   "source": [
+    "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."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "5713c295",
+   "metadata": {},
+   "source": [
+    "from IPython.display import Markdown\n",
+    "\n",
+    "if K >= max(S0*np.exp(np.dot(CHOLESKY,[grid_points[-1]]*2) + MU)):\n",
+    "    display(Markdown('<font color=\\'red\\'> K always greater than the maximum asset values. Stop the run, the payoff is 0</font>'))"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ecd5856d",
+   "metadata": {},
+   "source": [
+    "## Maximum Computation"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6d35363b",
+   "metadata": {},
+   "source": [
+    "### Precision utils"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "7794ce7c",
+   "metadata": {},
+   "source": [
+    "FRAC_PLACES = 2\n",
+    "\n",
+    "def round_factor(a):\n",
+    "    precision_factor = 2 ** FRAC_PLACES\n",
+    "    return round(a * precision_factor) / precision_factor\n",
+    "\n",
+    "\n",
+    "def floor_factor(a):\n",
+    "    precision_factor = 2 ** FRAC_PLACES\n",
+    "    return np.floor(a * precision_factor) / precision_factor"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ed89df19",
+   "metadata": {},
+   "source": [
+    "### Affine and maximum arithmetic definitions\n",
+    "\n",
+    "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:"
+   ]
+  },
+  {
+   "metadata": {},
+   "cell_type": "markdown",
+   "source": "  $\\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.",
+   "id": "74a02871489d9ec9"
+  },
+  {
+   "metadata": {},
+   "cell_type": "markdown",
+   "source": [
+    "$y_i$ is obtained through the dot product between the matrix $\\mathbf{L}$ and the standard multivariate Gaussian sample:\n",
+    "\n",
+    "$$    y_i = \\Delta x \\cdot \\sum_kl_{ik}d_k + x_{min} \\cdot \\sum_k l_{ik} $$\n",
+    "$\\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.\n",
+    "$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."
+   ],
+   "id": "6d7e2e6ce33a07ce"
+  },
+  {
+   "cell_type": "code",
+   "id": "4f6f55cb",
+   "metadata": {},
+   "source": [
+    "from functools import reduce\n",
+    "\n",
+    "from classiq import QNum, get_expression_numeric_attributes, qfunc, Output\n",
+    "from classiq.qmod.symbolic import max as qmax\n",
+    "\n",
+    "a = STEP_X / SCALING_FACTOR\n",
+    "b = np.log(S0[0]) + MU[0] + MIN_X * CHOLESKY[0].sum()\n",
+    "c = (\n",
+    "        SCALING_FACTOR\n",
+    "        * (\n",
+    "        np.log(S0[1])\n",
+    "        + MU[1]\n",
+    "        - (np.log(S0[0]) + MU[0])\n",
+    "        + MIN_X * sum(CHOLESKY[1] - CHOLESKY[0])\n",
+    "    )\n",
+    "        / STEP_X\n",
+    ")\n",
+    "c = round_factor(c)\n",
+    "\n",
+    "\n",
+    "def get_affine_formula(assets, i):\n",
+    "    return reduce(\n",
+    "        lambda x, y: x + y,\n",
+    "        [\n",
+    "            assets[j] * round_factor(SCALING_FACTOR * CHOLESKY[i, j])\n",
+    "            for j in range(NUM_ASSETS)\n",
+    "            if CHOLESKY[i, j]\n",
+    "        ],\n",
+    "    )\n",
+    "\n",
+    "\n",
+    "def calculate_max_reg_type():\n",
+    "    x1 = QNum(\"x1\", NUM_QUBITS, False, 0)\n",
+    "    x2 = QNum(\"x2\", NUM_QUBITS, False, 0)\n",
+    "    expr = qmax(get_affine_formula([x1, x2], 0), get_affine_formula([x1, x2], 1) + c)\n",
+    "    size_in_bits, sign, fraction_digits = get_expression_numeric_attributes(\n",
+    "        [x1, x2], expr\n",
+    "    )\n",
+    "    return size_in_bits, fraction_digits\n",
+    "\n",
+    "\n",
+    "MAX_NUM_QUBITS = calculate_max_reg_type()[0]\n",
+    "MAX_FRAC_PLACES = calculate_max_reg_type()[1]"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "code",
+   "id": "0074a98c",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "@qfunc\n",
+    "def affine_max(x1: QNum, x2: QNum, res: Output[QNum]):\n",
+    "    res |= qmax(get_affine_formula([x1, x2], 0), get_affine_formula([x1, x2], 1) + c)"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b37716e7",
+   "metadata": {
+    "tags": []
+   },
+   "source": [
+    "## Direct Method\n",
+    "The direct exponential amplitude loading encodes in $\\tilde{f}$ the following function:\n",
+    "\\begin{equation}\n",
+    "\\tilde{f}(x)=\n",
+    "    \\begin{cases}\n",
+    "        e^{-a\\hat{x}}, & \\text{if } \\frac{x}{2^P} \\geq \\frac{\\log(K) -b'}{b}\\\\\n",
+    "      Ke^{-(b'+ ax_{max})}, & \\text{if } \\frac{x}{2^P} < \\frac{\\log(K) -b'}{b}\n",
+    "    \\end{cases}\n",
+    "\\end{equation}\n",
+    "\n",
+    "where $\\hat{x}$ is the binary complement of $x$ ($\\hat{x}=x-x_{max}$) and $x_{max}=2^R-1$, the maximum value that can be stored in $|x\\rangle$ register. For loading $e^{-a\\hat{x}}$ the $|r\\rangle$ is initialized to all zeros and one controlled rotation for each qubit is performed as shown in Figure \\ref{fig:circuit}. The rotations angles are: $\\theta_i = 2\\arccos \\left({\\sqrt{e^{-a2^i}}}\\right)$. All the probabilities of getting a $|0\\rangle^{\\otimes{R}}$ in the $|r\\rangle$ are then collected by a multi-controlled X (MCX) gate and stored in the $|1\\rangle$ state of a target qubit."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "06717f7c",
+   "metadata": {},
+   "source": [
+    "from classiq import X, RY, apply_to_all, CReal, QArray, repeat, QBit, Constraints, synthesize, create_model, allocate, inplace_prepare_state, within_apply, control, bind, show\n",
+    "from classiq.qmod.symbolic import acos, asin, sqrt, exp\n",
+    "\n",
+    "@qfunc\n",
+    "def exponential_amplitude_loading(\n",
+    "    exp_rate: CReal, x: QArray[QBit], aux: QArray[QBit], res: QBit\n",
+    ") -> None:\n",
+    "    apply_to_all(X, x)\n",
+    "    repeat(\n",
+    "        x.len,\n",
+    "        lambda index: control(\n",
+    "            x[index],\n",
+    "            lambda: RY(2 * acos(1 / sqrt(exp(exp_rate * (2**index)))), aux[index])\n",
+    "        ),\n",
+    "    )\n",
+    "    apply_to_all(X, x)\n",
+    "\n",
+    "    aux_num = QNum(\"aux_num\", aux.len, False, 0)\n",
+    "    bind(aux, aux_num)\n",
+    "    control(aux_num == 0, lambda: X(res))\n",
+    "    bind(aux_num, aux)"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "code",
+   "id": "ae5249a2",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "def get_payoff_expression(x, size, fraction_digits):\n",
+    "    payoff = sqrt(\n",
+    "        qmax(\n",
+    "            S0[0]\n",
+    "            * exp(\n",
+    "                STEP_X / SCALING_FACTOR * (2 ** (size - fraction_digits)) * x\n",
+    "                + (MU[0] + MIN_X * CHOLESKY[0].sum())\n",
+    "            ),\n",
+    "            K,\n",
+    "        )\n",
+    "    )\n",
+    "    return payoff\n",
+    "\n",
+    "\n",
+    "def get_strike_price_theta_direct(x: QNum):\n",
+    "    x_max = 1 - 1 / (2**x.size)\n",
+    "    payoff_max = get_payoff_expression(x_max, x.size, x.fraction_digits)\n",
+    "\n",
+    "    return 2 * asin(np.sqrt(K) / payoff_max)\n",
+    "\n",
+    "\n",
+    "@qfunc\n",
+    "def direct_load_amplitudes(geq_reg: QNum, max_reg: QNum, aux_reg: QNum, ind_reg: QBit):\n",
+    "    exp_rate = (1 / (2**max_reg.fraction_digits)) * a\n",
+    "    control(\n",
+    "        geq_reg == 1,\n",
+    "        lambda: exponential_amplitude_loading(exp_rate, max_reg, aux_reg, ind_reg),\n",
+    "    )\n",
+    "    strike_price_theta = get_strike_price_theta_direct(max_reg)\n",
+    "    control(geq_reg == 0, lambda: RY(strike_price_theta, ind_reg))\n",
+    "\n",
+    "\n",
+    "@qfunc\n",
+    "def direct_payoff(max_reg: QNum, aux_reg: QNum, ind_reg: QBit):\n",
+    "    geq_reg = QBit(\"geq_reg\")\n",
+    "    within_apply(\n",
+    "        lambda: asset_geq_strike_price(max_reg, geq_reg),\n",
+    "        lambda: direct_load_amplitudes(geq_reg, max_reg, aux_reg, ind_reg),\n",
+    "    )\n",
+    "\n",
+    "@qfunc\n",
+    "def asset_geq_strike_price(\n",
+    "    x: QNum,\n",
+    "    res: Output[QBit],\n",
+    ") -> None:\n",
+    "    a = STEP_X / SCALING_FACTOR\n",
+    "    b = np.log(S0[0]) + MU[0] + MIN_X * CHOLESKY[0].sum()\n",
+    "    COMP_VALUE = (np.log(K) - b) / a\n",
+    "    res |= x > floor_factor(COMP_VALUE)\n",
+    "\n",
+    "@qfunc\n",
+    "def rainbow_direct(\n",
+    "    x1: QNum,  \n",
+    "    x2: QNum,\n",
+    "    aux_reg: QNum,\n",
+    "    ind_reg: QBit,\n",
+    ") -> None:\n",
+    "    inplace_prepare_state(probabilities, 0, x1)\n",
+    "    inplace_prepare_state(probabilities, 0, x2)\n",
+    "    max_out = QNum(\"max_out\")\n",
+    "    within_apply(\n",
+    "        lambda: affine_max(x1, x2, max_out),\n",
+    "        lambda: direct_payoff(max_out, aux_reg, ind_reg),\n",
+    "    )\n",
+    "\n",
+    "\n",
+    "@qfunc\n",
+    "def main(\n",
+    "    x1: Output[QNum],  \n",
+    "    x2: Output[QNum],\n",
+    "    aux_reg: Output[QNum],\n",
+    "    ind_reg: Output[QBit],\n",
+    ") -> None:\n",
+    "    allocate(MAX_NUM_QUBITS, aux_reg)\n",
+    "    allocate(NUM_QUBITS, x1)\n",
+    "    allocate(NUM_QUBITS, x2)\n",
+    "    allocate(1, ind_reg)\n",
+    "    rainbow_direct(x1, x2, aux_reg, ind_reg)\n",
+    "\n",
+    "constraints = Constraints(max_width=23)\n",
+    "qmod = create_model(main, constraints=constraints)\n",
+    "print(\"Starting synthesis\")\n",
+    "qprog = synthesize(qmod)\n",
+    "show(qprog)"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "markdown",
+   "id": "661a1859",
+   "metadata": {},
+   "source": [
+    "## IQAE algorithm"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "ecb49c61",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "from classiq import cfunc, Z\n",
+    "from classiq.qmod.builtins.classical_execution_primitives import iqae, save\n",
+    "\n",
+    "\n",
+    "@qfunc\n",
+    "def qmci_oracle(ind: QBit):\n",
+    "    Z(ind)\n",
+    "\n",
+    "@cfunc\n",
+    "def cmain():\n",
+    "    iqae_res = iqae(epsilon=0.05, alpha=0.1)\n",
+    "    save({\"iqae_res\": iqae_res})"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "code",
+   "id": "3dddfbea",
+   "metadata": {},
+   "source": [
+    "from classiq import CInt, QCallable, grover_operator, power\n",
+    "\n",
+    "@qfunc\n",
+    "def grover_algorithm(\n",
+    "    k: CInt,\n",
+    "    oracle_operand: QCallable[QArray[QBit]],\n",
+    "    sp_operand: QCallable[QArray[QBit]],\n",
+    "    x: QArray[QBit],\n",
+    "):\n",
+    "    sp_operand(x)\n",
+    "    power(k, lambda: grover_operator(oracle_operand, sp_operand, x))"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "code",
+   "id": "19a8b150",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "def get_main():\n",
+    "    @qfunc\n",
+    "    def main(\n",
+    "        k: CInt,\n",
+    "        ind_reg: Output[QBit],\n",
+    "    ) -> None:\n",
+    "        full_reg = QArray(\"full_reg\")\n",
+    "        allocate(2 * NUM_QUBITS + MAX_NUM_QUBITS + 1, full_reg)\n",
+    "        grover_algorithm(\n",
+    "            k,\n",
+    "            lambda x: qmci_oracle(x[x.len - 1]),\n",
+    "            lambda x: rainbow_direct(\n",
+    "                x[0:NUM_QUBITS],\n",
+    "                x[NUM_QUBITS : 2 * NUM_QUBITS],\n",
+    "                x[2 * NUM_QUBITS : 2 * NUM_QUBITS + MAX_NUM_QUBITS],\n",
+    "                x[x.len - 1],\n",
+    "            ),\n",
+    "            full_reg,\n",
+    "        )\n",
+    "        state_reg = QArray(\"state_reg\")\n",
+    "        bind(full_reg, [state_reg, ind_reg])\n",
+    "\n",
+    "    return main"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "code",
+   "id": "355f9bd4",
+   "metadata": {},
+   "source": [
+    "from classiq import execute\n",
+    "\n",
+    "def synthesize_and_execute(post_process):\n",
+    "    constraints = Constraints(max_width=25)\n",
+    "    qmod = create_model(\n",
+    "        get_main(),\n",
+    "        constraints=constraints,\n",
+    "        classical_execution_function=cmain,\n",
+    "    )\n",
+    "    print(\"Starting synthesis\")\n",
+    "    qprog = synthesize(qmod)\n",
+    "    show(qprog)\n",
+    "    print(\"Starting execution\")\n",
+    "    res = execute(qprog).result()\n",
+    "    iqae_res = res[0].value\n",
+    "    parsed_res = post_process(res[0].value)\n",
+    "\n",
+    "    return (qmod, qprog, iqae_res, parsed_res)"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "markdown",
+   "id": "16899f05",
+   "metadata": {},
+   "source": [
+    "## Post Process\n",
+    "We need to add to the post-processing function a term:\n",
+    "\\begin{equation}\n",
+    "\\begin{split}\n",
+    "&\\mathbb{E} \\left[\\max\\left(e^{b \\cdot z}, Ke^{-b'}\\right) \\right] e^{b'} - K \\\\\n",
+    "= &\\mathbb{E} \\left[\\max\\left(e^{-a\\hat{x}}, Ke^{-b'-ax_{max}}\\right) \\right]e^{b'+ ax_{max}} - K \n",
+    "\\end{split}\n",
+    "\\end{equation}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "514ff089",
+   "metadata": {},
+   "source": [
+    "import sympy\n",
+    "\n",
+    "def parse_result_direct(iqae_res):\n",
+    "    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}))\"\n",
+    "    payoff_func = sympy.lambdify(sympy.symbols(\"x\"), payoff_expression)\n",
+    "    payoff_max = payoff_func(1 - 1 / (2**MAX_NUM_QUBITS)) \n",
+    "\n",
+    "    option_value = iqae_res.estimation * (payoff_max**2) - K\n",
+    "    confidence_interval = np.array(iqae_res.confidence_interval) * (payoff_max**2) - K\n",
+    "    return (option_value, confidence_interval)"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "markdown",
+   "id": "10ab8023",
+   "metadata": {},
+   "source": [
+    "# Run Method"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "af4639dd",
+   "metadata": {},
+   "source": [
+    "qmod_direct, qprog_direct, iqae_res_direct, parsed_res_direct = synthesize_and_execute(\n",
+    "    parse_result_direct\n",
+    ")\n",
+    "result, conf_interval = parsed_res_direct\n",
+    "print(f\"raw iqae results: {iqae_res_direct.estimation} with confidence interval {iqae_res_direct.confidence_interval}\")\n",
+    "print(f\"option estimated value: {result} with confidence interval {conf_interval}\")"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "markdown",
+   "source": [
+    "## References\n",
+    "\n",
+    "<a name='QALROP'>[1]</a>: [Francesca Cibrario et al., Quantum Amplitude Loading for Rainbow Options Pricing. Preprint](https://arxiv.org/abs/2402.05574v2)\n"
+   ],
+   "id": "d0e80f6569d777d2"
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.14"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/research/rainbow_options/RainbowOptions-IntegrationMethod.ipynb b/research/rainbow_options/RainbowOptions-IntegrationMethod.ipynb
new file mode 100644
index 00000000..af9cb2c4
--- /dev/null
+++ b/research/rainbow_options/RainbowOptions-IntegrationMethod.ipynb
@@ -0,0 +1,642 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "80dc8829",
+   "metadata": {},
+   "source": [
+    "# Rainbow options with Integration\n",
+    "In this Notebook we will go through the implementation of the Integration Method for the rainbow option presented in [[1]](#QALROP)"
+   ]
+  },
+  {
+   "metadata": {},
+   "cell_type": "markdown",
+   "source": [
+    "### Introduction\n",
+    "\n",
+    "In finance, a crucial aspect of asset pricing pertains to derivatives.\n",
+    "Derivatives are contracts whose value is contingent upon another source, known as the underlying.\n",
+    "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).\n",
+    "This process relies on mathematical models, considering variables like current asset prices, time to expiration, volatility, and interest rates."
+   ],
+   "id": "bc041860b4197872"
+  },
+  {
+   "cell_type": "markdown",
+   "id": "01d7b531",
+   "metadata": {},
+   "source": [
+    "## Data Definitions\n",
+    "\n",
+    "The problem inputs are:\n",
+    "- `NUM_QUBITS`: the number of qubits representing an underlying asset\n",
+    "- `NUM_ASSETS`: the number of underlying assets\n",
+    "- `K`: the strike price\n",
+    "- `S0`: the arrays of underlying assets prices\n",
+    "- `dt`: the number of days to the maturity date\n",
+    "- `COV`: the covariance matrix that correlate the underlyings\n",
+    "- `MU_LOG_RET`: the array containing the mean of the log return of each underlyings\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "3ad37909",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "import numpy as np\n",
+    "import scipy\n",
+    "\n",
+    "NUM_QUBITS = 2\n",
+    "NUM_ASSETS = 2\n",
+    "\n",
+    "K = 190\n",
+    "S0 = [193.97, 189.12]\n",
+    "dt = 250\n",
+    "\n",
+    "COV = np.array([[0.000335, 0.000257], [0.000257, 0.000418]])\n",
+    "MU_LOG_RET = np.array([0.00050963, 0.00062552])"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "MU = MU_LOG_RET * dt\n",
+    "CHOLESKY = np.linalg.cholesky(COV) * np.sqrt(dt)\n",
+    "SCALING_FACTOR = 1 / CHOLESKY[0, 0]"
+   ],
+   "id": "1d608bccabd7ae78",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "markdown",
+   "id": "31eec2f4",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2024-04-01T09:40:28.838996Z",
+     "start_time": "2024-04-01T09:40:28.226590Z"
+    }
+   },
+   "source": [
+    "## Gaussian State preparation\n",
+    "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=$`NUM_QUBITS`) according to the formula $N=2^n$. \n",
+    "\n",
+    "$$\n",
+    "\\sum_{i=0}^{N-1} \\sqrt{p(w_i)}\\left|w_i\\right\\rangle \n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "69184d44",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "def gaussian_discretization(num_qubits, mu=0, sigma=1, stds_around_mean_to_include=3):\n",
+    "    lower = mu - stds_around_mean_to_include * sigma\n",
+    "    upper = mu + stds_around_mean_to_include * sigma\n",
+    "    num_of_bins = 2**num_qubits\n",
+    "    sample_points = np.linspace(lower, upper, num_of_bins + 1)\n",
+    "\n",
+    "    def single_gaussian(x: np.ndarray, _mu: float, _sigma: float) -> np.ndarray:\n",
+    "        cdf = scipy.stats.norm.cdf(x, loc=_mu, scale=_sigma)\n",
+    "        return cdf[1:] - cdf[0:-1]\n",
+    "\n",
+    "    non_normalized_pmf = (single_gaussian(sample_points, mu, sigma),)\n",
+    "    real_probs = non_normalized_pmf / np.sum(non_normalized_pmf)\n",
+    "    return sample_points[:-1], real_probs[0].tolist()\n",
+    "\n",
+    "\n",
+    "grid_points, probabilities = gaussian_discretization(NUM_QUBITS)\n",
+    "\n",
+    "STEP_X = grid_points[1] - grid_points[0]\n",
+    "MIN_X = grid_points[0]"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "markdown",
+   "id": "67a3a5bf",
+   "metadata": {},
+   "source": [
+    "### SANITY CHECK"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "86b78cd0",
+   "metadata": {},
+   "source": [
+    "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."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "5b51052a",
+   "metadata": {},
+   "source": [
+    "from IPython.display import Markdown\n",
+    "\n",
+    "if K >= max(S0*np.exp(np.dot(CHOLESKY,[grid_points[-1]]*2) + MU)):\n",
+    "    display(Markdown('<font color=\\'red\\'> K always greater than the maximum asset values. Stop the run, the payoff is 0</font>'))"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7d236ed2",
+   "metadata": {},
+   "source": [
+    "## Maximum Computation"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "86d93299",
+   "metadata": {},
+   "source": [
+    "### Precision utils"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "4b0bc081",
+   "metadata": {},
+   "source": [
+    "FRAC_PLACES = 2\n",
+    "\n",
+    "def round_factor(a):\n",
+    "    precision_factor = 2 ** FRAC_PLACES\n",
+    "    return round(a * precision_factor) / precision_factor\n",
+    "\n",
+    "\n",
+    "def floor_factor(a):\n",
+    "    precision_factor = 2 ** FRAC_PLACES\n",
+    "    return np.floor(a * precision_factor) / precision_factor"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f36b4e2e",
+   "metadata": {},
+   "source": [
+    "### Affine and maximum arithmetic definitions"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "dc03563a",
+   "metadata": {},
+   "source": [
+    "from functools import reduce\n",
+    "\n",
+    "from classiq import QNum, get_expression_numeric_attributes, qfunc, Output\n",
+    "from classiq.qmod.symbolic import max as qmax\n",
+    "\n",
+    "a = STEP_X / SCALING_FACTOR\n",
+    "b = np.log(S0[0]) + MU[0] + MIN_X * CHOLESKY[0].sum()\n",
+    "c = (\n",
+    "        SCALING_FACTOR\n",
+    "        * (\n",
+    "        np.log(S0[1])\n",
+    "        + MU[1]\n",
+    "        - (np.log(S0[0]) + MU[0])\n",
+    "        + MIN_X * sum(CHOLESKY[1] - CHOLESKY[0])\n",
+    "    )\n",
+    "        / STEP_X\n",
+    ")\n",
+    "c = round_factor(c)\n",
+    "\n",
+    "\n",
+    "def get_affine_formula(assets, i):\n",
+    "    return reduce(\n",
+    "        lambda x, y: x + y,\n",
+    "        [\n",
+    "            assets[j] * round_factor(SCALING_FACTOR * CHOLESKY[i, j])\n",
+    "            for j in range(NUM_ASSETS)\n",
+    "            if CHOLESKY[i, j]\n",
+    "        ],\n",
+    "    )\n",
+    "\n",
+    "def calculate_max_reg_type():\n",
+    "    x1 = QNum(\"x1\", NUM_QUBITS, False, 0)\n",
+    "    x2 = QNum(\"x2\", NUM_QUBITS, False, 0)\n",
+    "    expr = qmax(get_affine_formula([x1, x2], 0), get_affine_formula([x1, x2], 1) + c)\n",
+    "    size_in_bits, sign, fraction_digits = get_expression_numeric_attributes(\n",
+    "        [x1, x2], expr\n",
+    "    )\n",
+    "    return size_in_bits, fraction_digits\n",
+    "\n",
+    "\n",
+    "MAX_NUM_QUBITS = calculate_max_reg_type()[0]\n",
+    "MAX_FRAC_PLACES = calculate_max_reg_type()[1]"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "code",
+   "id": "1f613230",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "@qfunc\n",
+    "def affine_max(x1: QNum, x2: QNum, res: Output[QNum]):\n",
+    "    res |= qmax(get_affine_formula([x1, x2], 0), get_affine_formula([x1, x2], 1) + c)"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b90dc16e",
+   "metadata": {
+    "tags": []
+   },
+   "source": [
+    "## Integration Method\n",
+    "\n",
+    "The comparator collects the probabilities $g(r)$ of $|r\\rangle$ state until $|r\\rangle$ register is lower than $|x\\rangle$:\n",
+    "\\begin{equation}\n",
+    "\\begin{split}\n",
+    "&\\sum_{r=0}^{2^R-1}{\\sqrt{g(r)}}|x\\rangle|r\\rangle|r\\leq x\\rangle \\\\\n",
+    "= &|x\\rangle \\otimes \\left[ \\sum_{r=0}^{x}{\\sqrt{g(r)}} |r\\rangle |1\\rangle + \\sum_{r=x}^{2^R-1}{\\sqrt{g(r)}} |r\\rangle |0\\rangle \\right]\n",
+    "\\end{split}\n",
+    "\\end{equation}\n",
+    "Collecting the probability to have $r\\leq x$ we can define the function:\n",
+    "\\begin{equation}\n",
+    "\\tilde{h}(x)=\\sum_{r=0}^{x}g(r)\n",
+    "\\end{equation}\n",
+    "Evaluating the probability to get a $|1\\rangle$ results in $\\sum_{x = 0}^{2^R-1}{\\tilde{h}(x)}$.\n",
+    "To obtain a given function $\\tilde{h}$ a proper function $g(r)$ should be chosen.\n",
+    "The $g(r)$ for $r=0$ value must therefore be\n",
+    "$\n",
+    "g(0) = \\tilde{h}(0)\n",
+    "$\n",
+    "and for all the other $r$:\n",
+    "$$\n",
+    "g(r) = \\tilde{h}(r)-\\tilde{h}(r-1)\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "75122cba",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "from classiq import QBit, prepare_exponential_state, bind\n",
+    "\n",
+    "@qfunc\n",
+    "def integrator(x: QNum, ref: QNum, res: QBit) -> None:\n",
+    "    exp_rate = (1 / (2**x.fraction_digits)) * a\n",
+    "    prepare_exponential_state(-exp_rate, ref)\n",
+    "    x_uint = QNum(\"x_uint\", x.size, False, 0)\n",
+    "    bind(x, x_uint)\n",
+    "    res ^= x_uint >= ref\n",
+    "    bind(x_uint, x)"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "code",
+   "id": "95c12d95",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "from classiq import control, RY\n",
+    "from classiq.qmod.symbolic import asin, exp, sqrt\n",
+    "\n",
+    "def get_strike_price_theta_integration(x: QNum):\n",
+    "    exp_rate = (1 / (2**x.fraction_digits)) * a\n",
+    "    B = (exp((2**x.size) * exp_rate) - 1) / exp(exp_rate)\n",
+    "    A = 1 / exp(exp_rate)\n",
+    "    C = S0[0] * exp((MU[0] + MIN_X * CHOLESKY[0].sum()))\n",
+    "    return 2 * asin(sqrt((K - (C * A)) / (C * B)))\n",
+    "\n",
+    "@qfunc\n",
+    "def integration_load_amplitudes(\n",
+    "    geq_reg: QNum, max_reg: QNum, integrator_reg: QNum, ind_reg: QBit\n",
+    "):\n",
+    "    control(geq_reg == 1, lambda: integrator(max_reg, integrator_reg, ind_reg))\n",
+    "    strike_price_theta = get_strike_price_theta_integration(max_reg)\n",
+    "    control(geq_reg == 0, lambda: RY(strike_price_theta, ind_reg))"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "code",
+   "id": "c4ab6177",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "@qfunc\n",
+    "def asset_geq_strike_price(\n",
+    "    x: QNum,\n",
+    "    res: Output[QBit],\n",
+    ") -> None:\n",
+    "    a = STEP_X / SCALING_FACTOR\n",
+    "    b = np.log(S0[0]) + MU[0] + MIN_X * CHOLESKY[0].sum()\n",
+    "    COMP_VALUE = (np.log(K) - b) / a\n",
+    "    res |= x > floor_factor(COMP_VALUE)\n"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "code",
+   "id": "07c9e4f7",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "from classiq import within_apply\n",
+    "\n",
+    "@qfunc\n",
+    "def integration_payoff(max_reg: QNum, integrator_reg: QNum, ind_reg: QBit):\n",
+    "    geq_reg = QBit(\"geq_reg\")\n",
+    "    within_apply(\n",
+    "        lambda: asset_geq_strike_price(max_reg, geq_reg),\n",
+    "        lambda: integration_load_amplitudes(geq_reg, max_reg, integrator_reg, ind_reg),\n",
+    "    )"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "code",
+   "id": "c7870263",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "from classiq import Constraints, create_model, synthesize, show, allocate, allocate_num, inplace_prepare_state\n",
+    "\n",
+    "@qfunc\n",
+    "def rainbow_integration(\n",
+    "    x1: QNum, \n",
+    "    x2: QNum,\n",
+    "    integrator_reg: QNum,\n",
+    "    ind_reg: QBit,\n",
+    ") -> None:\n",
+    "    inplace_prepare_state(probabilities, 0, x1)\n",
+    "    inplace_prepare_state(probabilities, 0, x2)\n",
+    "    max_out = QNum(\"max_out\")\n",
+    "    within_apply(\n",
+    "        lambda: affine_max(x1, x2, max_out),\n",
+    "        lambda: integration_payoff(max_out, integrator_reg, ind_reg),\n",
+    "    )\n",
+    "\n",
+    "\n",
+    "@qfunc\n",
+    "def main(\n",
+    "    x1: Output[QNum], \n",
+    "    x2: Output[QNum],\n",
+    "    integrator_reg: Output[QNum],\n",
+    "    ind_reg: Output[QBit],\n",
+    ") -> None:\n",
+    "    allocate_num(MAX_NUM_QUBITS, False, MAX_FRAC_PLACES, integrator_reg)\n",
+    "    allocate(NUM_QUBITS, x1)\n",
+    "    allocate(NUM_QUBITS, x2)\n",
+    "    allocate(1, ind_reg)\n",
+    "    rainbow_integration(x1, x2, integrator_reg, ind_reg)\n",
+    "\n",
+    "\n",
+    "constraints = Constraints(max_width=23)\n",
+    "qmod = create_model(main, constraints=constraints)\n",
+    "print(\"Starting synthesis\")\n",
+    "qprog = synthesize(qmod)\n",
+    "show(qprog)"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "markdown",
+   "id": "507817ef",
+   "metadata": {},
+   "source": [
+    "## IQAE algorithm"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "7e9554eb",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "from classiq import cfunc, Z\n",
+    "from classiq.qmod.builtins.classical_execution_primitives import iqae, save\n",
+    "\n",
+    "\n",
+    "@qfunc\n",
+    "def qmci_oracle(ind: QBit):\n",
+    "    Z(ind)\n",
+    "\n",
+    "@cfunc\n",
+    "def cmain():\n",
+    "    iqae_res = iqae(epsilon=0.05, alpha=0.1)\n",
+    "    save({\"iqae_res\": iqae_res})"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "code",
+   "id": "2420716d",
+   "metadata": {},
+   "source": [
+    "from classiq import CInt, QCallable, grover_operator, power, QArray\n",
+    "\n",
+    "@qfunc\n",
+    "def grover_algorithm(\n",
+    "    k: CInt,\n",
+    "    oracle_operand: QCallable[QArray[QBit]],\n",
+    "    sp_operand: QCallable[QArray[QBit]],\n",
+    "    x: QArray[QBit],\n",
+    "):\n",
+    "    sp_operand(x)\n",
+    "    power(k, lambda: grover_operator(oracle_operand, sp_operand, x))"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "code",
+   "id": "0c8c1dbc",
+   "metadata": {
+    "collapsed": false
+   },
+   "source": [
+    "def get_main():\n",
+    "    @qfunc\n",
+    "    def main(\n",
+    "        k: CInt,\n",
+    "        ind_reg: Output[QBit],\n",
+    "    ) -> None:\n",
+    "        full_reg = QArray(\"full_reg\")\n",
+    "        allocate(2 * NUM_QUBITS + MAX_NUM_QUBITS + 1, full_reg)\n",
+    "        grover_algorithm(\n",
+    "            k,\n",
+    "            lambda x: qmci_oracle(x[x.len - 1]),\n",
+    "            lambda x: rainbow_integration(\n",
+    "                x[0:NUM_QUBITS],\n",
+    "                x[NUM_QUBITS : 2 * NUM_QUBITS],\n",
+    "                x[2 * NUM_QUBITS : 2 * NUM_QUBITS + MAX_NUM_QUBITS],\n",
+    "                x[x.len - 1],\n",
+    "            ),\n",
+    "            full_reg,\n",
+    "        )\n",
+    "        state_reg = QArray(\"state_reg\")\n",
+    "        bind(full_reg, [state_reg, ind_reg])\n",
+    "\n",
+    "    return main"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "code",
+   "id": "19deb769",
+   "metadata": {},
+   "source": [
+    "from classiq import execute\n",
+    "\n",
+    "def synthesize_and_execute(post_process):\n",
+    "    constraints = Constraints(max_width=25)\n",
+    "    qmod = create_model(\n",
+    "        get_main(),\n",
+    "        constraints=constraints,\n",
+    "        classical_execution_function=cmain,\n",
+    "    )\n",
+    "    print(\"Starting synthesis\")\n",
+    "    qprog = synthesize(qmod)\n",
+    "    show(qprog)\n",
+    "    print(\"Starting execution\")\n",
+    "    res = execute(qprog).result()\n",
+    "    iqae_res = res[0].value\n",
+    "    print(\"raw iqae results:\", iqae_res.estimation, iqae_res.confidence_interval)\n",
+    "    parsed_res = post_process(res[0].value)\n",
+    "\n",
+    "    return (qmod, qprog, iqae_res, parsed_res)"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "markdown",
+   "id": "bd024abf",
+   "metadata": {},
+   "source": [
+    "## Post Process\n",
+    "\n",
+    "We need to add to the post-processing function a term:\n",
+    "\n",
+    "\\begin{equation}\n",
+    "\\begin{split}\n",
+    "\\mathbb{E} \\left[\\max\\left(\\frac{e^{a(x+1)} - 1}{e^{a(x_{max} +1)}-1}c + \\frac{1}{e^a} , Ke^{-b'}\\right)\\right]  e^{b'} - K \\\\\n",
+    "=\\mathbb{E} \\left[\\max\\left(\\frac{e^{a(x+1)} - 1}{e^{a(x_{max} +1)}-1}, \\frac{Ke^{-b'}}{c} - \\frac{e^{-a}}{c}\\right)\\right]ce^{b'} + e^{b'}e^{-a} - K\n",
+    "\\end{split}\n",
+    "\\end{equation}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "f62ac3e6",
+   "metadata": {},
+   "source": [
+    "def parse_result_integration(iqae_res):\n",
+    "    exp_rate = (1 / (2**MAX_FRAC_PLACES)) * a\n",
+    "    B = (np.exp((2**MAX_NUM_QUBITS) * exp_rate) - 1) / np.exp(exp_rate)\n",
+    "    A = 1 / np.exp(exp_rate)\n",
+    "    C = S0[0] * np.exp((MU[0] + MIN_X * CHOLESKY[0].sum()))\n",
+    "\n",
+    "    option_value = (iqae_res.estimation * (C * B)) + (C * A) - K\n",
+    "    confidence_interval = (\n",
+    "        (np.array(iqae_res.confidence_interval) * (C * B)) + (C * A) - K\n",
+    "    )\n",
+    "    return (option_value, confidence_interval)"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "markdown",
+   "id": "697fa83f",
+   "metadata": {},
+   "source": "# Run method"
+  },
+  {
+   "cell_type": "code",
+   "id": "fde9fc59",
+   "metadata": {},
+   "source": [
+    "qmod_integration, qprog_integration, iqae_res_integration, parsed_res_integration = (\n",
+    "    synthesize_and_execute(parse_result_integration)\n",
+    ")\n",
+    "result, conf_interval = parsed_res_integration\n",
+    "print(f\"raw iqae results: {iqae_res_integration.estimation} with confidence interval {iqae_res_integration.confidence_interval}\")\n",
+    "print(f\"option estimated value: {result} with confidence interval {conf_interval}\")"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "markdown",
+   "source": [
+    "## References\n",
+    "\n",
+    "<a name='QALROP'>[1]</a>: [Francesca Cibrario et al., Quantum Amplitude Loading for Rainbow Options Pricing. Preprint](https://arxiv.org/abs/2402.05574v2)"
+   ],
+   "id": "f859b2b3ac750050"
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": "",
+   "id": "75c5a77d0cc56f98",
+   "outputs": [],
+   "execution_count": null
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "ipykernel"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.14"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}