The Sum and Difference of Two Lognormal Random Variables

Abstract

We have presented a new unified approach to model the dynamics of both the sum and difference of two correlated lognormal stochastic variables. By the Lie-Trotter operator splitting method, both the sum and difference are shown to follow a shifted lognormal stochastic process, and approximate probability distributions are determined in closed form. Illustrative numerical examples are presented to demonstrate the validity and accuracy of these approximate distributions. In terms of the approximate probability distributions, we have also obtained an analytical series expansion of the exact solutions, which can allow us to improve the approximation in a systematic manner. Moreover, we believe that this new approach can be extended to study both (1) the algebraic sum of N lognormals, and (2) the sum and difference of other correlated stochastic processes, for example, two correlated CEV processes, two correlated CIR processes, and two correlated lognormal processes with mean-reversion.

1. Introduction

“Given two correlated lognormal stochastic variables, what is the stochastic dynamics of the sum or difference of the two variables?”; or equivalently “What is the probability distribution of the sum or difference of two correlated lognormal stochastic variables?” The solution to this long-standing problem has wide applications in many fields such as telecommunication studies [1–6], financial modelling [7–9], actuarial science [10–12], biosciences [13], physics [14], and so forth. Although the lognormal distribution is well known in the literature [15, 16], yet almost nothing is known of the probability distribution of the sum or difference of two correlated lognormal variables. However, it is commonly agreed that the distribution of either the sum or difference is neither normal nor lognormal.

The aforesaid problem can be formulated as follows. Given two lognormal stochastic variables 𝑆₁ and 𝑆₂ obeying the following stochastic differential equations: 𝑑𝑆_𝑖𝑆_𝑖=𝜎_𝑖𝑑𝑍_𝑖,𝑖=1,2,(1.1) where 𝜎^2𝑖=Var(ln𝑆_𝑖), 𝑑𝑍_𝑖 denotes a standard Weiner process associated with 𝑆_𝑖, and the two Weiner processes are correlated as 𝑑𝑍₁𝑑𝑍₂=𝜌𝑑𝑡, the time evolution of the joint probability distribution function 𝑃(𝑆₁,𝑆₂,𝑡;𝑆₁₀,𝑆₂₀,𝑡₀) of the two correlated lognormal variables is governed by the backward Kolmogorov equation 𝜕𝜕𝑡₀+𝐿𝑃𝑆₁,𝑆₂,𝑡;𝑆₁₀,𝑆₂₀,𝑡₀=0for𝑡>𝑡₀,(1.2) where 𝐿=12𝜎²¹𝑆²¹⁰𝜕²𝜕𝑆²¹⁰+𝜌𝜎₁𝜎₂𝑆₁₀𝑆₂₀𝜕²𝜕𝑆₁₀𝜕𝑆₂₀+12𝜎²²𝑆²²⁰𝜕²𝜕𝑆²²⁰(1.3) subject to the boundary condition 𝑃𝑆₁,𝑆₂,𝑡;𝑆₁₀,𝑆₂₀,𝑡₀⟶𝑡=𝛿𝑆₁−𝑆₁₀𝛿𝑆₂−𝑆₂₀.(1.4) This joint probability distribution function tells us how probable the two lognormal variables assume the values 𝑆₁ and 𝑆₂ at time 𝑡>𝑡₀, provided that their values at 𝑡₀ are given by 𝑆₁₀ and 𝑆₂₀. Since 𝑃(𝑆₁,𝑆₂,𝑡;𝑆₁₀,𝑆₂₀,𝑡₀) is known in closed form as follows: 𝑃𝑆=₁,𝑆₂,𝑡;𝑆₁₀,𝑆1₂₀,𝑡₀2𝜋𝜎₁𝜎₂√1−𝜌²𝑆₁𝑆₂exp−ln𝑆¹−ln𝑆₁₀+(1/2)𝜎²¹𝑡−𝑡₀²×exp𝜌ln𝑆₁2𝜎²¹1−𝜌²𝑡−𝑡₀−ln𝑆₁₀+(1/2)𝜎²¹𝑡−𝑡₀ln𝑆₂−ln𝑆₂₀+(1/2)𝜎²²𝑡−𝑡₀×exp−ln𝑆₂𝜎₁𝜎₂1−𝜌²𝑡−𝑡₀−ln𝑆₂₀+(1/2)𝜎²²𝑡−𝑡₀²2𝜎²²1−𝜌²𝑡−𝑡₀,(1.5) the probability distribution of the sum or difference, namely 𝑆^±≡𝑆₁±𝑆₂, of the two correlated lognormal variables can be obtained by evaluating the integral 𝑃_±𝑆^±,𝑡;𝑆₁₀,𝑆₂₀,𝑡₀=^∞0𝑑𝑆₁𝑑𝑆₂𝑃𝑆₁,𝑆₂,𝑡;𝑆₁₀,𝑆₂₀,𝑡₀𝛿𝑆₁±𝑆₂−𝑆^±.(1.6) Despite that many methods have been developed to address the problem, a closed-form representation for the probability distribution of the sum or difference is still missing. Hence, we must resort to numerical methods to perform the integration. Nevertheless, the numerically exact solution does not provide any information about the stochastic dynamics of the sum or difference explicitly.

