Price European options using Monte Carlo integration

# %pip install fire-opal qiskit_ibm_runtime qiskit qiskit-finance scipy matplotlib numpy qctrl-visualizer

import matplotlib.pyplot as plt
import numpy as np
from dataclasses import dataclass
from scipy.integrate import quad
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.qasm2 import dumps
import qctrlvisualizer as qv
import fireopal as fo
from qctrlvisualizer import QCTRL_STYLE_COLORS
from matplotlib.axes import Axes
from matplotlib.figure import Figure
from qctrlvisualizer.utils import create_figure
from qctrlvisualizer.style import qctrl_style
from qiskit_finance.circuit.library import LogNormalDistribution


plt.style.use(qv.get_qctrl_style())

import warnings

warnings.filterwarnings("ignore")

# Authenticate your Q-CTRL account.
api_key = "YOUR_Q-CTRL_API_KEY"
fo.authenticate_qctrl_account(api_key=api_key)

token = "YOUR_IBM_CLOUD_API_KEY"
instance = "YOUR_IBM_CRN"

credentials = fo.credentials.make_credentials_for_ibm_cloud(
    token=token, instance=instance
)

# Replace with a backend available to your IBM Quantum instance.
backend_name = "desired_backend"

@dataclass
class EuropeanOptionPricing:
    spot_price: float
    volatility: float
    interest_rate: float
    maturity_years: float

    @property
    def mu(self) -> float:
        return (
            self.interest_rate - 0.5 * self.volatility**2
        ) * self.maturity_years + np.log(self.spot_price)

    @property
    def sigma(self) -> float:
        return self.volatility * np.sqrt(self.maturity_years)

    @property
    def mean(self) -> float:
        return np.exp(self.mu + self.sigma**2 / 2)

    @property
    def variance(self) -> float:
        return (np.exp(self.sigma**2) - 1) * np.exp(2 * self.mu + self.sigma**2)

    @property
    def stddev(self) -> float:
        return np.sqrt(self.variance)

    def get_high_low(self) -> tuple[float, float]:
        low = round(np.maximum(0, self.mean - 3 * self.stddev), 3)
        high = round(self.mean + 3 * self.stddev, 3)
        return high, low

    def get_strike_window(self) -> tuple[float, float]:
        high, low = self.get_high_low()

        strike_low = (high - low) / 2 + low
        strike_high = (high - low) / 4 + strike_low

        strike_window = (strike_low, strike_high)

        return strike_window

    def get_uncertainty_model(self, num_qubits: int) -> LogNormalDistribution:
        high, low = self.get_high_low()
        uncertainty_model = LogNormalDistribution(
            num_qubits=num_qubits, mu=self.mu, sigma=self.sigma**2, bounds=(low, high)
        )
        return uncertainty_model

    def get_exact_solution(self, return_error: bool = False):
        def lognormal_pdf(x, mu, sigma):
            if x <= 0:
                return 0
            return (1 / (x * sigma * np.sqrt(2 * np.pi))) * np.exp(
                -((np.log(x) - mu) ** 2) / (2 * sigma**2)
            )

        strike_window = self.get_strike_window()

        result, error = quad(
            lognormal_pdf,
            strike_window[0],
            strike_window[1],
            args=(self.mu, self.sigma),
        )
        if return_error:
            return result, error

        else:
            return result

    def get_exact_solution_with_discretized_values(self, num_qubits) -> float:
        uncertainty_model = self.get_uncertainty_model(num_qubits=num_qubits)
        low, high = self.get_strike_window()
        values = uncertainty_model.values
        probabilities = uncertainty_model.probabilities

        return np.sum(
            np.where((values >= low) & (values <= high), 1, 0) * probabilities
        ).item()

    def get_instance_parameters(self) -> dict:
        parameters = {}
        parameters["Spot price"] = self.spot_price
        parameters["Interest rate"] = self.interest_rate
        parameters["Maturity years"] = self.maturity_years
        parameters["Volatility"] = self.volatility
        return parameters


