Options pricing using a quantum Monte Carlo algorithm on Amazon Braket

Financial trading strategies involve calculating the value of financial instruments such as options, bonds, and interest rate derivatives. These calculations are typically governed by stochastic differential equations, often with no closed-form solutions. Although an exact solution may not be possible, an approximate, nearly optimal solution can still be found using Monte Carlo techniques. Monte Carlo estimation relies on taking random samples to simulate a vast number of possible scenarios. These scenarios are then averaged to evaluate the expected value of the underlying financial asset. While Monte Carlo methods are useful and ubiquitous, they are computationally expensive and time-consuming to run, preventing traders from estimating the values of financial instruments in near real-time.

In this post, we explore a quantum algorithm [1] for Monte Carlo estimation (QMC) for options pricing using Amazon Braket. In principle, QMC promises a quadratic speedup over its classical counterpart. This speedup, however, is out of reach for the current generation of quantum hardware, as QMC requires fault-tolerant devices with stringent requirements [2]. Nevertheless, QMC is a great example to illustrate how quantum algorithms can provide speedups over classical methods and to investigate the error scaling at a hands-on, small-scale example [3].

Pricing an option

In finance, a derivative is defined as a contract between two involved parties who can benefit from varying prices of an underlying asset, such as a stock or a bond. An option is a type of derivative that allows the holder the right to buy (call option) or sell (put option) the underlying asset at a pre-agreed price (the strike price) at or within a specified time frame (called the maturity time or exercise window). The holder only pays for the right to exercise the option, which allows the holder the choice of whether to exercise it at maturity time. If the holder chooses to exercise the option, the seller must buy or sell the underlying asset as agreed in the contract. Pricing an option refers to determining the amount the holder must pay to purchase the contract. The payoff function of an option, with S_T being the value of the underlying asset at time T and K being the strike price, is given as follows:

The task of pricing an option then becomes to determine the expected value at present time of the payoff function defined earlier. In simple cases, under a few assumptions, the Black-Scholes-Merton (BSM) model [4, 5] can be used to determine the value of the options, thus determining the appropriate price. The asset price S(t) is treated as a random variable X that follows a known or often assumed to follow a log-normal stochastic process. The BSM differential equation is given as follows:

where, X_t is the asset price at time t, a is called the “drift”, s is the volatility, and dW_t is the increment of an accompanying Brownian motion. In simple cases, with constant drift and volatility, equation (2) can be analytically solved. However, beyond simple cases, one needs to rely on running Monte Carlo sampling to solve equation (2) and determine the asset price.

Going forward, we will work with the example of an Asian call option which has the following payoff function:

where K is the strike price and A_T is the average asset price computed over the maturity time T. Here we will use the following geometric averaging method to compute A_T:

For constant drift and volatility, the A_T in eq. (4) can be analytically computed. Aside from geometric averaging, Asian payoff can also be computed using arithmetic averaging which cannot be analytically solved and one needs to use numerical methods. In this work, we will work with the geometric case for simplicity.

Classical Monte Carlo approach and computational challenges

In the Monte Carlo approach, instead of solving the associated SDE, one directly samples from the random process representing the underlying asset. To do this, one generates a large number of pricing trajectories over a period of time and estimates the value of the expected payoff by taking the weighted average over the trajectories. As a result, one obtains an estimate of the expected payoff function with associated error tolerance.

Let’s consider a random variable (X, x), where p(x) denotes the probability of the sample x of X and f(x) is a function for which we are interested in computing the average. Then the Monte Carlo approach evaluates the expected mean value using the following finite sum:

where the summation is performed over the number of Monte Carlo samples.

In options pricing, the random variable X represents the asset price, x_i becomes a sample of the asset price, representing one of the many possible financial scenarios, p(x_i) is the probability of each sample, and f(x_i) is the individual payoff (as shown in equation 3) for each discrete sample of the asset price. This way, the summation in equation (6) can be used to determine the average payoff function covering many possible financial scenarios represented by a number of Monte Carlo samples. Since this is still an approximation method, there is an associated estimation error which scales with the number of classical Monte Carlo (CMC) sample, N_CMC as follows:

