A wide variety of phenomena whose states change over time at random exhibit a number of time periods in which the states remain unchanged. These range across stock market data, diffusion of large molecules within cells, movement of bacteria, electricity markets, and more; see e.g. Chapter 1 of . For example, in finance, although the celebrated Black–Scholes–Merton model has been widely used to describe fluctuations of the prices of financial products, the model uses a geometric Brownian motion, which fails to capture the constant-price periods that many low-volume equities experience. This motivates the construction of time-changed stochastic processes exhibiting such constant periods. As a simple example, consider a Brownian motion composed with the inverse of an independent -stable subordinator (or inverse -stable subordinator for short). The resulting stochastic process , called a time-changed Brownian motion
, shows constant periods and has variance growing at the rate of. Since the stability index takes a value in the interval and the variance of the Brownian motion grows at the rate of , particles represented by the time-changed Brownian motion spread at a slower rate than the regular Brownian particles for . Because of this, the particular time-changed Brownian motion and its variants are often referred to as subdiffusion processes. Figure 2 gives a graphical comparison of sample paths of the inverse -stable subordinator and the time-changed Brownian motion with different values of . As the simulations show, we expect to observe longer constant periods for a smaller value of . The time-changed Brownian motion is not a Markov process and its densities satisfy the time-fractional heat equation where denotes the Caputo fractional derivative of order (see [22, 23]). For more details about properties of and its various extensions (including stochastic differential equations) as well as their continuous-time random walk counterparts, see e.g. [24, 29] and references therein.
Various methods for estimating the stability index of a stable distribution have been presented in the literature; see e.g. [1, 2, 5, 8, 9, 10, 11, 18, 21, 25]. For example, Hill’s estimator in  is popular and well-known for its robustness; it only assumes the underlying distribution has power law tails. In , the authors proposed a shift-invariant version of Hill’s estimator to increase the robustness. On the other hand, in , they revealed a connection between Hill’s estimator and a least-squares estimator obtained from a - transformation of the power law tails. The maximum likelihood estimator (MLE) was established in , overcoming the difficulty that the stable density does not have a closed form. The MLE in general requires a particular distributional form, so it is not as robust as Hill’s estimator; however, it uses all of the data, rather than only the largest order statistics based on which Hill’s estimator and its variants are constructed. In , another estimator using all of the data was constructed via the analysis of the -wright function together with the method of moments.
, the authors compared six known parameter estimation methods in the context of subdiffusions, where the lengths of the constant periods follow a one-sided (totally right-skewed)-stable distribution. This is the situation we also consider in this paper. However, we propose a new estimator for the stability index based on the recent development of numerical approximations of subdiffusion processes in [13, 14, 16, 17]. Our idea is to regard real data exhibiting constant periods as a numerically approximated path of some time-changed process and to use the method of moments to estimate . Note, however, that the method of moments does not apply to the stable distribution itself in the usual sense since it does not have the first moment. We instead apply the method to the “number” (rather than the “lengths”) of observed constant periods. Cahoy’s estimator in  is based on the logarithmic moments of the stable distribution and has a connection to the estimator to be proposed in this paper; see Remark 3.
Our estimator is especially suited for use on data sets with many observations over a long period of time (i.e. when both the number of observed paths and the time horizon are large) since it is asymptotically unbiased and consistent; see Remark 5(a). Moreover, it is robust in the sense that it only requires the underlying subordinator to have power law probability tails; see Remark 6. One major advantage of using our estimator as opposed to some other well-known estimators is that it allows one to estimate the entire range of values with reasonable accuracy and precision; see Remark 7(a).
2 Stable subordinators and their inverses
Throughout the paper, all stochastic processes are assumed to be defined on a probability space and denotes the expectation under . A stochastic process having independent and stationary increments is called a Lévy process. For the purpose of this paper, let us consider a one-dimensional Lévy process with increasing, càdlàg paths (right continuous paths with left limits) starting at , which is usually called a subordinator. The distribution of a subordinator is characterized by its Laplace transform
where , called the Laplace exponent of , is a Bernstein function on with . A -stable subordinator is a subordinator with for , where is a parameter called the (stability) index.
The probability tails of the -stable subordinator are very different from those of a Brownian motion . Namely, for , while shows the exponential decay
exhibits the slower, polynomial decay
where is the Gamma function (see p.37 and p.76 of ). Due to the presence of the heavy tail, does not have the th moment for , so the method of moments is not applicable in the usual sense. The stable subordinator has strictly increasing paths that go to as with infinitely many jumps, where the jump times are dense in (see Theorem 21.3 of ).
Define the inverse or first hitting time process of a -stable subordinator by
for . The process is called an inverse -stable subordinator. Since has strictly increasing paths starting at 0, has continuous, nondecreasing paths starting at 0. Figure 2 illustrates the inverse relationship between the sample paths of and . The jumps of correspond to the constant periods of , and the inverse relation implies
, the random variablepossesses finite exponential moments, i.e. for any (see e.g. [14, 15]). In particular, the following formula holds for the th moment of for any (see e.g. Proposition 5.6 of  or Example 3.1 of ) and plays a key role in deriving an estimator for the index :
3 Derivation of an Estimator and its Properties
Fix an equidistant step size and a time horizon . Simulate a sample path of the stable subordinator of index , which has independent and stationary increments, by setting and then following the rule for where ’s are i.i.d. random variables having the same distribution as (which can be generated via an algorithm presented in ). Stop this procedure upon finding the integer satisfying For each path of , the number exists as a finite number since as . Note that is an -valued random variable depending on the initially chosen constants and .
for . The sample paths of are nondecreasing step functions with constant jump size and the length of the th constant time interval given by . Indeed, whenever . In particular, each path of has a total of constant periods, and a.s.,
Note that even though and do not coincide, they are asymptotically equivalent as ; i.e. for fixed and ,
Based on this observation, we define a method-of-moments-like estimator for the parameter implicitly as a solution to the equation
where is the sample mean of a random sample of size from the distribution of the random variable . This is not the method of moments in the usual sense as we do not have the equality for a fixed time horizon ; however, the relation (7) makes the definition of our estimator reasonable when is sufficiently large. In the remainder of this section, we discuss how large should be as well as the properties of the estimator . First, as the following proposition shows, satisfying (8) uniquely exists as long as and , where is the Euler–Mascheroni constant.
Let . If , then the function in (6) is a smooth, strictly increasing bijection from to . Moreover, if , then is convex on .
The smoothness of the function follows by the smoothness of the Gamma function. Suppose and observe that
where is the digamma function defined by for . The digamma function and its derivative have the series representations and (see 6.3.16 and 6.4.10 in ). Since is strictly increasing, for any . Combining this with (9) yields for all . Since and , it follows that is a strictly increasing bijection from to .
Now, suppose . To prove the convexity of , observe that
Thus, the convexity follows upon verifying that for all . Since is increasing and is decreasing,
The latter is positive since , which completes the proof. ∎
In the remainder of the section, assume that . We have so far defined the method-of-moments-like estimator only conditionally on the event that . To define on the entire probability space, we formally set if and if . As a result, can be expressed as
This allows to take 0 and 1 even though the true parameter cannot. However, in practical situations where is very large, this will not be a serious issue since, as the following theorem shows, the probability can be made as small as we want by taking large enough.
Let . Let be the sample mean of a random sample of size from the distribution of satisfying (3). Then as ,
Moreover, there exists a constant not depending on , or such that for any ,
The construction of our estimator for the parameter is completely different from those of the existing estimators discussed in Section 1, except for the one proposed in . Namely, we analyze the distribution of , which represents the “number” of the observed constant periods minus 1, rather than the distribution of the “lengths” of the constant periods. On the other hand, in , the author derived an estimator using the logarithmic moment of the stable random variable together with the equality in distribution of and . (The latter equality is expressed as in that paper.) The estimator for is given for each fixed by , where denotes the sample mean corresponding to the theoretical mean . We can observe a connection of Cahoy’s estimator at with our method-of-moments-like estimator. Indeed, is calculated from observed values of the random variable , which can be approximated by with small; thus, is connected with the random variable arising in the method-of-moments-like estimation.
The following theorem concerns asymptotic properties of the method-of-moments-like estimator for large sample size (i.e. when many observations of the same subdiffusion process are available as in the real-life example given in Section 5). In particular, part (c) allows us to determine how large should be in order to achieve a target level for the asymptotic variance of , which is important in applications. Recall that the Gamma function is convex on and attains the minimum at .
If , then . If , then since is increasing. In both cases, a.s.
(b) By (a), it suffices to show that . Note that as due to (5). Note also that
Take large enough so that and . Then by (5), it follows that
Clearly, is a strictly increasing bijection from to with inverse , and for all . Thus,
from which the inequality follows, as desired.
are i.i.d. with finite second moment, so by the central limit theorem, aswhere is the theoretical variance of . We wish to apply the delta method to , but since does not exist and any integer-order derivative of vanishes on , the delta method fails if , even using higher order Taylor approximations.
To apply the delta method without any technical issues, we smooth out the function at as follows:
where with small is a smooth, strictly increasing, convex function such that with and denoting the derivatives from the left and right, respectively. Then is a smooth, strictly increasing function on , and hence, exists and is positive regardless of the value of . Therefore, for the modified estimator defined by by the delta method, as , and in particular,
Moreover, since is convex on and concave on when due to Proposition 1 (as the convexity of implies the concavity of ), it follows that
(a) (Asymptotic unbiasedness and consistency) Theorem 4(b) shows that , when regarded as an estimator indexed by both and , is asymptotically unbiased and consistent as the indices go to infinity.
(Robustness) Our estimation method can be applied to a more general subordinator whose Laplace exponent in (1) has the asymptotic behavior
(The special case when for all recovers a -stable subordinator.) Indeed, in that case, the Laplace transform of the function satisfies as (see e.g. Proposition 3.1 in ), and hence, by the Tauberian theorem (see e.g. Section 1.7 of ), as . In other words, equality (2) with approximately holds for large enough . Hence, when is large, it is reasonable to use the estimator defined in (10) to estimate the value of .
Note that the asymptotic condition (20) means the subordinator has power law probability tails. Indeed, by the proof of Lemma 3.4 in , for each fixed , which is asymptotically equivalent to as if and only if condition (20) holds, but by the Tauberian theorem, the latter is equivalent to the statement that as .
4 Numerical Comparison to Existing Methods
In this section, we use simulations in Matlab to compare the performance of the method-of-moments-like estimator defined in (10) (MOM-like estimator for short) to Cahoy’s estimator in , Hill’s estimator in , and the Meerschaert–Scheffler estimator in  (MS estimator for short). We chose the latter three estimators for comparison purposes since i) they are simple to apply, ii) Cahoy’s estimator is also constructed via the method of moments, iii) Hill’s estimator has been widely employed in practice, and iv) both Cahoy’s and MS estimators rely on all the data like the MOM-like estimator. We calculate Hill’s estimator based on the largest 10% of the data as that is common in practice.
Note that this section focuses on data that are realizations of the discretized inverse -stable subordinator but not on data coming from a time-changed process of the form . This is because under some assumptions on the outer process , a procedure for estimating for data following is indeed the same as the procedure for estimating for data following . We postpone a detailed discussion of this matter to Section 5, where we treat real data observable by means of a time-changed process.
We first generate paths of the discretized time change on a fixed time interval , calculate the sample mean for the observed paths, and obtain the MOM-like estimate via equation (10), where is determined as the number of constant periods minus 1 for the th path. On the other hand, as mentioned in Remark 3, Cahoy’s estimator requires realizations of the time change ; however, only the paths of the discretized time change are available here. Therefore, instead of appearing in Remark 3, we use to obtain an approximate version of Cahoy’s estimate. Using the same data set but after aggregating all the constant periods observed in the paths, we calculate Hill’s and MS estimates based on the “lengths” (rather than the “number”) of the constant periods via the formulas provided in .
Note that our simulation relies on the choice of three hyper-parameters — the step size , the final time , and the number of observed paths. So we simulate using various combinations of the hyper-parameters over different values of . The results in three particular cases with fixed to be 1 are summarized in the tables in Figure 3, which indicates that the MOM-like estimator i) improves its performance as and increase regardless of the value of , ii) is more accurate than Hill’s and MS estimators for large values of , and iii) performs better than the approximate version of Cahoy’s estimator for small values of .
Next, with still fixed, we take and , which provide the setting for a real data to be treated in Section 5, and repeat the above estimation procedure 100 times to produce 100 estimates based on each of the four different methods. We then calculate the sample mean and sample variance for each method. The results summarized in Figure 4 show that the MOM-like estimator gives a reasonably accurate mean with a very low variance for the entire range of values.
We now turn our attention to the effect of the value of on the MOM-like estimation. It is natural to expect that the accuracy of the MOM-like estimation increases as decreases, which is also indicated by the upper bound for the pathwise estimation error in (19). Here, we again record the sample means of estimates based on 100 repetitions of the MOM-like estimation with and fixed but this time with different values of . Figure 5 gives the simulation results with chosen from the four values , where the corresponding sample variances are all within 0.00025. It shows that when the true value is small, the MOM-like estimator performs better with a smaller value of ; however, for large , choosing a smaller value of from the particular four values contributes little to improving the estimation results. Unfortunately, we do not have a straightforward criterion for choosing the value of ; it depends on the accuracy level that one wants to achieve, and generally speaking, a small value is recommended when the data exhibits long constant periods (which implies the true value is small). On the other hand, we suggest taking with for any possible so that the quantity appearing in the fundamental estimates in (5) is guaranteed positive and hence provides a meaningful lower bound for . For the latter purpose, taking suffices since we always assume . In Section 5, where we deal with real data with and , we take since Figures 4 and 5 guarantee satisfactory performance of the MOM-like estimator for these values of and with .
(a) The above simulations may seem to suggest that, for realizations of ,
the MOM-like estimator outperforms the approximate version of Cahoy’s estimator when the true is small. However, Cahoy’s estimator has the following advantages: i) it does not involve a root finding algorithm and ii) it is valid for any fixed and comes with explicit formulas for confidence intervals (whereas the MOM-like estimator requires equation (
and comes with explicit formulas for confidence intervals (whereas the MOM-like estimator requires equation (8) to be solved with a certain algorithm and can only provide approximate confidence intervals for large via the normal approximation depending on small in the proof of Theorem 4(c)). On the other hand, when a decent computing environment is available and is very large as is often the case with real data, the MOM-like estimator may become a suitable option as it allows one to estimate the entire range of values with reasonable accuracy and precision.
(b) In  , the authors
introduced a modified version of the cumulative distribution function (CDF) for the lengths of constant periods. The modification accounts for the fact that the beginning and ending of a given constant period in real data may have actually occurred at time points when the data was not recorded. The latter is an important issue when parameter estimation is carried out based on the “lengths” of constant periods, while it becomes less of an issue with the MOM-like or Cahoy’s estimation since the “number” of constant periods
is not affected by the exact timing of the beginning and ending of each constant period. Being able to avoid a discussion of this subtle issue as well as the simple formula for finding the estimate
is an advantage of using the MOM-like or Cahoy’s estimation.
In the above discussion, we did not compare the modified CDF method to the other four estimation methods since (i) it requires much more computing power and (ii) it considers the setting for real data observed at pre-specified discrete time points (while we used simulated paths of
, the authors introduced a modified version of the cumulative distribution function (CDF) for the lengths of constant periods. The modification accounts for the fact that the beginning and ending of a given constant period in real data may have actually occurred at time points when the data was not recorded. The latter is an important issue when parameter estimation is carried out based on the “lengths” of constant periods, while it becomes less of an issue with the MOM-like or Cahoy’s estimation since the “number” of constant periods is not affected by the exact timing of the beginning and ending of each constant period. Being able to avoid a discussion of this subtle issue as well as the simple formula for finding the estimate is an advantage of using the MOM-like or Cahoy’s estimation. In the above discussion, we did not compare the modified CDF method to the other four estimation methods since (i) it requires much more computing power and (ii) it considers the setting for real data observed at pre-specified discrete time points (while we used simulated paths ofin which the beginning and ending of each constant period may occur at any time point in ).
5 Real-Life Application: Low-Volume Stock Modeling
This section illustrates how to apply the MOM-like estimator to real data collected from a stock market. Traditionally, low-volume stocks have been difficult to model since they often feature periods in which no trades are made. As a result, the price often has long constant periods throughout the trading day, and therefore, it is a good candidate to model using a time-changed process . However, since an observed data is always discrete, we regard it as a realization of the discretized process . In terms of the discretized time change , we assume that the underlying subordinator has power law probability tails with index even though the use of an exponentially tempered power law might be more appropriate, as pointed out e.g. in [20, 26]. (Note that our purpose here is simply to illustrate how to apply our estimation method to given data exhibiting constant periods.) On the other hand, is assumed to be a geometric Brownian motion with representation , where is a Brownian motion independent of the subordinator and the constants and are additional parameters to be estimated. We focus on the logarithmic stock price , where , and note that the independence assumption implies that is independent of .
Here, we use second-by-second trading data of the low-volume NASDAQ stock with ticker WSTL for the nine weeks starting from September 30, 2019 (consisting of trading days), obtained from Bloomberg Terminal. Figure 7 provides the observed path of the logarithmic price on a day, which clearly exhibits some constant periods. We assume that each of the observed paths is a realization of the discretized process with on the time interval , with being the number of seconds in a trading day. Recall from Figures 4 and 5 that the MOM-like estimator performs well for these hyper-parameters if the data is observed by means of .
By construction of , each path of takes values and has constant periods. We claim that the number of constant periods observed in a path of equals a.s. To see this, recall that is independent of and note that, given the information of the path of , if and only if for at least one . (For example, if and the remaining ’s are nonzero and distinct from one another, then .) Moreover,
and the random vectorgiven is non-degenerate multivariate Gaussian and hence has a continuous distribution, so
which yields a.s. Since is impossible, a.s. In other words, a sample path of is a step function with constant periods that are completely ascribed to the constant periods of the corresponding path of .
With this observation in mind, we estimate the three parameters , and as follows. If an observed path of constantly changes in value (i.e. the values at any two successive time points are distinct), then we consider the number of constant periods to be (i.e. ). In the other extreme case when an observed path stays constant over the entire time period, we set the number of constant periods to be (i.e. ). Due to our observation in the previous paragraph, the obtained number minus 1 for each of the days is considered a realization of the random variable , and consequently, we obtain a total of realizations of . The extremely large value makes the observed value of fall in the interval , and equation (8) gives as the MOM-like estimate for . (For comparison, Cahoy’s, Hill’s and MS estimates are 0.2950, 0.7501 and 0.6605, respectively. The discrepancy between the estimated values indicates the power law distribution may not be an appropriate model here.)
Next, following the idea presented in e.g. , we remove all the constant periods from each of the 44 observed paths and record all the jump sizes over the 44 days. Since each jump size is given in the form with , we regard the recorded jump sizes as a random sample drawn from , where . The obtained sample gives the sample mean
and sample standard deviation. Figure 7 provides a simulated path of the time-changed process