Constrained Wiener Processes and Their Financial Applications

Show more

1. Background

It is natural to model many financial variables as Wiener processes: examples may be found in asset returns, interest rates, bond yields as well as inflation and commodity prices. In fact, any stationary variable (and any variable integrated of order zero $I\left(0\right)$ ) is amenable to such modelling. The markets relating to such variables reflect their various characteristics in terms of long run behaviour, volatility and higher order moments.

Exotic derivatives are designed to exploit the higher order characteristics through their dependence on the evolution of variables over finite time periods. This is where the behaviour of extrema of Wiener processes becomes relevant and important.

There are well known results for the extrema for a one dimensional Wiener processes. For such a process, the extrema are either the maximum or minimum. But in practice financial variables cannot be viewed to operate in isolation, and it is the interplay between several variables that leads to complexity and richness in financial modelling.

For example we may wish to model:

• the behaviour of the asset returns involved in an overall portfolio;

• the interplay between bond yields and interest rates in a fixed interest market;

• the dependence of asset behaviour on economic variables such as inflation, employment and wage growth.

Multidimensional Wiener processes are the natural framework for such modelling. However the extrema processes in 2 or higher dimension are richer and more varied than in the one dimensional case. The purpose of this paper is to study the distribution of these extrema.

2. Heat Equation

Thee relationship between the Wiener process and the heat equation was first studied extensively by [1] , who approached the solution of a Dirichlet problem for the heat equation with methods using probabilistic expectations. In this paper, we do the reverse, that is we solve the probabilistic problem of the distribution of extrema by expressing it as a Dirichlet problem and then using elementary methods for partial differential equations (PDEs).

By focussing on the Wiener process, rather than the Black Scholes PDE (which are logically interrelated), we obtain simplifications in the analysis that obviate the practical issues found in [2] . Once the distributions of extrema are found, they may be applied directly and intuitively to the pricing of exotic options.

3. Characterisation of Extrema

The simplest way of viewing the role of extrema is in discrete time, with a random walk in D dimensions. Such a walk starts at $0$ and proceeds in steps of $\pm 1$ in one of the 2D possible directions, all with equal probability.

Thus after n periods there are ${\left(2D\right)}^{n}$ possible paths. After n steps a point ${x}_{n}\in {\mathbb{R}}^{D}$ is reached, and a path is denoted $x=\left(0\mathrm{,}{x}_{1}\mathrm{,}{x}_{2},\cdots \right)$

In general we wish to restrict the possible paths by imposing m linear constraints of the form ${a}_{i}^{\text{T}}x\le {y}_{i}$ where $i=1,2,\cdots ,m$ and ${a}_{i}\mathrm{,}{y}_{i}\in {\mathbb{R}}^{D}$ are specified vectors.

The set of admissible points is denoted $U=\left\{x|Ax\le y\right\}$ where $A={\left[{a}_{1}\mathrm{,}{a}_{2},\cdots \right]}^{\text{T}}$ and $y={\left[{y}_{1}\mathrm{,}{y}_{2},\cdots \right]}^{\text{T}}$ . We suppose that $A$ is of full row rank, so there is no redundancy in the constraints. Without loss of generality we may also take $\Vert {a}_{i}\Vert =1$ for all i. Clearly U is a closed convex set, though possibly infinite, and its boundary is denoted $\partial U=\left\{x|Ax=y\right\}$ .

Denote by ${f}_{n}\left(x\right)$ the probability of reaching $x\in {\mathbb{R}}^{D}$ after n steps under a path wholly contained in U. It is clear from the evolution of paths that

${f}_{n}\left(x\right)=\{\begin{array}{ll}1/2D{\displaystyle \underset{i}{\sum}}{f}_{n-1}\left(x\pm {e}_{i}\right)\hfill & x\in U\hfill \\ 0\hfill & x\notin U\hfill \end{array}$ (1)

where ${e}_{i}$ are the unit vectors in the ith direction in ${\mathbb{R}}^{D}$ . This recurrence relation defines ${f}_{n}\left(x\right)$ completely, starting with ${f}_{0}\left(0\right)=1$ .

It is well known that the limit as $n\to \infty $ of the random walk is the Wiener process in D dimensions [3] . We exploit this property to establish the characterisation of the constrained Wiener process in continuous time. However the random walk must be scaled for this to occur; as the standard random walk produces a variance of n/D after n steps, we need to employ a scaling factor of $\sqrt{n/D}$ in measuring the level attained after n steps.

For a given period t in continuous time, we thus consider nt steps to reach a level of $x\sqrt{n/D}$ . (The values of n must be chosen to give this meaning―that is, n/D must be a perfect square and nt must be an integer. However all this is possible since we consider only the behaviour as $n\to \infty $ .) Then the probability distribution ${\left(n/4D\right)}^{D/2}{f}_{nt}\left(x\sqrt{n/D}\right)$ might be considered to converge to a density, with the factor ${\left(n/4D\right)}^{D/2}$ being required to allow for the discrete spacing $\sqrt{n/D}/2$ across each dimension of $x$ . All levels, including the constraints $y$ need to be so scaled.

The following operators are used:

$\nabla =\frac{\partial}{\partial x}={\left[\frac{\partial}{\partial {x}_{1}}\mathrm{,}\frac{\partial}{\partial {x}_{2}},\cdots ,\frac{\partial}{\partial {x}_{m}}\right]}^{\text{T}}$

$div=tr(\nabla )=\frac{\partial}{\partial {x}_{1}}+\frac{\partial}{\partial {x}_{2}}+\cdots +\frac{\partial}{\partial {x}_{m}}$

$\Delta =tr\left({\nabla}^{2}\right)=\frac{{\partial}^{2}}{\partial {x}_{1}^{2}}+\frac{{\partial}^{2}}{\partial {x}_{2}^{2}}+\cdots +\frac{{\partial}^{2}}{\partial {x}_{D}^{2}}$

The discrete time problem in (1) converges to the heat PDE as follows.

Proposition 1. Suppose that the probability distribution converges to a continuous time density. Then that density is of the constrained Wiener process $\phi \left(t\mathrm{,}x\right)$ at time t, given by the solution of the PDE

${\phi}_{t}=1/2\Delta \phi $

with the boundary conditions $\phi =0$ on $\partial U$ .

Proof. Though the result is well known, its proof is very technical and not fully provided in the literature. Hence it is set out in Appendix A. The problem with the spatial boundary conditions is known as a Dirichlet problem.

Remark. The convergence for the one dimensional case is guaranteed by the

De Moivre-Laplace lemma, in which case ${f}_{n}\left(x\right)={2}^{-n}\left(\begin{array}{c}n\\ n+x/2\end{array}\right)$ and the limit is provided by

${f}_{n}\left(x\right)=\sqrt{2/\text{\pi}n}{\text{e}}^{-1/2{x}^{2}}+o\left(1/n\right)$

uniformly over compact sets of x. This leads directly to the Gaussian density for $\sqrt{n}/2{f}_{nt}\left(x\sqrt{n}\right)$ . In higher dimensions the limit may be derived by considering the Fourier transform

$E\left({\text{e}}^{i{\theta}^{\text{T}}x}\right)\mathrm{=}1/{2}^{D}{\left({\text{e}}^{i{\theta}_{1}{x}_{1}}+{\text{e}}^{-i{\theta}_{1}{x}_{1}}+{\text{e}}^{i{\theta}_{2}{x}_{2}}+{\text{e}}^{-i{\theta}_{2}{x}_{2}}+\cdots \right)}^{n}$

which is tantamount to proving the central limit theorem in D dimensions (Dym).

Remark. The appeal to the discrete time case as given above provides the intuition to the link with the heat equation. In fact a direct derivation is possible. If $\phi \left(t\mathrm{,}x\mathrm{,}y\right)$ denotes the constrained density under a D dimensional Wiener process, then it must satisfy the functional relationship

$\phi \left(s+t\mathrm{,}x\mathrm{,}y\right)={\displaystyle \int}\phi \left(t\mathrm{,}z,y\right)\phi \left(s\mathrm{,}x-z\mathrm{,}y-z\right)\text{d}z\text{\hspace{0.05em}}\text{\hspace{0.05em}}.$

Taking the later density $\phi \left(s\mathrm{,}x-z\mathrm{,}y-z\right)$ for a small time interval $s=\text{d}t$ leads to the heat equation.

Solving the Heat Equation Using Random Walks

The above result suggests that Dirichlet problems for the heat equation may be solved numerically with a sufficient number of steps of a random walk. This is an efficient alternative to finite difference methods for numerical solution.

An example is for the constraints

$\left[\begin{array}{c}{x}_{1}+{x}_{2}\le 1\\ {x}_{2}\le 0.5\end{array}\right]$

as set out in Figure 1.

In fact the constraints need not be linear; the same process applies to exponential constraints of the form

${\text{e}}^{{x}_{1}}+{\text{e}}^{{x}_{2}}\le 3$

as depicted in Figure 2.

Figure 1. Two constraints.

Figure 2. Exponential constraints.

4. Solution of the Heat Equation

The equation for the constrained Wiener process may be seen as a Dirichlet problem for the heat equation, with linear boundary conditions. This may be solved by transforming the boundary conditions into initial conditions (that is, at time $t=0$ ), using the technique known in the literature as Kelvin’s method of images (or heat poles [4] ).

The analogy with the heat equation is more than superficial. It is well known that the unconstrained heat equation ${\phi}_{t}=1/2\Delta \phi $ has as its unique solution the Gaussian distribution $\phi ={\text{e}}^{-1/2t{\left|x\right|}^{2}}{\left(2\text{\pi}t\right)}^{D/2}$ with a single heat source (i.e. initial condition) $\phi \left(\mathrm{0,}x\right)=\delta \left(x\right)$ where the right hand side is the Dirac delta function.

To formalise the solution, the density for the constrained Wiener process satisfies the Dirichlet (boundary) problem

$\partial \phi /\partial t=1/2\Delta \phi $ with boundary conditions $\phi =0$ on $\partial U$ (2)

4.1. Uniqueness

Proposition 2. If $\phi $ is a solution of the heat equation, then it is the unique solution.

Suppose that ${\phi}_{1}$ and ${\phi}_{2}$ are two such solutions, and let $\psi ={\phi}_{1}-{\phi}_{2}$ , which also satisfies the heat equation and is zero on the boundary. Then consider the energy integral $\underset{U}{\int}}{\psi}^{2}\text{d}x$ . This must be finite as $\psi $ is continuous and absolutely integrable, and has the time derivative

$\partial /\partial t{\displaystyle \underset{U}{\int}}{\psi}^{2}\text{d}x=2{\displaystyle \underset{U}{\int}}\psi \partial \psi /\partial t\text{d}x={\displaystyle \underset{U}{\int}}\psi \Delta \psi \text{d}x\mathrm{.}$

However we have the identity

$div\left(\nabla {\psi}^{2}\right)=2\psi \Delta \psi +2{\Vert \nabla \psi \Vert}^{2}$

so that

$\partial /\partial t{\displaystyle \underset{U}{\int}}{\psi}^{2}\text{d}x=-{\displaystyle \underset{U}{\int}}{\Vert \nabla \psi \Vert}^{2}\text{d}x+{\displaystyle \underset{\partial U}{\int}}\psi \nabla \psi \cdot \text{d}S$

where the surface integral on $\partial U$ follows from the divergence theorem and equals zero. Since both $\underset{U}{\int}}{\psi}^{2}\text{d}x$ and $\underset{U}{\int}}{\Vert \nabla \psi \Vert}^{2}\text{d}x$ are non-negative, we must have $\psi =0$ , giving only one solution to the heat equation.

4.2. An Initial Condition Problem

We will replace the boundary conditions under the Dirichlet problem with initial conditions that lead to the same solution. First we introduce the reflection operator about a hyperplane ${a}^{\text{T}}x=y$ . This is given by $T$ where

$Tx=\left(I-2a{a}^{\text{T}}\right)x+2ya$

It is clear that ${T}^{2}=I$ and leaves the hyperplane invariant. It is also clear that the transformation $T$ leaves the heat equation invariant, that is $\phi \left(t\mathrm{,}x\right)$ satisfies the heat equation, so too does $\phi \left(t\mathrm{,}Tx\right)$ . This follows immediately from

$\Delta \phi \left(t\mathrm{,}Tx\right)=tr\left[{\left(I-2a{a}^{\text{T}}\right)}^{\text{T}}{\nabla}^{2}\left(I-2a{a}^{\text{T}}\right)\phi \right]=\Delta \phi \left(t\mathrm{,}x\right)\mathrm{.}$

However the transformed equation has the initial condition $\phi \left(\mathrm{0,}Tx\right)=\delta \left(x\right)$ . We may apply the reflection for each of the constraints, denoting by ${T}_{i}$ that corresponding to the ith:

${T}_{i}x=\left(I-2{a}_{i}{a}_{i}^{\text{T}}\right)x+2{y}_{i}{a}_{i}$

The group generated by $\mathcal{T}=\left\{{T}_{i}\mathrm{|1}\le i\le m\right\}$ is in general a countably infinite Coxeter group with a complex structure. We use it to generate initial conditions under the following process.

Generating Process for Heat Sources/Sinks

Let ${S}_{0}=\left\{0\right\}$ be the initial singleton set. Then for $i=1,2,\cdots $ let ${S}_{i}=\mathcal{T}{S}_{i-1}$ after removing any points which have appeared previously in ${S}_{0}\cup {S}_{1}\cup \cdots \cup {S}_{i-2}$ . Then an alternative initial value problem can formulated as:

$\partial \phi /\partial t=1/2\Delta \phi $ with initial conditions $\phi \left(x\right)={\left(-1\right)}^{i}\delta \left(x-z\right)$ if $z\in {S}_{i}$ (3)

This method of construction ensures that:

• A point $x\in {S}_{i}$ for only one value of i;

• For any ${T}_{i}\in \mathcal{T}$ , the delta function at $x$ has the opposite sign to that at ${T}_{i}x$ .

By analogy with the heat equation, it is convenient to refer to the delta functions at time $t=0$ in 3 as (heat) sources where the corresponding coefficient is +1, and as (heat) sinks where the coefficient is −1.

The IC problem is clearly defined for the whole of ${\mathbb{R}}^{D}$ , and we then have our main result.

Proposition 3. If the solution of the Dirichlet problem in 2 can be extended to the whole space ${\mathbb{R}}^{D}$ , then the problem in 3 and the Dirichlet problem have the same solution, and it is unique.

Proof. We show that the solution to the IC problem, if it is well posed, satisfies the Dirichlet problem. Consider the operator ${T}_{i}$ . If $\phi \left(t\mathrm{,}x\right)$ solves the IC problem, then so too does $\phi \left(t\mathrm{,}{T}_{i}x\right)$ as the initial conditions are transformed, with the sign reversed. Thus $\phi \left(t\mathrm{,}x\right)$ and $-\phi \left(t\mathrm{,}{T}_{i}x\right)$ both solve the IC problem. But by uniqueness they must be equal.

On the hyperplane ${a}_{i}^{\text{T}}x={y}_{i}$ we have $x={T}_{i}x$ , and thus $\phi \left(t\mathrm{,}x\right)=-\phi \left(t\mathrm{,}x\right)$ , which implies that $\phi \left(t\mathrm{,}x\right)$ vanishes on every hyperplane, thereby satisfying the Dirichlet problem. By uniqueness of the solution, the two problems have the same solution in the domain U.

Conversely suppose that the Dirichlet problem can be extended to ${\mathbb{R}}^{D}$ . Then by uniqueness it must have the same solution as that of the IC problem, which involves the series of heat sources and sinks described in the above generating process.

The proposition extends the method of images used for the heat equation to many dimensions and types of constraint. However the images method is but one of several group theoretic methods, also known as similarity methods [5] . They are of practical importance as the solution to the IC problem in many cases can be written down as the summation of Gaussian densities, viz

$\phi \left(t\mathrm{,}x\right)=1/{\left(2\text{\pi}t\right)}^{D/2}{\displaystyle \underset{z\in {S}_{i}}{\sum}}{\left(-1\right)}^{i}{\text{e}}^{-{\Vert x-z\Vert}^{2}/2t}\mathrm{.}$

4.3. The Case D = 1

In the special case $D=1$ , we can have at most two constraints, being the maximum ${y}_{1}$ and minimum ${y}_{2}$ of the Wiener process, with ${a}_{1}=1,{a}_{2}=-1$ . In this case ${T}_{1}x=2{y}_{1}-x$ and ${T}_{2}x=2{y}_{2}-x$ . The initial conditions are therefore at the points shown below (letting ${z}_{1}=2{y}_{1},{z}_{2}=2{y}_{2}$ for simplicity)

$\begin{array}{ccc}\text{Resultofoperatoratlevel}& {T}_{1}& {T}_{2}\\ 0& {z}_{1}& {z}_{2}\\ 1& {z}_{1}-{z}_{2}& {z}_{2}-{z}_{1}\\ 2& {z}_{1}-\left({z}_{2}-{z}_{1}\right)& {z}_{2}-\left({z}_{1}-{z}_{2}\right)\\ 3& {z}_{1}-\left(2{z}_{2}-{z}_{1}\right)& {z}_{2}-\left(2{z}_{1}-{z}_{2}\right)\\ & \vdots & \vdots \end{array}$