This means that to achieve an estimation error e, the classical Monte Carlo sample grows as 1/e². Hence, to improve the accuracy of the estimation by a factor of 10, one needs to increase the number of classical Monte Carlo samples by a factor of 100, resulting in significant computational overhead. While pricing an individual asset might take only a few seconds, traders typically need to simulate a combination of a large number of assets influenced by different risk scenarios. This significantly increases the runtime and computational resource requirements if one wishes to perform computations with higher accuracy. Naturally, faster pricing estimates are highly desirable for more accurate and timely trading.

Speedup with Quantum Monte Carlo

It is theoretically possible to achieve a quadratic speedup in the runtime of Quantum Monte Carlo estimation, which uses the quantum amplitude estimation (QAE) algorithm. The essential idea is that the number of samples required in QAE to achieve a maximum estimation error of e, grows proportional to 1/e as opposed to the classical scaling of 1/e² [1,2,3]. This improves the estimation error scaling as follows:

where N_QMC refers to the number of Monte Carlo samples in the quantum algorithm. Hence, theoretically, the quantum approach would require quadratically fewer Monte Carlo samples to achieve the same level of estimation error. If a classical Monte Carlo computation requires 10,000 steps, the quantum algorithm would require only around 100 steps to achieve the same error scaling between accuracy and speed (up to a multiplicative factor). This potential speedup could enable up to 100 times improvement in speed and accuracy of options pricing, assuming equivalent Monte Carlo sample times for classical and quantum algorithms.

Whether a practical quantum computer could deliver on this speedup is another matter. It is important to note that these theoretical speedups are subject to certain assumptions and may not fully translate to practical implementations. For instance, options are more complex than the simple models considered here. Payoff functions of some of the options are path dependent and require considering the trajectory of the random variables in each Monte Carlo sample. Furthermore, quantum computers typically have slower clock speeds which, when considered, can push the problem size to much larger scale to achieve practical advantage. For a more detailed discussion of these caveats, see reference [2]. Nonetheless, the potential quantum advantage for Monte Carlo estimation remains promising, and exploring quantum-assisted workflows for financial modeling can provide valuable insights and opportunities for improving computational performance in the long term.

General Structure of the Quantum Monte Carlo algorithm

Figure 1: The unitary-circuit diagram of Quantum algorithm for Monte Carlo integration. Unitaries ‘A’ and ‘R’ are used to embed the underlying option. Repeated application of ‘Q’ controlled on ancilla qubits provide the quantum speedup (From reference [3]).

1. The first step is discretizing the values (x) of the random variable X (which represents the asset price). To denote N discrete values of the quantum state |x> we need n_d = log₂ (N) This is shown by the set of qubits in the top left circuit block in figure 1. Increasing n_d ultimately increases the number of discrete points used to represent the financial option which increases the discretization accuracy of the quantum algorithm. Note that this discretization accuracy is different from the accuracy of the Monte Carlo averaging dictated by e in equations 7 and 8.
2. Then the probabilities of the future prices of the asset are encoded into the amplitudes of a quantum state:This is achieved by the unitary operator shown in figure 1.
3. Payoff function in encoded in the form of angle q_x, so that sin(q_x) = f(x)^1/2 giving the quantum state to be:
   Encoding of the payoff function is obtained by applying the circuit denoted by unitary shown in figure 1.
4. The amplitude f(x)^1/2 is then encoded into amplitudes of n_pe ancilla qubits (shown by the lower half of the qubit in figure 1) to give the state:
   
5. The probability of measuring the previous state in |0> is the expectation value of the payoff function we are looking for.
6. Quantum Amplitude estimation [6] is then implemented using n_pe ancilla qubits through repeat application of unitary matrix Q (shown in figure 1) to estimate the amplitude of |0>. This repeated application of Q is equivalent to taking as many as N_QMC= 2^npe-1 ​Monte Carlo samples.

As discussed in reference [3], the estimation error for QMC variance scales as e µ 1/2^npe-1. Thus, applying unitary as many times as N_QMC = 2^npe-1 , the error scales as e µ 1/ N_QMC.This is the origin of the quadratic improvement in the error scaling. As a result, the runtime in quantum algorithm scales as 1/e which, compared to the runtime for the classical Monte Carlo algorithm (as 1/e²​), is quadratically faster in principle.

Options pricing using Quantum Monte Carlo

So far, we have discussed the theoretical aspects. Let’s now see how we can implement this on Amazon Braket and obtain numerical evidence for the expected speedup. The accompanying notebook demonstrates a quick start example of using Braket to price an Asian call option, showcasing the potential quadratic speedup of the quantum workflow compared to the classical approach.