In the lack of knowledge about the probability distribution of the sum or difference of two correlated lognormal variables, several analytical approximation methods which focus on finding a good approximation for the desired probability distribution have been proposed in the literature [1–6, 8, 17–27]. Essentially, these analytical approximations assume a specific distribution that the sum or difference of the two correlated lognormal variables follow, and then use a variety of methods to identify the parameters for that specific distribution. However, no mathematical justification for the specific distribution was apparently given. In spite of this shortcoming, these approximations attract considerable attention and have been extended to tackle the algebraic sums of 𝑁 correlated lognormal variables, too.

In this communication, we apply the Lie-Trotter operator splitting method [28] to derive an approximation for the dynamics of the sum or difference of two correlated lognormal variables. It is shown that both the sum and difference can be described by a shifted lognormal stochastic process. Approximate probability distributions of both the sum and difference of the lognormal variables are determined in closed form, and illustrative numerical examples are presented to demonstrate the accuracy of these approximate distributions. Unlike previous studies which treat the sum and difference in a separate manner, our proposed method thus provides a new unified approach to model the dynamics of both the sum and difference of two correlated lognormal stochastic variables. In addition, in terms of the approximate solutions, we are able to obtain an analytical series expansion of the exact solutions, which can allow us to improve the approximation systematically. Moreover, we believe that this new approach can be extended to study both (1) the algebraic sum of 𝑁 lognormals, and (2) the sum and difference of other correlated stochastic processes, for example, two correlated CEV processes, two correlated CIR processes, and two correlated lognormal processes with mean-reversion.

2. Lie-Trotter Operator Splitting Method

