Some Sequential Boundary Crossing Results for Geometric Brownian Motion and Their Applications in Financial Engineering

Abstract

This paper provides new explicit results for some boundary crossing distributions in a multidimensional geometric Brownian motion framework when the boundary is a piecewise constant function of time. Among their various possible applications, they enable accurate and efficient analytical valuation of a large number of option contracts traded in the financial markets belonging to the classes of barrier and look-back options.

1. Introduction

The joint law of the maximum (or the minimum) of a real-valued Brownian motion and its endpoint over a finite time interval is a central result in the study of Brownian motion, particularly with regard to the many applications of the theory in finance, medical imaging, robotics, and biology. It can be obtained as a consequence of the “reflection principle,” which derives from the strong Markov property of Brownian motion (Freedman [1]). Application of Girsanov’s theorem easily generalizes this seminal result to the case of a geometric Brownian motion (GBM), a frequently encountered diffusion process that is the building block of financial engineering. Alternatively, the law can be derived by a partial differential equation approach, using Kolmogorov’s equation for the transition density function of a diffusion process. The distribution of the first passage time by a one-dimensional GBM to a one-sided or a two-sided straight boundary then follows. A few cases where the boundary is curved have been handled (Barba Escribá [2]; Salminen [3]; Kunitomo and Ikeda [4]). For a comprehensive source of formulae, one may refer to Borodin and Salminen [5].

Various extensions to these results are often needed to solve practical engineering problems. In particular, one may look for joint distributions of the highest or lowest points hit over several time intervals, and one may deal with a multidimensional GBM. There are few known explicit results in these more general settings, partly because of the laborious analytical calculations involved, but also due to numerical obstacles: the rapidly increasing dimension of the boundary crossing problem leads to analytical solutions that are expressed in terms of functions that are hard to compute with accuracy. As far as one-sided boundaries are concerned, formulae have been published for the joint law of a sequence of maxima or minima of a one-dimensional GBM over several time intervals (Guillaume [6]), as well as for the joint law of the maxima or minima of a two-dimensional GBM over one time interval (Iyengar [7]; He et al. [8]). A formula for the joint law of the exit times of a one-dimensional GBM from two successive two-sided boundaries is also known (Guillaume [9]).

This paper focuses on a sequence of two one-sided straight boundaries conditional on two correlated GBMs, while the value of a third correlated GBM is taken into account at the endpoint of the time interval. The state space is thus three-dimensional. The choice of this particular distribution is motivated both by its usefulness in financial engineering applications and by the fact that it leads to tractable analytical solutions that can be computed with great accuracy and efficiency.

Section 2 presents the two main formulae of this paper. Section 3 deals with applications of these formulae and discusses their numerical implementation.

2. Main Formulae

This section contains the two main formulae of the paper. Let {𝑆₁(𝑡),𝑡≥0} and {𝑆₂(𝑡),𝑡≥0} be two geometric Brownian motions with constant drift coefficients, 𝛼₁ and 𝛼₂, respectively, under a probability measure ℙ, constant positive diffusion coefficients 𝜎₁ and 𝜎₂, respectively, and constant correlation coefficient 𝜌_1⋅2. In other words, the processes 𝑆₁(𝑡) and 𝑆₂(𝑡) evolve in time according to the following dynamics: 𝑑𝑑𝑆𝑆₁₂(𝑡)=𝛼₁𝑆₁(𝑡)𝑑𝑡+𝜎₁𝑆₁(𝑡)𝑑𝐵₁(𝑡),(𝑡)=𝛼₂𝑆₂(𝑡)𝑑𝑡+𝜎₂𝑆₂(𝑡)𝑑𝐵₂(𝑡),(2.1) where 𝐵₁(𝑡) and 𝐵₂(𝑡) are standard real-valued Brownian motions and 𝑑[𝐵₁,𝐵₂](𝑡)=𝜌_1⋅2𝑑𝑡. Let 𝐻₁, 𝐻₂, 𝐾₁, 𝐾₂, 𝐾₃ be positive real numbers and 𝑃_[⋅](𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝑡₁,𝑡₂) defined as one of the four following cumulative distribution functions, where 𝑡₂≥𝑡₁≥𝑡₀=0: (1)𝑃≜ℙ_[up-and-up]𝐻su₁p,_{0≤𝑡≤𝑡1}𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝑡₁,𝑡₂𝑆₁(𝑡)≤𝐻₁,𝑆₁𝑡₁≤𝐾₁,𝑆₂𝑡₁≤𝐾₂,sup_{𝑡1≤𝑡≤𝑡2}𝑆₂(𝑡)≤𝐻₂,𝑆₂𝑡₂≤𝐾₃,(2.2)(2)𝑃≜ℙin_{[down-and-down]}f_{0≤𝑡≤𝑡1}𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝑡₁,𝑡₂𝑆₁(𝑡)≥𝐻₁,𝑆₁𝑡₁≥𝐾₁,𝑆₂𝑡₁≥𝐾₂,inf_{𝑡1≤𝑡≤𝑡2}𝑆₂(𝑡)≥𝐻₂,𝑆₂𝑡₂≥𝐾₃,(2.3)(3)𝑃≜ℙ_{[down-and-up]}in𝐻f_{10≤𝑡≤𝑡1},𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝑡₁,𝑡₂𝑆₁(𝑡)≥𝐻₁,𝑆₁𝑡₁≥𝐾₁,𝑆₂𝑡₁≤𝐾₂,sup_{𝑡1≤𝑡≤𝑡2}𝑆₂(𝑡)≤𝐻₂,𝑆₂𝑡₂≤𝐾₃,(2.4)(4)𝑃≜ℙ_{[up-and-down]}su𝐻p_{10≤𝑡≤𝑡1},𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝑡₁,𝑡₂𝑆₁(𝑡)≤𝐻₁,𝑆₁𝑡₁≤𝐾₁,𝑆₂𝑡₁≥𝐾₂,inf_{𝑡1≤𝑡≤𝑡2}𝑆₂(𝑡)≥𝐻₂,𝑆₂𝑡₂≥𝐾₃.(2.5)

Let us introduce the following notations:𝑘₁=ln𝐾¹𝑆₁𝑡₀,𝑘²=ln𝐾²𝑆₂𝑡₀,𝑘³=ln𝐾³ℎ₁𝑆₂𝑡₀,=ln𝐻¹𝑆₁𝑡₀,ℎ²=ln𝐻²𝜇₁𝑆₂𝑡₀,=𝛼₁−𝜎²¹2,𝜇²=𝛼₂−𝜎²²2.(2.6) The next proposition provides an exact formula for the four above-mentioned cumulative distribution functions.

