Abstract

We present a method of deriving analytical solutions for a two-dimensional Black-Scholes-Merton equation. The method consists of three changes of variables in order to reduce the original partial differential equation (PDE) to a normal form and then solve it. Analytical solutions for two cases of option pricing on the minimum and maximum of two assets are derived using our method and are shown to agree with previously published results. The advantage of our solution procedure is the ability of splitting the original problem into several components in order to demonstrate some solution properties. The solutions of the two cases have a total of five components; each is a particular solution of the PDE itself. Due to the linearity of the two-dimensional Black-Scholes-Merton equation, any linear combination of these components constitutes another solution. Some other possible solutions as well as the solution properties are discussed.

1. Introduction

Black and Scholes [1] derived the solution to the value of a European-styled option on a stock under the assumptions that the stock follows a geometric Brownian motion. Starting with the geometric Brownian motion assumption, they derived a stochastic partial differential equation (PDE). With the key insight that the option could be combined with the underlying stock to create a hedged portfolio, they eliminated the stochastic portion of the PDE and derived the now well-known solution. To simplify the problem, they assumed the volatility of the stock and the discount rate were constant. Merton [2], besides discussing restrictions necessary for option pricing, also presents an alternate derivation of the Black-Scholes formula, allowing the interest rate to have a dynamic process rather than remaining constant. In addition, he describes the problem of pricing American-style options and derives explicit pricing formulas for both calls and puts, warrants, and “down-and-out” options.

Since the first result of the Black and Scholes [1] and Merton [2] solutions for European-style options, there have been many extensions, including early examples of Schwartz [3]; Cox, Ross, and Rubinstein (CRR) [4]; Rendleman and Bartter [5]; and Geske and Shastri [6]. Black and Scholes’ solution is also part of a recent study on bond options by Tomas and Yu [7]. Merton [2] provided an alternative derivation of the Black-Scholes formula that is valid under weaker assumptions and therefore more usable. Schwartz [3] applied finite difference methods to the option-pricing problem. CRR [4] and Rendleman and Bartter [5] study the binomial pricing of options as an alternative stochastic process to the Black-Scholes’ model. Geske and Shastri [6] examined some alternative numerical methods for valuing options, including Monte Carlo, binomial, and finite difference techniques. Recently, Tomas and Yu [7] studied call options on zero-coupon bonds, assuming a stochastic process for the price of the bond, rather than for interest rates in general. An asymptotic solution for the call option, with the leading order term being the Black and Scholes’ solution, was derived and studied.

There have also been extensions of option-pricing models to multiple dimensions. Early examples include Stulz [8], Johnson [9], Hull and White [10], Johnson and Shanno [11], and Gazizov and Ibragimov [12]. More recent examples include Villeneuve and Zanette [13], Grzelak et al. [14], Kim et al. [15], Sawangtong et al. [16], Naqeeb and Hussain [17], Wang et al. [18], and He and Lin [19, 20]. Stulz [8] examines the pricing of options on the minimum and maximum of two assets, where each asset follows the geometric Brownian motion process described in Black and Scholes. Johnson [9] extends the framework in Stulz to several assets. The original constant volatility assumption in the Black-Scholes/Merton models has been shown to fail to match market data. In order to address this, Hull and White [10] and Johnson and Shanno [11] consider the pricing of options on assets with stochastic volatility. They assume the stock price and its instantaneous variance are both stochastic. With this assumption, they each derive representative PDEs. Hull and White discuss the analytic solution to their PDE while Johnson and Shanno discuss the use of Monte Carlo simulation for the pricing. Since these original papers, there have been many other works examining nonconstant volatility. Grzelak et al. [14], for instance, provide an extension of stochastic volatility equity models with a stochastic Hull–White interest rate component. Very recently, He and Lin [19] modify the Heston-Hull-White model presented in Grzelak et al. [14] to capture the correlation between the stock price process and the interest rate while preserving analytical tractability.