It is observed that the probability distribution of the sum or difference of the two correlated lognormal variables, that is, 𝑃_±(𝑆^±,𝑡;𝑆₁₀,𝑆₂₀,𝑡₀), also satisfies the same backward Kolmogorov equation given in (1.2), but with a different boundary condition 𝑃_±𝑆^±,𝑡;𝑆₁₀,𝑆₂₀,𝑡₀⟶𝑡=𝛿𝑆₁₀±𝑆₂₀−𝑆^±.(2.1) To solve for 𝑃_±(𝑆^±,𝑡;𝑆₁₀,𝑆₂₀,𝑡₀), we first rewrite the backward Kolmogorov equation in terms of the new variables 𝑆^±0≡𝑆₁₀±𝑆₂₀ as 𝜕𝜕𝑡₀+𝐿₊+𝐿₀+𝐿₋𝑃_±𝑆^±,𝑡;𝑆⁺⁰,𝑆⁻⁰,𝑡₀=0,(2.2) where 𝐿₊=18𝜎²⁺𝑆⁺⁰²+2𝜎²¹−𝜎²²𝑆⁺⁰𝑆⁻⁰+𝜎²⁻𝑆⁻⁰²𝜕²𝜕𝑆⁰⁺²𝐿₀=14𝜎²¹−𝜎²²𝑆⁺⁰²+𝑆⁻⁰²+𝜎²¹+𝜎²²𝑆⁺⁰𝑆⁻⁰𝜕²𝜕𝑆⁺⁰𝜕𝑆⁻⁰𝐿₋=18𝜎²⁺𝑆⁻⁰²+2𝜎²¹−𝜎²²𝑆⁺⁰𝑆⁻⁰+𝜎²⁻𝑆⁺⁰²𝜕²𝜕𝑆⁰⁻²𝜎_±=𝜎²¹+𝜎²²±2𝜌𝜎₁𝜎₂.(2.3) The corresponding boundary condition now becomes 𝑃_±𝑆^±,𝑡;𝑆⁺⁰,𝑆⁻⁰,𝑡₀⟶𝑡=𝛿𝑆^±0−𝑆^±.(2.4) Accordingly, the formal solution of (2.2) is given by 𝑃_±𝑆^±,𝑡;𝑆⁺⁰,𝑆⁻⁰,𝑡₀=exp𝑡−𝑡₀𝐿₊+𝐿₀+𝐿₋𝛿𝑆^±0−𝑆^±.(2.5)

Since the exponential operator exp{(𝑡−𝑡₀)(𝐿₊+𝐿₀+𝐿₋)} is difficult to evaluate, we apply the Lie-Trotter operator splitting method [28] to approximate the operator by (see the appendix) 𝑂^LT_±=exp𝑡−𝑡₀𝐿_±exp𝑡−𝑡⁰𝐿₀+𝐿_∓,(2.6) and obtain an approximation to the formal solution 𝑃_±(𝑆^±,𝑡;𝑆⁺⁰,𝑆⁻⁰,𝑡₀), namely 𝑃^LT_±𝑆^±,𝑡;𝑆⁺⁰,𝑆⁻⁰,𝑡₀=𝑂^LT_±𝛿𝑆^±0−𝑆^±=exp𝑡−𝑡₀𝐿_±𝛿𝑆^±0−𝑆^±,(2.7) where the relation exp{(𝑡−𝑡₀)(𝐿₀+𝐿_∓)}𝛿(𝑆^±0−𝑆^±)=𝛿(𝑆^±0−𝑆^±) is utilized. For (𝑆⁻⁰/𝑆⁺⁰)²≪1, which is normally valid unless 𝑆₁₀ and 𝑆₂₀ are both close to zero, the operators 𝐿₊ and 𝐿₋ can be approximately expressed as 𝐿_±≅12𝜎^2±𝑆^0±2𝜕²𝜕𝑆^0±2(2.8) in terms of the two new variables: 𝑆⁺⁰=𝑆⁺⁰+𝜎²¹−𝜎²²𝜎²⁺𝑆⁻⁰𝑆₋₀,=𝑆⁻⁰+𝜎²⁻𝜎²¹−𝜎²²𝑆⁺⁰,(2.9) where 𝜎₊=𝜎₊/2and𝜎₋=(𝜎²¹−𝜎²²)/(2𝜎₋). Without loss of generality, we assume that 𝜎₁>𝜎₂. Obviously, both 𝑆⁺ and 𝑆⁻ are lognormal (LN) random variables defined by the stochastic differential equations 𝑑𝑆^±=𝜎_±𝑆^±𝑑𝑍_±,(2.10) and their closed-form probability distribution functions are given by 𝑓^LN𝑆^±,𝑡;𝑆^±0,𝑡₀=1𝑆^±2𝜋𝜎^2±𝑡−𝑡₀exp⎧⎪⎨⎪⎩−ln𝑆^±/𝑆^±0+(1/2)𝜎^2±𝑡−𝑡₀²2𝜎^2±𝑡−𝑡₀⎫⎪⎬⎪⎭(2.11) for 𝑡>𝑡₀. As a result, it can be inferred that within the Lie-Trotter splitting approximation both 𝑆⁺ and 𝑆⁻ are governed by a shifted lognormal process. It should be noted that for the Lie-Trotter splitting approximation to be valid, 𝜎^2±(𝑡−𝑡₀) needs to be small.