def visualize_option_pricing_problem(num_qubits: int, problem: EuropeanOptionPricing):
    uncertainty_model = problem.get_uncertainty_model(num_qubits=num_qubits)
    strike_window = problem.get_strike_window()

    # Plot distribution
    x = list(uncertainty_model.values)
    y = uncertainty_model.probabilities

    if len(x) > 20:
        for i in range(1, 10):
            if len(x[::i]) <= 20:
                x_ticks = x[::i]
                break
    else:
        x_ticks = x

    plt.plot(x, y)
    plt.xticks(x_ticks, size=15, rotation=90)
    plt.yticks(size=15)
    plt.xlabel("Spot price at maturity $S\_T$ (\$)", size=15)
    plt.ylabel("Probability", size=15)
    plt.ylim((0, max(uncertainty\_model.probabilities) * 1.1))
    plt.xlim((x[0], x[-1]))
    plt.tick\_params(axis="y", labelcolor=QCTRL\_STYLE\_COLORS[0])
    plt.vlines(
        x=x,
        ymax=max(uncertainty\_model.probabilities) * 1.1,
        ymin=0,
        color="lightgray",
        linewidth=1,
        zorder=1,
    )

    # Create payoff window
    y = [
        1 if value >= strike\_window[0] and value < strike\_window[-1] else 0
        for value in x
    ]
    first = y.index(1)
    last = len(y) - y[::-1].index(1) - 1
    y.insert(last + 1, 0)
    y.insert(first, 0)
    x.insert(last, x[last])
    x.insert(first, x[first])

    plt.twinx()
    plt.plot(x, y, "r-")
    plt.tick\_params(axis="y", labelcolor=QCTRL\_STYLE\_COLORS[1])
    plt.ylabel("Payoff", size=15)
    plt.xticks(x\_ticks, size=15, rotation=90)

    plt.yticks(size=15)
    plt.ylim((0, 1.1))
    plt.fill\_between(x, y, 0, alpha=0.2, color=QCTRL\_STYLE\_COLORS[1])
    plt.yticks([0, 1])
    plt.show()

num\_uncertainty\_qubits = 5  # Number of qubits to discretize the distribution
S = 2.0  # Initial spot price
vol = 0.6  # Annualized volatility of 60%
r = 0.01  # Annual interest rate of 1%
T = 10 / 365  # 10 days to maturity

problem = EuropeanOptionPricing(
    spot\_price=S, volatility=vol, interest\_rate=r, maturity\_years=T
)
visualize\_option\_pricing\_problem(num\_qubits=num\_uncertainty\_qubits, problem=problem)

def oracle(
    num\_qubits: int, strike\_window: tuple[float, float], high: float, low: float
) -> QuantumCircuit:
    """
    Creates the payoff marking circuit for our option pricing use case.

    Parameters
    ----------
    num\_qubits : int
        The number of data qubits.
    strike\_window : tuple[float, float]
        The low and high asset price for the option to pay out.
    high : float
        The maximum asset price to encode.
    low : float
        The minimum asset price to encode.

    Returns
    -------
    QuantumCircuit
        The quantum circuit to compute U|x>|y> = |x>|y + f(x) mod 2>,
        with n+1 qubits.
    """
    if strike\_window[1] < high:
        is\_window = True

    else:
        is\_window = False

    slope = (2**num\_qubits - 1) / (high - low)
    strike\_window = [
        int(np.ceil(slope * (strike\_price - low))) for strike\_price in strike\_window
    ]
    min\_strike\_bin = format(strike\_window[0], f"0{num\_qubits}b")[::-1]
    max\_strike\_bin = format(strike\_window[1], f"0{num\_qubits}b")[::-1]

    if is\_window:
        flip\_qubits = [
            i
            for i in range(len(max\_strike\_bin))
            if max\_strike\_bin[i] != min\_strike\_bin[i]
        ]
    else:
        flip\_qubits = []

    control\_qubits = [i for i in range(len(min\_strike\_bin)) if min\_strike\_bin[i] == "1"]

    qc = QuantumCircuit(num\_qubits + 1)
    if is\_window:
        qc.x(flip\_qubits)
    qc.mcx(control\_qubits + flip\_qubits, num\_qubits)
    if is\_window:
        qc.x(flip\_qubits)
    qc.name = "f(x)"
    return qc


def get\_preparation\_circuit(
    num\_qubits: int, problem: EuropeanOptionPricing
) -> QuantumCircuit:
    """
    Construct the circuit to encode the target probability, for the MCI function input.

    Parameters
    ----------
    num\_qubits: int
        Number of qubits to discretize the probability distribution.
    problem: EuropeanOptionPricing
        European option pricing problem.

    Returns
    -------
    QuantumCircuit
        State preparation circuit.
    """

    uncertainty\_model = problem.get\_uncertainty\_model(num\_qubits=num\_qubits)
    high, low = problem.get\_high\_low()
    strike\_window = problem.get\_strike\_window()

    spot\_price = QuantumRegister(num\_qubits, name="spot\_price")
    output\_qubit = QuantumRegister(1, name="output")
    circuit = QuantumCircuit(spot\_price, output\_qubit)

    circuit.compose(uncertainty\_model, qubits=spot\_price, inplace=True)
    circuit.barrier()
    circuit.compose(
        oracle(num\_qubits=num\_qubits, strike\_window=strike\_window, high=high, low=low),
        qubits=spot\_price[:] + output\_qubit[:],
        inplace=True,
    )
    circuit.barrier()

    return circuit

max\_iteration\_count = 20

preparation\_circuit = get\_preparation\_circuit(
    num\_qubits=num\_uncertainty\_qubits, problem=problem
)

