Pricing Product Options and Using Them to Complete Markets for Functions of Two Underlying Asset Prices ^†

Abstract

Options paying the product of put and/or call option payouts at different strikes on two underlying assets are observed to synthesize joint densities and replicate differentiable functions of two underlying asset prices. The pricing of such options is undertaken from three perspectives. The first perspective uses a geometric two-dimensional Brownian motion model. The second inverts two-dimensional characteristic functions. The third uses a bootstrapped physical measure to propose a risk charge minimizing hedge using options on the two underlying assets. The options are priced at the cost of the hedge plus the risk charge.

1. Introduction

Option markets considerably enhance the collection of functions of the stock price at the option maturity that may be purchased or sold. In the absence of options, the only functions one can access are straight lines with bonds delivering the intercept and stocks the slope. Options allow one to change the slope and access any twice differentiable function of the stock price as shown, for example, in Carr and Madan (1998).

Functions that have been priced this way include the logarithm of the S&P 500 index a month later in arriving at the VIX index. The CBOE skew prices the first three powers of the logarithm to build the Skew index. Others have used the procedure to study the risk neutral kurtosis and its effects (Bakshi et al. (2003)). It is also known from Breeden and Litzenberger (1978) that option prices synthesize risk neutral densities and permit the recovery of the density from the second derivative of option prices.

These results have natural extensions to two dimensions when one considers functions of a pair of stock prices at maturity. Functions differentiable four times may be replicated provided one trades product options. The payoff on a product option is the product of payoffs on two standard options be they puts or calls. Furthermore, joint densities may be synthesized from the prices of such product options. Consider, for example, a product call with strikes $K_1,K_2.$ For stock prices $S_1,S_2$ at maturity, the payoff or value at maturity T is

$$C(S_1,S_2,T;K_1,K_2)={\left(S_1-K_1\right)}^{+}{\left(S_2-K_2\right)}^{+}.$$

(1)

For a risk neutral density $f(S_1,S_2)$, the price of the product call $w(K_1,K_2)$ for an interest rate r is given by

$$w\left(K_1,K_2\right)=e^{-rT}{\int }_{K_1}^{\infty }{\int }_{K_2}^{\infty }\left(S_1-K_1\right)\left(S_2-K_2\right)f\left(S_1,S_2\right)d{S}_{1}d{S}_2.$$

Differentiation shows that

$$\frac{{\partial }^4}{\partial K_1^2\partial K_2^2}{e}^{rT}w(K_1,K_2)=f(K_1,K_2).$$

(2)

Other approaches to such recovery using basket option prices include Lipton (2001) and Carr and Laurence (2011).

Product options, therefore, play a critical role when one considers functions of pairs of prices. We, therefore, consider the pricing of such options. The first context is a simple generalization of the Black and Scholes (1973) and Madan (1973) geometric Brownian motion model to a bivariate normal density for the logarithm of the two stock prices. The product options are then priced using bivariate normal distribution functions and an assumed correlation value. The bivariate normal has normal marginals that do not fit the marginal risk neutral densities.

However, there is a large class of pure jump Lévy process models that fit option prices at each maturity. Among them is the three parameter variance gamma model of Madan and Seneta (1990), Madan et al. (1998). Recently, this model has been generalized to the four parameter bilateral gamma model of Küchler and Tappe (2008) that also fits marginal risk neutral densities. The bilateral gamma model is used to fit the two risk neutral marginals. Recent work reported in Madan and Schoutens (2020), Madan (2020) and Madan and Wang (2020a) formed multivariate densities consistent with prespecified bilateral gamma models with two additional dependency parameters. The multivariate bilateral gamma model has a closed form characteristic function.

The methods of Carr and Madan (1999) are extended to employ the two dimensional fast Fourier transform to price product options. In addition, methods reported in Madan and Wang (2020b) may be employed to price product options from the physical measures using the options in the two markets as hedge instruments.

The outline of the rest of the paper is as follows. Section 2 develops the two dimensional replication result. Section 3 presents the bivariate log normal pricing of product options. The two dimensional fast Fourier inversion of suitably modified product options is taken up in Section 4. The multivariate bilateral gamma model is described in Section 5. Sample computations and the effects on product option prices of changes in the dependency parameters of the multivariate bilateral gamma model are illustrated in Section 6. Physical measure pricing of product options is formulated and implemented in Section 7. Finally, Section 8 concludes.

2. Product Options and Functional Replication

Let $g(x,y)$ be a sufficiently smooth function of two variables. It may then be written as the integral of its second order cross partials as follows. For an arbitrary reference point $(x_0,y_0)$

$$g(x,y)=g(x_0,y)+g(x,y_0)-g(x_0,y_0)+{\int }_{x_0}^x{\int }_{y_0}^y\frac{{\partial }^{2}g(u,v)}{\partial x\partial y}dudv.$$

(3)

Now, we apply this result to the second cross partial itself to write

$$\begin{array}{lll}g(x,y)& =& g(x_0,y)+g(x,y_0)-g(x_0,y_0)\\ & & +{\int }_{x_0}^x{\int }_{y_0}^y\left[\begin{array}{c}\frac{{\partial }^{2}g(x_0,v)}{\partial x\partial y}+\frac{{\partial }^{2}g(u,y_0)}{\partial x\partial y}-\frac{{\partial }^{2}g(x_0,y_0)}{\partial x\partial y}\\ +{\int }_{x_0}^u{\int }_{y_0}^v\frac{{\partial }^{4}g(a,b)}{\partial x^2\partial y^2}dadb\end{array}\right]dudv\end{array}$$

(4)

$$\begin{array}{lll}& =& g(x_0,y)+g(x,y_0)-g(x_0,y_0)\\ & & +(x-x_0)\left[\frac{\partial g(x_0,y)}{\partial x}-\frac{\partial g(x_0,y_0)}{\partial x}\right]+(y-y_0)\left[\frac{\partial g(x,y_0)}{\partial y}-\frac{\partial g(x_0,y_0)}{\partial y}\right]\\ & & -\frac{{\partial }^{2}g(x_0,y_0)}{\partial x\partial y}(x-x_0)(y-y_0)\\ & & +{\int }_{x_0}^x{\int }_{y_0}^y{\int }_{x_0}^u{\int }_{y_0}^v\frac{{\partial }^{4}g(a,b)}{\partial x^2\partial y^2}dadbdudv\end{array}$$

(5)

The first row is a constant plus a function of x and a function of y that may be constructed using a bond and options on x and options on y. The second row is a product of x and options on y plus a product of y and options on x. The third row involves the product. The fourth order integral may be analyzed as follows. For $x>x_0$ and $y>y_0$, we write

$$\begin{array}{ccc}& & {\int }_{x_0}^x{\int }_{y_0}^y{\int }_{x_0}^u{\int }_{y_0}^v\frac{{\partial }^{4}g}{\partial x^2\partial y^2}\left(a,b\right)dadbdudv\hfill \end{array}$$

$$\begin{array}{ccc}& =& {\int }_{x_0}^x{\int }_{y_0}^{y}dadb{\int }_a^x{\int }_b^y\frac{{\partial }^{4}g}{\partial x^2\partial y^2}\left(a,b\right)dudv\hfill \end{array}$$

(6)

$$\begin{array}{ccc}& =& {\int }_{x_0}^x{\int }_{y_0}^{y}dadb\frac{{\partial }^{4}g}{\partial x^2\partial y^2}\left(a,b\right)\left(x-a\right)(y-b)\hfill \end{array}$$

(7)

$$\begin{array}{ccc}& =& {\int }_{x_0}^{\infty }{\int }_{y_0}^{\infty }dadb\frac{{\partial }^{4}g}{\partial x^2\partial y^2}\left(a,b\right){\left(x-a\right)}^{+}{(y-b)}^{+}\hfill \end{array}$$

(8)

In the positive orthant relative to the point $x_0,y_0$, one may replicate such functions with positions in product calls at strikes $(a,b).$ In the region $x<x_0$ and $y>y_0$, one may write

$$\begin{array}{ccc}& & {\int }_x^{x_0}{\int }_{y_0}^y{\int }_u^{x_0}{\int }_{y_0}^v\frac{{\partial }^{4}g}{\partial x^2\partial y^2}\left(a,b\right)dadbdudv\hfill \end{array}$$

$$\begin{array}{ccc}& =& {\int }_x^{x_0}{\int }_{y_0}^{y}dadb{\int }_x^a{\int }_b^y\frac{{\partial }^{4}g}{\partial x^2\partial y^2}\left(a,b\right)dudv\hfill \end{array}$$

(9)

$$\begin{array}{ccc}& =& {\int }_x^{x_0}{\int }_{y_0}^{y}dadb\frac{{\partial }^{4}g}{\partial x^2\partial y^2}\left(a,b\right)\left(a-x\right)(y-b)\hfill \end{array}$$

(10)

$$\begin{array}{ccc}& =& {\int }_0^{x_0}{\int }_{y_0}^{\infty }dadb\frac{{\partial }^{4}g}{\partial x^2\partial y^2}\left(a,b\right){\left(a-x\right)}^{+}{(y-b)}^{+}\hfill \end{array}$$

(11)

In the second orthant relative to $x_0,y_0$ products of puts on x and calls on y are employed. Similarly, in the third orthant, it is products of puts, and, in the fourth, it is calls on x and puts on $y.$ The general result may be stated as follows.

$$\begin{array}{ccc}\hfill g(x,y)& =& g(x_0,y)+g(x,y_0)-g(x_0,y_0)\hfill \\ & & +(x-x_0)\left[\frac{\partial g(x_0,y)}{\partial x}-\frac{\partial g(x_0,y_0)}{\partial x}\right]+(y-y_0)\left[\frac{\partial g(x,y_0)}{\partial y}-\frac{\partial g(x_0,y_0)}{\partial y}\right]\hfill \\ & & -\frac{{\partial }^{2}g(x_0,y_0)}{\partial x\partial y}(x-x_0)(y-y_0)\hfill \\ & & +{\int }_0^{x_0}{\int }_0^{y_0}dadb\frac{{\partial }^{4}g}{\partial x^2\partial y^2}(a,b){(a-x)}^{+}{(b-y)}^{+}\hfill \end{array}$$

$$\begin{array}{c}+{\int }_{{}_0}^{x_0}{\int }_{y_0}^{\infty }dadb\frac{{\partial }^{4}g}{\partial x^2\partial y^2}(a,b){(a-x)}^{+}{(y-b)}^{+}\hfill \end{array}$$

(12)

$$\begin{array}{c}+{\int }_{x_0}^{\infty }{\int }_0^{y_0}dadb\frac{{\partial }^{4}g}{\partial x^2\partial y^2}(a,b){(x-a)}^{+}{(b-y)}^{+}\hfill \end{array}$$

(13)

$$\begin{array}{c}+{\int }_{x_0}^{\infty }{\int }_{y_0}^{\infty }dadb\frac{{\partial }^{4}g}{\partial x^2\partial y^2}(a,b){(x-a)}^{+}{(y-b)}^{+}.\hfill \end{array}$$

(14)

A classic application of functional replication in one dimension is the pricing of the variance swap contract or the computation of the VIX index. Other applications include the computation of the CBOE skew index. In two dimensions, applications would include the replication of straddles and strangles on the spread of one asset price over another.

3. Multivariate Geometric Brownian Motion and Product Options

Suppose two assets, $S_1,S_2$, are driven by correlated Brownian motions with mean rates of return ${\mu }_1,{\mu }_2$ with volatilities ${\sigma }_1,{\sigma }_2$ and correlation $\rho$. For a pair of correlated standard Brownian motions $W_1,W_2$, the asset price dynamics under the true or physical measures are

$$\begin{array}{ccc}\hfill d{S}_1& =& {\mu }_1{S}_{1}dt+{\sigma }_1{S}_{1}d{W}_1\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill d{S}_2& =& {\mu }_2{S}_{2}dt+{\sigma }_2{S}_{2}d{W}_2\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill d{W}_{1}d{W}_2& =& \rho dt\hfill \end{array}$$

(17)

Consider a product call paying at maturity for strike $K_1,K_2$ the cash flow ${(S_1-K_1)}^{+}{(S_2-K_2)}^{+}$. Let the value of the option prior to maturity be $C(S_1,S_2,t)$ if, at time t, the stock prices are $S_1,S_2$. The differential of the call price is

$$\begin{array}{ccc}\hfill dC& =& C_{t}dt+C_{S_1}{\mu }_1{S}_{1}dt+C_{S_2}{\mu }_2{S}_{2}dt\hfill \end{array}$$

(18)

$$\begin{array}{ccc}& & +\frac{1}{2}{C}_{S_1{S}_1}{\sigma }_1^2{S}_1^{2}dt+\frac{1}{2}{C}_{S_2{S}_2}{\sigma }_2^2{S}_2^{2}dt+C_{S_1{S}_2}\rho {\sigma }_1{\sigma }_2{S}_1{S}_{2}dt\hfill \end{array}$$

(19)

$$\begin{array}{ccc}& & +C_{S_1}{\sigma }_1{S}_{1}d{W}_1+C_{S_2}{\sigma }_2{S}_{2}d{W}_2\hfill \end{array}$$

(20)

Consider a short call that is delta hedged with hedge returns

$$\begin{array}{ccc}& & C_S{\mu }_{1}Sdt+C_S{\sigma }_{1}Sd{W}_1\hfill \end{array}$$

(21)

$$\begin{array}{ccc}& & C_V{\mu }_{2}Vdt+C_V{\sigma }_{2}Vd{W}_2.\hfill \end{array}$$

(22)

Then, the hedged short call has a risk free return that must be the interest on value, or

$$-C_t-\frac{1}{2}{C}_{S_1{S}_1}{\sigma }_1^2{S}_1^2-\frac{1}{2}{C}_{S_2{S}_2}{\sigma }_2^2{S}_2{}^2-C_{S_1{S}_2}\rho {\sigma }_1{\sigma }_2{S}_1{S}_2=r\left(C_{S_1}{S}_2+C_{S_2}{S}_2-C\right).$$

(23)

Reversing time, the value function is a solution to the equation

$$C_T=r{S}_1{C}_{S_1}+r{S}_2{C}_{S_2}+\frac{1}{2}{C}_{S_1{S}_1}{\sigma }_1^2{S}_1^2+\frac{1}{2}{C}_{S_2{S}_2}{\sigma }_2^2{S}_2^2+C_{S_1{S}_2}\rho {\sigma }_1{\sigma }_2{S}_1{S}_2-rC$$

(24)

subject to the boundary condition

$$C(S,V,T)={(S_1-K_1)}^{+}{(S_2-K_2)}^{+}.$$

(25)

From the Feynmann–Kac relationship between partial differential equations and expectations, the value is given by

$$C(S_1,S_2,t)=e^{-rt}E\left[{\left(S_1{e}^{rt+{\sigma }_1\sqrt{t}{z}_1-{\sigma }_1^{2}t/2}-K_1\right)}^{+}{\left(S_2{e}^{rt+{\sigma }_2\sqrt{t}{z}_2-{\sigma }_2^{2}t/2}-K_2\right)}^{+}\right],$$

(26)

where $(z_1,z_2)$ has a standard bivariate normal distribution with correlation ρ. The product call option pricing formula may be developed on solving this integration problem. One may also assert Equation (26) directly by appealing to the existence of a single risk neutral measure; however, for completeness one should also display the explicit measure change. However, the explicit replication strategy is not as clear in such an approach.

The domain of integration is

$$\begin{array}{ccc}\hfill z_1& >& d=\frac{ln(K_1/S_1)}{{\sigma }_1\sqrt{t}}-\frac{r}{{\sigma }_1}\sqrt{t}+\frac{{\sigma }_1}{2}\sqrt{t}\hfill \end{array}$$

(27)

$$\begin{array}{ccc}\hfill z_2& >& e=\frac{ln(K_2/S_2)}{{\sigma }_2\sqrt{t}}-\frac{r}{{\sigma }_2}\sqrt{t}+\frac{{\sigma }_2}{2}\sqrt{t}\hfill \end{array}$$

(28)

The forward product call price $e^{rt}C(S_1,S_2,t)$ is given by the double integral

$${\int }_d^{\infty }{\int }_e^{\infty }\left(S_1{e}^{rt+{\sigma }_1\sqrt{t}{z}_1-{\sigma }_1^{2}t/2}-K_1\right)\left(S_2{e}^{rt+{\sigma }_2\sqrt{t}{z}_2-{\sigma }_2^{2}t/2}-K_2\right)f(z_1,z_2)d{z}_{1}d{z}_2$$

