Mixed-Stable Models: An Application to High-Frequency Financial Data

The paper extends the study of applying the mixed-stable models to the analysis of large sets of high-frequency financial data. The empirical data under review are the German DAX stock index yearly log-returns series. Mixed-stable models for 29 DAX companies are constructed employing efficient parallel algorithms for the processing of long-term data series. The adequacy of the modeling is verified with the empirical characteristic function goodness-of-fit test. We propose the smart-Δ method for the calculation of the α-stable probability density function. We study the impact of the accuracy of the computation of the probability density function and the accuracy of ML-optimization on the results of the modeling and processing time. The obtained mixed-stable parameter estimates can be used for the construction of the optimal asset portfolio.


Introduction
The increased availability of high-frequency data has caused a great interest in the research of this subject. The main applications belong to financial engineering, ranging from risk management to options hedging, transaction execution, portfolio optimization, and forecasting. Thus, Bailey and Steeley [1] compared forecasts of the volatility of the Australian dollar exchange rate to alternative measures of ex post volatility, using high-frequency data. Degiannakis and Filis [2] examined the importance of combining high-frequency financial information, along with the oil market fundamentals, to gain incremental forecasting accuracy for oil prices, showing that although the oil market fundamentals are helpful for long-run forecasting horizons, the combination of the latter with high-frequency financial data significantly improve oil price forecasts. Zhang and Wang [3] employed the MIDAS model and the high-frequency data of four stock market indices to forecast WTI and Brent crude oil prices at a lower frequency. The results indicated that high-frequency stock market indices have a certain advantage over the lower-frequency data in forecasting monthly crude oil prices, and the MIDAS model using high-frequency data proves superior to the ordinary model. Göncü and Yang compared variance-gamma and normal-inverse Gaussian distributions with the benchmark of generalized hyperbolic distribution in terms of their fit to the empirical distribution of Chinese high-frequency stock market index returns, showing that as the time scale of log-returns decrease, a normal-inverse Gaussian model consistently outperforms the variance-gamma model [4]. Remarkably, this result for the normal-inverse Gaussian model is consistent with findings of Belovas from the the same year [5], based on five Standard & Poor's stock market indices (covering the period of 10 years, 2006-2016) log-returns. Koopman et al. [6] investigated how dependence between high-frequency price changes of financial stocks (10-second frequency for 10 US financial stocks for January 2012 to Decemeber 2012) varies within the day. Schabek et al. [7] examined high-frequency log-returns of the Zagreb Stock Exchange CROBEX Index (from September 2017 to March 2018) to assess the reactions of the index to macroeconomic announcements regarding the Croatian economy within an ultrashort time periods. Tony Cai et al. [8] examined the estimation of a high-dimensional minimum-variance portfolio based on the high-frequency returns from S&P 100 Index constituents during the years 2003-2013. Huang and Gao [9] used high-frequency Bitcoin trading data (from 1 January 2012 to 12 August 2019) to explore the Bitcoin return predictability.
Ongoing COVID-19 turmoil has set a new trend in high-frequency studies. Thus, Ambros et al. [10], using 30 min tick returns, studied the impact of changes in the number of COVID-19 news on eight different stock markets during the initial two months of the coronavirus crisis, showing that COVID-19-related news impact stock market volatility positively in Europe, but less so in other markets. Yousaf and Ali [11] analyzed return and volatility transmission between major cryptocurrencies (Bitcoin, Ethereum and Litecoin) during the pre-COVID-19 and COVID-19 periods (it is noteworthy that the authors advised the investors to decrease their investments in Bitcoin).
Methodologies based on high-frequency data analysis also can be found in neural science and real-time network traffic management (cf. Kaklauskas [12]). A summary of the literature covering high-frequency and intra-daily data research is presented in [13] and the references therein.
The paper extends the study of applying the mixed-stable models to the analysis of large sets of high-frequency financial data; see [4,14,15]. In this research, we apply the parallel computing approach (cf. [16]) to the mixed-stable modeling of high-frequency data. We often observe many zero returns in the high-frequency return data in practice because the underlying asset price does not change at given short-time intervals. The mixed-stable model is well suited to cope with this specific feature.
We introduce the smart-∆ approach to the calculation of the α-stable probability density function and consider the impact of the accuracy of the computation of the probability density function and the accuracy of the maximum-likelihood optimization on the results of the modeling and processing time.
The mixed-stable model for financial data was first applied in [17]. The preliminary research [18] was dedicated to the analysis of empirical data covering one specific day. In the present study, we analyze high-frequency data for the whole year. This will drastically increase the size of data sets and computing time. We address this issue using more efficient numerical methods and parallel algorithms.
The paper is organized as follows. The first part is the introduction. Section 2 describes the real data used in our research, their aggregation and their transforms. Section 3 introduces our modeling methodology and smart-∆ method for calculating the α-stable probability distribution function. In Section 4, we discuss the impact of computations' and optimization' accuracy on the modeling and processing time and present the modeling results. The last section is devoted to the concluding remarks.

