Correlated Log-Normal Random Variables under a Multiscale Volatility Model

The study focuses on extending the fast mean-reversion volatility, which was developed by the author in a previous work, to the multiscale volatility model so that it can express a well-separated time scale. The leading-order term and ﬁ rst-order correction terms are analytically computed using the perturbation theory based on the Lie – Trotter operator splitting method. Finally, the study is concluded by deriving the numerical results that further validate the e ﬀ ectiveness of the model.


Introduction
A stochastic process provides a useful tool to analyze time series data and wide applications in many fields such as physics, finance, biotechnology, and telecommunication studies. The transition density function of a continuous-time process plays an important role in understanding and explaining the dynamics of the process. However, the transition density functions are unknown for general diffusion processes except for a few special cases (refer to Aït-Sahalia [1], Black and Scholes [2], Cox et al. [3], and Vasicek [4]). So finding analytical approximations to them is important as numerical methods, such as finite-difference method, Monte Carlo simulation, and Fourier inversion, because of much faster and precise at least under certain model parameter regime.
In particular, the multivariate log-normal distribution is a widely used stochastic model in social sciences. What is the probability of the sum or difference of log-normal random variables? The solution to this question has wide applications in many fields such as finance [5], physics [6], and actuarial science [7]. However, to the best of the author of this paper's knowledge, almost nothing is known about the question yet (see [8]).
Fat-tailed distribution and volatility clustering are widely utilized to yield stylized facts in the field of financial modeling. It is known that high volatility leads to high volatility and low volatility leads to low volatility. Generally, such phenomena are observed on different time scales and are characterized by the tendency of short-run volatility clustering (or fast mean-reversion) and long-run volatility clustering. In addition, these phenomena can only be delineated by multiscale models and not by single-scale models. Therefore, this study focuses on extending the fast meanreversion volatility (FMRV) model, which is a previously proposed model by the author, to multiscale volatility (MSV) model. This extension would aid the researchers and practitioners, such as Adrian and Rosenberg [9], Chernov et al. [10], Fouque et al. [11], and Gallant et al. [12], in obtaining financial derivatives, especially those with derived empirical evidence.
Based on the observation that the well-separated time scale exists in financial time series data, in this work, a slow varying process is incorporated into the previously derived results. Moreover, approximate transition density under the MSV model is derived by the perturbation theory based on the Lie-Trotter operator splitting method.

Problem Formulation.
We denote x 1 ± x 2 and x 10 ± x 20 as x ± and x ± 0 , respectively. The probability distribution of the sum or difference of the two correlated log-normal distributions can be obtained by calculating the integral where u is the joint probability distribution of the two lognormal random variables and δðx ± − x ± 0 Þ is the Dirac delta function. However, a closed-form representation for this probability distribution does not exist. Thus, we derive an analytic approximation of the transition density function of the sum or difference of the two correlated log-normal distributions (refer to [13]).
We start with a process ðX 1 , X 2 , Y, ZÞ describing the conditionally log-normal random variables X 1 and X 2 and two processes Y and Z to express a well-separated time scale. These two processes represent a fast scale factor and a slow scale factor of the volatility, respectively. The dynamics of the joint process ðX 1 , X 2 , Y, ZÞ are given by the following stochastic differential equations (SDEs) under a risk-neutral measure: where W 1 , W 2 , W Y , and W Z are standard Brownian motions correlated as follows: with Here, the functions f 1 : ℝ 2 ⟶ ℝ + and f 2 : ℝ 2 ⟶ ℝ + are assumed to be bounded, smooth, and strictly positive. The functions g and h satisfy the Lipschitz and growth conditions such that the corresponding SDE yields a unique strong solution. The parameters ε, δ, m, and ν are positive constants with the same order of ε ≈ δ ≪ 1. Notably, the process Y is an ergodic process such that its invariant distribution is given by the Gaussian probability distribution function as This provides an important averaging tool for the unobserved process Y that is well documented in [11]. Notation h·i is adopted for the expectation with respect to this invariant distribution.
By using the four-dimensional Feynman-Kac formula (cf. [14]), we obtain a singularly and regularly perturbed partial differential equation (PDE) problem given by where Here, L 0 is the infinitesimal generator of the Ornstein-Uhlenbeck (OU) process Y. L 1 consists of the mixed partial derivatives due to the correlations of the two Brownian motions W 1 and W Y and W 2 and W Y , respectively. L 2 corresponds to the operator of a generalized version of the twodimensional standard Brownian motion at the volatility levels f 1 ðy, zÞ and f 2 ðy, zÞ instead of constant volatilities σ 1 and σ 2 , respectively. K 1 includes the mixed partial derivatives due to the correlations of the two standard Brownian motions W 1 and W Z and W 2 and W Z , respectively. K 2 is the infinitesimal generator of the process Z. Finally, K 3 denotes the mixed partial derivative due to the correlation of the two standard Brownian motions W Y and W Z .

