# A Note on the Distribution of Multivariate Brownian Extrema

Authors: Marcos Escobar; Julio Hernandez
Journal: International Journal of Stochastic Analysis
 (2014)

---

## Abstract

This paper presents a closed-form solution for the joint probability of the endpoints and minimums of a multidimensional Wiener process for some correlation matrices. This is the only explicit expressions found in the literature for this joint probability. The analysis can only be carried out for special correlation structures as it is related to the fundamentals regions of irreducible spherical simplexes generated by
reflections and the link to the method of images. This joint distribution can be used in financial mathematics to obtain prices of credit or market related products in high dimension. The solution could be generalized to account for stochastic volatility and other stylized features of the financial markets.


---

## Body

## 1. Introduction

The paper finds closed-form expressions for the joint density/distribution function of the endpoints and extrema of an-dimensional Wiener process. This targeted function represents a density function in terms of the endpoints of the Wiener process while it represents a distribution function of the minimums of the processes underlying. The results found in this paper can be applied to processes that can be derived from a Wiener process using suitable transformations like log-normal processes and Ornstein-Uhlenbeck processes with suitable parameters. The problem of finding such joint density/distribution has attracted research since the late nineteenth century, see, for example, [1]. Closed-form solutions are known for the cases of one and two dimensions; see [2, 3] with either a minimum or a maximum or both. The problem was recently solved in three dimensions; see [4], via the method of images; therefore it only applies to a restricted set of correlation values.In this paper we generalize the results in [4] to any dimension n. The correlation structures for which a solution via the Image Method can be found depend on the dimension (n). In principle there could be as many as 5 different correlation structures, without permutations, in dimension 4 or as little as three in dimensions 3, 5, and higher than 8. This is derived using the fundamental regions for irreducible groups generated by reflections as presented in Table IV page 297 in [5] together with the finding in [6, 7].The paper is organized as follows. Section2 introduces the function of interest as the solution of a partial differential equation (PDE) with specific initial and boundary conditions. Then it simplifies this PDE to a heat equation. Section 3 provides the main results described per dimension. Section 4 comments on applications as well as possible extensions. Section 5 concludes.


## 2. Preliminaries and Simplification to Heat Equation