results\_scenario\_1 = fo.integrate\_monte\_carlo(
    problem={"preparation\_circuit": dumps(preparation\_circuit.decompose(reps=5))},
    integrator\_options={
        "max\_iteration\_count": max\_iteration\_count,
        "full\_output": True,
    },
    credentials=credentials,
    backend\_name=backend\_name,
)

results\_scenario\_1.status()

estimates = [
    result["estimate"]
    for result in results\_scenario\_1.result()["metadata"]["iterations"]
]
variance = [
    result["variance"]
    for result in results\_scenario\_1.result()["metadata"]["iterations"]
]
iteration\_count = [
    result["iteration\_count"]
    for result in results\_scenario\_1.result()["metadata"]["iterations"]
]
evaluation\_count = [
    result["evaluation\_count"]
    for result in results\_scenario\_1.result()["metadata"]["iterations"]
]

@qctrl\_style()
def plot\_estimation\_history(
    histories: dict,
    figure: Figure | None = None,
    known\_value: float | None = None,
    x\_axis\_log: bool = False,
) -> None:
    """
    Plot the history of estimates and confidence intervals.

    Parameters
    ----------
    histories : dict[str, Any]
        A dictionary where keys are labels for different integrators
        and values are tuples of evaluation\_counts, estimates, and variance.
    figure : Figure or None, optional
        A matplotlib Figure object for the plot. If None, a new figure is created,
        by default None.
    known\_value : float or None, optional
        The known value of the integral, if known. If provided, a horizontal
        dashed line will be drawn at this value, by default None.
    x\_axis\_log : bool, optional
        If True, the x-axis will be set to a logarithmic scale, by default False.
    """
    figure = create\_figure(figure)
    axes = figure.subplots(nrows=1, ncols=1)
    assert isinstance(axes, Axes)

    if known\_value is not None:
        axes.axhline(
            known\_value, linestyle="dashed", color="black", alpha=0.2, label="Known"
        )

    for label, (eval\_count, estimate, variance) in histories.items():
        std\_err = 1.96 * np.sqrt(variance)
        axes.fill\_between(eval\_count, estimate - std\_err, estimate + std\_err, alpha=0.1)
        axes.plot(eval\_count, estimate, "o-", label=label, markersize=2, linewidth=1)

    axes.legend()
    if x\_axis\_log:
        axes.set\_xscale("log")

    axes.set\_xlabel("Evaluation count")
    axes.set\_ylabel("Estimate")

plot\_estimation\_history(
    histories={"QMCI": (evaluation\_count, estimates, variance)},
    known\_value=problem.get\_exact\_solution\_with\_discretized\_values(
        num\_uncertainty\_qubits
    ),
)
plt.title(
    f'Data Qubits: {num\_uncertainty\_qubits+1} | 2Q gate counts: {preparation\_circuit.decompose(reps=5).count\_ops()["cx"]} \nmax iterations: {max\_iteration\_count} | backend: {backend\_name}'
)
plt.show()

high, low = problem.get_high_low()
strike_price = 2.05
num_uncertainty_qubits = 3
distribution = fo.make_lognormal_distribution(
    uncertainty=num_uncertainty_qubits,
    mu=problem.mu,
    sigma=problem.sigma**2,
    bounds=(low, high),
)
payoff_function = fo.make_european_call_delta_objective(
    bounds=(low, high), strike_price=strike_price
)


def get_exact_solution_delta_objective(
    num_qubits: int, problem: EuropeanOptionPricing, strike_price: float
):
    uncertainty_model = problem.get_uncertainty_model(num_qubits)
    return (
        np.exp(-problem.interest_rate * problem.maturity_years)
        * np.sum(
            uncertainty_model.probabilities[uncertainty_model.values > strike_price]
        ).item()
    )

results_scenario_2 = fo.integrate_monte_carlo(
    problem=fo.make_monte_carlo_problem(
        objective_function_parameters=payoff_function,
        distribution_parameters=distribution,
    ),
    credentials=credentials,
    backend_name=backend_name,
    integrator_options={
        "max_iteration_count": max_iteration_count,
        "full_output": True,
    },
)

estimates = [
    result["estimate"]
    for result in results_scenario_2.result()["metadata"]["iterations"]
]
variance = [
    result["variance"]
    for result in results_scenario_2.result()["metadata"]["iterations"]
]
iteration_count = [
    result["iteration_count"]
    for result in results_scenario_2.result()["metadata"]["iterations"]
]
evaluation_count = [
    result["evaluation_count"]
    for result in results_scenario_2.result()["metadata"]["iterations"]
]

plot_estimation_history(
    histories={"QMCI": (evaluation_count, estimates, variance)},
    known_value=get_exact_solution_delta_objective(
        num_qubits=num_uncertainty_qubits, problem=problem, strike_price=strike_price
    ),
)
plt.title(
    f"Data Qubits: {num_uncertainty_qubits*2} \nmax iterations: {max_iteration_count} | backend: {backend_name}"
)
plt.show()

from fireopal import print_package_versions

print_package_versions()