Abstract

This paper aims to efficiently implement the maximum likelihood estimator (MLE) for Hurst exponent, a vital parameter embedded in the process of fractional Brownian motion (FBM) or fractional Gaussian noise (FGN), via a combination of the Levinson algorithm and Cholesky decomposition. Many natural and biomedical signals can often be modeled as one of these two processes. It is necessary for users to estimate the Hurst exponent to differentiate one physical signal from another. Among all estimators for estimating the Hurst exponent, the maximum likelihood estimator (MLE) is optimal, whereas its computational cost is also the highest. Consequently, a faster but slightly less accurate estimator is often adopted. Analysis discovers that the combination of the Levinson algorithm and Cholesky decomposition can avoid storing any matrix and performing any matrix multiplication and thus save a great deal of computer memory and computational time. In addition, the first proposed MLE for the Hurst exponent was based on the assumptions that the mean is known as zero and the variance is unknown. In this paper, all four possible situations are considered: known mean, unknown mean, known variance, and unknown variance. Experimental results show that the MLE through efficiently implementing numerical computation can greatly enhance the computational performance.

1. Introduction

Signals of nature [1–6], medicine [7–14], business [15–18], and society [19–22] usually appear to be a strong long-term correlation. These signals can be differentiated by only one indicator, fractal dimension or Hurst exponent; therefore, many researchers are attracted to the study of how to estimate fractal dimension or Hurst exponent. In order to analyze the characteristics of fractal signals, users can determine the fractal dimension (). Among estimators, the box-counting technique [23–27] is a direct nonmodeling method. In general, engineers are fond of adopting indirect modeling methods like fractional Brownian motion (FBM) or fractional Gaussian noise (FGN), because they are more meaningful than direct nonmodeling methods. FBM or FGN first estimates the Hurst exponent (), a real number in (0, 1), and then calculates the fractal dimension via the relation [28]. The Hurst exponent is the only one parameter dominating the characteristics of FBM or FGN. FBM is a statistically self-similar nonstationary random process, which makes analysis difficult [29, 30], but the increment of FBM, FGN, is a strict-sense stationary process and has power spectral density (PSD) behaving asymptotically as [29, 31].

In real applications, signals must be sampled in advance; sampled FBM is called discrete-time fractional Brownian motion (DFBM), and sampled FGN is called discrete-time fractional Gaussian noise (DFGN), which has been proven to be a regular process [32]. Many natural and biomedical signals can be modeled as DFBM or DFGN [7–12, 33]. Among estimators, the maximum likelihood estimator (MLE) [29] provides the optimal accuracy; one of its approximate versions, called the Whittle estimator [34, 35], provides the second optimal accuracy. The aim of the Whittle estimator is to provide faster estimation with slight inaccuracy. Other quick versions include the variance method [29, 31, 36], moving-average (MA) method [33], and autoregressive (AR) method [32, 37].

Although the accuracy of estimating the Hurst exponent by the MLE is the best, it is easy to induce computational problems and enormous computational expenditure. For example, evaluating the inverse of an autocovariance matrix may be numerically unstable, especially when is close to 1. Under this situation, the autocovariance matrix almost becomes singular because the autocovariance matrix of DFGN changes very slowly [35]. This problem will cause computational inaccuracy, leading to wrong explanations for physical signals of interest. On the other hand, the cost often makes users hesitate to apply the MLE to quick response systems, and thus the theoretical value of the MLE is generally much higher than its practical applications.

When taking a closer look at the structure of the autocovariance matrix, a combination of the Levinson algorithm [38] and Cholesky decomposition [39] can solve computational problems and reduce computational cost. Accordingly, users will be encouraged to adopt the MLE even in the quick response situations, and then the MLE has a better opportunity to become the first choice in the future, especially when computer speed continues to be increased up to a certain level.

When the MLE was first proposed by Lundahl et al. [29], the analysis and evaluation of the MLE were based on the assumptions that the mean of DFGN is zero and the variance is unknown. It is only applicable to physical signals modeled as DFBM, but not suitable for the model of DFGN. When signals are modeled as DFGN, it is easy for users to obtain wrong estimation results, unless they take the sample mean out of the original signals beforehand. Therefore, it is necessary for practical signals to give a complete consideration of four possible cases, including known mean, unknown mean, known variance, and unknown variance; moreover, each unknown mean also considers two approaches for comparison: the sample mean and the mean estimated by MLE. In terms of the practical situation of a realization of physical signals, users can choose one case to estimate the Hurst exponent and then the fractal dimension.

The rest of this paper is organized as follows. Section 2 briefly describes mathematical preliminaries. Section 3 introduces practical considerations for the MLE. Section 4 shows how to implement the MLE in an efficient way. Section 5 discusses experimental results. Finally, Section 6 concludes the paper with some facts.