(29)

There are four terms and

$$e^{rt}C(S_1,S_2,t)=I_1-I_2-I_3+I_4$$

(30)

where

$$\begin{array}{ccc}\hfill I_1& =& S_1{S}_2{e}^{2rt}{\int }_d^{\infty }{\int }_e^{\infty }{e}^{{\sigma }_1\sqrt{t}{z}_1-{\sigma }_1^{2}t/2}{e}^{{\sigma }_2\sqrt{t}{z}_2-{\sigma }_2^{2}t/2}f(z_1,z_2)d{z}_{1}d{z}_2\hfill \end{array}$$

(31)

$$\begin{array}{ccc}\hfill I_2& =& K_1{S}_2{e}^{rt}{\int }_d^{\infty }{\int }_e^{\infty }{e}^{{\sigma }_2\sqrt{t}{z}_2-{\sigma }_2^{2}t/2}f(z_1,z_2)d{z}_{1}d{z}_2\hfill \end{array}$$

(32)

$$\begin{array}{ccc}\hfill I_3& =& K_2{S}_1{e}^{rt}{\int }_d^{\infty }{\int }_e^{\infty }{e}^{{\sigma }_1\sqrt{t}{z}_1-{\sigma }_1^{2}t/2}f(z_1,z_2)d{z}_{1}d{z}_2\hfill \end{array}$$

(33)

$$\begin{array}{ccc}\hfill I_4& =& K_1{K}_2{\int }_d^{\infty }{\int }_e^{\infty }f(z_1,z_2)d{z}_{1}d{z}_2\hfill \end{array}$$

(34)

Here, f is the bivariate normal density.

We then have that

$$\begin{array}{ccc}\hfill I_4& =& 1-\left(N\left(d\right)+N\left(e\right)-bvncdf(d,e,\rho )\right)\hfill \\ & =& {\int }_d^{\infty }\frac{1}{\sqrt{2\pi }}{e}^{-w^2/2}N\left(\frac{\rho w-e}{\sqrt{1-{\rho }^2}}\right)dw\hfill \end{array}$$

(35)

where $bvncdf$ is the bivariate normal cummulative distribution function.

For $I_3$, we write the joint density as

$$f(z_1,z_2)=\frac{1}{\sqrt{2\pi }}{e}^{-z_1^2/2}\frac{1}{\sqrt{2\pi (1-{\rho }^2)}}{e}^{-\frac{{(z_2-\rho z_1)}^2}{2(1-{\rho }^2)}}$$

(36)

It follows that

$$\begin{array}{ccc}\hfill I_3& =& K_2{S}_1{e}^{rt}{\int }_d^{\infty }{e}^{{\sigma }_1\sqrt{t}{z}_1-{\sigma }_1^{2}t/2}\frac{1}{\sqrt{2\pi }}{e}^{-z_1^2/2}N\left(\frac{\rho z_1-e}{\sqrt{1-{\rho }^2}}\right)\hfill \end{array}$$

(37)

$$\begin{array}{ccc}& =& K_2{S}_1{e}^{rt}{\int }_{d-{\sigma }_1\sqrt{t}}^{\infty }\frac{1}{\sqrt{2\pi }}{e}^{-w^2/2}N\left(\frac{\rho (w+{\sigma }_1\sqrt{t})-e}{\sqrt{1-{\rho }^2}}\right)dw\hfill \end{array}$$

(38)

Similarly,

$$I_2=K_1{S}_2{e}^{rt}{\int }_{e-{\sigma }_2\sqrt{t}}^{\infty }\frac{1}{\sqrt{2\pi }}{e}^{-w^2/2}N\left(\frac{\rho (w+{\sigma }_2\sqrt{t})-d}{\sqrt{1-{\rho }^2}}\right)dw$$

(39)

Finally,

$$\begin{array}{ccc}\hfill I_1& =& S_1{S}_2{e}^{2rt}{\int }_d^{\infty }{\int }_e^{\infty }{e}^{{\sigma }_1\sqrt{t}{z}_1-{\sigma }_1^{2}t/2}{e}^{{\sigma }_2\sqrt{t}{z}_2-{\sigma }_2^{2}t/2}\times \hfill \end{array}$$

(40)

$$\begin{array}{ccc}& & \frac{1}{\sqrt{2\pi }}{e}^{-z_1^2/2}\frac{1}{\sqrt{2\pi (1-{\rho }^2)}}{e}^{-\frac{{(z_2-\rho z_1)}^2}{2(1-{\rho }^2)}}\hfill \end{array}$$

(41)

$$\begin{array}{ccc}& =& S_1{S}_2{e}^{2rt}{\int }_d^{\infty }{e}^{{\sigma }_1\sqrt{t}{z}_1-{\sigma }_1^{2}t/2}\frac{1}{\sqrt{2\pi }}{e}^{-z_1^2/2}d{z}_1\times \hfill \end{array}$$

(42)

$$\begin{array}{ccc}& & {\int }_e^{\infty }{e}^{{\sigma }_2\sqrt{t}{z}_2-{\sigma }_2^{2}t/2}\frac{1}{\sqrt{2\pi (1-{\rho }^2)}}{e}^{-\frac{{(z_2-\rho z_1)}^2}{2(1-{\rho }^2)}}d{z}_2\hfill \end{array}$$

(43)

$$\begin{array}{ccc}& =& S_1{S}_2{e}^{2rt}{\int }_d^{\infty }{e}^{{\sigma }_1\sqrt{t}{z}_1-{\sigma }_1^{2}t/2}\frac{1}{\sqrt{2\pi }}{e}^{-z_1^2/2}d{z}_1\times \hfill \end{array}$$

(44)

$$\begin{array}{ccc}& & {\int }_{\frac{e-\rho z_1}{\sqrt{1-{\rho }^2}}}\frac{1}{\sqrt{2\pi }}{e}^{-w^2/2}{e}^{{\sigma }_2\sqrt{t}(\rho z_1+\sqrt{1-{\rho }^2}w)-{\sigma }_2^{2}t/2}dw\hfill \end{array}$$

(45)

$$\begin{array}{ccc}& =& S_1{S}_2{e}^{2rt}{\int }_d^{\infty }{e}^{{\sigma }_1\sqrt{t}{z}_1-{\sigma }_1^{2}t/2}\frac{1}{\sqrt{2\pi }}{e}^{-z_1^2/2}{e}^{{\sigma }_2\sqrt{t}\rho z_1}\times \hfill \end{array}$$

(46)

$$\begin{array}{ccc}& & {\int }_{\frac{e-\rho z_1}{\sqrt{1-{\rho }^2}}}\frac{1}{\sqrt{2\pi }}{e}^{-w^2/2+{\sigma }_2\sqrt{t}\sqrt{1-{\rho }^2}w--{\sigma }_2^{2}t/2}dw\hfill \end{array}$$

(47)

$$\begin{array}{ccc}& =& S_1{S}_2{e}^{2rt}{\int }_d^{\infty }{e}^{{\sigma }_1\sqrt{t}{z}_1-{\sigma }_1^{2}t/2}\frac{1}{\sqrt{2\pi }}{e}^{-z_1^2/2}{e}^{{\sigma }_2\sqrt{t}\rho z_1}\times \hfill \end{array}$$

(48)

$$\begin{array}{ccc}& & {\int }_{\frac{e-\rho z_1}{\sqrt{1-{\rho }^2}}}\frac{1}{\sqrt{2\pi }}{e}^{-{(w-{\sigma }_2\sqrt{t}\sqrt{1-{\rho }^2})}^2/2-{\rho }^2{\sigma }_2^{2}t/2}dw\hfill \end{array}$$

(49)

$$\begin{array}{ccc}& =& S_1{S}_2{e}^{2rt}{\int }_d^{\infty }{e}^{{\sigma }_1\sqrt{t}{z}_1-{\sigma }_1^{2}t/2}\times \hfill \end{array}$$

(50)

$$\begin{array}{ccc}& & \frac{1}{\sqrt{2\pi }}{e}^{-z_1^2/2}{e}^{{\sigma }_2\sqrt{t}\rho z_1-{\rho }^2{\sigma }_2^{2}t/2}N\left(\frac{\rho z_1+(1-{\rho }^2){\sigma }_2\sqrt{t}-e}{\sqrt{1-{\rho }^2}}\right)d{z}_1\hfill \end{array}$$

(51)