Alternatively, we can also approximate the operator 𝐿₋ by 𝐿₋≅12𝜎²⁻𝑅⁻⁰𝜕²𝜕𝑅⁰⁻²,(2.12)where𝜎₋=(𝜎²¹−𝜎²²)𝑆⁺⁰/2 and 𝑅⁻⁰=𝑆⁻⁰+12𝜎²¹𝜎−²⁻𝜎²²𝑆⁺⁰.(2.13) It is not difficult to recognize that 𝑅⁻ follows the square-root (SR) stochastic process defined by the stochastic differential equation 𝑑𝑅⁻=𝜎₋√𝑅⁻𝑑𝑍₋,(2.14) and has the closed-form probability distribution function 𝑓^SR𝑅⁻,𝑡;𝑅⁻⁰,𝑡₀=2𝜎²⁻𝑡−𝑡₀𝑅⁻⁰𝑅⁻exp−2𝑅⁻+𝑅⁻⁰𝜎²⁻𝑡−𝑡₀𝐼¹4√𝑅⁻𝑅⁻⁰𝜎²⁻𝑡−𝑡₀(2.15) for 𝑡>𝑡₀, where 𝐼₁(⋅) is the modified Bessel function of the first kind of order one. Accordingly, we have shown that within the Lie-Trotter splitting approximation, which requires 𝜎²⁻(𝑡−𝑡₀) to be small, 𝑆⁻ can be described by a shifted square-root process, too.

Moreover, in terms of the approximate solutions 𝑃^LT_±(𝑆^±,𝑡;𝑆⁺⁰,𝑆⁻⁰,𝑡₀), we can express the exact solutions 𝑃_±(𝑆^±,𝑡;𝑆⁺⁰,𝑆⁻⁰,𝑡₀) in the following form: 𝑃_±=𝑆^±,𝑡;𝑆⁺⁰,𝑆⁻⁰,𝑡₀+𝑃^LT_±𝑡𝑡0𝑆^±,𝑡;𝑆⁺⁰,𝑆⁻⁰,𝑡₀=1+𝑑𝜉e_𝑡𝑡0xp𝜉−𝑡₀𝐿_±𝐿⁰+𝐿_∓𝑃_±𝑆^±,𝑡;𝑆⁺⁰,𝑆⁻⁰,𝜉𝑑𝜉₁ℒ_±𝑡₀−𝜉₁+^𝑡𝑡₀𝑑𝜉₁^𝑡𝜉₁+_𝑡𝑡0𝑑𝜉₂ℒ_±𝑡₀−𝜉₁ℒ_±𝑡₀−𝜉₂𝑑𝜉₁^𝑡𝜉₁𝑑𝜉₂^𝑡𝜉₂+_𝑡𝑡0𝑑𝜉₃ℒ_±𝑡₀−𝜉₁ℒ_±𝑡₀−𝜉₂ℒ_±𝑡₀−𝜉₃𝑑𝜉₁^𝑡𝜉₁𝑑𝜉₂^𝑡𝜉₂𝑑𝜉₃^𝑡𝜉₃𝑑𝜉₄⋯𝑃^LT_±𝑆^±,𝑡;𝑆⁺⁰,𝑆⁻⁰,𝑡₀,(2.16) where ℒ_±(𝜏)==ex𝐿p₀−𝜏𝐿_±𝐿⁰+𝐿_∓exp𝜏𝐿^±+𝐿_∓+𝜏1!𝐿⁰+𝐿_∓,𝐿^±+𝜏²2!𝐿⁰+𝜏₃+𝐿_∓,𝐿^±,𝐿^±3!𝐿⁰+𝐿_∓,𝐿^±,𝐿^±,𝐿^±+⋯.(2.17) The integrals over the temporal variables {𝜉_𝑖;𝑖=1,2,3,…} can be evaluated analytically. If we keep terms up to the order of (𝑡−𝑡₀)², then 𝑃_±(𝑆^±,𝑡;𝑆⁺⁰,𝑆⁻⁰,𝑡₀) can be approximated by 𝑃_±𝑆^±,𝑡;𝑆⁺⁰,𝑆⁻⁰,𝑡₀=𝑃+^LT_±𝑡−𝑆^±𝑡₀,𝑡;𝑆⁺⁰,𝑆⁻⁰,𝑡₀+1𝐿₀+𝐿_∓𝑃^LT_±𝑆^±,𝑡;𝑆⁺⁰,𝑆⁻⁰,𝑡₀2𝑡−𝑡⁰²𝐿⁰+𝐿_∓²−𝐿₀+𝐿_∓,𝐿^±𝑃^LT_±𝑆^±,𝑡;𝑆⁺⁰,𝑆⁻⁰,𝑡₀.(2.18) This analytical series expansion can allow us to improve the approximate solutions systematically.