2. Mathematical Preliminaries

In this section, some related models are reviewed, including FBM, FGN, DFBM, and DFGN. For consistency, the notation is used to denote a continuous-time random process and a discrete-time random process.

FBM, represented by where belongs to a sample space , is a generalization of Brownian motion. For conciseness, the short notation is adopted in place of . According to Mandelbrot and Van Ness [31], FBM is formally defined by the following relations: where is the Hurst exponent with a value lying between 0 and 1 and the increments of FBM, , are zero mean, Gaussian, and independent increments of ordinary Brownian motion. Its symmetric form is described as follows: When equals 0.5, FBM becomes the ordinary Brownian motion. Unfortunately, FBM is a nonstationary process, whose Wigner-Ville spectrum (WVS) is given by the following expression [30]: In spite of FBM being a time-varying process, the increment of FBM is a stationary and self-similar process, called FGN.

In real applications, discrete data are collected; sampled data of FBM are expressed as , where is the sampling time. The increments of DFBM, called DFGN, are denoted by . DFGN is a normally distributed and stationary process with zero mean, whose autocorrelation functions (ACFs) are given by the following equation: where [29, 40]. The ACF, , behaves asymptotically as , [40].

3. Practical Considerations for the MLE

It is well known from the properties of DFGN that the probability density function (PDF) of DFGN can be expressed as follows [29]: where is the dataset and is the autocovariance matrix; that is, or , where is the ACF as expressed by (4).

In real applications, some physical signals can be modeled as either DFBM or DFGN. If the signals of interest are modeled as DFBM, their increment, DFGN, will not be affected by displacement. However, if the signals of interest are modeled as DFGN, signal displacement will result in a very severe error unless the displacement problem of signals is considered in advance. The reason for displacement may be modeling error, measurement error, inappropriate operation, or apparatus baseline calibrating error, and so forth. In order to avoid the error resulting from displacement, two approaches are considered to estimate displacement: one is to maximize PDF over the mean; the other is to simply take the sample mean out of signals. Considering that the PDF of DFGN has two explicit parameters, the mean and variance, each parameter may be known or unknown, and each unknown mean includes two approaches, all together, there are four cases covering six approaches.

3.1. Case 1: Known Mean (Displacement) and Known Variance

Under this case, there are no mean and variance necessary to be estimated before estimating the Hurst exponent. In theory, this is the best case since information about the mean and variance is provided. For convenience, the logarithm of PDF will be maximized instead of PDF, which produces the same result since logarithmic operation still preserves the monotonic property of a function. Without loss of generality, displacement is set to be 0. From (5), the logarithm of PDF is as follows:

Since constant terms and coefficients do not affect the maximum, a compact form is described as follows: where and is known to users.

3.2. Case 2: Known Mean (Displacement) and Unknown Variance

The case first proposed by Lundahl et al. [29] is like this. Likewise, displacement is assumed to be 0. The Hurst exponent can be estimated by using the following equation:

It is well known that the logarithm of PDF is expressed as follows: By maximizing the over , it follows that By substituting (11) into (10), the final function to be maximized is Likewise, the terms that do not affect maximization will be omitted, and thus a compact form is described as follows:

3.3. Case 3: Unknown Mean (Displacement) and Known Variance

Let measurement data be , where can be modeled as DFGN with zero mean and is a column vector with each element being constant ; that is, . The Hurst exponent can be estimated by using the following two approaches based on two estimators about .

Approach 1. First maximize the logarithm of PDF over by taking derivative with respect to , and then maximize the logarithm of the maximum PDF on estimated over , that is,
The unknown displacement of DFGN is assumed to be , and thus the PDF will be where . Therefore, (15) can be simplified as
First, maximize the over by taking derivative with respect to and the operation is equivalent to maximizing the . The estimator of is derived from the Appendix as follows: where and . It is easy to check that where and . Next, by substituting (17) into (16), the final function to be maximized is Likewise, the terms without affecting maximization are omitted, and thus a compact form is described as follows:

Approach 2. Use the sample mean to replace the previous estimator of . Other procedures are the same as the ones of Approach 1. The sample mean is the simplest method to estimate the mean, which is

3.4. Case 4: Unknown Mean (Displacement) and Unknown Variance

This case is the most general in real applications. Like Case 3, measurement data are assumed to be and the unknown variance is . Similarly, the Hurst exponent is estimated by using the following two approaches.

Approach 1. Like Case 3, how to estimate the Hurst exponent is described as follows:
First, maximize the over and by taking derivatives with respect to and , respectively, and then the estimators of and are derived as follows: Likewise, the terms that do not affect maximization are omitted, and thus a compact form is described as follows:

Approach 2. Use the sample mean to replace the previous estimator of . Other procedures are the same as the ones of Approach 1.

The final step of each case is to estimate the Hurst exponent, but it needs some tips. A direct maximization over is unfeasible because the Hurst exponent is an implicit parameter. Therefore, the golden section search [41] was adopted to find out the maxima of (7), (13), (20), and (24) in this paper.

4. Efficient Procedures for the MLE

In this section, the computational stability and efficiency of using the MLE for the Hurst exponent are studied. Since computing the inverse and determinant of an autocovariance matrix is sensitive to the data size, especially when is close to 1 [35]. Also, for a large dataset, storing the whole autocovariance matrix requires a large amount of computer memory. Therefore, a reliable and efficient procedure is necessary for estimating the Hurst exponent, especially when users use an ordinary computer with less memory and lower CPU speed. After carefully studying the structure of an autocovariance matrix, a combination of the Levinson algorithm and Cholesky decomposition can be applied to efficiently compute the inverse and determinant of an autocovariance matrix, and then the iterative structures of the two algorithms can be exploited to estimate the Hurst exponent without storing any matrix and performing any matrix multiplication. For convenience, some notations are listed below for further quotation: where are the predictor coefficients of order and are the prediction error powers of order . These coefficients are iteratively computed by the Levinson algorithm. It is worth noting that when numerically calculating is computed instead of because the term easily approaches to zero numerically as the data size grows larger. Thus, using the following equation to compute is essential:

Obviously, the Levinson algorithm and Cholesky decomposition can be used to save the time of computing the inverse and determinant of the autocovariance matrix. However, estimating the Hurst exponent by the currently mentioned structure of implementation still needs matrix computation and storing, which requires excessive computational time and computer memory. When taking a closer look at the term of (7) or (13), a magic and helpful form appears as follows: where

With (31), storing any matrix in the process of computation is no longer necessary, which is also a very efficient step.

In this paper, the golden section search was adopted to find out the maxima of (7), (13), (20), and (24). In the process of searching for each maximum, computing the inner terms of (7), (13), (20), or (24) is necessary, such as or .

In order to efficiently estimate the Hurst exponent by using (7) or (13), first compute , , by using the Levinson algorithm, then by using (26), by using (30) and (31), by using (29), and by using (28). The details of determining the Hurst exponent by using the MLE are described in the following procedure.

Procedure 1. By efficiently computing or , consider the following stpes:(1)initialize ;(2)compute and by using the Levinson algorithm;(3)compute by using (26);(4)compute ;(5)add a new element into the vector like (30);(6)perform ; if , then go to Step 2 or go to the next step;(7)compute using (29) and using (28);(8)compute for (7) or for (13).

Obviously, it is unnecessary to store any matrix and execute any matrix multiplication in a series of computation except for vector storing and multiplication. The efficient procedure not only saves computer memory but also storing time.

Next, an efficient procedure for computing (20) or (24) is considered. Each procedure considers two approaches to estimating the mean: the sample mean and the mean estimated by MLE. For the first approach, users first take the sample mean out of the original signals and then call the function evaluation of (7) or (13). For the second approach, users must use an efficient implementing procedure for (17) to estimate the mean. Based on the composition structure of the mean estimated by MLE, (17) can be decomposed into the following equation: where

When carefully observing the term , it can be decomposed into the following equation: where

In order to efficiently compute (20) or (24), first, compute , by using the Levinson algorithm, then, by using (26), by using (34), by using (36), by using (32), by using (37), by using (35), and by using (28). The details of determining the Hurst exponent by using the MLE are described in the following procedure.

Procedure 2. By efficiently computing or , consider the following steps:(1)initialize ;(2)compute and using the Levinson algorithm;(3)compute by using (26);(4)compute and ;(5)add a new element into the vector like (34) and like (36);(6)perform ; if , then go to Step 2 or go to the next step;(7)compute by using (32) and by using (37);(8)compute by using (35) and by using (28);(9)compute for (20) or for (24).

Similar to Procedure 1, it is not necessary to store any matrix and execute matrix multiplication in a series of computation except for vector storing and multiplication. Without considering the efficient computation of , users first implement the Levinson algorithm to obtain and and then store and , obtain by using , and finally take matrix computation for . The traditional procedure only saves the time of inverting matrix by using the Levinson algorithm but overlooks the potential efficacy of the predictor coefficients and prediction error powers iteratively generated. With this stable and efficient implementation, the practicability of the MLE will be greatly enhanced.

5. Results and Discussion