Data
In the previous research of intra-daily data from German DAX component stock returns ( [18]), we analyzed high-frequency data series of 29 stocks from DAX that represent just one business-active day of the year (17 August 2007). The present paper extends the research and deals with the whole year's intra-daily data (from 1 January 2007 to 27 December 2007; 251 days in total). The year 2007 is of particular interest to economists and financial analysts. Moreover, empirical data of this year comprise a valuable test case in creating and testing special models, because 2007 is the first year of the Global Financial Crisis of 2007-2008. Before the current COVID-19 turmoil, it was considered by many experts to have been the most severe financial crisis since the Great Depression, contributing significantly to the Eurozone crisis.
We aggregate raw inhomogeneous intra-daily data into equally-spaced homogeneous intra-daily time series. The aggregation is done with the previous-tick interpolation. A linear interpolation relies on future information, whereas the previous-tick interpolation is based on the information already known (cf. [19]).

Previous-Tick Interpolation
We denote times of raw inhomogeneous intra-daily series as {t i } and the corresponding prices as {P i }, where P i = P(t i ). The aggregated homogeneous high-frequency series {P j } is obtained at timest j = t 0 + j∆t with the step ∆t, where the index j identifies the regularly spaced sequence. By means of the previous-tick interpolation, we obtain that Having obtained equally-spaced price series, we can calculate the corresponding logarithmic returns series {X j }:

Models for Financial Data
Classical techniques in financial engineering heavily relied on the assumption that the random variables under investigation follow a normal distribution. However, time series observed in finance often deviate from the Gaussian model, exhibiting fat tails and asymmetry (cf. [18,19]). In such a situation, the classical approach's appropriateness for the modeling of returns is highly questionable. On the other hand, financial asset returns are the cumulative outcome of a vast number of pieces of information and individual decisions arriving almost continuously. Hence, in the presence of heavy tails, it is natural to assume that they are approximately governed by a non-Gaussian distribution (cf. [20]).
We can distinguish several fundamental reasons why models with the α-stable paradigm are used in financial engineering. The first is that stable random variables justify the generalized central limit theorem, which states that stable distributions are the only asymptotic distributions for adequately scaled and centered sums of independent identically distributed random variables (see [20]). The second one is that stable distributions are heavy-tailed (cf. [18,19]). All but one of the stable distributions have infinite variance, which implies that observations of a large magnitude can be expected and may dominate the sums of these random variables. It is not correct to treat these observations as outliers, since their exclusion takes away much of the significance of the original data; indeed, it is specifically these observations that may be of the most significant interest. The third reason is that stable distributions are asymmetric and leptokurtic [21]. Since stable distributions can accommodate the heavy tails and asymmetry, they ensure a perfect fit for empirical data. They are particularly valuable models for data sets covering extreme events such as market crashes or natural catastrophes. The fourth reason is that stable distributions are a more flexible tool compared to the normal distribution. As was pointed out by Cont, a parametric model to successfully reproduce specific empirical features of asset returns must have at least four parameters: a parameter describing the decay of the tails (stability index), an asymmetry parameter allowing the left and right tails to behave differently, a scale (volatility) parameter and ultimately a location parameter [22]. Our recent research on the comparison of models corroborates this opinion [15].