2
Advances in Mathematical Physics

Multiscale Analysis.
Since it is difficult to solve the PDE problem (5), we are interested in the following asymptotic expansions: such that u ε,δ ± is a series of the general term ε j/2 δ j/2 u i,j ± . Substituting the expansion (9) into (5) gives u ε,0 ± and u ε,1 ± , given by the solutions of the PDEs

Perturbation Theory Based on Lie-Trotter Operator Splitting Method
Theorem 1. Assume that the partial derivative of u i,j ± with respect to y does not grow as much as ð∂u i,j ± /∂yÞ~e y 2 /2 as y tends to infinity. Then, the leading-order term of the expansion (9) with j = 0 is independent of the fast-scale variable y, and it takes the following form: whereσ Here, σ 1 ðzÞ and σ 1 ðzÞ are defined by respectively, in terms of ϕ (the invariant distribution of Y).
Proof of Theorem 1. Applying the expansion with j = 0 to u 0,0 ± leads to By multiplying (15) with ε and then letting ε tend to zero, we obtain the ordinary differential equation (ODE) Reiterating the fact that the operator L 0 is the generator of the OU process Y, the solution u 0,0 ± of this ODE must be a constant with respect to the y variable because of the assumed growth condition: In addition, we have the PDE L 0 u 1,0 ± + L 1 u 0,0 ± = 0. Applying the y-independence of u 0,0 ± to this PDE reduces it to the ODE L 0 u 1,0 ± =0. Therefore, u 1,0 ± is independent of y because of the same reason as u 0,0 ± : Hitherto, the first two terms u 0,1 ± and u 1,0 ± did not depend on the current level y of the fast-scale volatility, which is driving the process Y.
To simplify the equation, one can continue to eliminate the terms of order 1, ffiffi ε p , ε, and so forth. From the order 1 terms of (15), we get L 0 u 2,0 ± + L 1 u 1,0 ± + L 2 u 0,0 ± = 0. Then, this PDE becomes due to the y-independence of u 1,0 ± . Equation (19) is a Poisson equation for u 2,0 ± with respect to the operator L 0 in the y variable. It is well known from the Fredholm alternative (solvability condition) that it has a solution only if the source term L 2 u 0,0 ± is centered with respect to the invariant distribution of Y; namely, 3 Advances in Mathematical Physics Then, u 0,0 ± solves the PDE (20).
Note that u 0,0 ± is identical to the transition density function of Lo's result such that the coefficients σ 1 and σ 2 are replaced by σ 1 ðzÞ and σ 2 ðzÞ, respectively. Furthermore, we employ the leading-order solution u 0,0 ± obtained above to derive the first-order correction terms u 1,0 ± and u 0,1 ± , respectively.

Theorem 2.
Assume that the partial derivative of u i,j ± with respect to y does not grow as much as ðu i,j ± /∂yÞ~e y 2 /2 as y tends to infinity. The correction term u 1,0 ± is independent of the y variable, and the first-order correction term satisfies the PDE where A = hL 1 L −1 0 ðL 2 − hL 2 iÞi. Then, one obtains the solution of the PDE where the constant parameters Here, the functions ϕ 1 , ϕ 2 , and ϕ 12 are defined by the solution of Proof of Theorem 2. The order ffiffi ε p terms in (15) lead to L 0 u 3,0 ± + L 1 u 2,0 ± + L 2 u 1,0 ± = 0 which is a Poisson equation for u 3,0 ± whose centering condition is given by Meanwhile, from (19) and the operator hL 2 i, we get for some function cðt, x 1 , x 2 , z 1 Þ that does not depend on the y variable. Substituting (27) into (26), a PDE for u 1,0 ± is obtained as This implies that u 1,0 ± is y-independent. Since we focus on the first-order correction terms of u 0,0 ± , we reconstruct (28) with respect toũ 1,0 ± as follows: To calculate the operatorÃ, we aim to derive Thus, the operatorÃ can be expressed as where V 1 ðzÞ, V 2 ðzÞ, V 12 ðzÞ, and V 21 ðzÞ are those given in Theorem 2. Subsequently, we obtain the result of Theorem 2 by direct computation.
Similarly, the correction term u 0,1 ± can be obtained.