In order to analyze four possible cases and compare their efficiency, the generating algorithm proposed by Lundahl et al. [29] was used to generate DFGN because the realizations produced by this algorithm possess fine correlation structure and long-term dependency.

For more convincing facts, a wider range of Hurst exponents and data sizes were considered, including , 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, and 0.99 (totally, 13 Hurst exponents), as well as , 256, 512, 1024, 2048, and 4096 (totally, 6 types of data sizes). For each data size, 100 realizations of white Gaussian noise were generated by a Gaussian random generator to form 100 realizations of DFGN for each Hurst exponent.

All estimations were performed with the same computing specifications: (1) hardware: a computer of Intel Core i7-2600 processor up to 3.40 GHz and a RAM of 8.00 GB (7.89 GB available); (2) operating system: Windows 7 Professional Service Pack 1; (3) programming software: MATLAB R2011b 64-bit (win64); (4) optimization algorithm: golden section search with threshold 0.0001, which takes 21 iterations and totally 22 function evaluations [42]. Table 1 shows the experimental results, each value representing the mean of mean-squared errors (MSEs) of 100 realizations over 13 Hurst exponents, simply denoted as mean mean-squared errors (MMSEs).

On the other hand, the function evaluation time is recorded and is used to compare with the implementation time of the traditional MLE for efficiency analysis. Table 2 lists the average time (in seconds) of 13 Hurst exponents spent by each approach in two executing procedures, with and without considering the computational efficiency, as well as their corresponding time ratio.

The best results are almost acquired among Case 1, and the second best results are among Case 3, and the worst results are among Case 4 from Table 1. This is reasonable because both mean and variance are known for Case 1, but both mean and variance are unknown for Case 4. It is worth noting that the accuracy of Case 3 is better than Case 2. This indicates that accuracy is more related to known variance than known mean. Generally speaking, a most practical situation is Case 4 with both unknown mean and variance; the situation of Case 1 is less likely practical. Table 1 suggests that using the sample mean instead of the mean estimated by MLE is also a reliable approach.

It is easy to see that the computational cost of the traditional MLE for the Hurst exponent needs , whereas two newly proposed procedures only need . In addition to lower computational complexity, storing data by vector instead of matrix also help raise the computational efficacy. Table 2 suggests that without matrix calculation, time saving is obvious, especially when the data size grows larger. For example, with the data size of 4096, the ratio of each proposed efficient procedure to the traditional one reaches up to 80 times. The ratio will be more tremendous especially for computers of limited resources. These results will contribute to the position of the MLE for estimating Hurst exponent.

6. Conclusions

In parameter estimation, both accuracy and efficiency are generally difficult to coexist. Accordingly, how to weigh the accuracy and efficiency before estimating parameters is usually a matter of a dilemma. The MLE for Hurst exponent is considered optimal in accuracy, whereas the computational cost of the MLE was once considered as tremendous, which hinders the MLE from being recommended to quick response systems. Fortunately, the Levinson algorithm and Cholesky decomposition can be combined to improve the computational efficiency, and further overcome the dilemma. On the other hand, a potential modeling problem of physical signals is also considered. The first proposed MLE for Hurst exponent only considered one case with given mean as zero, which is only suitable for signals of DFBM. However, many physical signals are like the model of DFGN, with nonzero means. Therefore, users must take the sample mean out of the original signal before using the MLE, or a direct computation will easily lead to a severely wrong result, further providing a wrong signal explanation. In order to extend the MLE for Hurst exponent to signals of DFGN, four possible cases are considered: known mean, unknown mean, known variance, and unknown variance. The experimental results show that the computational cost is largely reduced by a combination of Levinson algorithm and Cholesky decomposition. Moreover, numerical stability is also provided to help users avoid numerical mistakes due to negligence. After balancing inherent accuracy with boosted efficiency, the MLE might be the preferred option for estimating Hurst exponent in the near future. More importantly, this idea for efficiently implementing the MLE can be extended to other variants of the MLE for other fields, making real-time computation with best accuracy more possible.

Appendix

Proof of Case 3. In this appendix, maximizing with respect to is proved to be , where , , , and by using mathematical induction. Obviously, denotes the sum of all elements of matrix . For clarity, the subscript is used to emphasize the dependence on the data size during the procedure of proof. Under this situation, , where and .
For , the trivial case, it follows that By maximizing the above quantity, it follows that Therefore, it follows that , which is consistent with the equality , for , as desired. So, the proposition is true for . Next, assume that , for some integers ; that is, Finally, let ; then it follows that By maximizing the quantity with respect to , it follows that Therefore, it follows that as desired. Kay [43] also provides another derivation from a more general form of a linear model without considering the Hurst exponent .