Gazizov and Ibragimov [12] apply the Lie group theory to study these PDEs and provide examples for constructing exact (invariant) solutions. In a more recent example, Villeneuve and Zanette [13] examine the pricing of American-style options on two stocks that each follow the Black-Scholes framework, by adapting the alternating direction implicit algorithm developed by Peaceman and Rachford [21]. Kim et al. [15] apply a nonuniform grid finite difference approach where they add additional grid points near the nonsmooth portion of the payoff function to enhance accuracy and efficiency relative to a comparison Monte Carlo simulation. Sawangtong et al. [16] provide an analytic solution for the Black-Scholes equation with two assets based on the Liouville-Caputo fractional derivative. Naqeeb and Hussain [17] derive the multiasset Black-Scholes-Merton formulation and make links to the Hamilton-Jacobi equation of mechanics. Wang et al. [18] apply a finite difference method to the multidimensional fractional Black-Scholes model for the cases of one asset and three assets. He and Lin [20] examine exchange options, where the two assets each follow geometric Brownian motion. They introduce liquidity risk and changing economic conditions into the development of the model and modify the process by which the stocks are governed.

In this study, we revisit the two-dimensional problem in Stulz [8], in which there are two assets (stocks) both following the geometric Brownian motion process found in Black and Scholes [1]. Analytical solutions for two cases of option pricing on the minimum and maximum of two assets are derived using a three-step substitution method. The results are shown to agree with previously published results.

We derive the analytical solutions given in (11) and (11’) of Stulz [8], which are for a general case and a simpler case of exercise price , respectively. These solutions are derived in Stulz [8] by a statistical approach (see also Cox and Ross [22] and Harrison and Kreps [23]). Our method here utilizes applied mathematics techniques for solving the initial-boundary value problem (IBVP) of the PDE. In particular, we use the change of variables, including similarity variables, to convert the PDE to a simpler normal form and then solve it. The advantage of our solution procedure is the ability of splitting the problem into several components so that solution properties can be revealed. First, the original IBVP is converted into two or three sub-IBVPs, for the simpler or general cases, respectively. This demonstrates the linear superposition property of the problem. Secondly, in solving each sub-IBVP, we use three changes of variables, and each of these shows the solution properties such as geometric symmetry and scaling between variables. See the discussions below after Equation (23) and in Section 4.

Here is the IBVP from Stulz [8]: we look for a function, which satisfies the two-dimensional governing PDE: subject to the initial and boundary conditions:

The rest of the paper is organized as follows. Section 2 is devoted to the simpler case of where we derived the analytical solution (11’) of Stulz [8]. In Section 3, we study the general case and derived the analytical solution (11) of Stulz [8]. We then make some concluding remarks in Section 4.

2. Solution for a Simpler Case with the Exercise Price

We first study a simpler case of PDE (1) subject to the initial and boundary conditions (2) and (3) where This means that (1) and (3) are the same, but (2) becomes

Or

To solve IBVPs (1), (5), and (3), we notice that the problem is linear so that it can be split into two problems with , where satisfies (1), (3), and and satisfies (1), (3), and

Let us solve for first with the following three steps:

Step 1. We make a change of variables to eliminate the first three terms on the right-hand side of PDE (1), i.e.,

This is a first-order PDE, and the characteristic equations are

Integrating the first equation, we have or and the second equation, or where , and are constants. Therefore, a set of characteristics can be and Notice that with , a set of alternative characteristics can be (with ) and , and this will be used in (24) below. The set of characteristics ( and ) suggests the following change of variables:

Taking partial derivatives, we have

Now, we substitute these partial derivatives into PDE (1) to obtain where

Step 2. To convert PDE (12) to a normal form, we make a change of variables:

Taking partial derivatives, we have

Substituting these into (12), we find that

Step 3. Equation (15) can be solved by using similarity variables:

