Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
23,984 characters · 8 sections · 22 citation commands
Adaptive stable distribution and Hurst exponent by method of moments moving estimator for nonstationary time series
Keywords: nonstationary time series, stable distribution, Hurst exponent, adaptive models, methods od moments, heavy tails
Working with parametric probability distributions, e.g. alpha-stable (stable,stablebook) here, there is a difficult question of choosing its parameters - often done with static estimation: optimization of a single set of parameters, usually through maximization of some function: $F=\frac{1}{T}\sum_{t=1}^T f(\theta,x_t)$, especially log-likelihood in popular MLE (maximal likelihood estimation) using $f(\theta,x)=\ln(\rho_{\theta}(x))$, where $\rho_{\theta}(x)$ is PDF (probability distribution function) for the assumed parametric family. It uses all datapoints with the same $1/T$ contributions, what is appropriate for stationary time series.
However, real-life e.g. financial data is usually non-stationary, requiring some adaptive estimation of evolving time-dependent parameters e.g. $\theta_t=(\mu_t,\sigma_t,\alpha_t)$, visualized in Fig. (ref), (ref). In classical ARMA-ARCH garch approaches there is assumed some arbitrary dependence between parameters for succeeding times. To avoid their bias, there was proposed moving estimator approach adaptive where we shift estimator to use local past values $\{x_\tau\}_{\tau<t}$ with exponentially weakening weights for the old values:
for $\bar{\eta}\equiv 1-\eta$ learning rate, usually above 0.9. For the final evaluation, also to optimize $\eta$, there is used their average:
Such separate estimation for each position is computationally costly, e.g. we can minimize ((ref)) with the previous parameters as the starting point for optimization. It can be simplified especially for estimators being averages, which can be modified to moving estimator by just switching to exponential moving average (EMA). For example EPD (exponential power distribution) family with $\rho(x)\sim \exp(-|x|^\kappa)$ type density has $\sigma^\kappa = E[|x-\mu|^\kappa]$ MLE (maximum likelihood estimator) - changing average to EMA adaptive usually gives better log-likelihood than ARMA-ARCH, also here summarized in Fig. (ref).
However, usually MLE estimation is much more complex than just averaging, but moments can be easily extended to adaptive with moving averages, and allow to estimate parameters. Working with heavy tails, we need to use moments for low powers as higher become infinite, suggesting to use absolute central moments: $m_p= E[|x-\mu|^p]$ for one or a few powers, which for EMA becomes:
and previously was discussed mystudent for moving estimator of Student's t-distribution for all 3 parameters: $\mu,\sigma,\nu$.
Here this approach was adapted for another basic heavy-tailed distribution family: alpha-stable stable,stablebook, obtained e.g. by averaging i.i.d. random variables of not necessarily finite variance in the Generalized Central Limit Theorem genCTL, hence should be often universal for real-life date.
Stable distribution also leads to $H_q=q/\alpha$ generalized Hurst exponent multifractal for i.i.d time series, describing $\textrm{E}\left(|x_{t+\tau}-x_t|^q\right) \sim \tau^{q H_q}$ behavior, hence the proposed moving $\alpha$ estimator can be also used for adaptive estimation of Hurst exponent.
For dataset of 107 years Dow Jones Industrial Average (DJIA) daily log-returns and 10 years for 29 its recent companies, there was tested such adaptive estimation especially of $\sigma$, leading to essentially better log-likelihood than GARCH(1,1) garch, also than EPD adaptive.
However, obtained log-likelihood has turned out slightly worse than using adaptive Student's t-distribution mystudent, probably because it covers $\rho(x)\sim |x|^{-\nu-1}$ tails for arbitrarily large power, while stable distribution switches to Gaussian from 3rd power as variance becomes finite.
There was used 1900-2007 daily Dow Jones index\footnote{Source of DJIA time series: http://www.idvbook.com/teaching-aid/data-sets/the-dow-jones-industrial-average-data-set/}, working on $x_t=\ln(v_{t+1}/ v_t)$ sequence of daily log-returns.
Figure (ref) additionally contains such evaluation of log-returns for 29 out of 30 companies used for this index in September 2018. Daily prices for the last 10 years were downloaded from NASDAQ webpage (www.nasdaq.com) for all but DowDuPont (DWDP) - there were used daily close values for 2008-08-14 to 2018-08-14 period ($2518$ values) for the remaining 29 companies: 3M (MMM), American Express (AXP), Apple (AAPL), Boeing (BA), Caterpillar (CAT), Chevron (CVX), Cisco Systems (CSCO), Coca-Cola (KO), ExxonMobil (XOM), Goldman Sachs (GS), The Home Depot (HD), IBM (IBM), Intel (INTC), Johnson&Johnson (JNJ), JPMorgan Chase (JPM), McDonald's (MCD), Merck&Company (MRK), Microsoft (MSFT), Nike (NKE), Pfizer (PFE), Procter&Gampble (PG), Travelers (TRV), UnitedHealth Group (UNH), United Technologies (UTX), Verizon (VZ), Visa (V), Walmart (WMT), Walgreens Boots Alliance (WBA) and Walt Disney (DIS).
A distribution is referred as stable if a linear combination of two independent random variables from this distributions has the same distribution, up location $\mu$ and scale $\sigma$ parameters. They are also referred as L\'{e}vy-stable alpha distributions as probably were first studied by Paul L\'{e}vy stable.
Beside location $\mu$ and scale $\sigma$, they also have stability parameter $0<\alpha <2$ for $\rho(x)\sim 1/|x|^{\alpha+1}$ probability density of tails (heavy), or $\alpha=2$ being Gaussian distribution with $\rho(x)\sim \exp(-x^2)$ tails.
While generally there can be also discussed further asymmetry/skewness parameter $\beta\in [-1,1]$. For simplicity let us start without it: assume $\beta=0$, what leads to a simple characteristic function, with some PDF/CDF shown in Fig. (ref):
Characteristic function of sum of independent random variables is product of their characteristic functions: $\varphi_{X+Y}(t)=\varphi_{X}(t) \varphi_Y (t)$, making ((ref)) indeed a stable variable: for the same $\alpha$, location and scale parameter:
It suggests their universality for real-life data, especially through the Generalized Central Limit Theorem genCTL. It says that averaging sequence of independent, identically distributed variables leads to a stable distribution, generalizing the Cental Limit Theorem to infinite variance case.
For method of moments we will use absolute central moments: $E[|x-\mu|^p]$ for not necessarily integer power $p\in \mathbb{R}$. Using Wolfram Mathematica there was calculated moment formula as the below integral (for $\beta=0$), finite for $p<\alpha$:
Having a $\{x_t \}_{t=1..T}$ data sample, fixing $\alpha$ and using some $\mu$ estimator e.g. approximate $\hat{\mu}=T^{-1} \sum_t x_t$ as just mean, the above formula gives simple estimator of scale parameter $\sigma$:
The used $p$ has to be in $(-1,\alpha)$ range, where the possibility to use non-integer $p$ might be crucial for the $p<\alpha$ requirement.
Additionally, using various $p$ for such $\sigma$ estimation has different uncertainty depending on $\alpha$, as shown in Fig. (ref) - suggesting to optimize $p$ e.g. based on the used $\alpha$ range, or even modify $p$ dynamically. For $\alpha=2$ the optimal $p$ is $2$ variance estimation. For $\alpha\in (1,2)$ the optimal $p$ is $\approx (\alpha-1)/2$.\\
To estimate $\alpha$, a natural direct way is to divide such averages for two different powers $p_1,p_2$, removing $\sigma$ dependence:
Choosing some $p_1\neq p_2$, the $M_{\alpha p_1}/M_{\alpha p_2}$ is monotonous with $\alpha$ (examples in Fig. (ref)), we can e.g. put its behavior into a table and interpolate based on the averages to estimate $\alpha$, e.g. done as \verb"find"$\alpha$ in the code in Fig. (ref).
The finally found parameters might be worth tuning by some small modification, e.g. in Fig. (ref) there was used $\sigma m=1.05$ to multiply the final $\sigma$, allowing to slightly improve log-likelihood.
Above methods of moments can be easily adapted for moving estimator by just replacing averages with exponential moving averages - uniform weights with exponentially weakening.
For the center $\mu$ we can use just a basic adaptation below - it is optimal only for the Gaussian case ($\nu\to \infty$), hence generally it could be slightly improved. However, for the discussed data the gains were already nearly negligible.
The most crucial is $\sigma$ scale parameter adaptive estimation, as e.g. in ARCH family but in more agnostic way, here using ((ref)) formula for a chosen $p\in (-1,2)$ power ($p<\min_t(\alpha_t)$), this time with (central absolute) moments evolving in time:
Finally, for $\alpha$ estimation we can use ((ref)) formula for analogously updated moments for some 2 different powers $p_1, p_2$ and some $\eta_3$ learning rate.
Figure (ref) contains used Mathematica code for adaptation of all 3 parameters, with their evolution for DJIA time series. Manual tuning has lead to 3 different learning rates there: $\eta_1=0.002, \eta_2=0.03, \eta_3=0.006$ for correspondingly $\mu,\sigma,\alpha$ (much faster evolution of volatility/scale parameter $\sigma$).
Figure (ref) shows evaluation using fixed $\mu=0$ center and various fixed $\alpha$ for single MLE $\sigma$ parameter, or $\sigma$ adapted using ((ref)) estimation with $p=1$ power and $\eta_2=0.05$ learning rate, just slightly worse that for Student's t-distribution from mystudent.
The $\alpha$ evolution, unavailable in standard ARMA-ARCH approaches, evaluates local tail shapes, probability of potentially destabilizing extreme events - allowing to indeed interpret it as stability, complementing popular volatility evaluation similar to $\sigma$. For example it shows that in 1967-1983 there was nearly Gaussian distribution, also confirmed in mystudent. We can observe events like Black Monday leading to sudden drops of market stability. Having evolving evaluation of market stability brings a chance to understand various factors influencing it, and maybe try to control them.
The standard approach to include asymmetry for stable distribution is through $\beta\in[-1,1]$ skewness parameter - modifying the characteristic function for $\alpha \in (1,2]$ to:
It might be worth to include its estimation, what could be done e.g. by extending to EMA updated moments maintaining sign: $E[\textrm{sgn}(x)|x|^p]$. However, just fixing $\beta$ can already bring slight improvement of log-likelihood, like $\beta=-0.3$ in Fig. (ref).\\
Alternative general approach, e.g. mentioned for Student's t-distribution mystudent, and suggested by tail analysis in Fig. (ref), is just gluing two distributions of slightly different parameters - to be used for $(x_t-\mu_t)/\sigma_t$ online normalized variables.
For example, as $\alpha$ in Fig. (ref) turns out $\Delta\approx 0.05$ lower for $x<\mu$ negative tail, we can estimate $\alpha$ as for $x>\mu$ positive tail, and use $\alpha-\Delta$ for the negative tail. Working on $(x_t-\mu_t)/\sigma_t$ normalized values, we need probability density for normalized $\mu=0,\sigma=1$, or approximated and potentially different for left/right tails: $\sigma_l,\sigma_r \approx 1$. Finally, including continuity for $x=\mu$, we get such asymmetric PDF as:
for $\rho_0(\alpha)=\rho_{0,1,\alpha}(0)=\Gamma(1+1/\alpha)/\pi$ density in $\mu$. Unfortunately for DJIA data it did not allowed to improve log-likelihood.
There is popular estimation of multifractal Hurst exponent, for example as $\textrm{E}\left(|x_{t+\tau}-x_t|^q\right) \sim \tau^{q H_q}$ behavior for $H_q$ generalized Hurst exponent.
While its main motivation is evaluation of long-range behavior, for heavy tails already for i.i.d. variables it obtains e.g. $H_q=q/\alpha$ for $\alpha$ parameter stable distribution - allowing to estimate it in adaptive way through moving estimator of $\alpha$.
Moreover, the proposed adaptive estimation allows to include in considerations, or even remove the heavy tail dependence - allowing for compensated Hurst exponent indeed focused only on evaluation of long range behavior. For this purpose, we can for example normalize the variables by transforming with estimated CDF of stable or Student's t-distribution, and then inverted CDF of target distribution e.g. $N(0,1)$ normalized Gaussian: $$ x_t \to \textrm{CDF}^{-1}_{N(0,1)}(\textrm{CDF}_{\theta_t}(x_t))$$ then calculating Hurst exponent for such normalized variables.
This article extends the method of moments to adaptive case for EMA moving estimation of $(\mu_t,\sigma_t,\alpha_t)$ parameters of alpha-stable distribution analogously to previous articles for exponential power distribution adaptive and Student's t-distribution mystudent.
Figure (ref) suggests stable distribution is appropriate to model tails of financial data, however, the obtained log-likelihood was slightly worse than for Student's t-distribution, probably due to covering 1/polynomial tails only up to degree 3 or their differences visualized in Fig. (ref), suggesting to consider some hybrid between Student's t and stable distribution, what is planned for further work, e.g. just switching to Student for when evolving $\alpha$ exceeds 2, or maybe automatically optimizing shape of parametrized distribution based on the data.
Adaptive estimation of $\alpha$ distribution allows to continuously evaluate market stability, e.g. to understand and try to control it, react to crucial events. Also can be used as additional variable to improve predictions of various models.
Such found evolving distribution can be used as initial stage e.g. by normalization to nearly uniform distribution in $[0,1]$ by $x_t\to \textrm{CDF}_{\theta_t}(x_t)$ for further analysis e.g. of joint distributions, or long-range dependencies e.g. with Hurst exponent.
Examples of plans for further work: