The Lognormal Characteristic Function in Several Dimensions, with Application to Asian Options

Andrew P. Leung

Independent Researcher, Melbourne, Australia.

DOI: 10.4236/jmf.2020.103024                

Abstract

Many attempts have been made to derive a simple expression for the characteristic function of the lognormal distribution. This would be useful for computing the sum of lognormal variables, either with each other, or with other statistical variables. In this paper, we provide a simple formula for the characteristic function, which is exact, closed and computable. An extension to the sum of correlated lognormals, used in the pricing of Asian options, is a consequence of this approach.

Keywords

Lognormal, Characteristic Function, Multidimensional, Complex Functions, Asian Options

1. Introduction

The usefulness of finding a simple expression for the characteristic function (CF, a.k.a Fourier transform FT) of the lognormal distribution has been made clear in a variety of contexts: from insurance and finance to engineering to biology. The difficulty is that the CF has an essential singularity at 0, and therefore a branch cut must be made [1]. [2] was one of the first to address the density of total claims over an infinite time interval, which bypasses a need for the CF. The failure in finding a simple expression for the CF has led to the study of asymptotic approximations using orthogonal polynomials [3] [4] [5] or higher order functions [6] [7]. Truncated series expansions of Laplace transforms, which are a form of approximation, are entertained in [8]. A recent survey of approximation methods may be found in [9]. A critique of some of these papers is provided below in context.

Thus, we have the conclusions in [6], while citing [10]:

“A closed form expression of the Laplace transform of the lognormal distribution does not exist.”

and further in [7]:

“Recently, Asmussen et al. (2016a) reported a closed form approximation for the Laplace transform … that works quite well for small values of the $\sigma$ parameter and is asymptotically equivalent to the approximated Laplace transform.”

and from [2]:

“The usual method of evaluating the distribution function requires the computation of many convolutions of the conditional distribution of the amount of a claim given that a claim has occurred. When the expected number of claims is large, the computation can become unwieldy even with modern large scale electronic computers.”

The last quotation relates to convolutions, which are obviated if the CF of the claims amount distribution can be found.

In the context of Asian options, we find a similar view from [11]:

“Contrary to the case of geometric averaging, the distribution for the arithmetic average is not known. This implies that no exact closed form solution is available for European-style Asian options based on the arithmetic average. Because no analytical formula is available for European-style Asian options, it does not come as a surprise that an analytical solution to the similar American-style problem is not feasible either.”

Nonetheless, when viewed from complex function theory, the lognormal is a well behaved function; in fact it is entire (i.e. integral, with no finite singularities) in the complex plane. [10] exploits this property to derive a simple and computable formula for the CF for part of the complex plane. In this paper, we extend the formula to the whole complex plane. The closest approach, using complex function theory is that of [12], but even he resorts to asymptotic expansions.

In this paper, we provide and prove an exact expression for the lognormal CF, which can be easily verified numerically by the inverse Fourier transform. It is not the purpose of this paper to assess how it relates to the approximations, asymptotic or otherwise, in the literature. In fact, that is a task for the authors of those approximations.

2. Notation

For a variance parameter ${\sigma }^2$ and mean $\mu$ the pdf of the lognormal distribution may be written as

$p\left(x\right)=\frac{1}{\sqrt{2\pi {\sigma }^2}x}{\text{e}}^{-\frac{{\left(\ln x+\mu \right)}^2}{2{\sigma }^2}}$

(where $\mu =\frac{{\sigma }^2}{2}$, for x to have unit mean, which is assumed without loss of generality). Its FT is using $x={\text{e}}^{\zeta }$

$\hat{p}\left(z\right)=\frac{1}{\sqrt{2\pi {\sigma }^2}}{\displaystyle \underset{\zeta \in R}{\int }}{\text{e}}^{iz{\text{e}}^{\zeta }}{\text{e}}^{-\frac{{\left(\zeta +\mu \right)}^2}{2{\sigma }^2}}\text{d}\zeta$ (FT)

which is a difficult integral to evaluate numerically because of the presence of the oscillatory term ${\text{e}}^{iz{\text{e}}^{\zeta }}$, which grows exponentially for real z.

Remark 1 There is considerable ambiguity as to the meaning of “closed form” for a mathematical expression. In its strictest sense, it admits only arithmetic expressions. In its most general sense, it admits the operations of integration and differentiation, as many functions (e.g. the gamma function) need to be defined in terms of these operations. However the term precludes algorithms and numerical solutions of a system of equations.

3. Cauchy Integration

A simple approach, introduced by [10], resolves this difficulty. As the integrand in FT ${\text{e}}^{iz{\text{e}}^{\zeta }}{\text{e}}^{-\frac{{\left(\zeta +\mu \right)}^2}{2{\sigma }^2}}$ is an entire function of $\zeta$ in the complex plane, we can replace it with an integral along the contour $\zeta$ in the upper half plane bounded by

· the horizontal segments $x=\pm X$, and

· the vertical segments $y=0$ and $y=\frac{\pi }{2}$.

Alternatively we can use the contour in the lower half plane bounded by

· the horizontal segments $x=\pm X$, and

· the vertical segments $y=0$ and $-\frac{\pi }{2}$

This is depicted in Figure 1. Using the Cauchy theorem for these contours, it is possible to replace the integral in FT with integrals along the horizontal segments and allow $\left|X\right|\to \infty$.

This is justified since the lognormal is an entire function of $\zeta$, provided we show that the integrals on the vertical segments tend to zero as $\left|X\right|\to \infty$. We distinguish between contours in the upper and lower half planes, as is seen below.

For the upper half plane, the integral along the vertical segment for $0\le y\le \frac{\pi }{2}$, given X and z is:

$\frac{1}{\sqrt{2\pi {\sigma }^2}}{\displaystyle \underset{0\le y\le \frac{\pi }{2}}{\int }}\underset{\text{oscillatory}}{\underset{︸}{{\text{e}}^{iz{\text{e}}^{X+iy}}}}\underset{\text{dampening}}{\underset{︸}{{\text{e}}^{-\frac{{\left(X+iy+\mu \right)}^2}{2{\sigma }^2}}}}\text{d}y$ (1)

Now the dampening term in 1 has the real part:

$\Re {\left(X+\mu +iy\right)}^2=\frac{1}{2}\left[{\left(X+\mu +iy\right)}^2+{\left(X+\mu -iy\right)}^2\right]={\left(X+\mu \right)}^2-y^2.$

Thus the dampening term has no definitive influence on the convergence of the vertical integral as $X\to \pm \infty$. The oscillatory term can be dealt with as follows.

Figure 1. Contour integrals.

3.1. The Case $\Re \left(z\right)\ge 0$

The integral is over a finite range for y in the upper half plane and tends to 0 if $\Re \left(iz{\text{e}}^{iy}\right)<0$, as with $z=a+ib$ :

$\Re \left(iz{\text{e}}^{iy}\right)=\Re \left(a+ib\right)\left(i\cos y-\sin y\right)=-a\sin y-b\cos y.$

If $0\le y\le \frac{\pi }{2}$, so that $siny\mathrm{,}cosy\ge 0$, then $\Re \left(iz{\text{e}}^{iy}\right)\le 0$ for $a=\Re \left(z\right)\ge 0$, $b=\Im \left(z\right)\ge 0$.

3.2. The Case $\Re \left(z\right)\le 0$

The integral is over a finite range for y in the lower half plane and again tends to 0 if $\Re \left(iz{\text{e}}^{iy}\right)<0$. For $-\frac{\pi }{2}\le y\le 0$, $siny\le 0$ and $cosy\ge 0$, this applies to $\Re \left(z\right)\le 0$, $\Im \left(z\right)\ge 0$.

Thus for $\Im \left(z\right)\ge 0$ we have in the limit as $\left|X\right|\to \infty$, we derive an alternative expression for the lognormal CF:

$\hat{p}\left(z\right)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\{\begin{array}{ll}{\displaystyle \underset{\zeta \in R}{\int }}{\text{e}}^{-z{\text{e}}^{\zeta }}{\text{e}}^{-\frac{{\left(\zeta +\mu +\frac{\pi }{2}i\right)}^2}{2{\sigma }^2}}\text{d}\zeta \hfill & \text{if}\Re \left(z\right)\ge 0\hfill \\ {\displaystyle \underset{\zeta \in R}{\int }}{\text{e}}^{z{\text{e}}^{\zeta }}{\text{e}}^{-\frac{{\left(\zeta +\mu -\frac{\pi }{2}i\right)}^2}{2{\sigma }^2}}\text{d}\zeta \hfill & \text{if}\Re \left(z\right)\le 0\hfill \end{array}$ (2)

But the above expressions for $\hat{p}\left(z\right)$ may be analytically continued to $\Im \left(z\right)<0$ within the indicated ranges of $\Re \left(z\right)$, as the integrals converge, regardless of the sign of $\Im \left(z\right)$ ¹. It is also evident that $\hat{p}\left(z\right)$ then obeys the general equality for FTs of a real distribution:

$\overline{\stackrel{^}{p}\left(z\right)}=\hat{p}\left(-\bar{z}\right)\mathrm{.}$

For convenience, denote $\tilde{z}=-\text{sign}\left[\Re \left(z\right)\right]z$, so that 2 may be written succinctly as:

$\hat{p}\left(z\right)=\frac{1}{\sqrt{2\pi {\sigma }^2}}{\displaystyle \underset{\zeta \in R}{\int }}{\text{e}}^{\tilde{z}{\text{e}}^{\zeta }}{\text{e}}^{-\frac{{\left(\zeta +\mu +\frac{\pi }{2}i\right)}^2}{2{\sigma }^2}}\text{d}\zeta$ (3)

and so

$p\left(x\right)=\frac{1}{{\left(2\pi \right)}^{3/2}\sigma }{\displaystyle \iint {\text{e}}^{-izx}{\text{e}}^{\tilde{z}{\text{e}}^{\zeta }}{\text{e}}^{-\frac{{\left(\zeta +\mu +\frac{\pi }{2}i\right)}^2}{2{\sigma }^2}}\text{d}\zeta \text{d}z}$

Remark 2 It has been suggested (anonymously) that Equation (2) is equivalent to the formula in [13], written in our notation as

$\hat{p}\left(z\right)=\frac{1}{2\pi i}{\displaystyle \underset{k-i\infty }{\overset{k+i\infty }{\int }}}\Gamma \left(s\right){\text{e}}^{-\mu +\ln \left(z\right)-\text{sign}\left(z\right)\frac{i\pi }{2}+\frac{{\sigma }^2{s}^2}{2}}\text{d}s\text{for}z\in R\backslash \left\{0\right\}$ (4)

where k is any number $k>0$ and $\Gamma \left(s\right)$ is presumably the gamma function. This alleged equivalence is clearly false:

· The term for $\Gamma \left(s\right)$ does not appear in 2;

· The above expression is an integral over a vertical strip, unlike the horizontal strip in 2;

· Moreover the author does not exploit or attempt to compute the expression in 4; instead he resorts to asymptotic approximations.

Remark 3 The approach of [2], for Poisson distributed claim frequency, is to solve the following integral equation

$g\left(x\right)=\lambda {\text{e}}^{-\lambda }p\left(x\right)+\frac{\lambda }{x}{\displaystyle \underset{0}{\overset{x}{\int }}}yp\left(y\right)g\left(x-y\right)\text{d}y$ (5)

where $g\left(x\right)$ is the density of total claims x over a unit period and $\lambda$ a Poisson parameter. The equation was developed to avoid convolutions in computing $g\left(x\right)={\text{e}}^{-\lambda }{\displaystyle \underset{n=0}{\sum }}\frac{{\lambda }^n}{n!}{p}^{\ast n}$. In fact the FT has a simple form $\hat{g}\left(z\right)={\text{e}}^{\lambda \left(\hat{p}-1\right)}$, if we can derive the FT $\hat{p}\left(z\right)$, and use the inverse FT. However, given the form of $p\left(x\right)$, 5 has a very complicated kernel. It is claimed (anonymously) that this is “state of the art” for insurance purposes. However:

· the equation is an inhomogeneous Volterra integral equation of the 2nd kind, the most intractable possible. Since it applies to any claim amount density $p\left(x\right)$, it cannot be used to solve for $\hat{p}\left(x\right)$ unless other conditions are imposed.

· Generally it needs to be solved by discretization of the variable x and numerical quadrature of a system of equations, which renders the approach only slightly less problematic than the convolution approach;

· No solution, nor an approach to one, is offered, approximate or otherwise;

· It does not and cannot deal in any way with the CF of the lognormal, the theme of this paper, which is useful and important in fields other than insurance, and for a multitude of different purposes (e.g. electrokinetics).

4. Computation

There is clearly a branch cut along the imaginary axis (in fact, two branch cuts originating from 0 along the imaginary axis). By taking the inverse transform $\frac{1}{2\pi }{\displaystyle \int {\text{e}}^{-izx}\hat{p}\left(z\right)\text{d}z}$, we should recover $p\left(x\right)$. This is confirmed in Figure 2, where the p, the lognormal pdf (the black line) corresponds closely to the inverse FT for $\hat{p}$.

We calculate the FT of the lognormal average, where each $u_i$ is lognormally distributed, such as:

$v=\frac{u_1+u_2+\cdots +u_n}{n}=\frac{U}{n}.$ (6)

Taking averages allows us to compare the effect of convolution for various levels of n. The lognormal sum U has the CF $\hat{p}{\left(z\right)}^n$ and thus the average lognormal v has CF $\hat{p}{\left(z/n\right)}^n$. This simplifies greatly the procedures discussed in [5].

The FTs of $v\left(x\right)$ and that of $U\left(x\right)$ may be written and are related by

$v\left(x\right)=\frac{1}{2\pi }{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}}{\text{e}}^{-ixz}\hat{v}\left(z\right)\text{d}z$ (7)

$=\frac{1}{2\pi }{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}}{\text{e}}^{-ixz}\hat{p}{\left(z/n\right)}^n\text{d}z=\frac{n}{2\pi }{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}}{\text{e}}^{-inxz}\hat{p}{\left(z\right)}^n\text{d}z$ (8)

$=\frac{n}{{\left(2\pi \right)}^{n/2+1}{\sigma }^n}\{\begin{array}{cc}{\displaystyle \underset{0}{\overset{\infty }{\int }}}{\left[{\text{e}}^{-ixz}{\displaystyle \int }{\text{e}}^{-z{\text{e}}^{\zeta }}{\text{e}}^{-\frac{{\left(\zeta +\mu +\frac{\pi }{2}i\right)}^2}{2{\sigma }^2}}\text{d}\zeta \right]}^n\text{d}z& \text{if}\Re \left(z\right)\ge 0\\ {\displaystyle \underset{-\infty }{\overset{0}{\int }}}{\left[{\text{e}}^{-ixz}{\displaystyle \int }{\text{e}}^{z{\text{e}}^{\zeta }}{\text{e}}^{-\frac{{\left(\zeta +\mu -\frac{\pi }{2}i\right)}^2}{2{\sigma }^2}}\text{d}\zeta \right]}^n\text{d}z& \text{if}\Re \left(z\right)\le 0\end{array}$ (9)

so that

$v\left(x\right)=\frac{n}{2\pi }{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}}{\text{e}}^{-inxz}\hat{p}{\left(z\right)}^n\text{d}z$ (10)

$=\frac{2n}{{\left(2\pi \right)}^{n/2+1}{\sigma }^n}\Re {\displaystyle \underset{0}{\overset{\infty }{\int }}}{\left[{\text{e}}^{-ixz}{\displaystyle \int }{\text{e}}^{-z{\text{e}}^{\zeta }}{\text{e}}^{-\frac{{\left(\zeta +\mu +\frac{\pi }{2}i\right)}^2}{2{\sigma }^2}}\text{d}\zeta \right]}^n\text{d}z$ (11)

so that

$v\left(x\right)=\frac{n}{2\pi }{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}}{\text{e}}^{-inxz}\hat{p}{\left(z\right)}^n\text{d}z$ (12)

$=\frac{2n}{{\left(2\pi \right)}^{n/2+1}{\sigma }^n}\Re {\displaystyle \underset{0}{\overset{\infty }{\int }}}{\left[{\text{e}}^{-ixz}{\displaystyle \int }{\text{e}}^{-z{\text{e}}^{\zeta }}{\text{e}}^{-\frac{{\left(\zeta +\mu +\frac{\pi }{2}i\right)}^2}{2{\sigma }^2}}\text{d}\zeta \right]}^n\text{d}z$ (13)

Figure 2. Lognormal convolutions for various levels of n, with $\sigma =1$.

The charts of v for various values of n are shown in Figure 2 above. This provides a computable closed form expression for the lognormal CF and its convolutes. For each x, this requires only 2 integrations, in contrast to the integral Equation in (5).

Remark 4 The integrations used for Figure 2 were derived using Matlab and the Chebfun2 package. The form of the CF means that the number of convolutions does not hinder its efficiency. For the avoidance of doubt, Matlab and Chebfun do not offer solutions of 3, but merely facilitate integrations in 10, the Chebfun version being the more accurate and efficient. Integrations are inevitable, especially with the integral Equation in (5), and the task is to minimize the number required.

Remark 5 It is easy to see that the FT of a mixture of lognormals and other distributions, such as the normal, may be derived in the same way. Clearly the lognormals do not need to be identically distributed. Where the lognormals are correlated, it is possible to extend the one dimensional results to n dimensions, as is shown in the next section.

Remark 6 Where the sum for U in 6 has negative components, it is straightforward to allow for them. In this event U can assume negative values, so that the expression in 10 does not apply. To deal with this, consider the term

$\hat{p}\left(-z\right)=\frac{1}{\sqrt{2\pi {\sigma }^2}}{\displaystyle \underset{\zeta \in R}{\int }}{\text{e}}^{-iz{\text{e}}^{\zeta }}{\text{e}}^{-\frac{{\left(\zeta +\mu \right)}^2}{2{\sigma }^2}}\text{d}\zeta .$

Then 10 still applies, with the roles of $\Re \left(z\right)\ge 0$ and $\Re \left(z\right)\ge 0$ reversed. This may be regarded as a special case of the following section. This, and the next section, generalize the results in [14].

Remark 7 [9] also consider the technique of contour integration suggested by [10], with computation via Chebyshev polynomials. Their expression for the lognormal CF is

$\hat{p}\left(z\right)={\displaystyle \int }{\text{e}}^{-z\sqrt{2\pi }x}{\text{e}}^{-\frac{izx}{{\left(2\sigma \right)}^2}}{\text{e}}^{-\frac{1}{2}{x}^2}\text{d}x.$

· [10] applies only for $\Re \left(z\right)\ge 0$, and branch cuts at $z=0$ need to be made, so a simple integral would not suffice.

· The integral has a quadratic exponent, and is thus solvable analytically. This contradicts all the literature.

5. Multivariate Lognormal

For a vector $x$, the multivariate lognormal is the exponential of a multivariate normal of $x$, and thus has the form

$p\left(x\right)=\frac{{\text{e}}^{-\frac{1}{2}{\left(lnx\right)}^{\text{T}}{\Sigma }^{-1}\left(lnx\right)}}{\sqrt{{\left(2\pi \right)}^n\left|\Sigma \right|}{\displaystyle \prod }{x}_1{x}_2\cdots x_n}$

where $\Sigma$ is the covariance matrix of the elements of $x=\left[x_1,x_2,\cdots ,x_n\right]$, Since ${\Sigma }^{-1}$ is positive definite, it can be written as

${\Sigma }^{-1}=A^{\text{T}}A$

for some symmetric matrix $A$, so that

$p\left(x\right)=\frac{{\text{e}}^{-\frac{1}{2}{\Vert A\left(lnx\right)\Vert }^2}}{\sqrt{{\left(2\pi \right)}^n\left|\Sigma \right|}{\displaystyle \prod }{x}_1{x}_2\cdots x_n}$

and

$\hat{p}\left(z\right)=\frac{\left|\Sigma \right|}{\sqrt{{\left(2\pi \right)}^n}}{\displaystyle \int }\frac{1}{{\displaystyle \prod }{x}_1{x}_2\cdots x_n}{\text{e}}^{iz\cdot x}{\text{e}}^{\frac{1}{2}{\Vert A\left(lnx\right)\Vert }^2}\text{d}x$

and letting $x={\text{e}}^{\zeta }$ we get the simplification for $z\in C^n$ :

$\hat{p}\left(z\right)=\frac{\left|A\right|}{\sqrt{{\left(2\pi \right)}^n}}{\displaystyle \int }{\text{e}}^{iz\cdot {\text{e}}^{\zeta }}{\text{e}}^{-\frac{1}{2}{\Vert A\zeta \Vert }^2}\text{d}\zeta .$

here we have the same problem as in FT: where any component of $z=\left[z_1\mathrm{,}{z}_2,\cdots ,z_n\right]$ has a positive real part, there is an uncomfortable oscillation. For convenience, denote $\tilde{z}=-\text{sign}\left[\Re \left(z\right)\right]z$ for the multidimensional case, as suggested by 3.

Fortunately we can deal with this problem by applying contour integration, component by component, according to the sign of $\Re \left(z\right)$. There are 2ⁿ possible signs, so that there are a similar number of branches around $z=0$, and we have the following 3:

$\hat{p}\left(z\right)=\frac{\left|A\right|}{\sqrt{{\left(2\pi \right)}^n}}{\displaystyle \int }{\text{e}}^{-i\tilde{z}\cdot {\text{e}}^{\zeta }}{\text{e}}^{-\frac{1}{2}{\Vert A\zeta \Vert }^2}\text{d}\zeta .$

A further simplification can be made by writing $r=A\zeta$ so that

$\hat{p}\left(z\right)=\frac{1}{\sqrt{{\left(2\pi \right)}^n}}{\displaystyle \int }{\text{e}}^{i\tilde{z}\cdot {\text{e}}^{A^{-1}r}}{\text{e}}^{-\frac{1}{2}{\Vert r\Vert }^2}\text{d}r.$

Letting $A^{-1}=\left[a_1\mathrm{,}{a}_2,\cdots ,a_n\right]$, this becomes

$\begin{array}{c}\hat{p}\left(z\right)=\frac{1}{\sqrt{{\left(2\pi \right)}^n}}{\displaystyle \int }{\text{e}}^{i\tilde{z}\cdot {\text{e}}^{A^{-1}r}}{\text{e}}^{-\frac{1}{2}{\Vert r\Vert }^2}\text{d}r\\ =\frac{1}{\sqrt{{\left(2\pi \right)}^n}}{\displaystyle \int }{\text{e}}^{i\tilde{z}\cdot {\text{e}}^{a_1{r}_1}}{\text{e}}^{i\tilde{z}\cdot {\text{e}}^{a_2{r}_2}}\cdots {\text{e}}^{i\tilde{z}\cdot {\text{e}}^{a_{n}rn}}{\text{e}}^{-\frac{1}{2}{\Vert r\Vert }^2}\text{d}{r}_1\text{d}{r}_2\cdots \text{d}{r}_n\\ =\frac{1}{\sqrt{{\left(2\pi \right)}^n}}{\displaystyle \underset{j=1}{\overset{n}{\prod }}}\left[{\displaystyle \int }{\text{e}}^{i\tilde{z}\cdot {\text{e}}^{a_{j}r}}{\text{e}}^{-\frac{1}{2}{r}^2}\text{d}r\right]\mathrm{.}\end{array}$ (14)

since 14 is a separable function of $r$. For computation this requires n integrations for each value of $z$.

Once $\hat{p}\left(z\right)$ is calculated, the inverse FT may be used to compute $p\left(x\right)$ in the usual way. From there, the distribution of any function of $x$, including a weighted sum of the form $y=b^{\text{T}}x$ may be found. In particular, this may be found as follows.

If $q\left(y\right)$ denotes the density of the weighted sum of y, (which must be positive) then $Q\left(y\right)={\displaystyle \underset{y}{\overset{\infty }{\int }}}q\left(s\right)\text{d}s$ is its cdf, whose extremes are 0 and 1. Thus:

$Q\left(y\right)={\displaystyle \underset{0\le y\le b^{\text{T}}x}{\int }}p\left(x\right)\text{d}x$

and since $q\left(y\right)=-\frac{\partial Q}{\partial y}$ :

$\hat{q}\left(z\right)=-{\displaystyle \int }{\text{e}}^{iyz}\frac{\partial Q}{\partial y}\text{d}y=-{\text{e}}^{iyz}Q|+iz{\displaystyle \int }{\text{e}}^{iyz}Q\text{d}y=1+iz\hat{Q}(\; z\; )$

$\hat{q}\left(z\right)=\frac{1}{\sqrt{{\left(2\pi \right)}^n}}\frac{\partial }{\partial z}{\displaystyle \underset{j=1}{\overset{n}{\prod }}}\left[{\displaystyle \int }{\text{e}}^{i\tilde{z}b\cdot {\text{e}}^{a_{j}r}}{\text{e}}^{-\frac{1}{2}{r}^2}\text{d}r\right]\mathrm{.}$

and the density of y is $q=-\frac{\partial Q}{\partial y}$, and its FT is

$\begin{array}{c}\hat{q}=-{\displaystyle \int }{\text{e}}^{iyz}\frac{\partial Q}{\partial y}\text{d}y=-{{\text{e}}^{iyz}Q|}_0^{\infty }+{\displaystyle \int }iz{\text{e}}^{iyz}Q\text{d}y\\ =1+{\displaystyle \int }iz{\text{e}}^{iyz}Q\text{d}y=1-{\displaystyle \int }\left[{\text{e}}^{i\left(b^{\text{T}}x\right)z}-1\right]p\left(x\right)\text{d}x\\ ={\displaystyle \int }{\text{e}}^{i\left(b^{\text{T}}x\right)z}p\left(x\right)\text{d}x=p\left(zb\right)=\frac{1}{\sqrt{{\left(2\pi \right)}^n}}{\displaystyle \int }{\text{e}}^{i\widetilde{zb}\cdot {\text{e}}^{A^{-1}r}}{\text{e}}^{-\frac{1}{2}{\Vert r\Vert }^2}\text{d}r.\end{array}$

Taking the inverse FT, we get finally:

$\hat{p}\left(z\right)=\frac{1}{\sqrt{{\left(2\pi \right)}^n}}{\displaystyle \int }{\text{e}}^{i\tilde{z}\cdot {\text{e}}^{A^{-1}r}}{\text{e}}^{-\frac{1}{2}{\Vert r\Vert }^2}\text{d}r.$

$p\left(y\right)=\frac{1}{\sqrt{{\left(2\pi \right)}^n}}{\displaystyle \int }{\text{e}}^{iyz}{\text{e}}^{i\widetilde{zb}\cdot {\text{e}}^{A^{-1}r}}{\text{e}}^{-\frac{1}{2}{\Vert r\Vert }^2}\text{d}r\text{d}z\mathrm{.}$ (15)

and

$\begin{array}{c}\hat{Q}\left(z\right)={\displaystyle \underset{0\le y\le b^{\text{T}}x}{\int }}{\text{e}}^{iyz}p\left(x\right)\text{d}x\text{d}y\\ =\frac{1}{iz}{\displaystyle \int }\left[{\text{e}}^{i\left(b^{\text{T}}x\right)z}-1\right]p\left(x\right)\text{d}x\\ =\frac{1}{\sqrt{{\left(2\pi \right)}^n}}\frac{1}{iz}{\displaystyle \int }\left[{\text{e}}^{i\left(b^{\text{T}}x\right)z}\right]p\left(x\right)\text{d}x-\frac{1}{iz}\end{array}$

and so using the inverse FT:

$\begin{array}{c}Q\left(y\right)=\frac{1}{\sqrt{{\left(2\pi \right)}^{n+2}}}{\displaystyle \int }\frac{{\text{e}}^{-iyz}}{iz}\left[{\text{e}}^{i\left(b^{\text{T}}x\right)z}\right]p\left(x\right)\text{d}x\text{d}z-{\displaystyle \int }\frac{{\text{e}}^{-iyz}}{iz}\text{d}z\\ =\frac{1}{\sqrt{{\left(2\pi \right)}^{n+2}}}\left[{\displaystyle \int }1+\frac{{\text{e}}^{i\left(b^{\text{T}}x-y\right)z}}{iz}p\left(x\right)\text{d}x\text{d}z\right]\end{array}$ (16)

since the last improper integral becomes:

$-{\displaystyle \int }\frac{{\text{e}}^{-isz}}{iz}\text{d}z=-\frac{2}{\pi }\left[{\displaystyle \underset{0}{\overset{\infty }{\int }}}{\text{e}}^{-isz}-{\text{e}}^{+isz}\right]\frac{1}{z}\text{d}z=\frac{2}{\pi }{\displaystyle \underset{0}{\overset{\infty }{\int }}}\frac{\sin sz}{z}\text{d}z=1$

so finally

$Q\left(y\right)=\frac{1}{\sqrt{{\left(2\pi \right)}^n}}\left[{\displaystyle \int }{\text{e}}^{i\left(b^{\text{T}}x\right)\frac{z}{y}}p\left(x\right)\text{d}x-\frac{1}{2\pi }\right]\mathrm{.}$

6. Asian Options

Asian options may be priced according to the distribution of the multivariate lognormal. The Asian call option is payable as the excess (over a given strike K) of the (arithmetic) average of prices at predetermined times $T=\left[t_1\mathrm{,}{t}_2,\cdots ,t_n\right]$. Assuming that prices $P_t$ follow a standard unit Wiener process $W_t$ (and assuming for simplicity that the stock has zero drift, no dividends and unit variance),

we have $\frac{\text{d}{P}_t}{P_t}=\text{d}{W}_t$, and thus

$P_t={\text{e}}^{W_t-\frac{1}{2}t}.$

which is lognormal. Thus the covariance of the exponent at two times, s and t, is $min\left(s\mathrm{,}t\right)$, and the covariance matrix is given by:

${\Sigma }_{ij}=\min \left(t_i,t_j\right).$

(This does not depend on any drift that might be introduced into the stock price process). The distribution of the average price at the relevant times may then be determined from the above procedure. This contrasts with the approximations based on two moments [15]. It also satisfies the need to derive the density of a sum of correlated lognormals as in [16], who derive the density of the continuous average price as a reciprocal gamma.

Example

We take $T=\left[\mathrm{1,2,3,4}\right]/4$ over a unit period, for which

$\Sigma =\left[\begin{array}{cccc}0.25& 0.25& 0.25& 0.25\\ 0.25& 0.50& 0.50& 0.50\\ 0.25& 0.50& 0.75& 0.75\\ 0.25& 0.50& 0.75& 1.00\end{array}\right]$

This contradicts the assertion in [17] that the Asian option has a lower value than the European. Since the Wiener process enjoys independent increments, the value of $W_t$ is affected by the value of $W_s$ for $s<t$, and so averaging prices across a period does not necessarily reduce volatility relative to the terminal price.

$A=\left[\begin{array}{cccc}+0.11401& -0.17467& +0.15360& -0.06066\\ -0.21426& +0.07441& +0.18842& -0.13985\\ +0.28868& +0.28868& -0.00000& -0.28868\\ +0.32827& +0.61694& +0.83121& +0.94521\end{array}\right]$

$A^{-1}=\left[\begin{array}{cccc}+1.61073& -2.01175& +1.15470& +0.15838\\ -2.46778& +0.69867& +1.15470& +0.29765\\ +2.17013& +1.76910& -0.00000& +0.40102\\ -0.85705& -1.31308& -1.15470& +0.45603\end{array}\right]$

Since the price $P_t={\text{e}}^{W_t-\frac{1}{2}t}$ has a weighting of ${\text{e}}^{-\frac{1}{2}t}$ when applied to the lognormal, we take $b={\text{e}}^{-\frac{1}{2}T}/n$. The distribution and cumulative distribution of the average stock price $Q={\displaystyle \underset{i=1}{\overset{n}{\sum }}}{P}_{t_i}/n$ are depicted in Figure 3. And the cumulative distributions are depicted in Figure 4.

With a strike of $K=0.4$, the Asian or European options can be priced as

$C={\displaystyle \underset{x\ge K}{\int }}\left(x-K\right)p\left(x\right)\text{d}x$