Taking partial derivatives, we have

Substituting these into (15), we find that and this is a first-order ordinary differential equation in . Integrating once, we get where is an arbitrary constant. Integrating one more time, we have where is also an arbitrary constant.

With use of Equations (10), (13), and (16), initial condition (6) for becomes the following condition for :

The condition if implies in (20), and from if , we have so that where .

Combining (22), (16), (13), and (10), we obtain the solution for :

The above solution procedure for (23) has the following three advantages: (1)The change of variables in (10) utilizes the characteristics derived right before (10) to simplify the PDE, and it also allows us to show, in Section 4 below, the geometric symmetry between variables and in (23) and (34).(2)The change of variables in (13) uses a combination of translation and scaling of variables, i.e., and respectively. The translation eliminates the term in (12), and the scaling simplifies the diffusivity constant from in (12) to 1/2 in (15)(3)The change of variables in (16) demonstrates the similarity variable property of the diffusion Equation (15); i.e., the and dependence of is such that is actually a function of only. See, for example, Kevorkian [24].

Now, let us solve for using the similar three steps listed above:

Step 1. We use the alternative characteristics ( and ) derived right before (10) to make a change of variables similar to (10):

Taking partial derivatives, we have

Substituting these partial derivatives into PDE (1), we obtain where

Step 2. To convert PDE (26) to a normal form, we make a change of variables similar to (13):

Taking partial derivatives, we have

Substituting these into (26), we find that

Step 3. Equation (29) is the same as (15). Therefore, we use the same similarity variables as in (16), i.e., to obtain the same similarity solution as in (20) for : where and are arbitrary constants.

With the use of Equations (24), (27), and (30), initial condition (7) for becomes the following condition for :

The condition if implies in (31), and from if we have so that where is defined in (22).

Combining (33), (30), (27), and (24), we obtain the solution for :

The solution procedure for (34) has similar advantages as that for (23). See the discussions immediately after Equation (23).

Finally, we combine in (23) with in (34) to get where and This is the solution for the simpler case where , as given in (11’) of Stulz [8]. (Notice that the initial condition (5) is also satisfied when )

3. Solution for the General Case where

We now focus on the general case ( not zero) of the Stulz problem, i.e., PDE (1), with initial and boundary conditions (2) and (3). Notice that the initial condition (2) can be rewritten as

Similar to the simpler case in Section 2, IBVPs (1), (3), and (36) are linear, and they can be split into three problems with , where satisfies (1), (3), and satisfies (1), (3), and and satisfies (1), (3), and

Let us solve for first using the three-step method described in Section 2.

Step 1. We make a change of variables similar to (24):

Taking partial derivatives, we have

Substituting these partial derivatives into PDE (1), we obtain where

Step 2. To convert PDE (42) to a normal form, we make a change of variables similar to (27):

Taking partial derivatives, we have

Substituting these into (42), we find that

Step 3. Equation (45) can be solved by using similarity variables similar to that in (16):

Taking partial derivatives, we have

Substituting these into (45) and denoting constant we obtain

With use of Equations (40), (43), and (46), initial condition (37) for becomes the following condition for

We now show that the solution to Equation (49) subject to condition (50) is the bivariate cumulative standard normal distribution function, i.e.,

Taking partial derivatives in (51), we have

We rewrite the integral in (54) as so that it can be split into two parts with

We make a substitution, in This implies , and the integral can be carried out as

Comparing this with (53) and with (52), we obtain

By the simple symmetry between variables and in (51), we have

Adding (58) to (59), we find that Equation (49) is satisfied. Condition (50) can also be verified due to the fact that the in (51) is the bivariate cumulative standard normal distribution function.

Combining (51), (46), (43), and (40), we obtain the solution for as and this is the first term of solution (11) in Stulz [8]. Notice that there was a typo in the first term of solution (11) of Stulz [8]. The expression of the second input variable inside the function should be as given in (60), not as in the first term of solution (11) of Stulz [8]. There was another similar typo in the second term of solution (11) of Stulz [8]: the expression of the second input variable inside the function should be as given in (A.12) below for , not