Empirical Moments
Having processed yearly high-frequency returns data for 29 stocks at different time steps, we found that almost all data series are asymmetric (some typical examples may be found in Table 1), and the empirical kurtosis shows that density functions of the series are more peaked than Gaussian density functions. However, it should be pointed out that rather often, empirical data exhibit the stagnation effect; i.e., series contain numerous zero returns. This phenomenon is especially characteristic of young emerging markets with low-liquidity financial instruments [17] and high-frequency financial data [18]. We have examined the obtained returns series for the stagnation effect. In Table 2, we show max and min lengths of returns series with zeros removed as well as max and min percent of zeros depending on the level of aggregation. As we can see in Table 2, a strong stagnation effect (43 to 82 percent zeros at 10 s time step) manifests itself in most high-frequency series. To take the zero effect into account, we apply a generalized mixed-stable model.

α-Stable Distribution
The α-stable distribution is usually described by its characteristic function ϕ(t): where α ∈ (0, 2], β ∈ [−1, 1], σ > 0, µ ∈ R. Here, α is the stability index (in financial modeling, it is generally assumed that 1 < α ≤ 2), β is the skewness, µ is the location parameter and σ is the scale parameter. If σ = 1 and µ = 0, then the distribution is called standard stable. An overview of stable distributions properties can be found in [20]. The probability density function of stable laws cannot be expressed in elementary functions, except for a few cases: Levy, Cauchy and Gaussian distributions. The canonical representation (2) has one serious disadvantage. Characteristic functions have discontinuities at all points with α = 1, β = 0. Therefore, for numerical purposes, it is advisable to use Nolan's parametrization: This parametrization is a variant of Zolotarev's (M) parametrization, with the density and the distribution function jointly continuous in all the four parameters ( [20]). The location parameters of these two representations are related by The probability distribution function of representations (2) and (3) are related by Here, Θ = (α, β, µ, σ), and p 0 (x, α, β) is Nolan's standard stable probability distribution function with the integral representation where h(x, t; α, β) = xt + β tan πα 2 (t − t α ), α = 1, xt − βt 2 π ln t, α = 1.
A precise and fast calculation of stable densities is a nontrivial task (cf. [16,20]). To deal with the integral representation of the probability density function of α-stable distribution (5), we introduce the smart-∆ approach, replacing the improper integral in (5) by a definite integral with the upper integration bound ∆ = ∆(α, ε). Here, α is the stability index, and by ε we denote the error of the approximation. Details of the technique are explained in Section 3.2.
The problem of parameter estimation in stable modeling is hampered by the lack of closed form for stable density functions. Hence, many statistical methods depending on the probability density function's explicit form cannot be applied. Comparative studies (see [23]) corroborate that the most accurate method of estimation is the maximum likelihood method. However, it is the most time-consuming. The vector of stable parameters Θ = (α, β, µ, σ) can be estimated from the returns {X j } by maximizing the log-likelihood function In [16], we have studied the log-likelihood target function profiles for artificially generated stable distributed data. We have obtained that the log-likelihood target function is of an uniextremal nature, often with a very flat surface in the extremum's neighborhood.
In the present research, we have examined the log-likelihood target function for real financial data. Target function (6) was calculated for returns series Allianz SE of size 7385, which was obtained using the previous-tick interpolation (1) with the time step ∆t = 1000. The ML solution vector for this data set is α = 1.541470, β = 0.004055, µ = 0.000005, σ = 0.001241. Figure 1 shows 3D cuts of the target function, obtained by fixing pairs of parameters. As one can see, the target function with real data exhibits qualitatively the same behavior as with the artificially generated data (cf. [16]).  To optimize the log-likelihood function, we use the Nelder-Mead method. Although this method is not the fastest one, it is robust and does not require the calculation of derivatives (gradient or Hessian).

Stable Probability Density Function Calculation: The Smart-∆ Approach
The lack of analytical representation of the probability distribution function (with a few exceptions: Gaussian, Cauchy and Levy distributions) hampers the practical implementation of stable models. We can evaluate the stable probability density function replacing the improper integral (5) with an approximation I 0 , with a tail error Calculating I 0 with ε 2 precision and dropping I ∆ evaluated with the same accuracy yields ε joint accuracy for the probability density function. For α > 1, as is usually assumed in financial engineering, we have Hence, the roughest way of evaluating ∆ (see [16]) is The error of the approximation does not exceed ε/2. Noticing a relation of the bound of I ∆ to an upper incomplete gamma function Γ(s, x), we can evaluate ∆ = ∆(α, ε) in a more subtle way, as a root of a nonlinear equation We can calculate ∆ as follows. For x → ∞, we have Γ(a, x) ∼ x a−1 e −x , yielding an approximate equation Next, ∆ can be expressed in terms of the Lambert W function, i.e., Here, W(x) stands for the Lambert W function (y = W(x) ⇔ x = ye y ). For x ≥ 0, we can calculate the principal branch of the Lambert W function using Halley's method (cf. Corless et al. [24]). If x < 0 (i.e., α < 1), then we apply an algorithm, proposed by [25], to calculate the branch W −1 . Note that in order to proceed from a standard stable density to a stable density, we interchange in expression (9) π coefficients with πσ.

Mixed-Stable Distribution
The mixed-stable model was introduced to deal with the problem of daily zero returns ( [17]). The probability density function of a mixed-stable random variable is where p(x, Θ) is the probability density function (4) of a stable distribution (2) and δ(x) is the Dirac delta function. The coefficient r ∈ (0, 1) is the index of stagnation. The empirical cumulative distribution functions of data series with the stagnation effect exhibit jumps at x = 0. Model (10) enables us to accommodate these jumps. The likelihood function of the mixed-stable model (10) is given by where {x 1 , x 2 , . . . , x n−k } is a non-zero returns set, obtained by excluding k zero returns from the initial returns set {X 1 , X 2 , . . . , X n }. By optimizing the first factor (1 − r) n−k r k , we obtain r max = k/n. The optimization of the product is equivalent to the optimization of the likelihood function of the stable distribution of the non-zero returns set {x 1 , x 2 , . . . , x n−k }. Hence the optimal vector Θ max is estimated with non-zero returns via stable log-likelihood function (6).
Having parameters of the mixed-stable law estimated, we proceed with modeling adequacy testing. Since we have a discontinuous distribution function, classic methods for continuous distributions (e.g., Kolmogorov-Smirnov or Anderson-Darling tests) are unsuitable. Therefore, we apply a special goodness-of-fit test based on characteristic functions (note that characteristic functions are uniformly continuous on the entire space) proposed in [26].

Results and Discussion
Using the maximum-likelihood method (see Section 3), we have estimated parameters of mixed-stable models (10) for DAX financial data. Next, we studied the impact of the accuracy of the probability density function (4) calculation (ε pd f ) and the maximumlikelihood optimization (ε ML ) on the results of the modeling and processing time. We found that insufficient accuracy results in faulty outcomes. We illustrate the facts in Table 3, which contains Θ = (α, β, µ, σ) estimates for the Deutsche Post AG data series, taken with time step ∆t = 10 s. For every set of estimates, the corresponding processing time (PT) and the outcome of the Koutrouvelis goodness-of-fit test ( [26]) for the adequacy of the modeling (KT) are indicated (the significance level of the test is 5%). As we see, we need at most ε pd f = 10 −9 and ε ML = 10 −6 to achieve plausible results. Further increase in accuracy levels shows the convergence of obtained estimates. For higher accuracies, at least seven significant digits are not changing. Table 3. Dependence of mixed-stable estimates on the accuracy of the probability density function (4) calculation (ε pd f ) and the accuracy the maximum likelihood optimization (ε ML ); Deutsche Post AG returns data series with time step ∆t = 10 sec. Processing time (PT) (sec) and Koutrovelis test outcome (KT, 5% significance level).  Note that because of the flat surface of the ML-target function (see Figure 1), the initial guess selection is essential for fast and correct optimization. Thus, the global optimization (either deterministic or stochastic) is the natural way of solving this problem. However, it requires further in-depth investigation.
An increase in accuracies ε pd f and ε ML causes a surge in processing time (see Table 3). In Table 4, we show the optimization time for 29 DAX high-frequency returns series taken with the time step ∆t = 10 s. These results were obtained using parallel algorithms introduced with early-stage research ( [14]) with 64 processes. Note that when ε pd f ∼ ε ML , convergence is naturally very bad. Next, ε pd f and ε ML cannot be chosen independently, We must select ε pd f in relation to ε ML (e.g., ε pd f = 10 −12 is sufficient for ε ML = 10 −7 ; however, the convergence stagnates if we choose ε pd f = 10 −12 for ε ML = 10 −8 ). The number of stagnated series for every accuracy pair (ε pd f and ε ML ) is given in parentheses next to the corresponding overall processing time (see Table 4). Table 4. Dependence of 29 DAX returns series processing time (sec) with 64 processes (for time step ∆t = 10) on the accuracy of the probability density function (4) calculation (ε pd f ) and the accuracy of the maximum likelihood optimization (ε ML ). The number of stagnated series for every accuracy pair (ε pd f and ε ML ) is given in parentheses. Finally, we present the mixed-stable modeling results for 29 DAX companies. These results were obtained with ε ML = 10 −7 accuracy of the optimization method and ε pd f = 10 −12 accuracy of mixed-stable density function calculation. We must stress that we could obtain these results in a reasonable time using only parallel computations. Table 5 contains estimates of a stagnation parameter r, stable parameters Θ = (α, β, µ, σ) and the outcome of the Koutrouvelis test (zero stands for "rejected", unity stands for "not rejected").
The mixed-stable model, with a corresponding set of estimated parameters, was accepted for almost all DAX companies, justifying the application of mixed-stable models in high-frequency finance analysis (see Table 5). Note that considerable differences in corresponding parameter sets are obtained from yearly (cf. Table 5) and daily (cf. [18]) empirical data. This indicates the necessity of large data sets processing in forecasting and portfolio construction.
Examples of the use of the estimated mixed-stable parameters for the selection of the optimal asset portfolio have been provided in our preliminary comparative research [15]. This study compared stable and mixed-stable models with mixed diffusion-jump, the mixture of normals, scaled-t, logistic and normal-inverse Gaussian models and identified the mixed-stable one as the most adequate model for the data under analysis. The modeling results can be applied for the optimal portfolio selection employing two different strategies (considered in the research) with and without the relation coefficients matrices. Table 5. Maximum-likelihood estimates of mixed-stable parameters for 29 DAX returns series with time step ∆t = 10. The accuracy of the calculation of the probability density function ε pd f = 10 −12 and the accuracy of the maximum likelihood optimization ε ML = 10 −7 . KT stands for the outcome of the Koutrouvelis test (0, "rejected"; 1, "not rejected").

Conclusions
Having processed yearly high-frequency returns data for 29 German DAX companies with different time steps, we found that almost all data series are asymmetric. Moreover, the empirical kurtosis shows that density functions of the series are more peaked than Gaussian. We have noticed a stagnation effect in obtained high-frequency returns series. These factors lead us to the application of mixed-stable models. We introduced the smart-∆ upper integration bound to deal with the computationally demanding α-stable probability density function in ML-optimization. Parallel algorithms were employed to deal with sizeable yearly data sets.
We have studied the impact of pdf-computation accuracy and the accuracy of MLoptimization on the results of the modeling and processing time. We constructed mixedstable models for all 29 DAX companies. The adequacy of models was tested with Koutrouvelis goodness-of-fit test based on the empirical characteristic functions. Almost all models were accepted with corresponding sets of estimated parameters, justifying the mixed-stable modeling in high-frequency data analysis. Obtained parameter estimates can be used in the construction of the optimal asset portfolio.
The next research objective is to compare the presented results (concerning the first year of the Global Financial Crisis of 2007-2008) with the planned analysis of ongoing COVID-19 crisis data. Two more points to be brought to the attention are the robustness of empirical results for different periods and comparative studies with other markets.