Let(Ω,F,F,Q) be a filtered probability space on the domain Ω with sigma algebra F, filtration F={Ft}t≥0, and a probability measure Q on (Ω,F). Consider i=1,…,n, where n represents the dimension. We introduce the following notation:(1)
Yi(t): stochastic processes with constant drift αi and constant diffusion coefficient σi;(2)
Wi(t): Brownian motion processes with mean 0 and variance t;(3)
instantaneous correlation〈dWi(t),dWj(t)〉=ρijdt, i≠j; covariance matrix: Σ, correlation matrix: R=(ρij), and column vector drift α′=(α1,…,αn). The underlying process Y can be expressed through the following stochastic differential equation (SDE):
(1)dYit=αidt+σidWi;
in order to simplify the notation, we will assume Yi(0)=0.Our main objective is to find the joint density/distribution function for the minimum, denoted asYi_(t)≡min⁡0<s<t⁡Yi(s), and endpoints of Y(t); this function is defined below:
(2)py1,…,yn,m1,…,mn,tdy1⋯dyn=PY1_t>Y1t∈dy1,…,Ynt∈dyn,Y1_t>m1,…,Yn_t>mn,
where we also require Yi(0)>mi and therefore mi<0.In generalYi(t) could be seen as a function of a stochastic process Si such that a transformation Yi=f(Si) leads to a Wiener process. For example, Si could be a log-normal process (f(·)=ln⁡(·)); see [8] for the general type of transformations f and processes allowed.Whenever reasonable and in order to shorten the length of notation, we will writep(y,m,t)=p(y1,…,yn,m,t)=p(y1,…,yn,m1,…,mn,t). The same convention will be used for analogous intermediate functions.It is important to realize that if (2) is known then we could also derive the joint density/distribution of endpoints and maximums and all cases in between. To see this, first define the maximum of a process by Yi¯(t)≡max⁡0<s<t⁡Yi(s), and let Mi>0, i=1,…,k; then note that in the case with no drift
(3)PY1t∈dy1,…,Ynt∈dyn,Y1¯t<M1,…,Yk¯t<Mk,Yk+1_t>mk+1,…,Yn_t>mn=PY1*t∈dy1*,…,Ynt∈dyn,Y1*_t>m1,…,Yk*_t>mk,Yk+1_t>mk+1,…,Yn_t>mn=py1,…,yn,m1,…,mn,tdy1⋯dyn,
with Yi*=-Yi, mi=-Mi, yi*=-yi, i=1,…,k, and therefore correlation between Yi* and Yj, Corr(Yi*,Yj)=-Corr(Yi,Yj).It is known (see [9]) that the function p(·) above satisfies the following PDE, known as Fokker-Planck equation:
(4)∂p∂ty,t=-∑i=1n‍αi∂p∂yiy,t+12∑i=1n‍σi2∂2p∂yi2y,t+∑i<j‍σiσjρij∂2p∂yi∂yjy,t.
Since Yi0=0, i=1,…,n, the initial condition is
(5)py,m,t=0=∏i=1n‍δyi,
where δ(x) is the Dirac delta function with a spike at yi=0. The following n boundary conditions are added to match the extrema:
(6)py1,…,yi=mi,…,yn,m,t=0,i=1,…,n.
The meaning of the above constraints, in terms of the Wiener process, is clear; as soon as any of the component Yi equals mi we constraint the function p(·) to be equal 0. Now, noticing the initial conditions Yi(0)>mi and the continuity of the Wiener paths, we conclude that the constraints (6) are equivalent to enforcing that each path of Yi and therefore Yi_(t) remains above thebarrier levels mi before and up to time t. From the analytical point of view, that is, looking at the density p(y,m,t) as a function on Rn, the constraint (6) enforces that the support of the density is inside the set Ω={y∈Rn:y≥m}.Equations (4), (5), and (6) will be referred to as the PDE problem associated to the density (2). This function p(·) will be obtained for a set of special correlations by solving (4), (5), and (6) by means of the method of images; the reader is referred to [6, 7] for accounts of the method of images.To obtain the solution of the system (4), (5), and (6), we need to simplify (4) to a heat equation. The simplification is performed in two steps: first, we use conventional methods to modify the probability function to one associated to processes without drift, unit volatility, and barriers at zero; then, in a second step, we use a Cholesky transformation to generate a new set of independent processes, hence uncorrelated, leading to a heat equation with new boundary conditions.Let us define(7)d=Σ,ai=1dΣi,i=1,…,n,(8)b=a′Σa-αa,(9)py,m,t=qy,m,texp⁡a·y+bt,
where Σi is the covariance matrix with column i replaced by column vector α. It follows that
(10)∂q∂t=∑i=1n‍σi22∂2q∂yi2+∑i<j‍ρijσiσj∂2q∂yii∂yj,
with boundary conditions:
(11)qy1,…,yi=mi,…,yn,m,t=0,i=1,…,n,
and delta Dirac initial condition. We continue to simplify the above PDE, by eliminating the parameters σi and mi, i=1,…,n. Consider the following change of variables and transformation:
(12)xi=yi-miσi,i=1,…,n,hx,0,t=qy,m,t∏i=1n‍σi-1,
leading to
(13)∂h∂t=∑i=1n‍12∂2h∂xi2+∑i<j‍ρij∂2h∂xi∂xj,
with boundary conditions and initial condition,
(14)hx1,…,xi=0,…,xn,m,t=0,i=1,…,n,hx,0,t=0=∏i=1n‍δxi-xi0.
Here xi0=-mi/σi, i=1,…,n.The last step eliminates the correlations in matrixR; for this we perform a Cholesky decomposition.Proposition 1.
Let one define new variablesz as follows:
(15)z=xL-1,
where L is an upper triangular matrix in the Cholesky decomposition of R (R=L′L). Then the PDE associated to z satisfies
(16)∂h∂t=∑i=1n‍12∂2h∂zi2,
with boundary and initial conditions given by
(17)hHi,t=0,i=1,…,n,(18)hz,0,t=0=∏i=1n‍δzi-zi0,
where z0=(z10,…,zn0)=x0L-1 and Hi=z·Li=0, i=1,…,n, where Li represents column i of L.Proof.
This follows directly from noticing that the processes associated tozi are uncorrelated with zero drift.The variablesz are independent so (16), (18), and (17) represent the Fokker-Planck equation of an uncorrelated Brownian motion with zero drift and constrained to a region bounded by n hyperplane (Hi). Moreover the vector Li represents the normal vector to the hyperplane Hi, i=1,…,n.


## 3. Solving the Heat Equation