It is not hard to show that the points at the nth level are given by

$\begin{array}{ccc}\text{Level}& {T}_{1}& {T}_{2}\\ 2n& 2n\left({y}_{1}-{y}_{2}\right)& 2n\left({y}_{2}-{y}_{1}\right)\\ 2n-1& 2n{y}_{1}-2\left(n-1\right){y}_{2}& 2n{y}_{2}-2\left(n-1\right){y}_{1}\end{array}$

Hence the density is given by

$\phi \left(t\mathrm{,}x\right)=1/\sqrt{2\text{\pi}t}{\displaystyle \underset{n=-\infty}{\overset{\infty}{\sum}}}\left[{\text{e}}^{-{\left(x-2n{y}_{1}+2n{y}_{2}\right)}^{2}/2t}-{\text{e}}^{-{\left(x-2n{y}_{1}+2\left(n-1\right){y}_{2}\right)}^{2}/2t}\right]\mathrm{.}$ (4)

This result was originally derived in a stochastic setting by [6] , and indirectly in an option pricing setting in [2] . However from a much older historical perspective, this problem is identical to that of finding the temperature in a rod of length 2, where both ends at $x=\pm 1$ are held at zero temperature and with an initial unit source at $x=0$ :

$\partial \phi /\partial t=1/2{\partial}^{2}\phi /\partial {x}^{2}$

and

$\phi \left(t,\pm 1\right)=0$ , $\phi \left(0,x\right)=\delta (x)$

The classical solution [7] derived from Fourier series is