Advances in Mathematical Physics
Theorem 3. Assume that the partial derivative of u i,j ± with respect to y does not grow as much as ðu i,j ± /∂yÞ~e y 2 /2 as y tends to infinity. The first correction term u 0,1 ± is independent of the y variable, and the first-order correction term satisfies the PDE problem where which makes obtaining the solution to the PDE feasible where Proof of Theorem 3. To obtain another first-order correction term, it is necessary to consider another singular perturbation problem. Applying expansion (9) with j = 0 and j = 1 leads to By multiplying (36) with ε and then letting ε tend to zero, we find the first two leading-order terms as follows: u 0,1 ± (the solution of L 0 u 0,1 ± = 0) must be a constant with respect to the y variable because the operator L 0 is the generator of the OU process Y. Since K 3 has a derivative with respect to the y variable and u 0,0 ± is independent of y, we obtain Moreover, because each term of L 1 has a derivative with respect to y, L 1 u 0,1 ± = 0 holds. In light of the above observations, the equation reduces to which implies u 1,1 ± does not depend on the y variable. Hence, the two terms u 0,1 ± and u 1,1 ± do not depend on the current level y of the fast-scale volatility driving process Y: This makes the continuous elimination of order 1 terms, e.g., ffiffi ε p and ε, possible. For the order 1 term, we have As u 1,0 ± and u 1,1 ± are y-independent, the preceding PDE (42) can be expressed as This represents a Poisson equation for u 2,1 ± with respect to the operator L 0 in the y variable, which has a solution only if is centered with respect to the invariant distribution of Y.
Since u 0,0 ± and u 0,1 ± do not depend on the variable y, we have Then, we obtain the result of Theorem 3 by direct computation.
As a result, one can approximate u ε,δ to the summation of the leading-order term u 0,0 ± and the first correction terms u 1,0 ± and u 0,1 ± as follows: Note that the distributions of the sum and difference of correlated log-normal random variables under the MSV model are observed to follow a combination of shifted lognormal distribution and their mixed partial derivatives. In addition, all the original parameters are absorbed in the group parameters V 1 ðzÞ, V 2 ðzÞ, V 12 ðzÞ, V 21 ðzÞ, SðzÞ, and 5

Advances in Mathematical Physics
RðzÞ. The present level y of the hidden process Y that drives the fast time scale volatility need not be specified in the present approximation.

Numerical Experiment
In this section, we strengthen the effectiveness of the derived result (46) by showing numerical results.
Several parameter sets are obtained from the study of Ma except for new group parameters. Three main results can be observed from the left of Figure 1. First, the constant volatility is the leading-order term (shifted log-normal distribution). Second, the short-scale volatility is a combination of the leading-order term and the first correction term driven by the fast moving fluctuation (shifted log-normal distributions under the FMRV model). Finally, two different time scale volatility is a combination of the leading-order term and the first correction terms (shifted log-normal distributions under the MSV model). The right of Figure 1 depicts the errors calculated by subtracting the log-normal distributions under the FMRV model and the MSV model from the log-normal distribution, respectively. Additionally, we applied the same result to the difference in the probability density function in Figure 2. Numerical results delineate that the short-scale volatility and the multiscale volatility shift to the left (or to the right) from the constant volatility. The left and right skew scenarios are presented through the firstorder correction terms, and the major discrepancies appear around the peak of the probability density function. There is a small gap between the FMRV model and the MSV model; this implies that the prices of financial derivatives with shortterm maturity levels are ruled by a fast-scale volatility process. On the contrary, those of financial derivatives with long-term maturity levels are dominated by a slow-scale volatility process in the area of financial modeling. Finally, a simulation result shows that our result could explain a financial point of view. Note that the numerical result clearly outlines that this domain is sensitive to the choice of the involved parameters and provides considerable flexibility to the shape of the transition density.

Final Remarks
We incorporate a slow varying process into the previously published result of the author. A two-factor stochastic volatility model is formulated that comprises a slow varying  Advances in Mathematical Physics process representing one persistent factor for volatility and one ergodic process for displaying rapidly moving fluctuations. Using perturbation theory based on the Lie-Trotter operator splitting method, we compute the leading-order term, the first-order correction terms, and the effectiveness of the slow-scale volatility factor by outlining numerical results. This paper offers various possible directions for further development. This result can provide a very useful guide for credit risk management (see [13,[15][16][17]). We leave these issues as future research topics.

Data Availability
Our paper contains numerical experimental results, and values for these experiments are included in the paper. The data is freely available.

Conflicts of Interest
The author declares that there is no conflict of interests regarding the publication of this paper.