The method of images (MofI) is now utilized to solve the system (16), (18), and (17). The nature of the method of images is to replace the boundary conditions by a set of fictitious source points. Then the solution of the original equation satisfying the given boundary conditions reduces to that of finding the solution without boundary conditions at the source points. In the case of linear differential equations, the process of obtaining the final solution divides into three distinct steps (see [7]):(1)
checking that the differential equation is suitable for MofI and then solving this differential equation for a point source in an infinite medium, but with no boundary conditions except that of good behavior at infinity;(2)
checking that the region of interest is suitable for MofI and then finding the set of image source at each of the reflecting regions;(3)
summing the solution of(1) over the set of images obtained by step (2).In general, step(1) is the simplest one. The Laplace and heat equations are well known to fit into MofI (see, e.g., [6, 7]). The solution in step (1) is usually known in closed form. For example, in our case, the solution without boundary conditions is well known. Consider the following heat equation, for an arbitrary initial point b0=(b10,…,bn0):
(19)∂f∂t=∑i=1n‍12∂2f∂zi2,fz,0=∏i=1n‍δzi-bi0.
By means of an application of the Fourier transform, the solution of the above equation can be expressed as
(20)fz,t=12πtn/2exp⁡-12t∑i=1n‍zi-bi02.We will show next that the regions where MofI applies are connected to correlation matrices; therefore feasible regions imply feasible correlations. For the issue of whether the bounded region allows for the method of images, we rely on [6]. Note hyperplanes Hi pass through zero so they cut out the surface of an n-dimensional sphere centered at the vertex 0 in a convex polygonal domain also known as spherical simplexes. These spherical simplexes divide the sphere into symmetric parts where the method of images applies if and only if these simplexes are irreducible fundamental regions generated by reflections (see [5, 6]).The next result relates these regions and in particular the dihedral angle between these hyperplanes with the correlation matrixR for the multivariate process y(t). Let us denote by θij the dihedral angle between hyperplanes Hi and Hj with θ representing the matrix of angles. The next proposition relates these angles and the correlation matrix R.Proposition 2.
The relationship between the dihedral angle between hyperplanesHi and Hj, i,j=1,…,n, and the correlation matrix R satisfies -cos⁡(θ)=R.Proof.
The dihedral angle between hyperplanesHi and Hj can be obtained from the normal vectors to the plane as
(21)cos⁡π-θij=Li′·Lj′.
Here we use the fact thatLi=1 for all i. In matricial form we have
(22)-cos⁡θ=L′L=R.Next we use the results in the seminal work of [5] in which a complete list of fundamental regions for irreducible groups generated by reflections in the case of spherical simplexes is provided (Table IV, page 297). We use the notation in Coxeter for irreducible groups for spherical simplexes in dimension n: An, Bn, Cn, Dn, En, Fn, and Gn. In particular we extract the feasible dihedral angles for each of these spherical simplexes and the dimension where they apply. These results are given next.Proposition 3.
The type of correlation matricesR where the method of images can be used to solve the system (16), (18), and (17) is provided next. The total number of source points needed, N, for an explicit expression of the solution is also provided.(i)
For anyn>3 there are at least 3 cases of correlation matrices:
(23)An⟺Ri,j=-cos⁡π3,j=i+1,0,else,N=n+1!,Bn⟺Ri,j=-cos⁡π3,j=i+1,i=n-2,j=n,0,else,N=n!·2n-1,Cn⟺Ri,j=-cos⁡π3,j=i+1,-cos⁡π4,i=n-1,j=n,0,else,N=n!·2n.(ii)
Ifn=6,7 or 8 then there is one more case (4 total):
(24)En⟺Ri,j=-cos⁡π3,j=i+1,i≠4,0,i=4,j=5,-cos⁡π3,i=3,j=5,0,else.Forn=6 we have N=6!·72, n=7 implies N=9!·8, and for n=8 then N=10!·192.(iii)
Ifn=5 then there are only three cases: A5, B5, and C5.(iv)
Ifn=4 then there are a total of 5 cases:
(25)F4⟺Ri,j=-cos⁡π3,j=i+1,i≠2,-cos⁡π4,i=2,j=3,0,else,N=1152,G4⟺Ri,j=-cos⁡π3,j=i+1,i≠3,-cos⁡π5,i=3,j=4,0,else,N=1202.(v)
Ifn=3 then there are only three cases, A3  (N=24), C3  (N=48), and G3:
(26)G3⟺Ri,j=-cos⁡π3,j=i+1,i≠2,-cos⁡π5,i=2,j=3,0,else,N=120.(vi)
Ifn=2 then there are infinitely many cases (any p):
(27)D2⟺R1,2=-cos⁡πp,N=2p.Proof.
This follows from the dihedral angles provided in Table IV, page 297, in [5] together with Proposition 2. In this table, the regions are represented by graphs. Every node represents a bounding hyperplane, and the branches indicate pairs of hyperplanes inclined at angles π/p, p>2. If a value of p is not provided then a p=3 is understood. Perpendicular hyperplanes are represented by nodes not joined by a branch. Note that all hyperplanes should intercept; therefore most dihedral angles are π/2. The cases obtained in this proposition follow easily from reading these graphs out.In particular, for dimensions 2 and 3, the number of feasible cases and source points reproduces those found in [2, 4], respectively.Corollary 4.
Given a collection of source points inRn{bj}j=1N where bj=(b1j,…,bnj) and the sign associated to each point g(j)j=1N, where g(j)=0,1, then the solution to the system (16), (18), and (17) can be found as follows:
(28)hz,t=∑j=1N‍12πtn/2-1g(j)exp⁡-12t∑i=1n‍zi-bij2.Note that (28) is basically a linear combination of n-dimensional gaussian density functions with mean zero and (1/t)I covariance matrix. Substituting (28) into (16) and further into (9) would lead to the targeted joint density/distribution in (2) and this would be again a linear combination of n-dimensional gaussian densities but now with nonzero means (depending on the image sources) and covariance matrix (1/t)Σ as in [4].