where $p\left(x\right)$ is the density for either the Asian price 14, or it can be taken as the European price at expiration:

Figure 3. Probability densities.

Figure 4. Cumulative densities.

$p\left(x\right)=\frac{{\text{e}}^{-\frac{{\left(x+1\right)}^2}{2}}}{\sqrt{2\pi }}$

In numerical terms these are easily computed:

$\text{Asian}=0.8912$

$\text{European}=0.4485$

Evidently the Asian is the more valuable as the average price over a period of times has a greater likelihood of exceeding the strike than at any particular given time.

7. Conclusion

Contrary to much of what is asserted in the literature, the lognormal distribution can be comfortably accommodated within complex function theory, as given in 3. No approximations are required, except to the extent that numerical integrations over infinite domains may not be exactly precise. This facilitates the application of the lognormal to all its areas of relevance, including the pricing of Asian options.

NOTES

¹However only the result for real z is need to apply the inverse FT.

²https://www.chebfun.org; developed at Oxford University with use of Chebyshev polynomials.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Thorin, O. and Wikstad, N. (1974) Calculation of Ruin Probabilities When the Claim Distribution Is Lognormal. ASTIN Bulletin, 9, 241-246.
[2]	Panjer, H.H. (1981) Recursive Evaluation of a Family of Compound Distributions. ASTM Bulletin, 12, 22-26. https://doi.org/10.1017/S0515036100006796
[3]	Holgate, P. (1989) The Lognormal Characteristic Function. Communications in Statistics: Theory and Methods, 18, 4539-4548. https://doi.org/10.1080/03610928908830173
[4]	Chakrabarti, K.B. (2013) The Characteristic Function of a Lognormal Process Derived Using the Saddle Point Method. International Journal of Electronics & Communication Technology, 4, 12.
[5]	Asmussen, S., Goffard, P.O. and Laub, P.J. (2015) Orthonormal Polynomial Expan sions and Lognormal Sum Densities. Norberg Festschrift, 9, 54.
[6]	Asmussen, S., Jensen, J. and Rojas-Nandayapa, L. (2014) On the Laplace Transform of the Lognormal Distribution. Methodology and Computing in Applied Probability, 18, 441-458. https://doi.org/10.1007/s11009-014-9430-7
[7]	Furman, E., Hackmann, D. and Kuznetsov, A. (2019) On Log-Normal Convolutions: An Analytical-Numerical Method with Applications to Economic Capital Determination. Insurance Mathematics and Economics, 90, 120-134. https://doi.org/10.1016/j.insmatheco.2019.10.003
[8]	Rossberg, A.G. (1998) Laplace Transforms of Probability Distributions and Their Inversions Are Easy on Logarithmic Scales. Journal of Applied Probability, 45, 531-541. https://doi.org/10.1239/jap/1214950365
[9]	Wiegand, M. and Nadarajah, S. (2018) Approximation Methods for Lognormal Characteristic Functions. Journal of Statistical Computation and Simulation, 88, 3650-3663. https://doi.org/10.1080/00949655.2018.1530775
[10]	Gubner, J.A. (2006) A New Formula for Lognormal Characteristic Functions. IEEE Transactions on Vehicular Technology, 55, 747-773. https://doi.org/10.1109/TVT.2006.878610
[11]	Hansen, A.T. and Jørgensen, P.L. (2000) Analytical Valuation of American-Style Asian Options. Management Science, 46, 1116-1136. https://doi.org/10.1287/mnsc.46.8.1116.12027
[12]	Schroder, M. (2008) On Constructive Complex Analysis in Finance: Explicit Formulas for Asian Options. Quarterly Journal of Applied Mathematics, 66, 633-658. https://doi.org/10.1090/S0033-569X-08-01107-X
[13]	Miles, J. (2018) The Laplace Transform of the Lognormal Distribution. https://arxiv.org/pdf/1803.05878.pdf
[14]	Lo, C.F. (2012) The Sum and Difference of Two Lognormal Random Variables. Journal of Applied Mathematics, 2012, Article ID: 838397. https://doi.org/10.1155/2012/838397
[15]	Milevsky, M.A. and Posner, S.E. (1998) Asian Options, the Sum of Lognormals, and the Reciprocal Gamma Distribution. The Journal of Financial and Quantitative Analysis, 33, 409-422. https://doi.org/10.2307/2331102
[16]	Halliwell, L.J. (2012) The Lognormal Random Multivariate. Casualty Actuarial Society, Arlington.
[17]	Li, W.P. and Chen, S. (2016) Pricing and Hedging of Arithmetic Asian Options via the Edgeworth Series Expansion Approach. The Journal of Finance and Data Science, 2, 1-25. https://doi.org/10.1016/j.jfds.2016.01.001