The Python code snippet that follows defines the Asian option payoff function for a given spot price, strike price, interest rate, and volatility. In a full quantum implementation, the payoff function and the probability distributions will be represented by unitary matrices ‘A’ and ‘R’ as shown in figure 1 earlier. For demonstration purposes, we write a Python function for the classical payoff calculation and use the PennyLane quantum simulator for the quantum amplitude estimation which, as explained before, is the origin of quadratic quantum advantage. By simulating the algorithms in this blog, we can see the theoretically expected speedup in numerical experiments with the Asian option.

import pennylane as qml
import numpy as np
import scipy
import scipy.stats
import itertools
from time import time
import matplotlib.pyplot as plt
def asian_payoff_func (spot, strike, rate, volatility, time, x):
    Dt = time / periods
    payoffs = []
    for item in x:
        price = [spot]
        price.append(price[-1] * np.exp(volatility * item
                    + (rate - 0.5 * volatility ** 2) * Dt))
        price = np.expand_dims(np.array(price[1:]), axis=0)
        payoff = asian_call_payoffs(price, strike, spot)[0]
        payoffs.append(payoff)
    if np.max(payoffs) == 0:
        payoffs[0] = 1e-10
    return np.array(payoffs)*np.exp(- rate * time)

The time averaging can be performed as given in the following code snippet:

def asian_call_payoffs(paths, strike, spot, payoff_type = 'geometric'):

    spots = np.full((paths.shape[0], 1), spot)
    paths = np.append(spots, paths, axis=1)

    if payoff_type == 'geometric':
        means = scipy.stats.mstats.gmean(paths, axis=1)
    elif payoff_type == 'arithmetic':
        means = np.mean(paths, axis=1)
    else:
        raise Exception('Error: payoff_type must be either geometric or arithmetic')
    asian_payoffs = means - strike
    asian_payoffs[asian_payoffs < 0] = 0
    return asian_payoffs

Running the simulation experiment using the PennyLane Simulator

Now we are ready to test the algorithm. The quantum computing hardware available today is not yet capable of running this algorithm at scale and we would need a larger, less noisy quantum computer capable of fault tolerance and error correction to run this experiment on real hardware. In the meantime, we use classical simulators to simulate ideal quantum devices and understand the performance factors of these algorithms. Here we will use PennyLane’s QuantumMonteCarlo module and lightning.qubit simulator using the code snippet that follows:

def QMC_Asian_call(spot, strike, rate, volatility, time, n_disc, variance, cutoff_factor, n_pe):
    N_disc, N_pe = 2**n_disc, 2**n_pe
    x, p = discrete_normal(n_disc, T, cutoff_factor)
    asian_payoff = asian_payoff_func(spot, strike, rate, volatility, time, x)
    renormalization = max(asian_payoff)
    payoff_func = lambda i: asian_payoff[i]/renormalization
    target_wires = range(n_disc+1)
    estimation_wires = range(n_disc+1, n_disc+n_pe+1)
    
    dev=qml.device("lightning.qubit",wires=(n_disc+n_pe+1))
    @qml.qnode(dev)
    def circuit():
        qml.templates.QuantumMonteCarlo(
            p,
            payoff_func,
            target_wires=target_wires,
            estimation_wires=estimation_wires,
        )
        return qml.probs(estimation_wires)
    phase_estimated=np.argmax(circuit()[:int(N_pe/2)])/N_pe
    estimated_value_qmc = (1 - np.cos(np.pi * phase_estimated)) / 2 * renormalization
    
    return np.array(estimated_value_qmc)

The function QMC_Asian_call can be run by simply calling it as follows:

results = QMC_Asian_call (spot_price, strike_price, risk_free_rate,  volatility,  maturity_period, n_d, variance, cutoff_factor, n_p)

Note that in this example we picked PennyLane’s lightning.qubit simulator. We used the Hybrid Jobs feature of Amazon Braket to parallelize the runs with PennyLane [7] , running multiple benchmarking experiments and at the same time reduce the overall wait time.

Results

The performance of the algorithm covers three factors: error scaling comparing classical to quantum Monte Carlo for a baseline set of prices, error scaling across a range of discretization, and the wall-time to complete the experiments.