### 3.1. Finding Source Points

In this section we describe a quasi-analytical procedure to find all source points associated to a spherical simplex and therefore to a correlation matrix.Let us denote(i)
Li: normal vector to hyperplane i, with i=1,…,n;(ii)
x0: initial point inside the region;(iii)
N: number of source points;(iv)
Lik (Li,Li1,…,Lik-1): normal vector to an hyperplane created after k consecutive reflections, each across the hyperplane associated to the normal vector in the sequence (Li,Li1,…,Lik-1);(v)
x0(Li,Li1,…,Lik-1), sg(x0(Li,Li1,…,Lik-1)): source point and sign of it, created after k consecutive reflections, each across the hyperplane associated to the normal vector in the sequence (Li,Li1,…,Lik-1).We show next a method to create the image of a hyperplane, the image of a point reflected across a given hyperplane (passing through zero), and the sign of the new point as needed by the method of images; this uses standard concepts from geometry.(1)
New source point: reflecting a pointx0 across a hyperplane passing through zero defined by the normal vector Li leads to a new point with equation:
(29)x0(Li)=x0+-2x0·LiLi.(2)
New hyperplane: reflecting a hyperplane with normal vectorLj across a hyperplane with normal vector Li leads to a new hyperplane with normal vector:
(30)LjLi=aLi+bLj,a=-2c,b=1,c=Li·Lj.Note the normal vector to bothLi and Lj must also be normal to Lj(Li); this is why the latter is in the same hyperplane as the formers. The equation is obtained then using two facts: first the norm of Lj(Li) is one; therefore (a2+b2+2abc=1) and, second, Lj(Li)·(-Li)=L2·L1; therefore c=-a-bc.(3)
Sign of new point: assume original point has a positive sign “+” then(31)sgx0(Li,Li1,…,Lik-1)=-1k+1.The algorithms are based on doing reflections across all originals hyperplanes and then repeating the procedure for all new hyperplanes. This method would lead to repeated values; therefore we also have to check if there are duplicates (hyperplanes or source points). The method stops after the known number of different sourcesN has been detected but in principle few iterations should lead to a good approximation.


## 3.1. Finding Source Points

