Quantum computational finance: Monte Carlo pricing of financial derivatives

This work presents a quantum algorithm for the Monte Carlo pricing of financial derivatives. We show how the relevant probability distributions can be prepared in quantum superposition, the payoff functions can be implemented via quantum circuits, and the price of financial derivatives can be extracted via quantum measurements. We show how the amplitude estimation algorithm can be applied to achieve a quadratic quantum speedup in the number of steps required to obtain an estimate for the price with high confidence. This work provides a starting point for further research at the interface of quantum computing and finance.

Figure 1

The price of a stock in the Black-Scholes-Merton model behaves stochastically as a geometric Brownian motion, changing each time step according to a log-normal distribution. Here, five sample price evolutions (in dollars) of a single stock are plotted as a function of time (in days). The resultant distribution is log-normally distributed, with mean (dashed black line) and one standard deviation (dotted red lines) illustrated. Pricing an option requires estimating the expected value of a payoff function on the stock at various times. The parameters are initial price $S_0=\$3$, drift $\alpha =0.1$, and volatility $\sigma =0.25$.

Figure 2

Using amplitude estimation for the quantum Monte Carlo pricing of financial derivatives. (a) The $n+1$ qubit phase estimation unitary is written in terms of $\mathcal{F}:=\mathcal{R}(\mathcal{A}\otimes {\mathcal{I}}_2)$, and the simple rotation unitaries $\mathcal{Z}:={\mathcal{I}}_{2^{n+1}}-2|0^{n+1}\rangle \langle 0^{n+1}|$ and $\mathcal{V}:={\mathcal{I}}_{2^{n+1}}-2{\mathcal{I}}_{2^n}\otimes |1\rangle \langle 1|$. (b) A visualization of the action of $\mathcal{Q}:=\mathcal{U}\mathcal{S}$, with $\mathcal{S}=\mathcal{V}\mathcal{U}\mathcal{V}$ and $\mathcal{U}=\mathcal{F}\mathcal{Z}{\mathcal{F}}^{†}$, on an arbitrary state $|\psi \rangle$ (red) in the span of $|\chi \rangle$ and $\mathcal{V}|\chi \rangle$. First, the action of $-\mathcal{S}$ on $|\psi \rangle$ is to reflect along $\mathcal{V}|\chi \rangle$, resulting in the intermediate $-\mathcal{S}|\psi \rangle$ (amber). Then, $-\mathcal{U}$ acts on $-\mathcal{S}|\psi \rangle$ by reflecting along $|\chi \rangle$. The resultant state $\mathcal{Q}|\psi \rangle$ (green) has been rotated anticlockwise by an angle $2\theta$ in the hyperplane of $|\chi \rangle$ and $|{\chi }^{\perp }\rangle$. (c) The phase estimation circuit. Here, $\mathcal{A}$ encodes the randomness by preparing a superposition in $|x\rangle$, while $\mathcal{R}$ encodes the random variable into the $|1\rangle$ state of an ancilla qubit according to Eq. (17). The output after both steps is the multiqubit state $|\chi \rangle$. Amplitude estimation then proceeds by invoking phase estimation to encode the rotation angle $\theta$ in a register of quantum bits that are measured to obtain the estimate $\widehat{\theta }$. (d) For pricing a European call option, the superposition prepared by $\mathcal{A}$ [or equivalently $\mathcal{G}$ in Eq. (35)] is a discretization of the normal distribution in $x$ with a fixed cutoff (e.g., $c=4$), approximating the Brownian motion of the underlying asset. In this case, $\mathcal{R}$ encodes the call option payoff in Eq. (5).

Figure 3

Scaling of the error in classical and quantum MC methods [defined in Eq. (57)] plotted against number of MC steps for a European call option in log-log scale with $S_0=\$100,K=\$50,r=0.05,\sigma =0.2,T=1$, and $D=24$. Subscripts $C$ and $Q$ denote the errors from classical MC and quantum phase estimation, respectively. Evidently, the error for the quantum algorithm (with a fitted slope of ${\zeta }_{\mathrm{Q}}=-0.982$) scales almost quadratically faster than the classical MC method (which has ${\zeta }_{\mathrm{C}}=-0.5$). The theoretical upper bound on the error in quantum algorithm is shown by the solid green curve, which corresponds to ${\zeta }_{\mathrm{Q}}=-1$.

Figure 4

Ratio of the quantum to classical scaling exponents. The quantum scaling is obtained by fitting the simulations results to the power law in Eq. (58), while the classical scaling exponent is taken to be $-0.5$. Results are plotted with varying strike price $K$ (in dollars) and fixing other parameters to be the same as in Fig. 3. An almost quadratic speedup is obtained for all chosen values of $K$.