3. Illustrative Numerical Examples

In Figure 1 we plot the approximate closed-form probability distribution function of the sum 𝑆⁺ given in (2.11) for different values of the input parameters. We start with 𝑆₁₀=110, 𝑆₂₀=100, 𝜎₁=0.25, and 𝜎₂=0.15 in Figure 1(a). Then, in order to examine the effect of 𝑆₂₀, we decrease its value to 70 in Figure 1(b) and to 40 in Figure 1(c). In Figures 1(d), 1(e), and 1(f) we repeat the same investigation with a new set of values for 𝜎₁ and 𝜎₂, namely 𝜎₁=0.3 and 𝜎₂=0.2. Without loss of generality, we set 𝑡−𝑡₀=1 for simplicity. The (numerically) exact results which are obtained by numerical integrations are also included for comparison. It is clear that the proposed approximation can provide accurate estimates for the exact values. Moreover, to have a clearer picture of the accuracy, we plot the corresponding errors of the estimation in Figure 2. We can easily see that major discrepancies appear around the peak of the probability distribution function, and that the estimation deteriorates as the correlation parameter 𝜌 decreases from 0.5 to −0.5. It is also observed that the errors increase with the ratio 𝑆⁻⁰/𝑆⁺⁰ as expected but they seem to be not very sensitive to the changes in 𝜎₁ and 𝜎₂.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1 

Probability density versus 𝑆1+𝑆2: The solid lines denote the distributions of the approximate shifted lognormal process, and the dash lines show the exact results. (a) 𝑆10=110, 𝑆20=100, 𝜎1=0.25, and 𝜎2=0.15; (b) 𝑆10=110, 𝑆20=70, 𝜎1=0.25, and 𝜎2=0.15; (c) 𝑆10=110, 𝑆20=40, 𝜎1=0.25, and 𝜎2=0.15; (d) 𝑆10=110, 𝑆20=100, 𝜎1=0.3, and 𝜎2=0.2; (e) 𝑆10=110, 𝑆20=70, 𝜎1=0.3, and 𝜎2=0.2; (f) 𝑆10=110, 𝑆20=40, 𝜎1=0.3, and 𝜎2=0.2.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2 

Error versus 𝑆1+𝑆2: The error is calculated by subtracting the approximate estimate from the exact result. (a) 𝑆10=110, 𝑆20=100, 𝜎1=0.25, and 𝜎2=0.15; (b) 𝑆10=110, 𝑆20=70, 𝜎1=0.25, and 𝜎2=0.15; (c) 𝑆10=110, 𝑆20=40, 𝜎1=0.25, and 𝜎2=0.15; (d) 𝑆10=110, 𝑆20=100, 𝜎1=0.3, and 𝜎2=0.2; (e) 𝑆10=110, 𝑆20=70, 𝜎1=0.3, and 𝜎2=0.2; (f) 𝑆10=110, 𝑆20=40, 𝜎1=0.3, and 𝜎2=0.2.