In this section we describe a quasi-analytical procedure to find all source points associated to a spherical simplex and therefore to a correlation matrix.Let us denote(i)
Li: normal vector to hyperplane i, with i=1,…,n;(ii)
x0: initial point inside the region;(iii)
N: number of source points;(iv)
Lik (Li,Li1,…,Lik-1): normal vector to an hyperplane created after k consecutive reflections, each across the hyperplane associated to the normal vector in the sequence (Li,Li1,…,Lik-1);(v)
x0(Li,Li1,…,Lik-1), sg(x0(Li,Li1,…,Lik-1)): source point and sign of it, created after k consecutive reflections, each across the hyperplane associated to the normal vector in the sequence (Li,Li1,…,Lik-1).We show next a method to create the image of a hyperplane, the image of a point reflected across a given hyperplane (passing through zero), and the sign of the new point as needed by the method of images; this uses standard concepts from geometry.(1)
New source point: reflecting a pointx0 across a hyperplane passing through zero defined by the normal vector Li leads to a new point with equation:
(29)x0(Li)=x0+-2x0·LiLi.(2)
New hyperplane: reflecting a hyperplane with normal vectorLj across a hyperplane with normal vector Li leads to a new hyperplane with normal vector:
(30)LjLi=aLi+bLj,a=-2c,b=1,c=Li·Lj.Note the normal vector to bothLi and Lj must also be normal to Lj(Li); this is why the latter is in the same hyperplane as the formers. The equation is obtained then using two facts: first the norm of Lj(Li) is one; therefore (a2+b2+2abc=1) and, second, Lj(Li)·(-Li)=L2·L1; therefore c=-a-bc.(3)
Sign of new point: assume original point has a positive sign “+” then(31)sgx0(Li,Li1,…,Lik-1)=-1k+1.The algorithms are based on doing reflections across all originals hyperplanes and then repeating the procedure for all new hyperplanes. This method would lead to repeated values; therefore we also have to check if there are duplicates (hyperplanes or source points). The method stops after the known number of different sourcesN has been detected but in principle few iterations should lead to a good approximation.


## 4. A Comment on Applications

One of the key application fields for our findings is multidimensional financial derivatives (see [10] or [11] for an introduction to financial markets and problems). Market derivatives are contracts payable at future times, called maturity, deriving their value from the performance of potentially several underlying tradeable stocks. There exists an extensive family of financial derivatives in the market which depends on first passage time (barriers) of the underlying stock price processes. In general, as pointed out by [12], adding barriers is a convenient method for reducing the cost of a derivative.Two main families of such products are double lookback options (see [2]), in particular double barrier options (see [13]), and mountain range derivatives (see [14]). The latter are high dimensional (n≥3) and were created by Société Générale in 1998. Examples of these products are Altiplano, Annapurna, and Atlas. The payoff, at maturity T, of an Annapurna is of the form
(32)∏i=1n‍1Yi_(T)>log⁡Ki.HereY stands for the log stock price while K’s are prespecified strike prices. Therefore the product pays a dollar at time T if and only if all stocks remain above given strike prices (Ki, i=1,…,n) during the relevant period of the product, that is, (0,T]. The price of such a product, C(0), is the expected value of this payoff. Therefore the price is related to our function 1 in the following manner:
(33)C0=∫∫Ωpy,log⁡K,Tdy,
where Ω=y∣yi>log⁡Ki,i=1,…,n.In the past this expression could be evaluated either via Monte Carlo simulations or directly solving the PDE equations, which are highly time consuming and inaccurate approaches for dimensions higher than 3. This price, under the feasible correlations described in Proposition3, could be now found as a linear combination of n-dimensional Gaussian cumulative distribution functions.Another family of products benefiting from this work are credit derivatives, in particular collateralized debt obligation (CDO) and akth to default product (see [15]). These were products at the very heart of the financial crisis in 2008. A CDO has a payoff similar to that of a market product under proper default assumptions. For instance, if a default is assumed to be triggered by the company’s assets crossing its constant debt at any time prior to maturity of the companies’ debt, as proposed in the seminal work of [16], then the key element in the price of a CDO is the default of k prespecified companies simultaneously, which can be expressed as follows:
(34)E∏i=1k‍1Yi_(T)<log⁡Di=∫∫Ωkpy,log⁡D,Tdy,
where Ωk=y∣yi>log⁡Di,i=1,…,k, Y stands for the asset value of the company, D represents its constant debt, and n=k. This type of products was highly mispriced due, in part, to the mathematical complexity of handling multidimensions and first passage time. Simpler approaches like those assuming default only at maturity (avoiding first passage time) or noneconomical approaches like modelling default as an exogenously given process were more prone to oversight.


## 5. Conclusions and Possible Generalizations

This paper describes the correlation matrices for which a closed-form solution for the joint density/distribution of the endpoints and the minimum of a Wiener process can be found. The results are also applicable to other processes like log-normal. The general solution requires a detailed geometrical analysis of certain partitions of then-dimensional sphere; therefore the technique can only provide the closed-form solutions for a specific set of correlations. The resulting densities can be used in several applications, in particular analytical expressions for prices of financial products, therefore validating the accuracy of numerical simulations.The method developed in the present paper could also be extended to allow for stochastic volatility and random correlation. It can be applied to maximums and minimums combined, as long as one extreme per dimension is considered. Finally, the solution could be the basis for further approximations like those based on perturbation theory (see [17]) which currently work under the assumption of independence among the processes underlying.



---

*Source: 575270-2014-11-16.xml*