$\phi \left(t\mathrm{,}x\right)=-2/y{\displaystyle \underset{n\mathrm{=0}}{\overset{\infty}{\sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{\text{e}}^{-{\left(2n+1\right)}^{2}{\text{\pi}}^{2}/{y}^{2}t/2}sin\left(n\text{\pi}{y}_{1}/y\right)sin\left(n\text{\pi}x-{y}_{1}/y\right)$ (5)

where $y={y}_{2}-{y}_{1}$ . The form of this solution is more complex than the Gaussian form above, as it contains oscillatory elements, and is thus less preferred for calculation purposes. The equality of (4) and (5) is known as Jacobi’s identity [7] .

Where only a maximum is relevant, the density reduces to the well known result

$\phi \left(t\mathrm{,}x\right)=1/\sqrt{2\text{\pi}t}\left[{\text{e}}^{-{x}^{2}/2T}-{\text{e}}^{-{\left(2y-x\right)}^{2}/2t}\right]\mathrm{.}$ (6)

5. Analytic Solutions

The examples given above demonstrate that analytic solutions of the extremum problem can always be found for $D=1$ . In higher dimensions this is not always so. This may be because:

• the Dirichlet problem may not be extensible to the whole of ${\mathbb{R}}^{D}$ [in this case multiple sheets of ${\mathbb{R}}^{D}$ may be useful [4] ; or

• from a numerical viewpoint, the heat sources and sinks cluster within a limited domain in ${\mathbb{R}}^{D}$ and computational accuracy becomes an issue.

Where these issues arise, it is always possible to resort to the numerical procedure of solving the heat equation by using random walks as presented in section 3.1. In addition it sometimes happens that a slight variation in the constraints will lead to an analytic solution. The possibilities are illustrated in this section.

As a general rule, the most successful cases where the IC approach leads to an analytic solution are where the heat sources and sinks become infinitely dispersed under the Coxeter reflection group, as evident for the 2 extrema case for $D=1$ . The most serious difficulties are where the Coexter group is finite [8] , as this may lead to the Dirichlet problem being inextensible, in contrast to the assertions given in [9] .

Remark. For computational purposes, it is convenient to transform the affine reflections involved into linear reflections. This can be done by imbedding ${\mathbb{R}}^{D}$ into a hyperspace of projective space ${\mathbb{R}}^{D+1}$ where the last coordinate becomes 1. Thus the linear reflection corresponding to the affine reflection $Tx=\left(I-2a{a}^{\text{T}}\right)x+2ya$ becomes

$\stackrel{\u02dc}{T}\left[\begin{array}{c}x\\ 1\end{array}\right]=\left[\begin{array}{cc}I-2a{a}^{\text{T}}& 2ya\\ 0& 1\end{array}\right]\left[\begin{array}{c}x\\ 1\end{array}\right]$

The example below are all given for $D=2$ , as this should suffice to illustrate what is possible. In the figures below, the density of the contour indicate the level of the Wiener density function. Heat sources are indicated as red points, and heat sinks as blue. The caption of each figure indicates the constraints imposed.

5.1. The Case m = 1

An example of a single constraint for $D=2$ is shown in Figure 3.

For this possibility, all densities are well behaved and are similar to the case for $D=1$ .

5.2. The Case m = 2

An example of two constraints for $D=2$ is shown in Figure 4.

Most of the cases are well behaved. However the constraints

$\left[\begin{array}{c}{x}_{1}+{x}_{2}\le 1\\ {x}_{2}\le 0.5\end{array}\right]$

lead to a finite Coxeter group, and it is easy to show that the additional constraint ${x}_{1}\ge 0.5$ is necessarily implied by the other two. This is an example of an inextensible Dirichlet problem. Adding this additional constraint, however, leads to a well posed Dirichlet solution as shown in Figure 5.

Figure 3. A single constraint.

Figure 4. Two constraints.

Figure 5. An invisible constraint.

However the original problem, without the additional constraint, might be approximated by varying the constraint levels to $\u03f5={10}^{-4}$ as shown in Figure 6.

Which can be compared with Figure 1.

It also of interest to approximate non-linear constraints with a set of linear ones. Here is an example which is relevant to evaluating exotic options as discussed in Section 7, and shown in Figure 7.

Figure 6. A slightly varied situation with two constraints.

Figure 7. An approximation to exponential constraints.

5.3. The Case m ≥ 3

Examples of multiple constraints are shown in Figure 8 and Figure 9.

No difficulties are evident, except that the larger the number of constraints, the greater the possibility of clustering and therefore of computational problems. Slight variations in the constraints may resolve these difficulties, as in Figure 6.

Figure 8. Three constraints.

Figure 9. Four constraints.

6. Time of First Breach

The density $\phi \left(t\mathrm{,}x\right)$ clearly depends on the level of the extrema y and is thus a cumulative distribution with respect to the extrema; this may be recognized by expressing it as $\phi \left(t\mathrm{,}x\mathrm{,}y\right)$ , with ${y}_{i}=max\left({a}_{i}^{\text{T}}x\right)$ . The joint density $\Phi \left(t\mathrm{,}x\mathrm{,}y\right)$ of the Wiener process $x$ and its extrema $y$ may then be found from

$\Phi \left(t\mathrm{,}x\mathrm{,}y\right)=\frac{\partial}{\partial y}\phi \left(t\mathrm{,}x\mathrm{,}y\right)$ .

In particular, if k denotes a specific extremum, the density of the process which obeys all the constraints apart from k, which is breached at some time, is:

$\Phi \left(t\mathrm{,}x\mathrm{,}y\backslash {y}_{k}\right)-\Phi \left(t\mathrm{,}x\mathrm{,}y\right)$

where the symbol\denotes exclusion of the relevant constraint. It is not difficult to incorporate time of breach of extremum k in this analysis. This is useful in pricing options (such as of the Parisian type), where the timing of breach and re-entry are relevant.

For this purpose, consider the situation where the time of first breach of extremum k occurs after time $\tau $ . Then the process must satisfy the Dirichlet condition until time $\tau $ at a level $z\le y$ satisfying all extrema, and breach extremum k after time $\tau $ . The density of the process at level $x$ is therefore

$\Phi \left(\tau \mathrm{,}z\mathrm{,}y\right)\Phi \left(t-\tau \mathrm{,}x-z,y\backslash {y}_{k}\right)\mathrm{.}$

Summation over all possible times $0\le \tau \le t$ and levels $z\le y$ therefore provides the density of meeting the extrema, apart from k and breaching k after time $\tau $ but before time t. This can be obviously generalized to breaches of multiple extrema at specified times. This provides a straightforward generalization of the result in [10] for $D=1,m=1$ .

7. Application to Option Pricing

We now conclude with our initial motivation for considering extrema of Wiener processes―option pricing. This is based, not on asset prices or variables being modelled as Wiener processes, but more commonly as asset returns being such

$\text{d}P/P=\mu \text{d}t+\Sigma \text{d}x$

where $x$ is the outcome of a Wiener process as before. This results in

$P={\text{e}}^{\left(\mu -1/2{\sigma}^{2}\right)t+\Sigma x}$

where ${\sigma}^{2}=tr\left({\Sigma}^{2}\right)$ . The exponent in the above price is a Wiener process with drift $\mu -1/2{\sigma}^{2}$ .

The general expression for call options with strike $K$ on the assets is

$C=E\left({P}_{t}-K\mathrm{|}{P}_{s}\in U\mathrm{,0}\le s\le t\right)$

The exponent appearing in the price is of the form $\left(\mu -1/2{\sigma}^{2}\right)t+\Sigma x$ . For $D\ge 2$ . Two difficulties arise.

First, the exponent contains the drift term $\left(\mu -1/2{\sigma}^{2}\right)t$ , which means that the constraints themselves contain drift. This can easily be dealt with using Girsanov’s theorem, or alternatively (and equivalently) by including drift in the heat equation. This is detailed in Appendix B.

Second, linear constraints on asset prices $P$ do not translate exactly into linear constraints for the Wiener process. This is because linear constraints for prices of the form ${a}^{\text{T}}P\le y$ imply a non-linear constraint on the exponents, that is ${a}^{\text{T}}{\text{e}}^{\left(\mu -1/2{\sigma}^{2}\right)t+\Sigma x}\le y$ . This is dealt with next.

7.1. Non-Linear Constraints for Wiener Processes

Where the time period involved does not allow the exponential of a Wiener process to be approximated, two approaches can be considered.

One of these is direct―to solve the heat equation with non-linear constraints. The random walk of Section 3.1 shows this can be handled numerically and efficiently.

The second approach is to approximate the non-linear constraints with an envelope of linear ones. The example of Section 3.1 shows how this might be feasible computationally.

7.2. Example

[2] use a similar technique to that employed above for finding analytic solutions for exotic option pricing for $D=1$ . The technique is also called `the method of images’. However the reason for this is less obvious than that of this paper as the authors consider, not the PDE for the Wiener process, but that for the option price V under the Black Scholes PDE, namely

$\partial V/\partial \tau +1/2{\sigma}^{2}{S}^{2}{\partial}^{2}V/\partial {S}^{2}+rS\partial V/\partial S-rV=0$

where r and ${\sigma}^{2}$ relate to the characteristics of the stock price S and $\tau =T$ is the time of expiration.

Its is straightforward to show this PDE can be reduced to the heat equation ${u}_{t}=1/2{u}_{xx}$ for $D=1$ using a transformation of the form

$S=K{\text{e}}^{x}$

$t=\left(T-\tau \right){\sigma}^{2}$

$V=Ku{\text{e}}^{-\left({\alpha}^{2}/4+\alpha +1\right)t-\alpha x/2}$

$\alpha =2r/{\sigma}^{2}-1.$

This can be accomplished in many ways by choosing the value of the constant K; this possibility leads to invariance of the PDE under the images involved.

Consider the example where the option has a payoff of $f\left(P\right)$ at time $t=T$ provided the spot price does not exceed the ceiling H. The spot price can be solved to give $P=S{\text{e}}^{\lambda t+\sigma x}$ where $\lambda =r-1/2{\sigma}^{2}$ and S is its initial price. The value of the payoff $f\left(P\right)$ at time T is

$C={\displaystyle \underset{x\le y\le lnH/\sigma}{\int}}f\left(S{\text{e}}^{\lambda T+\sigma x}\right)u\left(x\mathrm{,}y\right)\text{d}x\text{d}y$

The integral $\underset{y\le lnH/\sigma}{\int}}u\left(x\mathrm{,}y\right)\text{d}y$ is just the density of the Wiener process where

the maximum is $lnH/\sigma $ . From 6 this is given $1/\sqrt{2\text{\pi}t}\left[{\text{e}}^{-{x}^{2}/2T}-{\text{e}}^{-{\left(2lnH/\sigma -x\right)}^{2}/2T}\right]$ . Hence the expression for C becomes

$C=1/\sqrt{2\text{\pi}t}{\displaystyle \underset{x\le y\le \mathrm{ln}H/\sigma}{\int}}f\left(S{\text{e}}^{\lambda T+\sigma x}\right)\left[{\text{e}}^{--{x}^{2}/2T}-{\text{e}}^{-{\left(2lnH/\sigma -x\right)}^{2}/2T}\right]\text{d}x\text{d}y$

Now bring in a further transformation to remove drift: let $\lambda t+\sigma x=\sigma z$ for a new Wiener process z. The transformation may be written $z=x+\lambda /\sigma t$ , so that the change of measure under Girsanov’s theorem from x to z is ${\text{e}}^{\lambda /\sigma z-12{\left(\lambda /\sigma \right)}^{2}T}$ .

The above integral may then be written as

$\begin{array}{c}C={\text{e}}^{-rT}/\sqrt{2\text{\pi}T}[{\displaystyle \underset{x\le \mathrm{ln}H/\sigma}{\int}}f\left(S{\text{e}}^{\sigma z}\right){\text{e}}^{--{x}^{2}/2T}{\text{e}}^{\lambda /\sigma z-1/2{\left(\lambda /\sigma \right)}^{2}T}\text{d}z\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{\displaystyle \underset{x\le \mathrm{ln}H/\sigma}{\int}}f\left(S{\text{e}}^{\sigma z}\right){\text{e}}^{-{\left(2\mathrm{ln}H/\sigma -z\right)}^{2}/2T}{\text{e}}^{\lambda /\sigma z-1/2{\left(\lambda /\sigma \right)}^{2}T}\text{d}z]\\ ={\text{e}}^{-rT}/\sqrt{2\text{\pi}T}[{\displaystyle \underset{x\le \mathrm{ln}H/\sigma}{\int}}f\left(S{\text{e}}^{\sigma z}\right){\text{e}}^{-{\left(z-\lambda /\sigma T\right)}^{2}/2T}\text{d}z\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{H}^{2\lambda /{\sigma}^{2}}{\displaystyle \underset{x\le \mathrm{ln}H/\sigma}{\int}}f\left(S{\text{e}}^{\sigma z}\right){\text{e}}^{-{\left(z-2\mathrm{ln}H/\sigma -\lambda /\sigma T\right)}^{2}/2T}\text{d}z]\\ ={\text{e}}^{-rT}/\sqrt{2\text{\pi}T}[{\displaystyle \underset{x\le \mathrm{ln}H/\sigma}{\int}}f\left(S{\text{e}}^{\sigma z}\right){\text{e}}^{-{\left(z-\lambda /\sigma T\right)}^{2}/2T}\text{d}z\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{H}^{2\lambda /{\sigma}^{2}}{\displaystyle \underset{x\le \mathrm{ln}H/\sigma}{\int}}f\left(S{H}^{2}{\text{e}}^{\sigma z}\right){\text{e}}^{-{\left(z-\lambda /\sigma T\right)}^{2}/2T}\text{d}z]\end{array}$

Thus the value of the path-restricted call option is thus given by the difference between two path-unrestricted options. This result is consistent with [2] , who refer to the second integral as an image of the first. Clearly similar examples involving both extrema can be found.

8. The Case of Steel

Another application of a multi-dimensional Wiener process is in the case of commodities, e.g. steel (s). The major inputs are iron ore (i), coal (c), labor and energy. The first two inputs and the output may be measured in USD. These outputs may be analyzed by the Augmented Dickey Fuller test in isolation to be stationary, i.e. $I\left(0\right)$ , except possibly for the labor and energy inputs. For simplicity, we assume that the full process is stationary; time trends can be be dealt with by introducing drift as in Appendix B. The simplification is that commodities do not have an intrinsic return, compared to financial assets.

Data was provided for 600 observations during the period 1/26/2016 to 8/8/2016 [11] . The observed means and the covariance matrix of $SIC=\left(s,i,c\right)$ are as follows:

$\stackrel{\xaf}{SIC}=\left[\begin{array}{c}495.243\\ 66.560\\ 83.314\end{array}\right]$

and

$\Sigma =\left[\begin{array}{ccc}14,776& 696& 1888\\ 696& 111& 109\\ 1888& 109& 375\end{array}\right].$

The commodities can be modelled as a three dimensional Wiener process $SIC=Ax$ with $x$ being a 3 dimensional Wiener process, and with $A$ being given by the Cholesky decomposition of $A$ :

$A=\left[\begin{array}{ccc}0.043& 1.595& 121.544\\ -8.025& -3.596& 5.776\\ 2.621& -11.038& 15.681\end{array}\right]$

For an example of an exotic option, consider a down-and-out option, where SIC is subject to the particular floors over a one year period as follows:

$Ax+\stackrel{\xaf}{SIC}\le \left[\begin{array}{c}505.243\\ 76.560\\ 93.314\end{array}\right]\mathrm{.}$ (7)

Since the commodities are stationary in their own right, the price of a call option at 500 is given by

$Call={\displaystyle \int}\stackrel{\u02dc}{f}\left[s\right]max\left(\mathrm{0,}s-500\right)\text{d}s$

where $\stackrel{\u02dc}{f}\left[s\right]$ is the density of the steel price if the variable SIC obeys the floors, or equivalently the constraints 7. Thus it remains to assess the density $f\left(x\right)$ under the constraints, and when s exceeds 500 in its payoff.

The poles to be assessed under the constrained process are the sources:

$\mathcal{P}=\left\{\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right]\right\}$

and the sinks^{1}

$\mathcal{Q}=\left\{\left[\begin{array}{c}0.6129\\ -113.7736\\ 37.1549\end{array}\right]\mathrm{,}\left[\begin{array}{c}1.5069\\ 0.0716\\ 0.1944\end{array}\right]\right\}$

Thus the density of x is

$f\left(x\right)=\frac{1}{{\left(2\text{\pi}\right)}^{3/2}}\left[{\displaystyle \underset{p\in \mathcal{P}}{\sum}}\text{\hspace{0.05em}}{\text{e}}^{-\frac{{\left(x-p\right)}^{2}}{2}}-{\displaystyle \underset{q\in \mathcal{Q}}{\sum}}\text{\hspace{0.05em}}{\text{e}}^{-\frac{{\left(x-q\right)}^{2}}{2}}\right]\mathrm{.}$

^{1}Note that poles differing by more than 10^{−9} have been eliminated on the basis of computational inaccuracy.

From this expression the the density of $s={e}_{1}^{\text{T}}Ax$ can be computed. where ${e}_{1}^{\text{T}}=\left[\mathrm{1,0,0}\right]$ . If its steel (and its labor and energy inputs) were subject to time

trends, then drift as in Appendix B can be employed using a factor ${\text{e}}^{\kappa -\frac{1}{2}{\kappa}^{2}}$ for the trend $\kappa $ in steel s.

We now compute the density of $SIC=Ax$ . Since s has an historic mean of 495, we require only that ${e}_{1}^{\text{T}}Ax\ge 10$ . The value of the barrier option is then

$\frac{1}{{\left(2\text{\pi}\right)}^{3/2}}{\displaystyle \int}f\left(x\right)max\left[\mathrm{0,}{e}_{1}^{\text{T}}Ax-10\right]\text{d}x=\mathrm{0.991.}$

It is more convenient to evaluate this integral numerically in the result above.

In contrast the option price without constraints is

$\begin{array}{l}\frac{1}{\sqrt{2\text{\pi}14776}}{\displaystyle \underset{s>500}{\int}}\left[s-15.6224\right]{\text{e}}^{-\frac{{\left(s-495.2\right)}^{2}}{2\times 14776}}\text{d}s\\ =\sqrt{\frac{14776}{\text{\pi}}}{\text{e}}^{-\frac{{\left(500-495.2\right)}^{2}}{2\times 14776}}-\left(500-495.2\right)\Phi \left(1-\frac{495.2}{\sqrt{14776}}\right)=\mathrm{68.523.}\end{array}$

Thus the imposition of floors on the inputs during the one year process has a strong effect in reducing the variability of the steel price, and thus narrowing the option price. Of course more extreme examples are possible.

9. Conclusion

Though the connection between the Wiener process and the heat equation has long been known, there has been little research to exploit it. The context of constrained processes provides a natural setting for deriving computational, and even analytic, solutions for these problems. The power of group theoretic methods, embodied in Kelvin’s method of images, is only touched on in this paper, but it is evident they have a much greater role to play in mathematical finance.

Appendix A: Convergence of the Random Walk to the Wiener Process

The recurrence relation in (1) after 2 steps may be written as follows

${f}_{nt+2}\left(x\sqrt{n/D}\right)\mathrm{=}1/{\left(2D\right)}^{2}{\displaystyle \underset{i\mathrm{,}j}{\sum}}\text{\hspace{0.05em}}{f}_{nt}\left(x\sqrt{n/D}\pm {e}_{i}\pm {e}_{j}\right)$

or

$\begin{array}{l}{f}_{nt+2}\left(x\sqrt{n/D}\right)-{f}_{nt}\left(x\sqrt{n/D}\right)\\ =1/{\left(2D\right)}^{2}\left[{\displaystyle \underset{i\ne j}{\sum}}{f}_{nt}\left(x\sqrt{n/D}\pm {e}_{i}\pm {e}_{j}\right)-{f}_{nt}\left(x\right)+{\displaystyle \underset{i}{\sum}}{f}_{nt}\left(x\sqrt{n/D}\pm 2{e}_{i}\right)-{f}_{nt}\left(x\right)\right]\end{array}$

Denoting ${\psi}_{t}\left(x\right)={f}_{nt}\left(x\sqrt{n/D}\right)$ and $h=\sqrt{D/n}$ this may be expressed as

$\begin{array}{l}{\psi}_{t+2{h}^{2}/D}\left(x\right)-{\psi}_{t}\left(x\right)/2{h}^{2}/D\\ =1/16D\left[{\displaystyle \underset{i\ne j}{\sum}}\text{\hspace{0.05em}}{\psi}_{t}\left(x\pm h{e}_{i}\pm h{e}_{j}\right)-{\psi}_{t}\left(x\right)/{h}^{2}+{\displaystyle \underset{i}{\sum}}\text{\hspace{0.05em}}{\psi}_{t}\left(x\pm 2h{e}_{i}\right)-{\psi}_{t}\left(x\right)/{h}^{2}\right]\end{array}$ (8)

The Taylor series expansions are:

${\psi}_{t}\left(x\pm h{e}_{i}\pm h{e}_{j}\right)-{\psi}_{t}\left(x\right)=\pm h{e}_{i}^{\text{T}}\Delta {\psi}_{t}\pm h{e}_{j}^{\text{T}}\Delta {\psi}_{t}+1/2{h}^{2}{e}_{i}^{\text{T}}{\Delta}^{2}{\psi}_{t}{e}_{j}+O(h3)$

${\psi}_{t}\left(x\pm 2h{e}_{i}\right)-{\psi}_{t}\left(x\right)=\pm 2h{e}_{i}^{\text{T}}\Delta {\psi}_{t}+2{h}^{2}{e}_{i}^{\text{T}}{\Delta}^{2}{\psi}_{t}{e}_{i}+O(h3)$

There are $2D\left(2D-2\right)$ of the first of the above relations, and 2D of the second. In the above summations, the aggregate coefficients of ${e}_{i}^{\text{T}}\Delta {\psi}_{t}$ and ${e}_{i}^{\text{T}}{\Delta}^{2}{\psi}_{t}{e}_{j}$ for $i\ne j$ are zero, as for any term $x\pm h{e}_{i}\pm h{e}_{j}$ there is the opposing term $x\mp h{e}_{i}\mp h{e}_{j}$ . The only non-zero terms are ${e}_{i}^{\text{T}}{\Delta}^{2}{\psi}_{t}{e}_{i}$ , each with total coefficient $8\left(D-1\right)+8=8D$ . Thus as $n\to \infty $ , or equivalently $h\to 0$ , we have the limit $\phi =li{m}_{h\to 0}{\psi}_{t}\left(x\right)$ and

$\partial \phi /\partial t=1/2tr\left({\nabla}^{2}\phi \right)=1/2\Delta \phi $

which is the heat equation in D dimensions. The constraints $Ax\le y$ are similarly scaled as $Ax\sqrt{n/D}\le y\sqrt{n/D}$ .

Appendix B: Allowance for Drift in Wiener Processes

The stochastic approach to introducing (or removing) drift into Wiener processes is via Girsanov’s theorem, which allows the use of a Radon Nikodỳm change of measure for the process. Since we are dealing directly with densities under the PDE approach to such processes, it is not surprising that a more direct treatment is possible under this approach.

This can be achieved by considering the heat equation with linear drift, namely

${v}_{t}+{\kappa}^{\text{T}}\nabla v=1/2\Delta v$

where $\kappa \in {\mathbb{R}}^{D}$ is a specific drift vector. Thus if $u\left(t\mathrm{,}x\right)$ is a solution of the usual heat equation ${u}_{t}=1/2\Delta u$ , then $v\left(t\mathrm{,}x\right)=u\left(t\mathrm{,}x+\kappa t\right)$ is the solution to that with drift. It also follows that if u is constrained to vanish on surfaces of the form ${a}^{\text{T}}x=y$ , then must v vanish on the drift-adjusted surfaces ${a}^{\text{T}}\left(x-\kappa t\right)=y$ . Thus we can solve problems where the constraints themselves drift with time.

Fortunately there is an easy adjustment that allows us to connect the solutions of heat equations with or without drift. If u is a solution of the driftless equation, then consider the function

$v={\text{e}}^{-{\Vert x-\kappa t\Vert}^{2}/2t}{\text{e}}^{-{\Vert x\Vert}^{2}/2t}u={\text{e}}^{{x}^{\text{T}}\kappa -1/2{\Vert \kappa \Vert}^{2}t}u\mathrm{.}$

We have

${v}_{t}={\text{e}}^{{x}^{\text{T}}\kappa -1/2{\Vert \kappa \Vert}^{2}t}\left[-1/2{\Vert \kappa \Vert}^{2}u+{u}_{t}\right]$

${v}_{x}={\text{e}}^{{x}^{\text{T}}\kappa -1/2{\Vert \kappa \Vert}^{2}t}\left[\kappa u+{u}_{x}\right]$

${v}_{xx}={\text{e}}^{{x}^{\text{T}}\kappa -1/2{\Vert \kappa \Vert}^{2}t}\left[{\Vert \kappa \Vert}^{2}u+2{\kappa}^{\text{T}}{u}_{x}+{u}_{xx}\right]$

and thus

${v}_{t}+{\kappa}^{\text{T}}\nabla v\mathrm{=}1/2\Delta v$

satisfies the heat equation with drift. The factor ${\text{e}}^{{x}^{\text{T}}\kappa -1/2{\Vert \kappa \Vert}^{2}t}$ is precisely the change of measure prescribed by Girsanov’s theorem, but here it appears as a transformation for linear drift. Clearly this result may be generalized to non-constant drift.

References

[1] Doob, J.L. (1955) A Probability Approach to the Heat Equation. Transactions of the American Mathematical Society, 216-280. https://doi.org/10.1090/S0002-9947-1955-0079376-0

[2] Buchen, P. and Konstandatos, O. (2009) A New Approach to Pricing Double-Barrier Options. Applied Mathematical Finance, 16, 497-515. https://doi.org/10.1080/13504860903075480

[3] Marchetti, D. and da Silva, R. (1999) Brownian Motion Limit of Random Walks in Symmetric Nonhomogeneous Media. Brazilian Journal of Physics.

[4] Sommerfeld, A. (1949) Partial Differential Equations. Academic Press.

[5] Bluman, G.W. and Cole, J. (1969) The General Similarity Solution of the Heat Equation. Journal of Mathematics and Mechanics, 18.

[6] Freedman, D. (1983) Brownian Motion and Diffusion. Springer Verlag.

https://doi.org/10.1007/978-1-4615-6574-1

[7] Dym, H. and McKean, H.P. (1972) Fourier Series and Integrals. Academic Press.

[8] Coxeter, H. (1935) The Complete Enumeration of Finite Groups of the Form. Journal of the London Mathematical Society, s1-10.

[9] Muirhead, S. (2011) Pricing Multi-Asset Barrier Options Using the Generalised Reflection Principle. Master’s Thesis, Mathematics and Statistics, University of Melbourne.

[10] Shepp, L.A. (1979) The Joint Density of the Maximum and Its Location for a Wiener Process with Drift. Journal of Applied Probability, 16.

[11] Trading Economics Group. Steel.

https://tradingeconomics.com/analytics/features.aspx?source=/commodity/steel