Next, we apply the same sequence of analysis to the two approximate closed-form probability distribution functions of the difference 𝑆⁻ given in (2.11) and (2.15). Similar observations about the accuracy of the proposed approximation can be made for the difference 𝑆⁻, too (see Figures 3 and 4). However, contrary to the case of 𝑆⁺, the estimation performs better for positive correlation. Of the two different approximation schemes for the 𝑆⁻, the shifted LN process seems to have a comparatively better performance than the shifted SR process, as evidenced by the numerical results.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3 

Probability density versus 𝑆1−𝑆2: the dash lines denote the distributions of the approximate shifted lognormal process, the dotted lines indicate the distributions of the approximate shifted square-root process, and the solid lines show the exact results. (a) 𝑆10=110, 𝑆20=100, 𝜎1=0.25, and 𝜎2=0.15; (b) 𝑆10=110, 𝑆20=70, 𝜎1=0.25, and 𝜎2=0.15; (c) 𝑆10=110, 𝑆20=40, 𝜎1=0.25, and 𝜎2=0.15; (d) 𝑆10=110, 𝑆20=100, 𝜎1=0.3, and 𝜎2=0.2; (e) 𝑆10=110, 𝑆20=70, 𝜎1=0.3, and 𝜎2=0.2; (f) 𝑆10=110, 𝑆20=40, 𝜎1=0.3, and 𝜎2=0.2.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4 

Error versus 𝑆1−𝑆2: the error is calculated by subtracting the approximate estimate from the exact result. The dash lines denote the errors of the approximate shifted square-root process, and the solid lines show the errors of the approximate shifted lognormal process. (a) 𝑆10=110, 𝑆20=100, 𝜎1=0.25, and 𝜎2=0.15; (b) 𝑆10=110, 𝑆20=70, 𝜎1=0.25, and 𝜎2=0.15; (c) 𝑆10=110, 𝑆20=40, 𝜎1=0.25, and 𝜎2=0.15; (d) 𝑆10=110, 𝑆20=100, 𝜎1=0.3, and 𝜎2=0.2; (e) 𝑆10=110, 𝑆20=70, 𝜎1=0.3, and 𝜎2=0.2; (f) 𝑆10=110, 𝑆20=40, 𝜎1=0.3, and 𝜎2=0.2.

4. Conclusion

In this paper we have presented a new unified approach to model the dynamics of both the sum and difference of two correlated lognormal stochastic variables. By the Lie-Trotter operator splitting method, both the sum and difference are shown to follow a shifted lognormal stochastic process, and approximate probability distributions are determined in closed form. Illustrative numerical examples are presented to demonstrate the validity and accuracy of these approximate distributions. In terms of the approximate probability distributions, we have also obtained an analytical series expansion of the exact solutions, which can allow us to improve the approximation in a systematic manner. Moreover, we believe that this new approach can be extended to study both (1) the algebraic sum of 𝑁 lognormals, and (2) the sum and difference of other correlated stochastic processes, for example, two correlated CEV processes, two correlated CIR processes, and two correlated lognormal processes with mean-reversion.

Appendix

Lie-Trotter Splitting Approximation

Suppose that one needs to exponentiate an operator 𝐶 which can be split into two different parts, namely 𝐴 and 𝐵. For simplicity, let us assume that 𝐶=𝐴+𝐵, where the exponential operator exp(𝐶) is difficult to evaluate but exp(𝐴) and exp(𝐵) are either solvable or easy to deal with. Under such circumstances, the exponential operator exp(𝜀𝐶), with 𝜀 being a small parameter, can be approximated by the Lie-Trotter splitting formula [28]: exp𝜀𝐶=exp𝜀𝐴exp𝜀𝐵+𝑂𝜀².(A.1) This can be seen as the approximation to the solution at 𝑡=𝜀 of the equation 𝑑𝑌/𝑑𝑡=(𝐴+𝐵)𝑌 by a composition of the exact solutions of the equations 𝑑𝑌/𝑑𝑡=𝐴𝑌 and 𝑑𝑌/𝑑𝑡=𝐵𝑌 at time 𝑡=𝜀.