$$\begin{array}{ccc}& =& S_1{S}_2{e}^{2rt+\rho {\sigma }_1{\sigma }_{2}t}{\int }_{d-{\sigma }_1\sqrt{t}-{\sigma }_2\sqrt{t}\rho }^{\infty }\frac{1}{\sqrt{2\pi }}{e}^{-w^2/2}\times \hfill \end{array}$$

(52)

$$\begin{array}{ccc}& & N\left(\frac{\rho w+\rho {\sigma }_1\sqrt{t}+{\sigma }_2\sqrt{t}-e}{\sqrt{1-{\rho }^2}}\right)\hfill \end{array}$$

(53)

We, thus, obtain

$$\begin{array}{ccc}\hfill I_1& =& S_1{S}_2{e}^{2rt+\rho {\sigma }_1{\sigma }_{2}t}(1-\left(N\left(d_1\right)+N\left(e_1\right)-bvncdf(d_1,e_1,\rho )\right))\hfill \end{array}$$

(54)

$$\begin{array}{ccc}\hfill I_2& =& K_1{S}_2{e}^{rt}(1-\left(N\left(d_2\right)+N\left(e_2\right)-bvncdf(d_2,e_2,\rho \right)))\hfill \end{array}$$

(55)

$$\begin{array}{ccc}\hfill I_3& =& K_2{S}_1{e}^{rt}\left(1-\left(N\left(d_3\right)+N\left(e_3\right)-bvncdf(d_3,e_3,\rho \right))\right)\hfill \end{array}$$

(56)