The solution procedure for (60) has similar advantages as that for (23). See the discussions immediately after Equation (23).

The solution procedures for and are similar to that for using the three-step procedure of change of the variables given above. These are presented in Appendixes A and B, respectively, and omitted here for brevity.

Finally, by adding , , and together, from (60), (A.12), and (B.12), respectively, we arrive at solution (11) of Stulz [8].

Many of Stulz’ solutions have been programmed in MATLAB; see Stulz Model [25] in the References. We refer readers to the link there for example graphs of the solutions.

4. Concluding Remarks

We use a three-step substitution method to derive solution (11’) of Stulz [8] in Section 2 and solution (11) of Stulz [8] in Section 3. Notice that there is a typo in the first two terms of solution (11) of Stulz [8] (see the discussion after Equation (60)).

In the derivation process, we find five particular solutions, in (23), in (34), in (60), in (A.12), and in (B.12), for the two-dimensional Black-Scholes-Merton Equation (1). In fact, there are even simpler particular solutions of Equation (1), for examples, is a solution, is also a solution, and is also a solution. As given in Haug [26], there are a variety of exotic options on two risky assets, and our method may be helpful for studying their solutions. For example, for an exchange-one-asset-for-another option and the call option on the maximum of two assets, we can similarly split the problem and solve them using our three-step substitution method. The call option on the maximum of two assets, as a second example, was discussed in Stulz [8] using parity relationships. Our method here rederives the solution as given in Haug [26]. These two examples are presented in Appendix C (Two Additional Examples of Call Options on Two Risky Assets), which are omitted here for brevity. As a result, we additionally find five particular solutions of Equation (1), i.e., in (C.7), in (C.9), in (C.21), in (C.22), and in (C.23). In fact, is also a solution. Altogether in this paper, we have 14 particular solutions of Equation (1). Because (1) is linear, any linear combination of these solutions is also a solution. Therefore, these particular solutions are the building blocks for constructing solutions of various IBVPs. These as well as our solution technique could be used to explore the solutions of other exotic option-pricing problems on two assets.

Some of these solutions we obtained thus far have good geometric symmetry between variables. For example, solutions and are symmetric in variables and i.e., and Using the identity we can show that solutions in (23) and in (34) are also symmetric in variables and , i.e., and Finally, solutions in (60), in (A.12), in (B.12), in (C.21), in (C.22), and in (C.23) all have geometric symmetry with respect to their respective transformation variables and as show in Equation (51).

Our method may not be applied directly to American options with free boundaries or the case with stochastic volatility. For American options, even if we can reduce the partial differential equation to a simpler form, its analytical solution with the free boundaries is still difficult to find. For the case with stochastic volatility, while our method might work well with the option price dependence on the stock value, the price dependence on the variance is hard to study using our current method. Either of these directions requires further investigations and could be future research projects.

Appendix

A. Derivation of Solution for

Here, we derive the solution, , to the PDE (1) subject to conditions (3) and (38), using the three-step method described in Sections 2 and 3.

Step 1. We make a change of variables similar to (10) and (40):

Taking partial derivatives, we have

Substituting these partial derivatives into PDE (1), we obtain where

Step 2. To convert PDE (A.3) to a normal form, we make a change of variables similar to (43):

Taking partial derivatives, we have

Substituting these into (A.3), we find that

Step 3. Equation (A.6) is the same as (45), just with a different constant, and it can be solved by using the similarity variables given in (46), i.e.,

Taking partial derivatives, substituting into (A.6), and denoting constant we obtain

With the use of Equations (A.1), (A.4), and (A.7), the initial condition (38) for becomes the following condition for

Since PDE (A.9) and condition (A.10) are the same as PDE (49) and condition (50), we have solution (51) for , i.e.,