Figure 2 that follows shows the comparison of error scaling obtained from classical and quantum algorithms for Monte Carlo estimation for Spot price $100, Strike price $100, risk-free rate 0.05 %, volatility of 0.2 and maturity time one year with two time-intervals for averaging. We quantify the overall scaling by fitting the error to a function y = A x^z, where  is the coefficient and z is the scaling exponent. For classical Monte Carlo the scaling exponent is known to be z_C =-0.5.  The fitting reveals that for QMC, z_Q =-0.986, giving rise to a nearly quadratic speedup ofz_Q /z_c = 1.972.

Figure 2: Scaling of error estimation in classical and quantum algorithms for Monte Carlo integration applied to Asian Call option with Spot price $100, Strike price $100, risk free rate 0.05 %, Volatility of 0.2 and maturity time one year with two time-intervals for averaging. The discretization qubits was fixed at nd=6.

The exact scaling depends on a variety of factors, but we can still show that it remains roughly quadratic by varying other factors of the algorithm. In the preceding figure 2, we fixed the discretization qubit at n_d =6. Let’s now consider a variety of value for n_d and explore the theoretical quadratic speedup quantified by z_Q /z_c. In the following figure 3, we repeat the Asian call option estimation for a range of values of n_d and plot the ratio of quantum to classical scaling exponents z_Q /z_c. A nearly quadratic speedup is achieved for larger values of n_d, whereas smaller values result in a less-than-quadratic speedup. This is expected as lower values of n_d can cause higher errors due to discretization which may in turn impact the overall speedup of QMC.

Figure 3: Ratio of quantum to classical error scaling for a range of discretization qubits nd. An almost quadratic speedup is obtained for all strike prices nd.

Figure 4: Scaling of simulation time plotted against varying phase estimation qubits (npe) for five different values of discretization qubits (nd).

Let us now compare the runtime of the simulations performed here. Figure 4 shows the computational wall-time plotted against n_pe for different values of n_d​​. The wall-time grows exponentially with increasing n_pe  for each n_d ​​. Also, for a fixed n_pe ​​, a higher n_d​​  results in longer simulation wall-time.

Conclusion

In this blog post, we have demonstrated how a quantum algorithm for Monte Carlo can be run on Amazon Braket to estimate the price of an example Asian call option. The simulation experiments performed here demonstrate the theoretically predicted quadratic advantage of quantum Monte Carlo over its classical counterpart in estimation error scaling. This approach can in principle be applied to other derivative financial instrument pricing which are usually estimated using classical Monte Carlo methods.

We used PennyLane’s QuantumMonteCarlo module to simulate the QMC algorithm on PennyLane’s lightning.qubit  simulator. Amazon Braket provides adequate capabilities to further push the limits of these simulations, such as:

• Use Hybrid Jobs to parallelize PennyLane simulations,
• Use larger AWS computational instances to accommodate higher values of n_d and n_pe, thereby curbing discretization error and increasing number of QMC samples, and
• Use GPU instances with lightning.gpu, to speed up each simulation.

References

[1] Montanaro Ashley,Quantum speedup of Monte Carlo methods, Proc. R. Soc A.471: 20150301.http://doi.org/10.1098/rspa.2015.0301

[2] Alexander Dalzel et.al. AWS Center for Quantum Computing. “Quantum algorithms: A survey of applications and end-to-end complexities.” arXiv:2310.03011 (2023). https://arxiv.org/abs/2310.03011.

[3] Patrick Rebentrost, Brajesh Gupt, and Thomas R. Bromle, “Quantum computational finance: Monte Carlo pricing of financial derivatives” Phys. Rev. A 98, 022321 – Published 20 August 2018.

[4] F. Black and M. Scholes, “The Pricing of Options and Corporate Liabilities” Journal of Political Economy 81, 637 (1973).

[5] R. C. Merton, “Theory of Rational Option Pricing”, The Bell Journal of Economics and Management Science 4, 141 (1973).

[6] Gilles Brassard, Peter Hoyer, Michele Mosca, Alain Tapp, Quantum Computation and Quantum Information, Samuel J. Lomonaco, Jr. (editor), AMS Contemporary Mathematics, 305:53-74, 2002

[7] Matthew Beach, “Getting started with Amazon Braket Hybrid Jobs”, https://pennylane.ai/qml/demos/getting_started_with_hybrid_jobs/