$$\begin{array}{ccc}\hfill I_4& =& K_1{K}_2(1-(N\left(d\right)+N\left(e\right)-bvncdf(d,e,\rho ))\hfill \end{array}$$

(57)

$$\begin{array}{ccc}\hfill d_1& =& d-{\sigma }_1\sqrt{t}-\rho {\sigma }_2\sqrt{t}\hfill \end{array}$$

(58)

$$\begin{array}{ccc}\hfill e_1& =& e-{\sigma }_2\sqrt{t}-\rho {\sigma }_1\sqrt{t}\hfill \end{array}$$

(59)

$$\begin{array}{ccc}\hfill d_2& =& d-\rho {\sigma }_2\sqrt{t}\hfill \end{array}$$

(60)

$$\begin{array}{ccc}\hfill e_2& =& e-{\sigma }_2\sqrt{t}\hfill \end{array}$$

(61)

$$\begin{array}{ccc}\hfill d_3& =& d-{\sigma }_1\sqrt{t}\hfill \end{array}$$

(62)

$$\begin{array}{ccc}\hfill e_3& =& e-\rho {\sigma }_1\sqrt{t}\hfill \end{array}$$

(63)

$$\begin{array}{ccc}\hfill d& =& \frac{ln(K_1/S_1)}{{\sigma }_1\sqrt{t}}-\frac{r}{{\sigma }_1}\sqrt{t}+\frac{{\sigma }_1}{2}\sqrt{t}\hfill \end{array}$$

(64)

$$\begin{array}{ccc}\hfill e& =& \frac{ln(K_2/S_2)}{{\sigma }_2\sqrt{t}}-\frac{r}{{\sigma }_2}\sqrt{t}+\frac{{\sigma }_2}{2}\sqrt{t}\hfill \end{array}$$

(65)

The option price is

$$\begin{array}{ccc}\hfill C(S_1,S_2,t)& =& S_1{S}_2{e}^{rt+\rho {\sigma }_1{\sigma }_{2}t}\left(1-\left(N\left(d_1\right)+N\left(e_1\right)-bvncdf(d_1,e_1)\right)\right)\hfill \\ & & -K_1{S}_2\left(1-\left(N\left(d_2\right)+N\left(e_2\right)-bvncdf(d_2,e_2)\right)\right)\hfill \\ & & -K_2{S}_1\left(1-\left(N\left(d_3\right)+N\left(e_3\right)-bvncdf(d_3,e_3)\right)\right)\hfill \\ & & +K_1{K}_2{e}^{-rt}(1-\left(N\left(d\right)+N\left(e\right)-bvncdf(d,e)\right))\hfill \\ \hfill d& =& \frac{ln(K_1/S_1)}{{\sigma }_1\sqrt{t}}-\frac{r}{{\sigma }_1}\sqrt{t}+\frac{{\sigma }_1}{2}\sqrt{t}\hfill \\ \hfill e& =& \frac{ln(K_2/S_2)}{{\sigma }_2\sqrt{t}}-\frac{r}{{\sigma }_2}\sqrt{t}+\frac{{\sigma }_2}{2}\sqrt{t}\hfill \\ \hfill d_1& =& d-{\sigma }_1\sqrt{t}-\rho {\sigma }_2\sqrt{t}\hfill \\ \hfill e_1& =& e-{\sigma }_2\sqrt{t}-\rho {\sigma }_1\sqrt{t}\hfill \\ \hfill d_2& =& d-\rho {\sigma }_2\sqrt{t}\hfill \\ \hfill e_2& =& e-{\sigma }_2\sqrt{t}\hfill \\ \hfill d_3& =& d-{\sigma }_1\sqrt{t}\hfill \\ \hfill e_3& =& e-\rho {\sigma }_1\sqrt{t}\hfill \end{array}$$

(66)

Consider now the product of a put times a call with value

$$e^{-rt}E\left[{\left(K_1-S_1{e}^{rt+{\sigma }_1\sqrt{t}{z}_1-{\sigma }_1^{2}t/2}\right)}^{+}{\left(S_2{e}^{rt+{\sigma }_2\sqrt{t}{z}_2-{\sigma }_2^{2}t/2}-K_2\right)}^{+}\right]$$

(67)

The domain of integration is

$$\begin{array}{ccc}\hfill z_1& <& d=\frac{ln(K/S)}{{\sigma }_1\sqrt{t}}-\frac{r}{{\sigma }_1}\sqrt{t}+\frac{{\sigma }_1}{2}\sqrt{t}\hfill \end{array}$$

(68)

$$\begin{array}{ccc}\hfill z_2& >& e\hfill \end{array}$$

(69)

The result is of the form

$$I_2-I_1+I_3-I_4$$

(70)

and the integral $I_2$ is

$$\begin{array}{ccc}& & K_1{S}_2{e}^{rt}{\int }_{-\infty }^d{\int }_e^{\infty }{e}^{{\sigma }_2\sqrt{t}{z}_2-{\sigma }_2^{2}t/2}f(z_1,z_2)d{z}_{1}d{z}_2\hfill \\ & =& {\int }_{e-{\sigma }_2\sqrt{t}}^{\infty }\frac{1}{\sqrt{2\pi }}{e}^{-w^2/2}N\left(\frac{d-\rho {\sigma }_2\sqrt{t}-\rho w}{\sqrt{1-{\rho }^2}}\right)dw.\hfill \end{array}$$

There are four cases denoted $cc,cp,pc,$ and $pp$ for whether the option is a product of two calls, a call and a put, a put and a call, or two puts, respectively. For two calls and two puts, the forward price is $I_1-I_2-I_3+I_4$, and, for cases $cp,pc$, the forward price is given by $I_2-I_1+I_3-I_4.$ In each of the integrals, the probability is a bivariate normal probability of the appropriate quadrant defined by the point $(d,e)$.

4. Fast Fourier Transform for Product Options

For a joint density $f(x,y)$ of a pair of log returns, define the product call price by

$$w(a,b)={\int }_a^{\infty }{\int }_b^{\infty }\left(e^x-e^a\right)\left(e^y-e^b\right)f(x,y)dxdy$$

(71)

for a joint density $f(x,y).$

Further suppose the joint characteristic function $\varphi (u,v)$ of the joint density is available in closed form. By definition,

$$\varphi (u,v)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{e}^{iux+ivy}f(x,y)dxdy.$$

(72)

Following the methods of Carr and Madan (1999), define

$$\begin{array}{ccc}\hfill \gamma (a,b)& =& {\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{e}^{\alpha a+\beta b+iua+ivb}w(a,b)dadb\hfill \end{array}$$

(73)

$$\begin{array}{ccc}& =& {\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{e}^{i(u-i\alpha )a+i(v-i\beta )b}\times \hfill \end{array}$$

(74)

$$\begin{array}{ccc}& & {\int }_a^{\infty }{\int }_b^{\infty }\left(e^x-e^a\right)\left(e^y-e^b\right)f(x,y)dxdydadb\hfill \end{array}$$

(75)

$$\begin{array}{ccc}& =& {\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)dxdy\times \hfill \end{array}$$

(76)

$$\begin{array}{ccc}& & {\int }_{-\infty }^x{\int }_{-\infty }^y{e}^{i(u-i\alpha )a+i(v-i\beta )}\left(e^x-e^a\right)\left(e^y-e^b\right)dadb\hfill \end{array}$$

(77)

$$\begin{array}{ccc}& =& {\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)dxdy\times {\int }_{-\infty }^x{\int }_{-\infty }^{y}dadb\hfill \end{array}$$

(78)

$$\begin{array}{ccc}& & \left(\begin{array}{c}{e}^{i(u-i\alpha )a+i(v-i\beta )+x+y}-e^{i(u-i(\alpha +1\left)\right)a+i(v-i\beta )b+y}-\\ e^{i(u-i\alpha )a+i(v-i(\beta +1\left)\right)b+x}+e^{i(u-i(\alpha +1\left)\right)a+i(v-i(\beta +1\left)\right)b}\end{array}\right)\hfill \end{array}$$

(79)

$$\begin{array}{ccc}& =& \varphi \left(u-i(\alpha +1),v-i(\beta +1)\right)\times \hfill \end{array}$$

(80)

$$\begin{array}{ccc}& & \left(\begin{array}{c}\frac{-1}{(u-i\alpha )(v-i\beta )}+\frac{1}{(u-i(\alpha +1\left)\right)(v-i\beta )}+\\ \frac{1}{(u-i\alpha )(v-i(\beta +1\left)\right)}-\frac{1}{(u-i(\alpha +1\left)\right)(v-i(\beta +1\left)\right)}\end{array}\right)\hfill \end{array}$$

(81)

$$\begin{array}{ccc}& =& \frac{\varphi (u-i(\alpha +1),v-i(\beta +1\left)\right)}{\left({\left(u-i\alpha \right)}^2-i(u-i\alpha )\right)\left({\left(v-i\beta \right)}^2-i(v-i\beta )\right)}\hfill \end{array}$$

(82)

It follows that

$$w(a,b)=\frac{1}{4{\pi }^2}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{e}^{-iua-ivb}\frac{e^{-\alpha a-\beta b}\varphi (u-i(\alpha +1),v-i(\beta +1))}{\left({\left(u-i\alpha \right)}^2-i(u-i\alpha )\right)\left({\left(v-i\beta \right)}^2-i(v-i\beta )\right)}dudv$$

(83)

Hence, by two-dimensional fast Fourier inversion, we may price product options given the joint characteristic function.

For the product of put and a call defined for $\alpha >1,\beta >0$

$$\begin{array}{ccc}\hfill {\gamma }_{pc}(a,b)& =& {\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{e}^{i(u+i\alpha )a}{e}^{i(v-i\beta )b}{w}_{pc}(a,b)dadb\hfill \end{array}$$

(84)

$$\begin{array}{ccc}& =& {\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\int }_{-\infty }^a{\int }_b^{\infty }{e}^{i(u+i\alpha )a}{e}^{i(v-i\beta )}\times \hfill \end{array}$$

(85)

$$\begin{array}{ccc}& & \left(e^a-e^x\right)\left(e^y-e^b\right)f(x,y)dxdydadb\hfill \end{array}$$

(86)

$$\begin{array}{ccc}& =& {\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)dxdy{\int }_x^{\infty }{\int }_{-\infty }^y\times \hfill \end{array}$$

(87)

$$\begin{array}{ccc}& & e^{i(u+i\alpha )a}{e}^{i(v-i\beta )b}\left(e^a-e^x\right)\left(e^y-e^b\right)dadb\hfill \end{array}$$

(88)

$$\begin{array}{ccc}& =& \varphi \left(u+i(\alpha -1),v-i(\beta +1)\right)\times \hfill \end{array}$$

(89)

$$\begin{array}{ccc}& & \left[\begin{array}{c}\frac{1}{(u+i(\alpha -1\left)\right)(v-i\beta )}-\frac{1}{(u+i(\alpha -1\left)\right)(v-i(\beta +1\left)\right)}\\ -\frac{1}{(u+i\alpha )(v-i\beta )}+\frac{1}{(u+i\alpha )(v-i(\beta +1\left)\right)}\end{array}\right]\hfill \end{array}$$

(90)

$$\begin{array}{ccc}& =& \frac{\varphi \left(u+i(\alpha -1),v-i(\beta +1)\right)}{({(u+i\alpha )}^2-i(u+i\alpha ))({(v-i\beta )}^2-i(v-i\beta ))}\hfill \end{array}$$

(91)

Similarly, we will have

$$\begin{array}{ccc}\hfill {\gamma }_{cp}(a,b)& =& \frac{\varphi (u-i(\alpha +1),v+i(\beta -1\left)\right)}{({(u-i\alpha )}^2-i(u-i\alpha ))({(u+i\beta )}^2-i(u+i\beta ))}\hfill \end{array}$$

(92)

$$\begin{array}{ccc}\hfill {\gamma }_{pp}(a,b)& =& \frac{\varphi (u+i(\alpha -1),v+i(\beta -1\left)\right)}{({(u+i\alpha )}^2-i(u+i\alpha ))({(u+i\beta )}^2-i(u+i\beta ))}.\hfill \end{array}$$

(93)

It follows that

$$\begin{array}{ccc}\hfill w_{cc}(a,b)& =& \frac{1}{4{\pi }^2}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{e}^{-iua-ivb}\times \hfill \end{array}$$

(94)

$$\begin{array}{ccc}& & \frac{e^{-\alpha a-\beta b}\varphi (u-i(\alpha +1),v-i(\beta +1))}{\left({\left(u-i\alpha \right)}^2-i(u-i\alpha )\right)\left({\left(v-i\beta \right)}^2-i(v-i\beta )\right)}dadb\hfill \end{array}$$

(95)

$$\begin{array}{ccc}\hfill w_{pc}(a,b)& =& \frac{1}{4{\pi }^2}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{e}^{-iua-ivb}\times \hfill \end{array}$$

(96)

$$\begin{array}{ccc}& & \frac{e^{\alpha a-\beta b}\varphi \left(u+i(\alpha -1),v-i(\beta +1)\right)}{({(u+i\alpha )}^2-i(u+i\alpha ))({(v-i\beta )}^2-i(v-i\beta ))}dadb\hfill \end{array}$$

(97)

$$\begin{array}{ccc}\hfill w_{cp}(a,b)& =& \frac{1}{4{\pi }^2}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{e}^{-iua-ivb}\times \hfill \end{array}$$

(98)

$$\begin{array}{ccc}& & \frac{e^{-\alpha a+\beta b}\varphi (u-i(\alpha +1),v+i(\beta -1))}{({(u-i\alpha )}^2-i(u-i\alpha ))({(u+i\beta )}^2-i(u+i\beta ))}dadb\hfill \end{array}$$

(99)

$$\begin{array}{ccc}\hfill w_{pp}(a,b)& =& \frac{1}{4{\pi }^2}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{e}^{-iua-ivb}\times \hfill \end{array}$$

(100)

$$\begin{array}{ccc}& & \frac{e^{\alpha a+\beta b}\varphi (u+i(\alpha -1),v+i(\beta -1))}{({(u+i\alpha )}^2-i(u+i\alpha ))({(u+i\beta )}^2-i(u+i\beta ))}dadb\hfill \end{array}$$

(101)

For puts, the values of $\alpha ,\beta$ need to be above unity.

The inversion may be accomplished using a two dimensional fast Fourier Transform. The fast Fourier transform Y of an $m\times n$ matrix X evaluates the double sum

$$Y_{p+1,q+1}=\sum _{j=0}^{m-1}\sum _{k=0}^{n-1}{\omega }_m^{jp}{\omega }_n^{kq}{X}_{j+1,k+1},p=0,\cdots m-1,q=0,\cdots n-1$$

(102)

where

$$\begin{array}{ccc}\hfill {\omega }_m& =& e^{-2\pi i/m}\hfill \end{array}$$

(103)

$$\begin{array}{ccc}\hfill {\omega }_n& =& e^{-2\pi i/n}\hfill \end{array}$$

(104)

The integral to be evaluated is

$${\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{e}^{-iua}{e}^{-ivb}g(u,v)dudv$$

(105)

We may approximate by

$${\int }_{-M}^M{\int }_{-N}^N{e}^{-iua-ivb}g(u,v)dudv$$

(106)

We discretize with spacing ${\lambda }_u,{\lambda }_v$

$$\begin{array}{ccc}\hfill u_j& =& -M+{\lambda }_{u}j\hfill \end{array}$$

(107)

$$\begin{array}{ccc}\hfill v_k& =& -N+{\lambda }_{v}k\hfill \end{array}$$

(108)

We also discretize in strike space with spacing ${\eta }_a,{\eta }_b$ by

$$\begin{array}{ccc}\hfill a_p& =& -A+{\eta }_{a}p\hfill \end{array}$$

(109)

$$\begin{array}{ccc}\hfill b_q& =& -B+{\eta }_{b}q\hfill \end{array}$$

(110)

to approximate the double integral by the double sum

$$\begin{array}{ccc}\hfill W(a_p,b_q)& =& \sum _{j=0}^{m-1}\sum _{k=0}^{n-1}{e}^{-i(-M+{\lambda }_{u}j)(-A+{\eta }_{a}p)}{e}^{-i(-N+{\lambda }_{v}k)(-B+{\eta }_{b}q)}\hfill \end{array}$$

(111)

$$\begin{array}{ccc}& & \times g(-M+{\lambda }_{u}j,-N+{\lambda }_{v}k){\lambda }_u{\lambda }_v\hfill \end{array}$$

(112)

$$\begin{array}{ccc}& & e^{-iM(-A+{\eta }_{a}p)}W(a_p,b_q)e^{-iN((-B+{\eta }_{b}q)}\hfill \end{array}$$

(113)

$$\begin{array}{ccc}& =& \sum _{j=0}^{m-1}\sum _{k=0}^{n-1}{e}^{-i{\lambda }_u{\eta }_{a}jp}{e}^{-i{\lambda }_v{\eta }_{b}kq}\times \hfill \end{array}$$

(114)

$$\begin{array}{ccc}& & \left(\begin{array}{c}{e}^{iA(-M+{\lambda }_{u}j)}\\ \times g(-M+{\lambda }_{u}j,-N+{\lambda }_{v}k)\\ \times e^{iB(-N+{\lambda }_{v}k)}\end{array}\right){\lambda }_u{\lambda }_v\hfill \end{array}$$

(115)

Now, we set

$$\begin{array}{ccc}\hfill {\lambda }_u{\eta }_a& =& \frac{2\pi }{m}\hfill \end{array}$$

(116)

$$\begin{array}{ccc}\hfill {\lambda }_v{\eta }_b& =& \frac{2\pi }{n}\hfill \end{array}$$

(117)

Define

$$X_{j+1,k+1}=e^{iA(-M+{\lambda }_{u}j)}g(-M+{\lambda }_{u}j,-N+{\lambda }_{v}k)e^{iB(-N+{\lambda }_{v}k)}{\lambda }_u{\lambda }_v$$

(118)

Then, define Y via Equation (102). The option prices are obtained as

$$W\left(a_p,b_q\right)=e^{iM(-A+{\eta }_{a}p)}{Y}_{p+1,q+1}{e}^{iN((-B+{\eta }_{b}q)}.$$

(119)

5. Multivariate Bilateral Gamma

An application of the two dimensional fast Fourier transform on pricing product options requires a closed form multivariate characteristic function. The marginals of this joint characteristic function should be capable of calibrating risk neutral densities observed in the separate option markets for a given maturity. There are a number of pure jump Lévy processes that will fit marginal option prices.

Among them are the three parameter variance gamma model of Madan and Seneta (1990), Madan et al. (1998) and its four parameter extension to the bilateral gamma model proposed in Küchler and Tappe (2008). Recently, Madan and Schoutens (2020), Madan (2020), and Madan and Wang (2020a) have investigated multivariate models with closed form characteristic functions consistent with arbitrary marginal bilateral gamma models. Here, we present a summary of the required characteristic functions.

Consider first the bilateral gamma process for the marginals. Let $({\gamma }_p\left(t\right),t>0)$, $({\gamma }_n\left(t\right),t>0)$ be two independent standard gamma processes with unit mean and variance rates. Now, we introduce four parameters $b_p,c_p,b_n,c_n$ representing the scale and speed of positive and negative moves, respectively. The bilateral gamma process $X_{BG}\left(t\right)$ is given by

$$X_{BG}\left(t\right)=b_p{\gamma }_p\left(c_{p}t\right)-b_n{\gamma }_n\left(c_{n}t\right).$$

(120)

The process is a Lévy process of independent and identically distributed increments. It is also a pure jump process with the characteristic function

$$E\left[expiu{X}_{BG}\left(t\right)\right]={\left(\frac{1}{1-iu{b}_p}\right)}^{c_{p}t}{\left(\frac{1}{1+iu{b}_n}\right)}^{c_{n}t}.$$

(121)

The characteristic function has a Lévy–Khintchine decomposition in terms of a Lévy density $k\left(x\right)$ that announces the arrival rate of jumps of size x as follows.

$$E\left[expiu{X}_{BG}\left(t\right)\right]=exp\left({\int }_{-\infty }^{\infty }\left(e^{iux}-1\right)k\left(x\right)dx\right).$$

(122)

The specific form for the Lévy density is

$$k\left(x\right)=c_p\frac{exp\left(-\frac{x}{b_p}\right)}{\left|x\right|}{\mathbf{1}}_{x>0}+c_n\frac{exp\left(-\frac{\left|x\right|}{b_n}\right)}{\left|x\right|}{\mathbf{1}}_{x<0}.$$

(123)

The density for any horizon may be obtained by Fourier inversion of the characteristic function using the Fast Fourier Transform. The density also has a closed form in terms of the Whittaker W function (Küchler and Tappe (2008)). We denote this function by $W(a,b)$. For $x>0$, the density is given by

$$\begin{array}{ccc}\hfill f_{BG}\left(x\right)& =& {\left(\frac{1}{b_p}\right)}^{c_p}{\left(\frac{1}{b_n}\right)}^{c_n}\left(\frac{1}{{\left(1/b_p+1/b_n\right)}^{(c_p+c_n)/2}\Gamma \left(c_p\right)}\right)\times \hfill \\ & & x^{\left(c_p+c_n\right)/2-1}exp\left(-x/2(1/b_p+1/b_n)\right)\times \hfill \\ & & W\left(\frac{c_p-c_n}{2},\frac{c_p+c_n-1}{2},\left(\frac{1}{b_p}+\frac{1}{b_n}\right)x\right).\hfill \end{array}$$

(124)

For $x<0$, the roles of $b_p,c_p$ and $b_n,c_n$ are reversed.

The multivariate bilateral gamma model is presented in Madan and Schoutens (2020). The presentation follows Buchmann et al. (2019) where a multivariate Lévy process is constructed having a Lévy density with full support on ${\mathbb{R}}^n-\left\{0\right\}$ that is also consistent with prespecified variance gamma marginals displaying different kurtosis levels for each component. The construction in Madan and Schoutens (2020) extends that of Buchmann et al. (2019) to attain consistency with prespecified bilateral gamma marginals.

Let the marginal distributions be bilateral gamma with the parameters $b_{pi},c_{pi},b_{ni},$ and $c_{ni}$ for component i. The multivariate bilateral gamma model is the sum of a multivariate variance gamma plus independent bilateral gamma shocks. The multivariate variance gamma is a multivariate normal with the mean vector θ and covariance Σ, both of which are scaled by a single gamma variate g with a unit mean and variance ν. The multivariate variance gamma random vector $X_{MVG}$ may then be written as

$$X_{MVG}=\theta g+\sqrt{g}\sqrt{diag(\Sigma )}Z$$

(125)

where the vector $Z=(Z_i,i=1,\cdots ,n)$ is multivariate normal with zero means, unit variances, and correlation matrix C. The covariance matrix is

$$\Sigma =\sqrt{diag\left(\Sigma \right)}C\sqrt{diag(\Sigma )}.$$

(126)

In addition, there are independent bilateral gamma variates $Y_i$ with the parameters ${\tilde{b}}_{pi},{\tilde{c}}_{pi},{\tilde{b}}_{ni},$ and ${\tilde{c}}_{ni}$ and with $Y={(Y_1,\cdots ,Y_n)}^{\prime }$

$$X_{MBG}=X_{MVG}+Y.$$

(127)

The only dependency parameters are the correlation matrix C and the variance ν of the gamma variate g. All other parameters are derived from the parameters of the marginal processes. Specifically, we have that

$$\begin{array}{ccc}\hfill {\theta }_i& =& \frac{b_{pi}-b_{ni}}{\nu }\hfill \end{array}$$

(128)

$$\begin{array}{ccc}\hfill {\Sigma }_{ii}& =& \frac{2b_{pi}{b}_{ni}}{\nu }\hfill \end{array}$$

(129)

$$\begin{array}{ccc}\hfill {\tilde{b}}_{pi}& =& b_{pi}\hfill \end{array}$$

(130)

$$\begin{array}{ccc}\hfill {\tilde{b}}_{ni}& =& b_{ni}\hfill \end{array}$$

(131)

$$\begin{array}{ccc}\hfill {\tilde{c}}_{pi}& =& c_{pi}-\frac{1}{\nu }\hfill \end{array}$$

(132)

$$\begin{array}{ccc}\hfill {\tilde{c}}_{ni}& =& c_{ni}-\frac{1}{\nu }\hfill \end{array}$$

(133)

The parameter ν has a lower bound of the reciprocal of the minimum of all the marginal speed parameters ${min}_i(min(c_{pi},c_{ni}))$.

The characteristic function of the multivariate bilateral gamma vector is given by

$$\begin{array}{ccc}\hfill {\varphi }_{MBG}\left(u\right)& =& E\left[exp\left(i{u}^{\prime }{X}_{MBG}\right)\right]\hfill \\ & =& {\left(\frac{1}{1-i{u}^{\prime }\theta \nu +\frac{\nu }{2}{u}^{\prime }\Sigma u}\right)}^{\frac{1}{\nu }}\prod _j{\left(\frac{1}{1-i{u}_j{b}_{pj}}\right)}^{c_{pj}-1/\nu }{\left(\frac{1}{1+iu{b}_{nj}}\right)}^{c_{nj}-1/\nu }.\hfill \end{array}$$

(134)

The multivariate arrival rates or Lévy density $k\left(x\right)$ for the multivariate bilateral gamma model may be specified as follows.

$$k\left(x\right)=\tilde{m}\left(x\right)+\sum _{j=1}^n{k}_j\left(x_j\right)\prod _{\begin{array}{c}l\ne j\\ l=1\end{array}}^M{\mathbf{1}}_{x_l=0},$$

(135)

where

$$\begin{array}{ccc}\hfill \tilde{m}\left(z\right)& =& \frac{exp\left({\theta }^T{\Sigma }^{-1}x\right)}{\nu {\left(2\pi \right)}^{n/2-1}\sqrt{|\Sigma |}\sqrt{x^T{\Sigma }^{-1}x}}\times \hfill \\ & & exp\left(-\sqrt{\left({\theta }^T{\Sigma }^{-1}\theta +\frac{2}{\nu }\right)\left(x^T{\Sigma }^{-1}x\right)}\right)\hfill \end{array}$$

(136)

$$\begin{array}{ccc}\hfill k_j\left(x_j\right)& =& \frac{c_{nj-1/\nu }}{|x_j|}exp\left(-|x_j|/b_{nj}\right){\mathbf{1}}_{x_j<0}\hfill \end{array}$$

(137)

$$\begin{array}{ccc}& & +\frac{c_{pj}-1/\nu }{x_j}exp\left(-x_j/b_{pj}\right){\mathbf{1}}_{x_j>0}.\hfill \end{array}$$

(138)

For option pricing, one must incorporate the marginal convexity corrections to match the forward prices for the two assets. This is the characteristic function ${\varphi }_{BGMP}\left(u\right)$, and, in the two dimensional case, it is given by

$$\begin{array}{ccc}\hfill {\varphi }_{BGMP}\left(u\right)& =& {\varphi }_{MBG}\left(u\right)exp(i{u}_1{\omega }_1+i{u}_2{\omega }_2)\hfill \end{array}$$

(139)

$$\begin{array}{ccc}\hfill {\omega }_1& =& ln\left(S_1\right)+(r-q_1)t-ln{\varphi }_{MBG}\left((-i,0)\right)\hfill \end{array}$$

(140)

$$\begin{array}{ccc}\hfill {\omega }_2& =& ln\left(S_2\right)+(r-q_2)t-ln{\varphi }_{MBG}\left((0,-1)\right).\hfill \end{array}$$

(141)

6. Sample Computations Using the Multivariate Bilateral Gamma Model

For an example on product option pricing, we consider product options on JPM and SPY, both of which have many options trading. On 12 December 2019 for maturities below three months and moneyness, measured by the absolute value of the logarithm of strike to the spot price, below 0.3, there were 226 options trading on JPM and 1582 options on SPY. The number of days to maturity of traded maturities closest to a month was 29 days. The maturity of 29 days had 36 strikes for JPM and 74 strikes for SPY. These options may be employed to determine the marginal risk neutral distribution on JPM and SPY for the 29 day maturity. The marginal bilateral gamma parameters are presented in

Table 1. Marginal BG parameters for the 29 day maturity on 12 December 2019.

Asset	${\mathit{b}}_{\mathit{p}}$	${\mathit{c}}_{\mathit{p}}$	${\mathit{b}}_{\mathit{n}}$	${\mathit{c}}_{\mathit{n}}$	$\mathit{\omega }$
JPM	0.0241	18.2249	0.0398	16.8810	0.2142
SPY	0.0270	2.2784	0.0509	6.6806	0.2693

.

Figure 1 presents a graph of the observed and model option prices for the 29 day maturity on 12 December 2019 for the two assets.

The two dependency parameters were set as follows. The value of ν = 0.4828 was ten percent above the lower bound and the correlation was set at 0.6. For this parameter setting and the spot for both assets set to 100 with strikes ranging between 80 and 120 Figure 2 presents four graphs for product option prices computed by two dimensional fast Fourier inversion of the Fourier transform of modified product option prices as per Equations (82), (91), (92), and (93).

With a view toward studying the effects of dependency parameters on product option prices, we considered four product options with strikes five percent out of the money. The options were a product of two calls, a call and a put, a put and a call, and two puts. There were five settings for correlation from −0.8 to 0.8 in steps of 0.4. There were four settings for the percentage excess of ν above its lower bound of 0.25, 0.5, 0.75, and 1.0. These results are presented in

Table 2. Product of two calls.

	Correlation
ν	−0.8	−0.4	0	0.4	0.8
0.25	55.47	64.67	75.08	86.75	99.72
0.5	58.35	66.15	74.85	84.48	95.08
0.75	60.44	67.21	74.69	82.89	91.84
1.0	62.04	68.02	74.57	81.71	89.46

,

Table 3. Product of a call and a put.

	Correlation
ν	−0.8	−0.4	0	0.4	0.8
0.25	98.06	86.33	75.49	65.55	56.52
0.5	94.82	85.13	76.11	67.76	60.08
0.75	92.52	84.28	76.55	69.35	62.68
1.0	90.82	83.64	76.88	70.55	64.65

,

Table 4. Product of a put and a call.

	Correlation
ν	−0.8	−0.4	0	0.4	0.8
0.25	90.90	80.32	70.53	61.53	53.35
0.5	87.80	79.06	70.92	63.37	56.41
0.75	85.61	78.18	71.20	64.69	58.65
1.0	83.98	77.51	71.41	65.69	60.35

and

Table 5. Product of two puts.

	Correlation
ν	−0.8	−0.4	0	0.4	0.8
0.25	62.59	71.39	80.77	90.72	101.3
0.5	64.47	71.86	79.68	87.92	96.60
0.75	65.82	72.19	78.89	85.94	93.31
1.0	66.85	72.44	78.31	84.45	90.87

, for the four option types.

For the product of two calls, the option price rises with correlation. The effect of ν is positive for negative correlation, negative for positive correlation, and flat for zero correlation.

For the product of a call and a put, the option price falls with correlation. The effect of ν is negative for negative correlation, positive for positive correlation, and flat for zero correlation.

For the product of a put and a call, the option price falls with correlation. The effect of ν is negative for negative correlation, positive for positive correlation, and flat for zero correlation.

For the product of two puts, the option price rises with correlation. The effect of ν is positive for negative correlation, negative for positive correlation, and flat for zero correlation.

7. Pricing Product Options Using the Physical Measure

For product options, one has the ability to partially hedge the risk using options on the two underlying assets. The prices of options on these assets are informative on the terms at which the product options could trade. However, as the hedge is not perfect, there is residual risk. Recently, Madan and Wang (2020b) reported on an investigation of option pricing at the cost of a hedge plus a charge for residual risk to levels of risk acceptability. For acceptable risk defined by concave distortions of probability as proposed in Cherny and Madan (2009) following Artzner et al. (1999) and Kusuoka (2001), the residual risk charge is just a distorted expectation of the simulated residual risk taken at an appropriate stress level for the distortion employed. For further details, we recommend, apart from the cited papers, Madan and Schoutens (2016). Here, we use the distortion minmaxvar with stress parameter $\gamma$ and definition

$${\Psi }^{\left(\gamma \right)}\left(u\right)=1-{\left(1-u^{\frac{1}{1+\gamma }}\right)}^{1+\gamma }.$$

For the physical measure, we bootstrap pairs of returns from the data immediately prior to the pricing date for the option maturity of interest. Bilateral gamma marginals are fit to this data. The dependency parameters of the multivariate bilateral gamma model are estimated by matching the two dimensional empirical characteristic function to the theoretical counterpart.

The estimated parameters are then used to simulate a hundred thousand readings on pairs of returns. On the simulated space, a hedge is determined by determining the funds needed to cover the residual risk to a level of risk acceptability. This magnitude is the residual risk charge. The cost of the hedge is then obtained from the market prices of the hedging instruments. The product option is priced at the cost of the hedge plus the residual risk charge. This procedure is illustrated on a set of strikes for all four types of product options on data for XLE and XLP as on 12 December 2019.

The first step is the formulation of joint returns on the two assets to the option maturity. We work with an option maturity of a month or 21 business days. However, an exact 21 days may not be a relevant way of forming joint returns as some months may be longer and others shorter in terms of economic time. We, therefore, take the number of days to be random with a gamma distribution with a mean of 21 days. For the shape parameter, we take the value of 2 that places the mode at half the mean.

The number of days n is simulated as one plus the integer part of the gamma variate. For the random number of days, we randomly select a start date within the last thousand days and an end date at the start date plus n. The returns on the two assets between the start and end dates delivers a single reading on joint returns. The procedure is repeated 500 times to draw a sample of 500 monthly joint returns.

In modeling the joint returns, we follow Madan and Wang (2020b). First we lock in a regression of the first asset return on the second and for the pair XLE and XLP, we obtain the results

$$r_{1n}=-0.0018+0.4006r_{2n}+u_n$$

The $t-statistics$ for this regression were −0.92, and 6.38. The $R-square$ was 6.22 percent.

We then model $r_{2n}$ and $u_n$ to be multivariate bilateral gamma. The marginal bilateral gamma parameters for $r_2$ and u are displayed in

Table 6. Marginal bilateral gamma parameters.

	${\mathit{b}}_{\mathit{p}}$	${\mathit{c}}_{\mathit{p}}$	${\mathit{b}}_{\mathit{n}}$	${\mathit{c}}_{\mathit{n}}$
$r_2$	0.0134	2.7719	0.0111	2.7719
u	0.0361	1.0546	0.0359	1.0547

.

The dependency parameters $\nu$ and $\rho$ were estimated by matching the empirical joint characteristic functions at the values 0.9482 and −0.4769, respectively. The drifts in the Brownian motions to be time changed were 0.0025 and 0.0003, respectively, for $r_2$ and $u$. The corresponding volatilities were 0.0177 and 0.0523. Setting the initial price levels for the two assets as they were at 31 December 2019 at 60.03 and 62.95 for XLE and XLP, respectively, we simulated 100,000 readings for the two asset prices a month later from this joint law for the two returns.

Consider first a product call at the strikes 61.21 and 64 for XLE and XLP, respectively. The payout on the product call is

$$c(S_1,S_2)={\left(S_1-61.21\right)}^{+}{(S_2-64)}^{+}.$$

By taking positions in the assets and options on the two assets one may access option payouts from the markets of $g\left(S_1\right)$ and $h\left(S_2\right).$ These functions may be built up from payouts to the assets and option on the assets as a function of the positions taken. One may construct matrices G and H that evaluate the payout on the assets and the options at the 100,000 readings on the two assets. G and H are then 100,000 times $n_1,n_2$ respectively, where $n_1,n_2$ are the number of assets in the two markets. With positions given by vectors $a_1,a_2$ of dimensions $n_1,n_2$, we have a hundred thousand readings on hedge cash flows as

$$\begin{array}{ccc}\hfill g\left(S_1\right)& =& G{a}_1\hfill \\ \hfill h\left(S_2\right)& =& H{a}_2.\hfill \end{array}$$

Evaluating the product option payout at the hundred thousand pairs of outcomes the product option payout is also a vector $c(S_1,S_2)$ of length a hundred thousand.

We now propose the construction of the hedge or the positions $a_1,a_2$ in the two sets of assets. We took for hedging assets, the two underlying assets and options at moneyness multiples of a percentage point between plus or minus 10 percent to get 19 and 12 hedging strikes on XLE and XLP, respectively. The post hedge residual cash flow is

$$R(S_1,S_2)=c(S_1,S_2)-g\left(S_1\right)-h\left(S_2\right).$$

The residual is the post hedge shortfall on the product option liability and it is to be valued at its ask or upper price, which is the risk charge for holding the residual risk.

The lower price for a risk using distorted expectations is given by its expectation evaluated under a probability distortion. The distorted probability distribution for a risk X with distribution function $F_X\left(x\right)$ is $\Psi \left(F_X\left(x\right)\right)$ for a concave distortion $\Psi \left(u\right),$$0\le u\le 1.$ The upper price is just the negative of distortion expectation for the negated risk. Defining the distorted expectation by

$$\mathcal{E}\left(X\right)={\int }_{-\infty }^{\infty }xd\Psi \left(F_X\left(x\right)\right)$$

the upper price or risk charge is

$$rc=-\mathcal{E}(-X).$$

Positions in the hedging assets are found by minimizing the residual risk charge $rc$.

The distortion employed, as already stated, is minmaxvar, and we need to pick the stress level $\gamma$ for the distortion. As distorted expectations of hedge fund returns turn negative at stress levels of 0.75, indicating acceptability at the margin for this stress level, we employed the stress level of $\gamma =0.75$. In constructing hedges, it is imperative to keep the hedging assets neutral, and this was done by centering their cash flows using the sample mean of the payouts for each asset. The minimal risk charge was 1.9899. Figure 3 displays the hedge cash flows accessed on the two assets. The cost of the hedge was 1.2843, and the price for the product call was 3.2743.

With a view to reporting further on the quality of the hedge, we present the percentiles for the hedge cash flows, inclusive of the risk charge, when the payout on the product call is zero. We also present, in

Table 7. Short fall and zero payout percentiles.

	1	5	10	25	50	75	90	95	99
shortfall	−20.78	−11.06	−7.24	−2.81	−0.71	−0.22	0.11	0.30	1.03
zero payout	0.22	0.23	0.23	0.23	0.24	0.70	3.47	8.34	20.95

, the percentiles of the shortfall on the product option payout, again inclusive of the risk charge.

Table 8. Other product options.

Type	K1	K2	rc	hc	POP
PC	57.21	63	1.2009	1.4320	2.6329
PP	60.21	61	0.7096	2.9561	3.6657
CP	62.21	62.5	1.1374	0.2	1.3374

Three types are PC, PP, and CP for products of a put and a call, two puts, and a call and a put. rc, hc, and POP are the risk charge, hedge cost, and the product option price, respectively.

presents the strikes, risk charges, hedge costs, and product option prices for other types of product options hedged using risk charge minimizing hedges.

8. Conclusions

Product options paying the product of put and or call options payouts at different strikes on two underlying assets were shown to be useful in synthesizing joint densities and replicating sufficiently differentiable functions of two underlying assets. Theory for pricing such options was presented and implemented from three perspectives. The first employed assumptions on asset log returns being jointly driven by correlated Brownian motion. The second used the two dimensional inversion of joint characteristic functions. The third was based on using the physical measure and proposing a risk charge minimizing hedge-using options on the two underlying assets. The options were then priced at the cost of the hedge plus the risk charge for the residual or unhedged risk.