Proposition 2.1. Let 𝑁₃[⋅,⋅,⋅;𝜃_1⋅2,𝜃_1⋅3,𝜃_2⋅3] denote the joint trivariate cumulative distribution function of three standard normal random variables 𝑋₁, 𝑋₂, 𝑋₃, where 𝜃_𝑎⋅𝑏 is the correlation coefficient between 𝑋_𝑎 and 𝑋_𝑏, (𝑎,𝑏)∈{1,2,3}.
Then, for the up-and-up and the down-and-down distributions, one has P_[⋅]=𝐻₁𝑁,𝐻₂₃,𝐾₁,𝐾₂,𝐾₃,𝑡₁,𝑡₂𝜆𝑘¹−𝜇₁𝑡₁𝜎₁√𝑡₁,𝜆𝑘²−𝜇₂𝑡₁𝜎₂√𝑡₁,𝜆𝑘³−𝜇₂𝑡₂𝜎₂√𝑡₂;𝜌^1⋅2,𝑡₁𝑡₂𝜌_1⋅2,𝑡₁−exp2𝜇₂𝑡₂ℎ₂𝜎²²𝑁³⎡⎢⎢⎢⎢⎢⎣𝜆𝑘₁−𝜇₁𝑡₁𝜎₁√𝑡₁+𝜌^1⋅22𝜇₂√𝑡₁𝜎₂,𝜆𝑘²+𝜇₂𝑡₁𝜎₂√𝑡₁,𝜆𝑘³−2ℎ₂−𝜇₂𝑡₂𝜌_1⋅2𝜎₂√𝑡₂;,−𝜌_1⋅2𝑡₁𝑡₂,−𝑡₁𝑡₂−exp2𝜇₁⎤⎥⎥⎥⎥⎥⎦ℎ₁×𝑁₃𝜎²¹⎡⎢⎢⎢⎢⎢⎣𝜆𝑘₁−2ℎ₁−𝜇₁𝑡₁𝜎₁√𝑡₁,𝜆𝑘²−𝜇₂𝑡₁𝜎₂√𝑡₁−𝜌_1⋅22ℎ₁𝜎₁√𝑡₁,𝜆𝑘³−𝜇₂𝑡₂𝜎₂√𝑡₂−𝜌_1⋅22ℎ₁𝜌_1⋅2𝜎₁√𝑡₂;,𝜌_1⋅2𝑡₁𝑡₂,𝑡₁𝑡₂+exp2𝜇₁⎤⎥⎥⎥⎥⎥⎦𝜎²¹−4𝜇²𝜌_1⋅2𝜎₁𝜎₂ℎ¹+2𝜇²ℎ₂×𝑁₃𝜎²²⎡⎢⎢⎢⎢⎢⎣𝜆𝑘₁−2ℎ₁−𝜇₁𝑡₁𝜎₁√𝑡₁+𝜌^1⋅22𝜇₂√𝑡₁𝜎₂,𝜆𝑘²+𝜇₂𝑡₁𝜎₂√𝑡₁−𝜌_1⋅22ℎ₁𝜆𝑘₃𝜎₁√𝑡₁,−2ℎ₂−𝜇₂𝑡₂𝜎₂√𝑡₂+𝜌_1⋅22ℎ₁𝜎₁√𝑡₂;𝜌^1⋅2,−𝜌_1⋅2𝑡₁𝑡₂,−𝑡₁𝑡₂⎤⎥⎥⎥⎥⎥⎦,(2.7) where 𝜆=𝜆−=11ifif𝑃𝑃_[⋅][⋅]𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝑡₁,𝑡₂=𝑃_[up-and-up]𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝑡₁,𝑡₂,𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝑡₁,𝑡₂=𝑃_{[down-and-down]}𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝑡₁,𝑡₂.(2.8) And, for the up-and-down and the down-and-up distributions, one has 𝑃_[⋅]=𝐻₁𝑁,𝐻₂₃,𝐾₁,𝐾₂,𝐾₃,𝑡₁,𝑡₂𝜆𝑘¹−𝜇₁𝑡₁𝜎₁√𝑡₁,𝜆−𝑘²+𝜇₂𝑡₁𝜎₂√𝑡₁,𝜆−𝑘³+𝜇₂𝑡₂𝜎₂√𝑡₂;−𝜌^1⋅2,−𝑡₁𝑡₂𝜌_1⋅2,𝑡₁−exp2𝜇₂𝑡₂ℎ₂×𝑁₃𝜎²²⎡⎢⎢⎢⎢⎢⎣𝜆𝑘₁−𝜇₁𝑡₁𝜎₁√𝑡₁+𝜌^1⋅22𝜇₂√𝑡₁𝜎₂,𝜆−𝑘²−𝜇₂𝑡₁𝜎₂√𝑡₁,𝜆−𝑘³+2ℎ₂+𝜇₂𝑡₂−𝜌_1⋅2𝜎₂√𝑡₂;,𝑡₁𝑡₂𝜌_1⋅2,−𝑡₁𝑡₂−exp2𝜇₁⎤⎥⎥⎥⎥⎥⎦ℎ₁×𝑁₃𝜎²¹⎡⎢⎢⎢⎢⎢⎣𝜆𝑘₁−2ℎ₁−𝜇₁𝑡₁𝜎₁√𝑡₁,𝜆−𝑘²+𝜇₂𝑡₁𝜎₂√𝑡₁+𝜌_1⋅22ℎ₁𝜎₁√𝑡₁,𝜆−𝑘³+𝜇₂𝑡₂𝜎₂√𝑡₂+𝜌_1⋅22ℎ₁−𝜌_1⋅2𝜎₁√𝑡₂;,−𝑡₁𝑡₂𝜌_1⋅2,𝑡₁𝑡₂+exp2𝜇₁⎤⎥⎥⎥⎥⎥⎦𝜎²¹−4𝜇²𝜌_1⋅2𝜎₁𝜎₂ℎ¹+2𝜇²ℎ₂×𝑁₃𝜎²²⎡⎢⎢⎢⎢⎢⎣𝜆𝑘₁−2ℎ₁−𝜇₁𝑡₁𝜎₁√𝑡₁+𝜌^1⋅22𝜇₂√𝑡₁𝜎₂,𝜆−𝑘²−𝜇₂𝑡₁𝜎₂√𝑡₁+𝜌_1⋅22ℎ₁𝜆−𝑘₃𝜎₁√𝑡₁,+2ℎ₂+𝜇₂𝑡₂𝜎₂√𝑡₂−𝜌_1⋅22ℎ₁𝜎₁√𝑡₂;−𝜌^1⋅2,𝑡₁𝑡₂𝜌_1⋅2,−𝑡₁𝑡₂⎤⎥⎥⎥⎥⎥⎦,(2.9) where 𝜆𝜆==−11ifif𝑃𝑃_[⋅][⋅]𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝑡₁,𝑡₂=𝑃_{[up-and-down]}𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝑡₁,𝑡₂,𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝑡₁,𝑡₂=𝑃_{[down-and-up]}𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝑡₁,𝑡₂.(2.10)

Corollary of Proposition 2.1
The four following joint cumulative distribution functions, that will be useful in Section 3, are deduced from Proposition 2.1:(1)ℙsup_{0≤𝑡≤𝑡1}𝑆₁(𝑡)≤𝐻₁,𝑆₁𝑡₁≤𝐾₁,𝑆₂𝑡₁≤𝐾₂,sup_{𝑡1≤𝑡≤𝑡2}𝑆₂(𝑡)≤𝐻₂,𝑆₂𝑡₂≥𝐾₃=𝑃_[up-and-up]𝐻₁,𝐻₂,𝐾₁,𝐾₂,∞−𝑃_[up-and-up]𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,(2.11)(2)ℙinf_{0≤𝑡≤𝑡1}𝑆₁(𝑡)≥𝐻₁,𝑆₁𝑡₁≥𝐾₁,𝑆₂𝑡₁≥𝐾₂,inf_{𝑡1≤𝑡≤𝑡2}𝑆₂(𝑡)≥𝐻₂,𝑆₂𝑡₂≤𝐾₃=𝑃_{[down-and-down]}𝐻₁,𝐻₂,𝐾₁,𝐾₂,−∞−𝑃_{[down-and-down]}𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,(2.12)(3)ℙinf_{0≤𝑡≤𝑡1}𝑆₁(𝑡)≥𝐻₁,𝑆₁𝑡₁≥𝐾₁,𝑆₂𝑡₁≤𝐾₂,sup_{𝑡1≤𝑡≤𝑡2}𝑆₂(𝑡)≤𝐻₂,𝑆₂𝑡₂≥𝐾₃=𝑃_{[down-and-up]}𝐻₁,𝐻₂,𝐾₁,𝐾₂,∞−𝑃_{[down-and-up]}𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,(2.13)(4)ℙsup_{0≤𝑡≤𝑡1}𝑆₁(𝑡)≤𝐻₁,𝑆₁𝑡₁≤𝐾₁,𝑆₂𝑡₁≥𝐾₂,inf_{𝑡1≤𝑡≤𝑡2}𝑆₂(𝑡)≥𝐻₂,𝑆₂𝑡₂≤𝐾₃=𝑃_{[up-and-down]}𝐻₁,𝐻₂,𝐾₁,𝐾₂,−∞−𝑃_{[up-and-down]}𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃.(2.14)

Proof of Proposition 2.1 is provided in the Appendices A, B, and C.

The numerical implementation of Proposition 2.1 and its corollary is easy using Genz’s algorithm for the computation of trivariate normal cumulative distribution functions (Genz [10]).

In the next Proposition, we introduce a third correlated geometric Brownian motion that will serve as the endpoint of the joint distribution, and we show that this can still be analytically valued. Let {𝑆₃(𝑡),𝑡≥0} be a geometric Brownian motion driven by the following dynamics under the initial probability measure ℙ:𝑑𝑆₃(𝑡)=𝛼₃𝑆₃(𝑡)𝑑𝑡+𝜎₃𝑆₃(𝑡)𝑑𝐵₃(𝑡),(2.15) where 𝛼₃ and 𝜎₃ are real constants (𝜎₃>0), 𝐵₃(𝑡) is a real-valued standard Brownian motion and correlations are as follows:𝑑𝐵₁,𝐵₃(𝑡)=𝜌_1⋅3𝑑𝑡,𝑑𝐵₂,𝐵₃(𝑡)=𝜌_2⋅3𝑑𝑡.(2.16) Let 𝐾₄ be a positive real number and 𝑃_[⋅](𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄,𝑡₁,𝑡₂,𝑡₃) defined as one of the four following joint cumulative distribution functions, where 𝑡₃≥𝑡₂: (1)𝑃≜ℙ_[up-and-up]su𝐻p_{10≤𝑡≤𝑡1},𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄,𝑡₁,𝑡₂,𝑡₃𝑆₁(𝑡)≤𝐻₁,𝑆₁𝑡₁≤𝐾₁,𝑆₂𝑡₁≤𝐾₂,sup_{𝑡1≤𝑡≤𝑡2}𝑆₂(𝑡)≤𝐻₂,𝑆₂𝑡₂≤𝐾₃,𝑆₃𝑡₃≤𝐾₄(2,.17)(2)𝑃≜ℙinf_{[down-and-down]0≤𝑡≤𝑡1}𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄,𝑡₁,𝑡₂,𝑡₃𝑆₁(𝑡)≥𝐻₁,𝑆₁𝑡₁≥𝐾₁,𝑆₂𝑡₁≥𝐾₂,inf_{𝑡1≤𝑡≤𝑡2}𝑆₂(𝑡)≥𝐻₂,𝑆₂𝑡₂≥𝐾₃,𝑆₃𝑡₃≥𝐾₄(2,.18)(3)𝑃≜ℙin_{[down-and-up]}f_{0≤𝑡≤𝑡1}𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄,𝑡₁,𝑡₂,𝑡₃𝑆₁(𝑡)≥𝐻₁,𝑆₁𝑡₁≥𝐾₁,𝑆₂𝑡₁≤𝐾₂,sup_{𝑡1≤𝑡≤𝑡2}𝑆₂(𝑡)≤𝐻₂,𝑆₂𝑡₂≤𝐾₃,𝑆₃𝑡₃≤𝐾₄(2,.19)(4)𝑃≜ℙsu_{[up-and-down]}p_{0≤𝑡≤𝑡1}𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄,𝑡₁,𝑡₂,𝑡₃𝑆₁(𝑡)≤𝐻₁,𝑆₁𝑡₁≤𝐾₁,𝑆₂𝑡₁≥𝐾₂,inf_{𝑡1≤𝑡≤𝑡2}𝑆₂(𝑡)≥𝐻₂,𝑆₂𝑡₂≥𝐾₃,𝑆₃𝑡₃≥𝐾₄(2..20)

Let us introduce the following new notations:𝑘₄=ln𝐾⁴𝑆₃𝑡₀,𝜇³=𝛼₃−𝜎²³2.(2.21) The next proposition provides an exact formula for the four above-mentioned joint cumulative distribution functions.

