Model of continuous random cascade processes in financial markets

This article present a continuous cascade model of volatility formulated as a stochastic differential equation. Two independent Brownian motions are introduced as random sources triggering the volatility cascade. One multiplicatively combines with volatility; the other does so additively. Assuming that the latter acts perturbatively on the system, then the model parameters are estimated by application to an actual stock price time series. Numerical calculation of the Fokker--Planck equation derived from the stochastic differential equation is conducted using the estimated values of parameters. The results reproduce the pdf of the empirical volatility, the multifractality of the time series, and other empirical facts.


Introduction
In financial time series, past coarse-grained measures of volatility correlate better to future fine-scale volatility than the reverse process does. Such a causal structure of financial time series was first reported by Müller et al. [1]. Since then, the causal structure between time scales, the flow of information from a long-term to a short-term scale, was investigated empirically in financial markets; it has been supported by multiple studies [2,3] as a stylized fact of financial time series [4]. The asymmetric flow of information resembles an energy cascade found in conditions of turbulence. In a developed turbulent flow, the energy cascades from the macroscopic spatial scale, where energy is injected from the outside, to the microscopic spatial scale, where energy is dissipated as heat [5,6,7,8,9]. Gashghaie et al. investigated details of the self-similar transformation rule of the probability density function of price fluctuations and the nonlinear scaling law of the structure function (n-th moment of fluctuations), signifying the multifractality of the time series, in their study of the time series of foreign exchange. They pointed out the similarity of price changes in the financial time series to the velocity difference between two spatial points in turbulence [10,11]. The intermittency in turbulence is a phenomenon by which a laminar flow is interrupted by irregular bursts that occur suddenly. Such intermittency, which is frequently encountered in heterogeneous complex systems, is well known in financial markets as volatility clustering [4,12]. Intermittency at each time scale produces a characteristic hierarchical structure designated as multifractality [8,9].
In the developed turbulence, the process by which mechanically generated vortices on a macro scale deform and destabilize according to the Navier-Stokes equation and then split into smaller vortices is regarded as an energy cascade. A similar idea of modeling multifractal time series by a recursive random multiplication process from a coarse-grained scale to a microscopic scale has offered an attractive means of describing financial time series [13,14]. Chen et al. verified the statistics of multiplier factors in the random multiplication process of turbulent flow by empirical studies using measured data and numerical experiments of Navier-Stokes equations [15]. Results show that the multiplier factors connecting two adjacent layers follow a Cauchy distribution in which all moments diverge, and show that they are not independent. They show strongly negative correlation between the multiplier factors of adjacent layers. The authors verified the statistics of multipliers calculated backward from actual stock price fluctuations, finding a Cauchy distribution of multiplier factors and also the strongly negative correlation between the multiplier factors in financial markets. Results show that the discrete cascade model using the random multiplication process did not reproduce the statistical property of the multiplier factors. Therefore, as an alternative model, a discrete random multiplicative cascade process with additional additive stochastic process [16,17,18], or a model formulated as the Fokker-Planck equation considering the cascade process as a continuous Markov process [19,20,21,22,23] was proposed. Those models have been applied to stock market or foreign exchange market data, yielding empirical results including the statistics of multipliers.
This study examines a continuous cascade model of volatility formulated as a stochastic differential equation including two independent modes of Brownian motion: one has multiplicative coupling with volatility; the other has additive coupling as in the discrete random multiplicative cascade process with additional additive stochastic processes described above. The model parameters are estimated by its application to the stock price time series. Numerical calculation of the Fokker-Planck equation derived from the stochastic differential equation is conducted using the estimated values of parameters resulting in successful reproduction of the pdf of the empirical volatility and the multifractality of the time series.
where function ψ is designated as the analyzing wavelet. When using the delta function ψ(t) = δ(t − 1) − δ(t) as the analyzing wavelet, the wavelet transform is exactly the logarithmic return of the period s.
Here we use the second derivative of the Gaussian functions as In general, by using the n-th derivative of the function having asymptotic fast decay as the analyzing wavelet, one can remove the local trend of m-th order (m ≤ n − 1) because the function is orthogonal to m-th order polynomials. For the second derivative of the Gaussian functions, the linear trends of Z(t) with scale s have been eliminated in the wavelet transform W ψ Z[u, s].
The quantity used herein is the absolute value of the wavelet transform x(λ) = |W ψ Z[t 0 , s(λ)]| for arbitrary t 0 , where we use the variable λ = log L/s. Quantity x(λ) is thought to be a generalization of empirical volatility, whereas wavelet transform W ψ Z[u, s] is exactly the absolute value of logarithmic return when we use ψ(t) = δ(t − 1) − δ(t).
The following stochastic equation is used to start.
The power law behavior of the q-th moment E[x(λ) q ] (q-th structure function) as a function of scale s is proved by the solution (5) as the following.
We introduce an additional additive stochastic process as we have done into the discrete cascade model. We first consider the following stochastic differential equation.
The equation is produced on the assumption that Brownian motions dB M (λ) and dB A (λ) are mutually independent. The first two terms correspond to equation (4). The origin of those random sources triggering volatility cascade in financial markets remains unclear.
To solve the stochastic differential equation (7), we consider the following stochastic differential equation of which is the same as (4). Using the solution of (8) the solution of (7) is expressed as shown below: 4
(2) When considering the three scales λ 1 < λ 2 < λ 3 (s 1 > s 2 > s 3 ), the adjacent multipliers W (λ 2 , λ 1 ) = x(λ 2 )/x(λ 1 ) and W (λ 3 , show strongly negative correlation. Here we show property (1) and infer the existence of correlation between adjacent multipliers under some reasonable approximations. Parameter σ M is an important model parameter for the signal to have multifractality. As presented in a later section, in spite of the importance, the value of the parameter σ 2 M is tiny, about 0.02 − −0.03 in stock markets, irrespective of the stock issue. To specifically examined the role of additional stochastic processes, we investigate the 0-th order approximation of small σ M . When setting σ M = 0, the solution (10) becomes Therefore, the difference ∆x(λ, λ 0 ) = x(λ) − x(λ 0 ) follows a normal distribution N ( λ0 λ a A (u)e γM (u−λ0) du, λ λ0 (b A (u)) 2 e 2γM (u−λ0) du). If one simply assumes that x(λ 0 ) follows a normal distribution, then the ratio ∆x(λ 0 , λ 1 )/x(λ 0 ) of two independent stochastic variables following normal distributions follows a Cauchy distribution. So x(λ 1 )/x(λ 0 ) is the same.

Relation to discrete random cascade model
Assuming that ∆λ is sufficiently small, then when we use the following approximation of Ito's stochastic integration [24] as we obtain the discrete random cascade equation as The conditional expectation value of the square of x(λ + ∆λ), as the function of x 2 (λ), (14) shows that deviation of the quadratic curve from the origin results from the parameter b A (λ), as demonstrated from an empirical study in [18].

Constraint condition from the pdf of x(λ)
A remarkable feature of the probability density function (pdf) of the quantity Fig. 1 for the data examined in this study (see also Fig. 10 for the pdf of x(λ)). It indicates the constraint condition as Derivation of the constraint condition (15) is given in the Appendix below.
The additional additive stochastic process in the model (7) is expected to be a small perturbation to the basic model (4) to avoid violating multifractality. We also impose the following condition for all scales s: The power law scaling shown in Fig. 1, and condition (16) show the following constraint condition: Inserting (18) into (15), we also have the equation 2.1.5 Appendix: Derivation of (15) We introduce some notation for simplification of the description: From equation (13), we have We also have Because of the coincidence of the expected value and the standard deviation, we have E 2 (λ) = 2E 2 1 (λ) and E 2 (λ + dλ) = 2E 2 1 (λ + dλ). Inserting those equalities, and using approximation e σ 2 M dλ = 1, we have the constraint condition (15).

Fokker-Planck equation
We can derive the Fokker-Planck equation for the stochastic process {x(λ)} expressed by the stochastic differential equation (7) as the following, which is the master equation that the density of the transition probability p(x, λ|x 0 , λ 0 ) follows [24].
Therein, the functions D 1 (λ, x) and D 2 (λ, x) are defined as The k-th moment of the change δx(λ) = x(λ + δλ) − x(λ) induced by the infinitesimal scale transformation δλ is derived as shown below.
Therein, we used the identity p(y, λ|x, λ) = δ(y − x). Coefficients D 1 (λ, x) and D 2 (λ, x) show a relation to the first and second moments of δx(λ) in the following way.
Coefficients D k are designated as Kramers-Moyal coefficients [24,25]. We use equation (23) to estimate the function a A (λ) and b A (λ) and parameters γ M and σ M . To validate the model (7), it is necessary to confirm vanishing of the k-th moments for 3 ≤ k in the limit of δλ → 0. Renner et al. proposed almost identical equations with (20) in the literature [20,21], in which they deal with the price change itself as an analogy of the velocity difference in turbulence [19]. They derived a Fokker-Planck equation as a result of their empirical studies using Kramers-Moyal expansion of the Chapman-Kolmogorov equation, regarding the process as a Markovian process.

Data
We analyze the normalized average of the logarithmic stock prices of the constituent issues of FTSE 100 index listed on the London Stock Exchange for November 2007 through January 2009, which includes the Lehman shock of 15 September 2008 and the market crash of 8 October 2008.

Data processing
First, we calculate the average deseasonalized return of each issue δZ i (t) = log(S i (t)) − log(S i (t − δt)), which describes the average change of the portfolio as where µ i and σ i respectively denote the average and the standard deviation of δZ i and where N F represents the number of constituent stock issues (stocks). The constituents of FTSE100 Index are updated frequently. We selected N F = 111 stocks that remained listed on the London Stock Exchange throughout the period. Here, we set δt = 1 (min) and examine the 1-min log-return. We excluded the overnight price change and specifically examine the intraday evolutions of returns. To remove the effect of intraday U shape patterns of market activity from the time-series, the return was divided by the standard deviation of the corresponding time of the day for each issue i. Then we cumulate δZ(t) to obtain the path of process Z(kδt) (k = 1, . . . , L) ( Fig. 2(A)) as 3 Results

Multifractal analysis
First, we analyze the multifractal properties of the path Z(t) using an approach with wavelet-based multifractal formalism proposed by Muzy, Bacry, and Arneodo [26,27]. Initially, we define two mathematical terms. The Hölder exponent α(x 0 ) of a function f (x) at x 0 is defined as the largest exponent such that there exists an nth-order polynomial P n (x) and constant C that satisfy for x in a neighborhood of x 0 , characterizing the regularity of the function f (x) at x 0 . The singular spectrum D(α) is the Hausdorff dimension of the set where the Hölder exponent is equal to α, as For multifractal paths, the Hölder exponent α is distributed in a range; for paths of the Brownian motion, which are fractal, D(0.5) = 1 and D(α) = 0 for α = 0.5. Muzy, Bacry and Arneodo proposed the wavelet transform modulus maxima (WTMM) method based on continuous wavelet transform of function to calculate the singular spectrum D(α). We briefly sketch the WTMM method in the Appendix below. We calculate the partition function Z(q, s) of the q-th moment of wavelet coefficients using equation (29) for the path of our data. Results are presented in Fig. 2(B). The partition function Z(q, s) for each order q shows power law behavior in the range of scales s/L < 2 −5 . Exponents τ (q) are derived by the equation (30). Figure 2(C) shows that it is a convex function of q. Those results underscore the multifractality of the data path. The singular spectrum D(α) derived as the Legendre transformation of the function τ (q) by equation (31) is a convex function that has compact support [0.25, 0.79] taking the peak at α = 0.53, as shown in Fig. 2(D).

Appendix: WTMM method
The WTMM method builds a partition function from the modulus maxima of the wavelet transform defined at each scale s as the local maxima of |W ψ [f ](x, s)| regarded as a function of x. Those maxima mutually connect across scales and form ridge lines designated as maxima lines. The set L(s 0 ) is the set of all the maxima lines l which satisfy The partition function is defined by the maxima lines as Assuming power-law behavior of the partition function one can define the exponents τ (q). The singular spectrum D(α) can be computed using the Legendre transform of τ (q):

Parameter estimations
a A and γ M Parameters a A (λ) and γ M are estimated by taking the limit is fitted by a linear function as where dλ = λ 1 − λ 2 . As shown in Fig. 3(a), the first moment is well fitted by a linear function. Fitting of this kind is applied to various λ 1 = log(L/s1) and λ 2 = log(L/s2) combinations (Fig. 3(B)). Taking the limit dλ → 0 (ds/s = (s 2 − s 1 )/s 2 → 0), one obtains a A (λ 2 ) = lim dλ→0 a A (λ 2 , dλ)(a A (s 2 ) = lim ds/s→0 a A (s 2 , ds/s)) and γ M = lim dλ→0 γ M (λ 2 , dλ)(γ M = lim ds/s→0 γ M (s 2 , ds/s)). Fig. 4(A) presents examples of a A (s 2 , ds/s) and nonlinear fittings by the function log(a A (s, ds/s)) = a + b(ds/s) + c(ds/s) 2 . We estimate a A (s) by a A (s) = exp(a) for each line. The result is presented in Fig. 4(B). The solid line is the least-squares fit a A (s) to a power law function as log(a A (s)) = −1.50(0.41) + 0.58(0.11) log s, where the standard errors are in parentheses. The estimated exponent 0.58(0.11) is consistent with the constraint condition (18) where the standard error is the value in the parenthesis.
b A and σ M Similarly, we estimate parameters b A and σ M by taking the limit is fitted by a quadratic function (a regression against x 2 ) as As shown in Fig. 6(a), the second moment is well fitted by a quadratic function. Fitting of this kind is applied to various λ 1 and λ 2 combinations (Fig. 6(B)). Taking the limit dλ → 0, one obtains b 2 where the standard errors are in parentheses. The estimated exponent 1.26(0.13) is slightly higher than the constraint condition (18)(b 2 A (s) ∼ s). However, it is acceptable with accuracy. By a similar extrapolation log(σ 2 M (s, ds/s)) = a + b(ds/s) + c(ds/s) 2 , we estimate σ 2 M (s) = exp(a). Figure 8(A) presents an example of σ 2 M (s 2 , ds/s). and an estimate σ 2 M (s) by σ 2 M (s) = exp(a) for each line. The result is shown in Fig. 8(B). We estimate parameter σ 2 M by the average value weighted by the reciprocals of the standard errors.
Therein the standard error is in the parenthesis.

Numerical calculation of the Fokker-Planck equation
We confirmed that estimation of the parameter a A (λ) and b A (λ) by the E((x 1 − x 2 ) k |x 2 = x)/(λ 1 − λ 2 (k = 1, 2) is consistent with the constraint condition (18) with accuracy. If one imposes the other constraint (19), then the parameters have the following functional form of a A (λ(s)) = ǫs 0.5 b 2 A (λ(s)) = 2.27ǫs, where ǫ is a small parameter. The consistent range of ǫ found by estimation (33) and (36) is 0.15 ≤ ǫ ≤ 0.34. To fix parameter γ M and σ M , we use the empirical value of the scaling exponent τ (q), which is fitted by the quadratic function τ (q) = −1 + 0.52q − 0.013q 2 (see Fig. 2(C)). One can derive τ (q) = −1+(γ M + 1 2 σ 2 M )q− 1 2 σ 2 M q 2 for the basic model (4) without additional stochastic processes. Again using the assumption of slight perturbation, then from the coefficients of the quadratic function, the parameters γ M and σ M are expected to exist respectively in the neighborhood of γ M = 0.51 and σ M = 0.026. Next we try the value of the parameters γ M = 0.51, σ M = 0.026 and ǫ = 0.16 for numerical calculation of the Fokker-Planck equation. Results are presented in Fig. 10. The initial pdf of the numerical calculation represented by the dashed line was based on the measured pdf on scale s = 128(min). In the initial values, the fine fluctuation was smoothed using a spline function with the rationale that small fluctuations in the measured pdf are attributable to the finiteness of the number of observations. The tails are extrapolated using a power function with index −4.9 which is obtained empirically. For time evolution, the fourth-order explicit Runge-Kutta method was used. The solid line is the calculation result obtained using the estimated value of the parameters γ M = 0.51, σ M = 0.026 and ǫ = 0.16. The dotted line is the result obtained when ǫ = 0. The difference between the two was very small. The results closely matched the actual pdf.
Using results of the numerical calculation of the pdf obtained at each scale, we calculate the scaling exponent τ (q) as shown below.
The result is presented in Fig. 11. No difference exists between the two numerical calculation results. Both curves are convex upward, indicating multifractal properties. Comparison with measured values is also good. These results, when combined with consideration of the statistics of multipliers given in 2.1.2, underscore the effectiveness of the continuous cascade model (7) with additive stochastic processes proposed.

Discussion
The random cascade model has evolved as a model of developed turbulence. The original model, in which the stochastic process that connects each layer of the spatial scale is an independent random multiplication process, contradicts results obtained through empirical research. Therefore, an improved discrete random multiplicative cascade model with additional additive stochastic process was proposed along with a model formulated as a Fokker-Planck equation by considering cascade processes as a continuous Markov process. Moreover, those models have been applied to data analysis of the stock market and the foreign exchange market, where they have been successful. Herein, we propose a continuous cascade model formulated as a stochastic differential equation of volatility including two independent modes of Brownian motion: one has multiplicative coupling with volatility; the other has additive coupling, as in an improved discrete cascade model for the stock market, with effectiveness clarified by results of earlier research [18]. The model parameters were estimated by application to a stock price time series. The Fokker-Planck equation was derived from the stochastic differential equation as a master equation with the transition probability density function of volatility. Furthermore, the model parameters were estimated by its application to the average stock price time series made from FTSE100 constituents listed on the London Stock Exchange. At that time, as an alternative variable of volatility, the wavelet transform coefficient with the second derivative of the Gaussian function as an analyzing wavelet was used. Numerical calculation of the Fokker-Planck equation was conducted using the estimated parameter values. The results reported herein faithfully reproduce the results of an earlier empirical study. This model includes information about neither the time axis nor the sign of the price fluctuation, which are necessary for a model of price fluctuations. The actual stock market exhibits well known properties that break symmetry with respect to the time axis, such as the causal structure from long-term to short-term scale volatility described first in the Introduction and price-volatility correlation (Leverage effect) [4,12]. Therefore the extension of the random cascade model to encompass these phenomena remains as a subject for future work. [