Combining (A.11), (A.7), (A.4), and (A.1), we obtain the solution for as and this is the second term of the solution (11) in Stulz [8].

B. Derivation of Solution for

Here, we derive the solution, , to the PDE (1) subject to conditions (3) and (39), using the three-step method described in Sections 2 and 3.

Step 1. We make a change of variables similar to (10) and (40):

Taking partial derivatives, we have

Substituting these partial derivatives into PDE (1), we obtain

Step 2. To convert PDE (B.3) to a normal form, we make a change of variables similar to (43):

Taking partial derivatives, we have

Substituting these into (B.3), we find that

Step 3. Equation (B.6) is the same equation as (45), just with a different constant, and it can be solved by using the similarity variables given in (46), i.e.,

Taking partial derivatives, substituting into (B.6), and denoting constant we obtain

With use of Equations (B.1), (B.4), and (B.7), initial condition (39) for becomes the following condition for

Since PDE (B.9) and condition (B.10) are the same as PDE (49) and condition (50), we have solution (51) for , i.e.,

Combining (B.11), (B.7), (B.4), and (B.1), we obtain the solution for as where and are defined in (60) and (A.12), respectively. This solution given in (B.12) is the third term of solution (11) in Stulz [8].

C. Two Additional Examples of Call Options on Two Risky Assets

C.1. Exchange-One-Asset-for-Another Options

The IBVP for an exchange-one-asset-for-another option consists of PDE (1) with boundary condition (3) and initial condition where and are the quantity of assets and , respectively (see Margrabe [27]). Notice that the initial condition (C.1) can be rewritten as

Similar to the case in Section 2, we split IBVPs (1), (3), and (C.2) into two subproblems with , where satisfies (1), (3), and and satisfies (1), (3), and

Let us summarize the solution steps for using the three-step method described in Section 2. Similar to (10), we first use the change of variables: and we find that PDE (1) becomes (12). As shown in Section 2, Equation (12) can be converted into (15) by the change of variables (13). Using the similarity variable (16), Equation (15) can be solved to obtain (20). With use of (C.5), (13), and (16), initial condition (C.3) becomes as follows:

Therefore, the solution for is given by

Now, let us summarize the solution steps for using the three-step method described in Section 2. Similar to (24), we first use the change of variables and we find that PDE (1) becomes (26). As shown in Section 2, Equation (26) can be converted into (29) by the change of variables (27). Using the similarity variable (30), Equation (29) can be solved to obtain (31). With use of (C.8), (27), and (30), initial condition (C.4) becomes (C.6). Therefore, the solution for is given by

Combining (C.7) and (C.9) with , we obtain the solution to the option to exchange one asset for another as

C.2. Call Option on the Maximum of Two Assets

The IBVP for a call option on the maximum of two assets consists of PDE (1) with boundary condition (3) and initial condition

Notice that the initial condition (C.11) can be rewritten as

Similar to the case in Section 3, IBVPs (1), (3), and (C.12) can be split as , where satisfies (1), (3), and satisfies (1), (3), and and satisfies (1), (3), and

Let us summarize the solution steps for using the three-step method described in Section 3: (1)Similar to (40), we first use the change of variables:and PDE (1) becomes where (2)Similar to (43), we then use the change of variables:and (C.17) becomes

Finally, we use the change of variables (46) to obtain (49) with

With use of (C.16), (C.18), and (46), initial condition (C.13) becomes as follows:

Therefore, the solution for is given by

A similar procedure as above (see also the derivation in Appendix A) gives as

A procedure similar to that in Appendix B gives as where and are given in (60) and (A.12), respectively.

Finally, we use to obtain the solution to the call option on the maximum of two assets, as given in Section 5.7 of Haug [26].

Data Availability

No underlying data was collected or produced in this study.

Conflicts of Interest

The authors declare that they have no conflicts of interest.