Proposition 2.2. Let the real function Φ[𝑏₁,𝑏₂,𝑏₃,𝑏₄;𝜃_1⋅2,𝜃_1⋅3,𝜃_1⋅4,𝜃_2⋅3,𝜃_3⋅4], where {𝑏₁,𝑏₂,𝑏₃,𝑏₄}∈ℝ⁴ and each of the real numbers 𝜃_1⋅2,𝜃_1⋅3,𝜃_1⋅4,𝜃_2⋅3,𝜃_3⋅4 is included in ]−1,1[, be defined as follows: Φ𝑏₁=,𝑏₂,𝑏₃,_{𝑏1𝑥1=−∞}𝑏₄;𝜃_1⋅2,𝜃_1⋅3,𝜃_1⋅4,𝜃_2⋅3,𝜃_3⋅4^𝑏2_𝑥2=−∞^𝑏3_𝑥3=−∞1𝜙_2∣1𝜙_3∣2√8𝜋³exp−𝑥²²2−𝑥¹−𝜃_1⋅2𝑥₂²2𝜙^22∣1−𝑥³−𝜃_2⋅3𝑥₂²×𝑁⎡⎢⎢⎢⎢⎣𝑏₄2𝜙^23∣2−𝜃_1⋅4𝑥₁−𝜃_3⋅4∣1𝑥³−𝜃_1⋅3𝑥₁/1−𝜃^21⋅3𝜙_4∣1⋅3⎤⎥⎥⎥⎥⎦𝑑𝑥₃𝑑𝑥₂𝑑𝑥₁,(2.22) where 𝑁[⋅] is the standard normal cumulative distribution function and the following definitions apply: 𝜙_2∣1=1−𝜃^21⋅2,𝜙_3∣2=1−𝜃^22⋅3,𝜃_3⋅4∣1=𝜃^3⋅4−𝜃_1⋅3𝜃_1⋅4𝜙_3∣1,𝜙_4∣1⋅3=1−𝜃^21⋅4−𝜃^23⋅4∣1.(2.23) Then, for the up-and-up and the down-and-down distributions, one has P_[⋅]=𝐻₁Φ,𝐻⎡⎢⎢⎢⎢⎢⎣₂𝜆,𝐾₁𝑘,₁𝐾₂,𝐾₃,𝐾₄,𝑡₁,𝑡₂,𝑡₃−𝜇₁𝑡₁𝜎₁√𝑡₁,𝜆𝑘²−𝜇₂𝑡₁𝜎₂√𝑡₁,𝜆𝑘³−𝜇₂𝑡₂𝜎₂√𝑡₂,𝜆𝑘⁴−𝜇₃𝑡₃𝜃_1⋅2𝜎₃√𝑡₃;=𝜌_1⋅2,𝜃_1⋅3=𝑡₁𝑡₂𝜌_1⋅2,𝜃_1⋅4=𝑡₁𝑡₃𝜌_1⋅3,𝜃_2⋅3=𝑡₁𝑡₂,𝜃_3⋅4=𝑡₂𝑡₃𝜌_2⋅3−exp2𝜇₂⎤⎥⎥⎥⎥⎥⎦ℎ₂×Φ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣𝜆𝑘𝜎²²₁−𝜇₁𝑡₁𝜎₁√𝑡₁+𝜌^1⋅22𝜇₂𝑡₁𝜎₂√𝑡₁,𝜆𝑘²+𝜇₂𝑡₁𝜎₂√𝑡₁,𝜆𝑘³−2ℎ₂−𝜇₂𝑡₂𝜆⎛⎜⎜⎜⎝𝑘₄𝜎₂√𝑡₂,−𝜇₃𝑡₃𝜎₃√𝑡₃−𝜃^3⋅4−𝜃_1⋅3𝜃_1⋅41−𝜃^21⋅3𝜎2₂√ℎ²𝑡₂+𝜃_1⋅3𝜌_1⋅22𝜇₂𝑡₁𝜎₂√𝑡₁+𝜃^1⋅4𝜌_1⋅22𝜇₂𝑡₁𝜎₂√𝑡₁𝜃_1⋅2⎞⎟⎟⎟⎠;=𝜌_1⋅2,𝜃_1⋅3=𝑡₁𝑡₂𝜌_1⋅2,𝜃_1⋅4=𝑡₁𝑡₃𝜌_1⋅3,𝜃_2⋅3=−𝑡₁𝑡₂,𝜃_3⋅4=𝑡₂𝑡₃𝜌_2⋅3−exp2𝜇₁⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ℎ₁×Φ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣𝜆𝑘𝜎²¹₁−2ℎ₁−𝜇₁𝑡₁𝜎₁√𝑡₁,𝜆𝑘²−𝜇₂𝑡₁𝜎₂√𝑡₁−𝜌_1⋅22ℎ₁𝜎₁√𝑡₁,𝜆𝑘³−𝜇₂𝑡₂𝜎₂√𝑡₂−𝜌_1⋅22ℎ₁𝜆⎛⎜⎜⎜⎝𝑘₄𝜎₁√𝑡₂,−𝜇₃𝑡₃𝜎₃√𝑡₃−𝜃_1⋅42ℎ₁𝜎₁√𝑡₁−𝜃^3⋅4−𝜃_1⋅3𝜃_1⋅41−𝜃^21⋅3𝜌^1⋅22ℎ₁𝜎₁√𝑡₂−𝜃_1⋅32ℎ₁𝜃_1⋅2𝜎₁√𝑡₁⎞⎟⎟⎟⎠;=𝜌_1⋅2,𝜃_1⋅3=𝑡₁𝑡₂𝜌_1⋅2,𝜃_1⋅4=𝑡₁𝑡₃𝜌_1⋅3,𝜃_2⋅3=𝑡₁𝑡₂,𝜃_3⋅4=𝑡₂𝑡₃𝜌_2⋅3+exp2𝜇₁⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦𝜎²¹−4𝜇²𝜌_1⋅2𝜎₁𝜎₂ℎ¹+2𝜇²ℎ₂×Φ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣𝜆𝑘₁𝜎²²−2ℎ₁−𝜇₁𝑡₁𝜎₁√𝑡₁+𝜌^1⋅22𝜇₂𝑡₁𝜎₂√𝑡₁,𝜆𝑘²+𝜇₂𝑡₁𝜎₂√𝑡₁−𝜌_1⋅22ℎ₁𝜆𝑘₃𝜎₁√𝑡₁,−2ℎ₂−𝜇₂𝑡₂𝜎₂√𝑡₂+𝜌_1⋅22ℎ₁𝜆⎛⎜⎜⎜⎜⎜⎝𝑘₄𝜎₁√𝑡₂,−𝜇₃𝑡₃𝜎₃√𝑡₃+𝜃_1⋅4𝜌^1⋅22𝜇₂√𝑡₁𝜎₂−2ℎ¹−𝜃_3⋅4𝜎₁√𝑡₁−𝜃_1⋅3𝜃_1⋅41−𝜃^21⋅3𝜎2₂√ℎ²𝑡₂−𝜌_1⋅22ℎ₁𝜎₁√𝑡₂+𝜃_1⋅3𝜌^1⋅22𝜇₂√𝑡₁𝜎₂−2ℎ¹𝜃𝜎₁√𝑡₁_1⋅2⎞⎟⎟⎟⎟⎟⎠;=𝜌_1⋅2,𝜃_1⋅3=𝑡₁𝑡₂𝜌_1⋅2,𝜃_1⋅4=𝑡₁𝑡₃𝜌_1⋅3,𝜃_2⋅3=−𝑡₁𝑡₂,𝜃_3⋅4=𝑡₂𝑡₃𝜌_2⋅3⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦,(2.24) where 𝜆=𝜆−=11ifif𝑃𝑃_[⋅][⋅]𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄,𝑡₁,𝑡₂,𝑡₃=𝑃_[up-and-up]𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄,𝑡₁,𝑡₂,𝑡₃,𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄,𝑡₁,𝑡₂,𝑡₃=𝑃_{[down-and-down]}𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄,𝑡₁,𝑡₂,𝑡₃.(2.25) And, for the up-and-down and the down-and-up distributions, one has 𝑃_[.]=𝐻₁Φ,𝐻⎡⎢⎢⎢⎢⎢⎣𝜆₂,𝐾₁−,𝑘𝐾₂₁,𝐾₃,𝐾₄,𝑡₁,𝑡₂,𝑡₃+𝜇₁𝑡₁𝜎₁√𝑡₁,𝜆𝑘²−𝜇₂𝑡₁𝜎₂√𝑡₁,𝜆𝑘³−𝜇₂𝑡₂𝜎₂√𝑡₂,𝜆𝑘⁴−𝜇₃𝑡₃𝜃_1⋅2𝜎₃√𝑡₃;=−𝜌_1⋅2,𝜃_1⋅3=−𝑡₁𝑡₂𝜌_1⋅2,𝜃_1⋅4=−𝑡₁𝑡₃𝜌_1⋅3,𝜃_2⋅3=𝑡₁𝑡₂,𝜃_3⋅4=𝑡₂𝑡₃𝜌_2⋅3−exp2𝜇₂⎤⎥⎥⎥⎥⎥⎦ℎ₂×Φ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣𝜆−𝜎𝑘²²₁+𝜇₁𝑡₁𝜎₁√𝑡₁+𝜌^1⋅22𝜇₂𝑡₁𝜎₂√𝑡₁,𝜆𝑘²+𝜇₂𝑡₁𝜎₂√𝑡₁,𝜆𝑘³−2ℎ₂−𝜇₂𝑡₂𝜆⎛⎜⎜⎜⎝𝑘₄𝜎₂√𝑡₂,−𝜇₃𝑡₃𝜎₃√𝑡₃−𝜃^3⋅4−𝜃_1⋅3𝜃_1⋅41−𝜃^21⋅3𝜎2₂√ℎ²𝑡₂+𝜃_1⋅3𝜌_1⋅22𝜇₂√𝑡₁𝜎₂+𝜃^1⋅4𝜌_1⋅22𝜇₂√𝑡₁𝜎₂𝜃_1⋅2⎞⎟⎟⎟⎠;=−𝜌_1⋅2,𝜃_1⋅3=𝑡₁𝑡₂𝜌_1⋅2,𝜃_1⋅4=−𝑡₁𝑡₃𝜌_1⋅3,𝜃_2⋅3=−𝑡₁𝑡₂,𝜃_3⋅4=𝑡₂𝑡₃𝜌_2⋅3−exp2𝜇₁⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ℎ₁×Φ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣𝜆−𝜎𝑘²¹₁+2ℎ₁+𝜇₁𝑡₁𝜎₁√𝑡₁,𝜆𝑘²−𝜇₂𝑡₁𝜎₂√𝑡₁−𝜌_1⋅22ℎ₁𝜎₁√𝑡₁,𝜆𝑘³−𝜇₂𝑡₂𝜎₂√𝑡₂−𝜌_1⋅22ℎ₁𝜆⎛⎜⎜⎜⎝𝑘₄𝜎₁√𝑡₂,−𝜇₃𝑡₃𝜎₃√𝑡₃−𝜃_1⋅42ℎ₁𝜎₁√𝑡₁−𝜃^3⋅4−𝜃_1⋅3𝜃_1⋅41−𝜃^21⋅3𝜌^1⋅22ℎ₁𝜎₁√𝑡₂−𝜃_1⋅32ℎ₁𝜃_1⋅2𝜎₁√𝑡₁⎞⎟⎟⎟⎠;=−𝜌_1⋅2,𝜃_1⋅3=−𝑡₁𝑡₂𝜌_1⋅2,𝜃_1⋅4=−𝑡₁𝑡₃𝜌_1⋅3,𝜃_2⋅3=𝑡₁𝑡₂,𝜃_3⋅4=𝑡₂𝑡₃𝜌_2⋅3+exp2𝜇₁⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦𝜎²¹−4𝜇²𝜌_1⋅2𝜎₁𝜎₂ℎ¹+2𝜇²ℎ₂×Φ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣𝜆−𝑘₁𝜎²²+2ℎ₁+𝜇₁𝑡₁𝜎₁√𝑡₁+𝜌^1⋅22𝜇₂√𝑡₁𝜎₂,𝜆𝑘²+𝜇₂𝑡₁𝜎₂√𝑡₁−𝜌_1⋅22ℎ₁𝜆𝑘₃𝜎₁√𝑡₁,−2ℎ₂−𝜇₂𝑡₂𝜎₂√𝑡₂+𝜌_1⋅22ℎ₁𝜆⎛⎜⎜⎜⎜⎜⎝𝑘₄𝜎₁√𝑡₂,−𝜇₃𝑡₃𝜎₃√𝑡₃+𝜃_1⋅4𝜌^1⋅22𝜇₂√𝑡₁𝜎₂−2ℎ¹−𝜃_3⋅4𝜎₁√𝑡₁−𝜃_1⋅3𝜃_1⋅41−𝜃^21⋅3𝜎2₂√ℎ²𝑡₂−𝜌_1⋅22ℎ₁𝜎₁√𝑡₂+𝜃_1⋅3𝜌^1⋅22𝜇₂√𝑡₁𝜎₂−2ℎ¹𝜃𝜎₁√𝑡₁_1⋅2⎞⎟⎟⎟⎟⎟⎠;=−𝜌_1⋅2,𝜃_1⋅3=𝑡₁𝑡₂𝜌_1⋅2,𝜃_1⋅4=−𝑡₁𝑡₃𝜌_1⋅3,𝜃_2⋅3=−𝑡₁𝑡₂,𝜃_3⋅4=𝑡₂𝑡₃𝜌_2⋅3⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦,(2.26) where 𝜆𝜆==−11ifif𝑃𝑃_[⋅][⋅]𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄,𝑡₁,𝑡₂,𝑡₃=𝑃_{[down-and-up]}𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄,𝑡₁,𝑡₂,𝑡₃,𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄,𝑡₁,𝑡₂,𝑡₃=𝑃_{[up-and-down]}𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄,𝑡₁,𝑡₂,𝑡₃.(2.27)

Corollary of Proposition 2.2
The four following cumulative distribution functions, that will be useful in Section 3, are deduced from Proposition 2.2:(1)ℙsup_{0≤𝑡≤𝑡1}𝑆₁(𝑡)≤𝐻₁,𝑆₁𝑡₁≤𝐾₁,𝑆₂𝑡₁≤𝐾₂,sup_{𝑡1≤𝑡≤𝑡2}𝑆₂(𝑡)≤𝐻₂,𝑆₂𝑡₂≤𝐾₃,𝑆₃𝑡₃≥𝐾₄=𝑃_[up-and-up]𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,∞−𝑃_[up-and-up]𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄,(2.28)(2)ℙinf_{0≤𝑡≤𝑡1}𝑆₁(𝑡)≥𝐻₁,𝑆₁𝑡₁≥𝐾₁,𝑆₂𝑡₁≥𝐾₂,inf_{𝑡1≤𝑡≤𝑡2}𝑆₂(𝑡)≥𝐻₂,𝑆₂𝑡₂≥𝐾₃,𝑆₃𝑡₃≤𝐾₄=𝑃_{[down-and-down]}𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,−∞−𝑃_{[down-and-down]}𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄,(2.29)(3)ℙinf_{0≤𝑡≤𝑡1}𝑆₁(𝑡)≥𝐻₁,𝑆₁𝑡₁≥𝐾₁,𝑆₂𝑡₁≤𝐾₂,sup_{𝑡1≤𝑡≤𝑡2}𝑆₂(𝑡)≤𝐻₂,𝑆₂𝑡₂≤𝐾₃,𝑆₃𝑡₃≥𝐾₄=𝑃_{[down-and-up]}𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,∞−𝑃_{[down-and-up]}𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄,(2.30)(4)ℙsup_{0≤𝑡≤𝑡1}𝑆₁(𝑡)≤𝐻₁,𝑆₁𝑡₁≤𝐾₁,𝑆₂𝑡₁≥𝐾₂,inf_{𝑡1≤𝑡≤𝑡2}𝑆₂(𝑡)≥𝐻₂,𝑆₂𝑡₂≥𝐾₃,𝑆₃𝑡₃≤𝐾₄=𝑃_{[up-and-down]}𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,−∞−𝑃_{[up-and-down]}𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄.(2.31)

Proof of Proposition 2.2 is provided in the Appendices A, B, and C.

Practical use of Proposition 2.2 and its corollary in engineering applications depend on the accuracy and the efficiency with which the function Φ can be numerically implemented. This question is dealt with in Section 3.

3. Applications and Numerical Implementation

This Section deals with applications of the results provided in Section 1 in financial engineering. Indeed, Propositions 2.1 and 2.2 can be used as the building blocks for the valuation and the risk management of a large number of option contracts. The main class of instruments under consideration will be barrier options. In their standard form, the latter are contracts whereby the holder is entitled (but not obligated) to buy (call option) or sell (put option) an asset at a prespecified future date (the option expiry) at a prespecified price (the strike price), on condition that the asset price has not (knock-out type) or has (knock-in type) crossed a specific upward level (called the up-barrier) or downward level (called the down-barrier) at any given time from valuation date to expiry. In standard option pricing theory, asset prices are modelled as geometric Brownian motions, so that boundary crossing distributions such as those provided in Section 1 are of immediate use in financial engineering.

It must be emphasized that the market for barrier options is huge. They are the most actively traded class of nonstandard options (usually referred to as exotic options). Many variations on the standard barrier option payoff have been designed to match investors’ demand more closely. One of them is the partial-time barrier option, which specifies that the barrier level is monitored during only a fraction of the option lifespan. The basic contracts were studied by Heynen and Kat [11], and more general contracts were valued by Armstrong [12] and Guillaume [13]. Another variation is the step barrier option (Guillaume [6, 9]), whereby the barrier evolves as a step function of time. It is also common to encounter outside barrier options, whose payoff is a function of two asset prices: one of them is compared with the strike price at expiry to determine the moneyness of the option, while the other one is monitored up until expiry to check whether the barrier level has been crossed. Outside barrier options were originally valued by Heynen and Kat [14], and they were further studied by Kwok et al. [15] and Wong and Kwok [16].

These different features: partial time barrier, step barrier, and outside barrier, often combine to allow for increased flexibility. But then, performing analytical valuation becomes more and more involved and practitioners have to turn to numerical methods that are slow and relatively inaccurate. This is when the results provided in Section 1 become valuable. Indeed, they enable to price a large number of barrier option contracts sharing the three above-mentioned innovative features (partial time, step, and outside barrier), that is, options based on two or even three assets that may knock out depending on whether the underlying assets move within a sequence of prespecified ranges of prices over all or part of the option lifespan. More precisely, Proposition 2.1 and its corollary are the building blocks for the valuation of sequential two-asset knock-out calls and puts whose four pay-off structures are defined as follows:(1)(𝐾−𝑆^𝑡(2)₂)⁺𝕀_{{sup10≤𝑡≤𝑡𝑆}^𝑡(1)_{<𝐻1,sup𝑡12≤𝑡≤𝑡𝑆}^𝑡(2)_<𝐻2} for a two-asset up-and-up knock-out put (substitute (𝑆^𝑡(2)₂−𝐾)⁺ for (𝐾−𝑆^𝑡(2)₂)⁺ for an up-and-up knock-out call),(2)±(𝑆^𝑡(2)₂−𝐾)⁺𝕀_{{inf10≤𝑡≤𝑡𝑆}^𝑡(1)_{>𝐻1,inf𝑡12≤𝑡≤𝑡𝑆}^𝑡(2)_>𝐻2} for a two-asset down-and-down knock-out call or put, (3)±(𝑆^𝑡(2)₂−𝐾)⁺𝕀_{{sup10≤𝑡≤𝑡𝑆}^𝑡(1)_{<𝐻1,inf𝑡12≤𝑡≤𝑡𝑆}^𝑡(2)_>𝐻2} for a two-asset up-and-down knock-out call or put,(4)±(𝐾−𝑆^𝑡(2)₂)⁺𝕀_{{inf10≤𝑡≤𝑡𝑆}^𝑡(1)_{>𝐻1,sup𝑡12≤𝑡≤𝑡𝑆}^𝑡(2)_<𝐻2} for a two-asset down-and-up knock-out put or call,

where:(i)𝕀_{} is the indicator function, (ii){𝑆₁(𝑡),𝑡≥0} and {𝑆₂(𝑡),𝑡≥0} are the price processes of two risky assets,(iii)𝐾 is the strike price and 𝑡₂ is the option expiry, (iv)𝐻₁ is a knock-out barrier monitored with respect to {𝑆₁(𝑡),𝑡≥0} over the time interval [𝑡₀,𝑡₁],(v)𝐻₂ is a knock-out barrier monitored with respect to {𝑆₂(𝑡),𝑡≥0} over the time interval [𝑡₁,𝑡₂].

The no-arbitrage price of these eight different option contracts is now provided by the following proposition.

Proposition 3.1. Let two geometric Brownian motions {𝑆₁(𝑡),𝑡≥0} and {𝑆₂(𝑡),𝑡≥0} be the price processes of two risky assets with constant volatilities 𝜎₁ and 𝜎₂, respectively, and constant correlation coefficient 𝜌_1⋅2. Let 𝛿₁ and 𝛿₂ be two constant pay-out rates on assets 𝑆₁ and 𝑆₂, respectively, and 𝑟 be the constant risk-free interest rate. Let 𝜐₁=𝑟−𝛿₁+𝜎₁𝜎₂𝜌_1⋅2 and 𝜐₂=𝑟−𝛿₂+𝜎²².
The value, at time 𝑡₀=0, of a two-asset sequential knock-out option with strike price 𝐾, with knock-out barriers 𝐻₁ and 𝐻₂ monitored with respect to 𝑆₁ over the time interval [𝑡₀,𝑡₁] and 𝑆₂ over the time interval [𝑡₁,𝑡₂], respectively, and with maturity 𝑡₂, 𝑡₂≥𝑡₁≥𝑡₀, is given by 𝑉𝑆₁𝑡₀,𝑆₂𝑡₀,𝑡₁,𝑡₂,𝐾,𝐻₁,𝐻₂=±𝐾exp−𝑟𝑡₂𝑃₁−𝑆₂𝑡₀exp−𝛿₂𝑡₂𝑃₂,(3.1) where the following specifications hold. (1) If the option is an up-and-up put, ±=1, and 𝑃₁ is obtained by substituting 𝛼₁=𝑟−𝛿₁ and 𝛼₂=𝑟−𝛿₂ in 𝑃_[up-and-up](𝐻₁,𝐻₂,𝐻₁,𝐻₂,𝐾,𝑡₁,𝑡₂) in Proposition 2.1, and 𝑃₂ is obtained by substituting 𝛼₁=𝜐₁ and 𝛼₂=𝜐₂ in 𝑃_[up-and-up](𝐻₁,𝐻₂,𝐻₁,𝐻₂,𝐾,𝑡₁,𝑡₂) in Proposition 2.1. (2) If the option is a down-and-down call,±=−1, and 𝑃₁ is obtained by substituting 𝛼₁=𝑟−𝛿₁ and 𝛼₂=𝑟−𝛿₂ in 𝑃_{[down-and-down]}(𝐻₁,𝐻₂,𝐻₁,𝐻₂,𝐾,𝑡₁,𝑡₂) in Proposition 2.1, and 𝑃₂ is obtained by substituting 𝛼₁=𝜐₁ and 𝛼₂=𝜐₂ in 𝑃_{[down-and-down]}(𝐻₁,𝐻₂,𝐻₁,𝐻₂,𝐾,𝑡₁,𝑡₂) in Proposition 2.1. (3) If the option is an up-and-down call,±=−1, and 𝑃₁ is obtained by substituting 𝛼₁=𝑟−𝛿₁ and 𝛼₂=𝑟−𝛿₂ in 𝑃_{[up-and-down]}(𝐻₁,𝐻₂,𝐻₁,𝐻₂,𝐾,𝑡₁,𝑡₂) in Proposition 2.1 and 𝑃₂ is obtained by substituting 𝛼₁=𝜐₁ and 𝛼₂=𝜐₂ in 𝑃_{[up-and-down]}(𝐻₁,𝐻₂,𝐻₁,𝐻₂,𝐾,𝑡₁,𝑡₂) in Proposition 2.1. (4) If the option is a down-and-up put, ±=1, and 𝑃₁ is obtained by substituting 𝛼₁=𝑟−𝛿₁ and 𝛼₂=𝑟−𝛿₂ in 𝑃_{[down-and-up]}(𝐻₁,𝐻₂,𝐻₁,𝐻₂,𝐾,𝑡₁,𝑡₂) in Proposition 2.1 and 𝑃₂ is obtained by substituting 𝛼₁=𝜐₁ and 𝛼₂=𝜐₂ in 𝑃_{[down-and-up]}(𝐻₁,𝐻₂,𝐻₁,𝐻₂,𝐾,𝑡₁,𝑡₂) in Proposition 2.1. By using Corollary of Proposition 2.1, one can value up-and-up and down-and-up call options, as well as down-and-down and up-and-down put options.

A sketch of proof of Proposition 3.1 is provided in the Appendices A, B, and C.

Just as Proposition 2.1 enables to value two-asset knock-out calls and puts, Proposition 2.2 can be used to value three-asset knock-out calls and puts. More specifically, if {𝑆₃(𝑡),𝑡≥0} is the price process of a third risky asset and 𝑡₃ is the new option expiry, with 𝑡₃≥𝑡₂≥𝑡₁≥𝑡₀, then the four following pay-off structures can be analytically tackled:(1)±(𝐾−𝑆^𝑡(3)₃)⁺𝕀_{{sup10≤𝑡≤𝑡𝑆}^𝑡(1)_{<𝐻1,sup𝑡12≤𝑡≤𝑡𝑆}^𝑡(2)_<𝐻2} for a three-asset up-and-up knock-out put or call,(2)±(𝑆^𝑡(3)₃−𝐾)⁺𝕀_{{inf10≤𝑡≤𝑡𝑆}^𝑡(1)_{>𝐻1,inf𝑡12≤𝑡≤𝑡𝑆}^𝑡(2)_>𝐻2} for a three-asset down-and-down knock-out call or put, (3)±(𝑆^𝑡(3)₃−𝐾)⁺𝕀_{{sup10≤𝑡≤𝑡𝑆}^𝑡(1)_{<𝐻1,inf𝑡12≤𝑡≤𝑡𝑆}^𝑡(2)_>𝐻2} for a three-asset up-and-down knock-out call or put,(4)±(𝐾−𝑆^𝑡(3)₃)⁺𝕀_{{inf10≤𝑡≤𝑡𝑆}^𝑡(1)_{>𝐻1,sup𝑡12≤𝑡≤𝑡𝑆}^𝑡(2)_<𝐻2} for a three-asset down-and-up knock-out put or call.

The no-arbitrage prices of the eight types of option under consideration are now provided by the following proposition.

Proposition 3.2. Let three geometric Brownian motions {𝑆₁(𝑡),𝑡≥0}, {𝑆₂(𝑡),𝑡≥0}, and {𝑆₃(𝑡),𝑡≥0} be the price processes of three risky assets with constant volatilities 𝜎₁, 𝜎₂, and 𝜎₃, respectively, and constant pairwise correlation coefficients 𝜌_𝑎⋅𝑏, (𝑎,𝑏)∈{1,2,3}, at any fixed time 𝑡≥0. Let 𝛿₁, 𝛿₂, and 𝛿₃ be three constant pay-out rates on assets 𝑆₁, 𝑆₂, and 𝑆₃, respectively, and 𝑟 be the constant risk-free interest rate. Let 𝜐₁=𝑟−𝛿₁+𝜎₁𝜎₃𝜌_1⋅3, 𝜐₂=𝑟−𝛿₂+𝜎₂𝜎₃𝜌_2⋅3, and 𝜐₃=𝑟−𝛿₃+𝜎²³.
The value, at time 𝑡₀=0, of a three-asset sequential knock-out option with strike price 𝐾, with knock-out barriers 𝐻₁ and 𝐻₂ monitored with respect to 𝑆₁ over the time interval [𝑡₀,𝑡₁] and 𝑆₂ over the time interval [𝑡₁,𝑡₂], respectively, with maturity 𝑡₃, 𝑡₃≥𝑡₂≥𝑡₁≥𝑡₀, is given by 𝑉𝑆₁𝑡₀,𝑆₂𝑡₀,𝑆₃𝑡₀,𝑡₁,𝑡₂,𝑡₃,𝐾,𝐻₁,𝐻₂=±𝐾exp−𝑟𝑡₃𝑃₁−𝑆₃𝑡₀exp−𝛿₃𝑡₃𝑃₂,(3.2) where the following specifications hold. (1) If the option is an up-and-up put, ±=1, and 𝑃₁ is obtained by substituting 𝛼₁=𝑟−𝛿₁, 𝛼₂=𝑟−𝛿₂ and 𝛼₃=𝑟−𝛿₃ in 𝑃_[up-and-up](𝐻₁,𝐻₂,𝐻₁,𝐻₂,𝐻₂,𝐾,𝑡₁,𝑡₂,𝑡₃) in Proposition 2.2, and 𝑃₂ is obtained by substituting 𝛼₁=𝜐₁, 𝛼₂=𝜐₂ and 𝛼₃=𝜐₃ in 𝑃_[up-and-up](𝐻₁,𝐻₂,𝐻₁,𝐻₂,𝐻₂,𝐾,𝑡₁,𝑡₂,𝑡₃) in Proposition 2.2. (2) If the option is a down-and-down call,±=−1, and 𝑃₁ is obtained by substituting 𝛼₁=𝑟−𝛿₁, 𝛼₂=𝑟−𝛿₂ and 𝛼₃=𝑟−𝛿₃ in 𝑃_{[down-and-down]}(𝐻₁,𝐻₂,𝐻₁,𝐻₂,𝐻₂,𝐾,𝑡₁,𝑡₂,𝑡₃) in Proposition 2.2, and 𝑃₂ is obtained by substituting 𝛼₁=𝜐₁, 𝛼₂=𝜐₂ and 𝛼₃=𝜐₃ in 𝑃_{[down-and-down]}(𝐻₁,𝐻₂,𝐻₁,𝐻₂,𝐻₂,𝐾,𝑡₁,𝑡₂,𝑡₃) in Proposition 2.2. (3) If the option is an up-and-down call, ±=−1, and 𝑃₁ is obtained by substituting 𝛼₁=𝑟−𝛿₁, 𝛼₂=𝑟−𝛿₂ and 𝛼₃=𝑟−𝛿₃ in 𝑃_{[up-and-down]}(𝐻₁,𝐻₂,𝐻₁,𝐻₂,𝐻₂,𝐾,𝑡₁,𝑡₂,𝑡₃) in Proposition 2.2 and 𝑃₂ is obtained by substituting 𝛼₁=𝜐₁, 𝛼₂=𝜐₂ and 𝛼₃=𝜐₃ in 𝑃_{[up-and-down]}(𝐻₁,𝐻₂,𝐻₁,𝐻₂,𝐻₂,𝐾,𝑡₁,𝑡₂,𝑡₃) in Proposition 2.2. (4) If the option is a down-and-up put,±=1, and 𝑃₁ is obtained by substituting 𝛼₁=𝑟−𝛿₁, 𝛼₂=𝑟−𝛿₂ and 𝛼₃=𝑟−𝛿₃ in 𝑃_{[down-and-up]}(𝐻₁,𝐻₂,𝐻₁,𝐻₂,𝐻₂,𝐾,𝑡₁,𝑡₂,𝑡₃) in Proposition 2.2 and 𝑃₂ is obtained by substituting 𝛼₁=𝜐₁, 𝛼₂=𝜐₂ and 𝛼₃=𝜐₃ in 𝑃_{[down-and-up]}(𝐻₁,𝐻₂,𝐻₁,𝐻₂,𝐻₂,𝐾,𝑡₁,𝑡₂,𝑡₃) in Proposition 2.2.By using Corollary of Proposition 2.2, one can value up-and-up and down-and-up call options, as well as down-and-down and up-and-down put options.

A sketch of proof of Proposition 3.2 is provided in the Appendices A, B, and C.

It is easy to notice that Proposition 3.2 nests Proposition 3.1 under the following three conditions:(1)𝜐₁=𝑟−𝛿₁+𝜎₁𝜎₂𝜌_1⋅2 and 𝜐₂=𝑟−𝛿₂+𝜎²²,(2)take 𝑃_[⋅](𝐻₁,𝐻₂,𝐻₁,𝐻₂,𝐾,𝑈,𝑡₁,𝑡₂,𝑡₃) with 𝑈→∞ if you are dealing with a put and with 𝑈→−∞ if you are dealing with a call,(3)compute ±[𝐾exp(−𝑟𝑡₂)𝑃₁−𝑆₂(𝑡₀)exp(−𝛿₂𝑡₂)𝑃₂].

It is useful, though, to state Proposition 3.1, as it uses more standard functions that can be numerically computed with known high accuracy and thus can serve as a benchmark to test the reliability of the numerical implementation of Proposition 3.2.

It can also be noticed that Proposition 3.2 enables to value sequential two-asset knock-out options having the property that the moneyness of the option at expiry is tested with respect to asset 𝑆₁ instead of asset 𝑆₂ as in Proposition 3.1. This is achieved by implementing Proposition 3.2 with the following specifications: (1)𝜐₁=𝑟−𝛿₁+𝜎²¹ and 𝜐₂=𝑟−𝛿₂+𝜎₁𝜎₂𝜌_1⋅2,(2)𝑡₃→𝑡₂, 𝜎₃→𝜎₁, 𝜌_1⋅3→1 and 𝜌_2⋅3→𝜌_1⋅2, (3)compute ±[𝐾exp(−𝑟𝑡₂)𝑃₁−𝑆₁(𝑡₀)exp(−𝛿₁𝑡₁)𝑃₂].

Similarly, Proposition 3.2 can be used to value sequential two-asset knock-out options with expiry 𝑡₃, so that, overall, this formula alone suffices to compute 40 different types of sequential multiasset knock-out options.

Proposition 3.2 implies the numerical implementation of Proposition 2.2. The function Φ in Proposition 2.2 involves three-dimensional numerical integration. Given the fact that the integrand is very smooth, a classical Gauss-Legendre quadrature can be used, along with a simple adaptive rule increasing the number of points in case the required accuracy is not reached after one iteration, starting with an 8-point quadrature. To test the efficiency and the accuracy of this numerical implementation, four series of tests have been implemented. In the first three categories of tests, the option parameters are specified so that the option can be valued by a known analytical formula involving numerical integration of lower dimension. More specifically, tests of category 1 reduce a three-asset up-and-up knock-out put to a two-asset up-and-up knock-out put by letting 𝐾₄→∞ in 𝑃_[up-and-up](𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄,𝑡₁,𝑡₂,𝑡₃); the benchmark is then Proposition 3.1, that involves the computation of trivariate normal cumulative distribution functions. The latter can be obtained with high precision by Genz’s algorithm (Genz [10]), as mentioned earlier.

Tests of category 2 reduce the three-asset up-and-up knock-out put to an outside knock-out put by letting 𝐾₁=𝐻₁, 𝐻₂→∞, 𝐾₃→∞, 𝐾₄→∞, and 𝑡₃→𝑡₂→𝑡₁ in 𝑃_[up-and-up](𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄,𝑡₁,𝑡₂,𝑡₃); the benchmark, then, is the formula for an outside barrier option as given by Heynen and Kat [14], that involves the computation of bivariate normal cumulative distribution functions. Again, the latter can be obtained with high precision using Genz [10].

Tests of category 3 reduce the three-asset up-and-up knock-out put to a single asset up-and-out put by letting 𝐻₂→∞, 𝐾₂=𝐾₁, 𝐾₃→∞, 𝐾₄→∞, 𝑡₃→𝑡₂→𝑡₁, 𝜎₂=𝜎₁ and 𝜌_1⋅2→1 in 𝑃_[up-and-up](𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄,𝑡₁,𝑡₂,𝑡₃); the benchmark, then, is the formula for a basic barrier option as given, for example, by Haug [17], that involves the computation of univariate normal cumulative distribution functions.

In the last set of tests, the three-asset up-and-up knock-out put cannot be valued as a simpler contract, so that no reduction of the dimension of numerical integration is possible. Our goal, then, is to find a robust numerical valuation method and study its convergence pattern with the results obtained by implementing Proposition 3.2. Neither using trees nor numerically solving a partial differential equation stands as the most suitable methods in this case, in view of the double difficulty of having three correlated diffusion processes and path-dependent features in the option payoff. On another hand, crude Monte Carlo simulation would be very poor both in terms of accuracy and efficiency, due to the required time discretization. That is why the most suitable numerical method in this case is probably a Brownian Bridge Monte Carlo approximation (Glasserman [18]), as it involves only 6 random draws of the asset prices for each simulation and avoids a discretization bias. This is indeed the method commonly used by market practitioners when they are faced with a pricing problem involving both multiasset and path-dependent features.

For each category of tests, a total of 10,000 different option prices have been computed with randomly drawn volatility, correlation, and expiry parameters. Let us now sum up the main findings of these numerical experiments. First, the average computational time when applying Proposition 3.2 on an ordinary personal computer (P7350 2 GHz) is 0.6 seconds, which is fast. The longest recorded computational time was 1.4 second. Speed could obviously be greatly increased by involving more hardware resources. In absolute terms, Proposition 3.2 can thus be considered as an efficient way of computing option prices. In relative terms, Proposition 3.2 is even extremely efficient. Indeed, 10,000,000 simulations are required to obtain a modest 10⁻⁴ digits convergence of the Brownian Bridge Monte Carlo approximation with Proposition 3.2, which amounts to a computational time of over 4 minutes! Such slowness becomes a serious issue in real-time trading or when large portfolios of options are valued. As far as accuracy is concerned, the results obtained using Proposition 3.2 always matched to at least 10⁻⁸ digits the analytical benchmarks in almost all the numerical experiments conducted. The very rare exceptions were observed when the absolute value of one of the correlation coefficients was above 99.7%. It must be noticed that in this couple of cases, the correlation matrix was found to be almost singular. Besides, risky assets with such levels of correlation are never encountered in the markets, so that it does not seem worthwhile to us to search for more complicated integration rules that would be better suited to these negligible subregions of integration. For all practical purposes, the numerical accuracy of our simple implementation of Proposition 3.2 can thus be deemed as very satisfactory. The results did not vary significantly between tests of categories 1, 2, and 3, which means that an increase in the number of option parameters that are given limit values does not entail instability in the quadrature rule.

To close this section, it is worth mentioning that the results provided in Section 1 can also apply to other popular derivatives contracts. Among them is the look-back option, which allows investors to sell at the highest and buy at the lowest over a time interval, thus optimizing market exit and market entry on a given asset. In its basic form, a fixed strike look-back call pays off (sup_{0≤𝑡≤𝑇}𝑆(𝑡)−𝐾)⁺ at expiry 𝑇. Many variations on this standard structure have been designed. One of them is the double look-back call, written on two assets 𝑆₁ and 𝑆₂, which provides investors with the following payoff at expiry 𝑇: (𝑎₁sup_{0≤𝑡≤𝑇}𝑆₁(𝑡)−𝑎₂sup_{0≤𝑡≤𝑇}𝑆₂(𝑡)−𝐾)⁺, where 𝑎₁ and 𝑎₂ are two positive weights. The valuation of this contract has been studied by He et al. [8]. A sequential double look-back call is a variation on the double look-back call. Its payoff at expiry 𝑡₂ writes (𝑎₁sup_{0≤𝑡≤𝑡1}𝑆₁(𝑡)−𝑎₂sup_{𝑡1≤𝑡≤𝑡2}𝑆₂(𝑡)−𝐾)⁺. Drawing on Proposition 2.1, the no-arbitrage price of this option is equal toexp−𝑟𝑡₂^∞ℎ₁₌₀^∞ℎ_2=−∞𝑎₁𝑆₁𝑡₀expℎ₁−𝑎₂𝑆₂𝑡₀expℎ₂−𝐾⁺𝜕²𝑃_[up-and-up]𝐻₁,𝐻₂,𝐻₁,𝐻₂,𝐻₂𝜕ℎ₁𝜕ℎ₂𝑑ℎ₂𝑑ℎ₁,(3.3) where 𝛼₁=𝑟−𝛿₁ and 𝛼₂=𝑟−𝛿₂ in Proposition 2.1.

The partial derivative in (3.3) can be determined by analytical differentiation of 𝑃_[up-and-up](𝐻₁,𝐻₂,𝐻₁,𝐻₂,𝐻₂), but this is relatively cumbersome. Given the simplicity of the computation of Proposition 2.1, it is faster and easier to compute those partial derivatives numerically. One can notice that the valuation of sequential double look-back options is only semianalytical since it still requires numerical double integration. For a given degree of the quadrature rule, it will not be as accurate as the valuation of two-asset and three-asset sequential knock-out options because of the inaccuracies arising from the numerical differentiation of Proposition 2.1.

4. Conclusion

This paper provides new explicit solutions for some sequential boundary crossing problems in a multidimensional geometric Brownian motion framework. Among their various possible applications, they form the basis for accurate and fast valuation of a large number of option contracts traded in the markets. The methods presented here could be used to solve higher-dimensional problems, especially the introduction of a third boundary monitored with respect to the third process in the model. However, the resulting formulae would become quite cumbersome and, more importantly, the quality of the numerical integration involved should have to be compared with a Brownian Bridge Monte Carlo approximation. Another possible extension would be to study the analytical tractability of introducing curved boundaries, but it should be pointed out that this form of boundary is rarely encountered in the options markets.

Appendices

A. Proof of Proposition 2.1

Proof is outlined for 𝑃_[up-and-up](𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝑡₁,𝑡₂) only. The other distributions are obtained in a similar manner.

The first step is to take the logarithms of 𝑆₁(𝑡₁), 𝑆₂(𝑡₁), and 𝑆₂(𝑡₂), in order to move from lognormal random variables to normal random variables. The sought cumulative distribution becomesℙsup_{0≤𝑡≤𝑡1}𝑋₁(𝑡)≤ℎ₁,𝑋₁𝑡₁≤𝑘₁,𝑋₂𝑡₁≤𝑘₂,sup_{𝑡1≤𝑡≤𝑡2}𝑋₂(𝑡)≤ℎ₂,𝑋₂𝑡₂≤𝑘₃,(A.1) where𝑋₁(𝑡)=ln𝑆¹(𝑡)𝑆₁𝑡₀,𝑋²(𝑡)=ln𝑆²(𝑡)𝑆₂𝑡₀.(A.2) For a fixed 𝑡≥0, 𝑋₁(𝑡) and 𝑋₂(𝑡) are absolutely continuous random variables, hence one can express (A.1) as the following triple integral:^𝑘1_𝑥1=−∞^𝑘2_𝑥2=−∞^𝑘3_𝑥3=−∞ℙ⎛⎜⎜⎜⎝sup^{0≤𝑡≤𝑡1}sup_{𝑡1≤𝑡≤𝑡2}𝑋₁(𝑡)≤ℎ₁,𝑋₁𝑡₁∈𝑑𝑥₁,𝑋₂𝑡₁∈𝑑𝑥₂,𝑋₂(𝑡)≤ℎ₂,𝑋₂𝑡₂∈𝑑𝑥₃⎞⎟⎟⎟⎠𝑑𝑥₃𝑑𝑥₂𝑑𝑥₁.(A.3) By conditioning and using the Markov property of processes 𝑋₁(𝑡) and 𝑋₂(𝑡), (A.3) becomes^𝑘1_𝑥1=−∞^𝑘2_𝑥2=−∞^𝑘3_𝑥3=−∞𝑃₁𝑥₁𝑃₂𝑥₁,𝑥₂𝑃₃𝑥₂,𝑥₃𝑑𝑥₃𝑑𝑥₂𝑑𝑥₁,(A.4) where𝑃₁𝑥₁=ℙsup_{0≤𝑡≤𝑡1}𝑋₁(𝑡)≤ℎ₁,𝑋₁𝑡₁∈𝑑𝑥₁𝑃₂,𝑃₃𝑥₁,𝑥₂=ℙ𝑋₂𝑡₁∈𝑑𝑥₂∣𝑋₁𝑡₁∈𝑑𝑥₁,𝑥₂,𝑥₃=ℙsup_{𝑡1≤𝑡≤𝑡2}𝑋₂(𝑡)≤ℎ₂,𝑋₂𝑡₂∈𝑑𝑥₃∣𝑋₂𝑡₁∈𝑑𝑥₂.(A.5) The functions 𝑃₁(𝑥₁) and 𝑃₂(𝑥₁,𝑥₂) are classical results. To obtain the function 𝑃₃(𝑥₂,𝑥₃), the following formula is used (Guillaume [6]):ℙ𝑋₂𝑡₁≤𝑥₂,sup_{𝑡1≤𝑡≤𝑡2}𝑋₂(𝑡)≤ℎ₂,𝑋₂𝑡₂≤𝑥₃=𝑁₂𝑥²−𝜇₂𝑡₁𝜎₂√𝑡₁,𝑥³−𝜇₂𝑡₂𝜎₂√𝑡₂;𝑡₁𝑡₂−exp2𝜎𝜇²²²ℎ₂𝑁²𝑥²+𝜇₂𝑡₁𝜎₂√𝑡₁,𝑥³−2ℎ₂−𝜇₂𝑡₂𝜎₂√𝑡₂;−𝑡₁𝑡₂,(A.6) where 𝑁₂[⋅,⋅;𝜌] is the bivariate standard normal cumulative distribution function with correlation coefficient 𝜌. Differentiating (A.6) twice with respect to variables 𝑥₂ and 𝑥₃, and applying the definition of conditional probability yields𝑃₃𝑥₂,𝑥₃=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝exp−(1/2)𝑥³−𝑥₂−𝜇₂𝑡₂−𝑡₁/𝜎₂√𝑡₂−𝑡₁²−exp2𝜇₂𝜎₂2𝜋𝑡₂−𝑡₁/𝜎²²ℎ₂−𝑥₂−(1/2)𝑥₃+𝑥₂−2ℎ₂−𝜇₂𝑡₂−𝑡₁/𝜎₂√𝑡₂−𝑡₁²𝜎₂2𝜋𝑡₂−𝑡₁⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠.(A.7) Then, the most cumbersome part of the proof consists in performing the necessary calculations to solve (A.4) in closed form and find the sum of trivariate standard normal cumulative distribution functions given by Proposition 2.1.

B. Proof of Proposition 2.2

Proof is outlined for 𝑃_[up-and-up](𝐻₁,𝐻₂,𝐾₁,𝐾₂,𝐾₃,𝐾₄,𝑡₁,𝑡₂,𝑡₃) only. The other distributions are obtained in a similar manner.

The following lemma is instrumental in proving Proposition 2.2.

Lemma B.1. Let 𝑋₁, 𝑋₂, and 𝑋₃ be three standard normal random variables defined on a probability space with measure ℙ. Denote by 𝜃_𝑎⋅𝑏 the correlation coefficient between 𝑋_𝑎 and 𝑋_𝑏, (𝑎,𝑏)∈{1,2,3}. Then, the joint density function of (𝑋₁,𝑋₂,𝑋₃) is given by ℙ𝑋=₁∈ex𝑑p𝑥₁,−𝑋𝑥₂₂₁∈𝑑𝑥₂,𝑋₃∈𝑑𝑥₃2−2𝜙1^22∣1𝑥₂−𝜃_1⋅2𝑥₁²−12𝜙^23∣1⋅2𝑥³−𝜃_1⋅3𝑥₁−𝜃_2⋅3∣1𝑥²−𝜃_1⋅2𝑥₁×(2𝜋)_3/2𝜙_2∣1²𝜙_2∣1𝜙_3∣1⋅2⁻¹,(B.1) where 𝜙_2|1=1−𝜃^21⋅2 is the standard deviation of 𝑋₂ conditional on 𝑋₁, 𝜃_2⋅3|1=(𝜃_2⋅3−𝜃_1⋅2𝜃_1⋅3)/𝜙_2|1 is the partial correlation between 𝑋₂ and 𝑋₃ conditional on 𝑋₁, and 𝜙_3|1⋅2=1−𝜃^21⋅3−𝜃^22⋅3|1 is the standard deviation of 𝑋₃ conditional on 𝑋₁ and on 𝑋₂.

Proof of Lemma B.1 is simple. It is known that 𝑋₂ is equal in law to 𝜃_1⋅2𝑋₁+𝜙_2|1𝑌₂, where 𝑌₂ is a standard normal random variable uncorrelated with 𝑋₁. By applying Cholesky’s decomposition to the correlation matrix associated with (𝑋₁,𝑋₂,𝑋₃), one also finds that 𝑋₃ admits the following representation:𝑋₃=𝜃_1⋅3𝑋₁+𝜃_2⋅3∣1𝑌₂+𝜙_3∣1⋅2𝑌₃⟷𝑋₃=𝜃_1⋅3𝑋₁+𝜃_2⋅3∣1𝑋₂−𝜃_1⋅2𝑋₁𝜙_2∣1+𝜙_3∣1⋅2𝑌₃,(B.2) where 𝑌₃ is a standard normal random variable uncorrelated with 𝑋₁ and 𝑌₂.

Thus, given 𝑋₁, 𝑋₂, and 𝑋₃ is normally distributed with mean 𝜃_1⋅3𝑋₁+𝜃_2⋅3|1((𝑋₂−𝜃_1⋅2𝑋₁)/𝜙_2|1) and standard deviation 𝜙_3|1⋅2. Hence, the density function of 𝑋₃ conditional on 𝑋₁ and 𝑋₂ writesℙ𝑋=₃∈ex𝑑p𝑥₃−∣𝑋₁1∈𝑑𝑥₁,𝑋₂∈𝑑𝑥₂2𝜙^23∣1⋅2𝑥³−𝜃_1⋅3𝑥₁−𝜃_2⋅3∣1𝑥²−𝜃_1⋅2𝑥₁𝜙_2∣1²𝜙^3∣1⋅2√2𝜋⁻¹,(B.3) which suffices to prove Lemma B.1.

Next, we proceed to outline proof of Proposition 2.2. By taking the logarithms of the asset prices at the times 𝑡₁, 𝑡₂, 𝑡₃, and integrating with respect to these variables, the sought cumulative distribution function becomes ^𝑘1_𝑥1=−∞^𝑘2_𝑥2=−∞^𝑘3_𝑥3=−∞^𝑘4_𝑥4=−∞ℙ⎛⎜⎜⎜⎝sup^{0≤𝑡≤𝑡1}sup_{𝑡1≤𝑡≤𝑡2}𝑋₁(𝑡)≤ℎ₁,𝑋₁𝑡₁∈𝑑𝑥₁,𝑋₂𝑡₁∈𝑑𝑥₂,𝑋₂(𝑡)≤ℎ₂,𝑋₂𝑡₂∈𝑑𝑥₃,𝑋₃𝑡₃∈𝑑𝑥₄⎞⎟⎟⎟⎠𝑑𝑥₄𝑑𝑥₃𝑑𝑥₂𝑑𝑥₁,(B.4) where 𝑋₃(𝑡)=ln(𝑆₃(𝑡)/𝑆₃(𝑡₀)) and the other terms are defined in the proof of Proposition 2.1. Then, by conditioning and using the Markov property of the processes 𝑋₁(𝑡), 𝑋₂(𝑡) and 𝑋₃(𝑡), (B.4) becomes^𝑘1_𝑥1=−∞^𝑘2_𝑥2=−∞^𝑘3_𝑥3=−∞^𝑘4_𝑥4=−∞𝑃₁𝑥₁𝑃₂𝑥₁,𝑥₂𝑃₃𝑥₂,𝑥₃𝑃₄𝑥₁,𝑥₃,𝑥₄𝑑𝑥₄𝑑𝑥₃𝑑𝑥₂𝑑𝑥₁,(B.5) where𝑃₁𝑥₁=ℙsup_{0≤𝑡≤𝑡1}𝑋₁(𝑡)≤ℎ₁,𝑋₁𝑡₁∈𝑑𝑥₁𝑃₂,𝑃₃𝑥₁,𝑥₂=ℙ𝑋₂𝑡₁∈𝑑𝑥₂∣𝑋₁𝑡₁∈𝑑𝑥₁,𝑥₂,𝑥₃=ℙsup_{𝑡1≤𝑡≤𝑡2}𝑋₂(𝑡)≤ℎ₂,𝑋₂𝑡₂∈𝑑𝑥₃∣𝑋₂𝑡₁∈𝑑𝑥₂𝑃₄,𝑥₁,𝑥₃,𝑥₄=ℙ𝑋₃𝑡₃∈𝑑𝑥₄∣𝑋₂𝑡₂∈𝑑𝑥₃,𝑋₁𝑡₁∈𝑑𝑥₁.(B.6) The functions 𝑃₁(𝑥₁), 𝑃₂(𝑥₁,𝑥₂), and 𝑃₃(𝑥₂,𝑥₃) are dealt with in the proof of Proposition 2.1. The function 𝑃₄(𝑥₁,𝑥₃,𝑥₄) is now derived. Since the processes {𝑆₁(𝑡),𝑡≥0}, {𝑆₂(𝑡),𝑡≥0} and {𝑆₃(𝑡),𝑡≥0} are geometric Brownian motions with constant drift and diffusion coefficients, the random variables 𝑋₁(𝑡₁), 𝑋₂(𝑡₂), and 𝑋₃(𝑡₃) are normally distributed. Elementary calculations yield their cross-correlations. Applying Cholesky’s decomposition to 𝑋₂(𝑡₂) and 𝑋₃(𝑡₃), one can see that any linear combination 𝑎₁𝑋₁(𝑡₁)+𝑎₂𝑋₂(𝑡₂)+𝑎₃𝑋₃(𝑡₃) for any set of real numbers {𝑎₁,𝑎₂,𝑎₃} is a sum of three uncorrelated normal random variables; hence, it is normally distributed. Therefore, using Fréchet’s characterization of multivariate normality (Fréchet [19]), the joint density function 𝑃₄(𝑥₁,𝑥₃,𝑥₄) is trivariate normal, and Lemma B.1 yields 𝑃₄𝑥₁,𝑥₃,𝑥₄=exp⎛⎜⎜⎝−12𝜙^24∣1⋅3𝜑³𝑥₄−𝜃_1⋅4𝜑₁𝑥₁−𝜃_3⋅4∣1𝜑²𝑥₃−𝜃_1⋅3𝜑₁𝑥₁𝜙_3∣1²×𝜙_4∣1⋅3⎞⎟⎟⎠𝜙_3∣1𝜎₁𝜎₂𝜎₃𝑡₁𝑡₂𝑡₃8𝜋³⁻¹,(B.7) where the notations are defined as follows:𝜑₁𝑥₁=𝑥¹−𝜇₁𝑡₁𝜎₁√𝑡₁,𝜑₂𝑥₃=𝑥³−𝜇₂𝑡₂𝜎₂√𝑡₂,𝜑₃𝑥₄=𝑥⁴−𝜇₃𝑡₃𝜎₃√𝑡₃𝜃_1⋅3,=𝑡₁𝑡₂𝜌_1⋅2,𝜃_1⋅4=𝑡₁𝑡₃𝜌_1⋅3,𝜃_3⋅4=𝑡₂𝑡₃𝜙_3∣1𝜌_2⋅3,=1−𝜃^21⋅3,𝜃_3⋅4∣1=𝜃^3⋅4−𝜃_1⋅3𝜃_1⋅4𝜙_3∣1,𝜙_4∣1⋅3=1−𝜃^21⋅4−𝜃^23⋅4∣1.(B.8) Then, the most cumbersome part of the proof consists in performing the necessary calculations to solve (B.5) in closed form and find the sum of functions given by Proposition 2.2.

C. Sketch of Proof of Propositions 3.1 and 3.2

Once endowed with Proposition 2.1, proof of Proposition 3.1 is straightforward, as an application of option pricing theory (Harrison and Pliska [20]) and Girsanov’s theorem in two dimensions. The substitutions 𝛼₁=𝜐₁ and 𝛼₂=𝜐₂ in Proposition 2.1 derive from the introduction of a new measure, which is equivalent to the risk-neutral measure, under which the market numeraire is asset 𝑆₂. If we denote by {𝐵₁(𝑡),𝑡≥0} and {𝐵₂(𝑡),𝑡≥0} the two correlated Brownian motions driving {𝑆₁(𝑡),𝑡≥0} and {𝑆₂(𝑡),𝑡≥0} under the risk-neutral measure, then the processes {𝐵₁(𝑡)−𝜎₂𝜌_1⋅2𝑡,𝑡≥0} and {𝐵₂(𝑡)−𝜎₂1−𝜌^21⋅2𝑡,𝑡≥0} are independent standard Brownian motions under the new measure.

Similarly, once endowed with Proposition 2.2, proof of Proposition 3.2 is straightforward, as an application of option pricing theory and Girsanov’s theorem in three dimensions. The substitutions 𝛼₁=𝜐₁, 𝛼₂=𝜐₂, and 𝛼₃=𝜐₃ in Proposition 2.2 derive from the introduction of a new measure, which is equivalent to the risk-neutral measure, under which the market numeraire is asset 𝑆₃. If we denote by {𝐵₁(𝑡),𝑡≥0}, {𝐵₂(𝑡),𝑡≥0}, and {𝐵₃(𝑡),𝑡≥0} the three correlated Brownian motions driving {𝑆₁(𝑡),𝑡≥0}, {𝑆₂(𝑡),𝑡≥0}, and {𝑆₃(𝑡),𝑡≥0} under the risk-neutral measure, then the three following processes:⎧⎪⎨⎪⎩𝐵𝐵₁(𝑡)−𝜎₃𝜌_1⋅3𝑡,𝑡≥0₂,(𝑡)−𝜎₃⎛⎜⎜⎜⎝𝜌_2⋅3−𝜌_1⋅2𝜌_1⋅31−𝜌^21⋅2⎧⎪⎨⎪⎩𝐵⎞⎟⎟⎟⎠𝑡,𝑡≥0₃⎫⎪⎬⎪⎭,(𝑡)−𝜎₃𝑡⎷1−𝜌^21⋅3−𝜌^2⋅3−𝜌_1⋅2𝜌_1⋅3²1−𝜌^21⋅2,𝑡≥0⎫⎪⎬⎪⎭(C.1) are independent standard Brownian motions under the new measure.