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.
29,303 characters · 8 sections · 20 citation commands
Adaptive Student's t-distribution with method of moments moving estimator for nonstationary time series
Keywords: nonstationary time series, Student's t-distribution, adaptive models, methods od moments, heavy tails
Choosing a parametric family of probability distributions, e.g. Student's t-distribution here, there is usually focus on intuitively static estimation: optimization of a single set of parameters $\theta$ for the entire dataset, usually through maximization of some evaluation like $F=\frac{1}{T}\sum_{t=1}^T f(\theta,x_t)$. For example 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. This way all datapoints have equal $1/T$ contributions, what seems a perfect choice for stationary time series.
In contrast, real life time series are often non-stationary, suggesting to use adaptive estimation adaptive instead - with evolving parameters, like $\theta_t=(\mu_t,\sigma_t,\nu_t)$ in Fig. (ref) for Student's t-distribution we will focus on. Moving estimator for each time $t$ will separately optimize $\theta_t$ parameters based on the previous values $\{x_\tau\}_{\tau<t}$ with weakening weights, to finally optimize:
A natural approach to estimate $\theta_t$ is optimizing analogous function $F_t$: using only the past values $\{x_\tau\}_{\tau<t}$, with exponentially weakening weights to get local behavior:
for $\bar{\eta}\in (0,1)$ learning rate usually above 0.9, also define $\eta=1-\bar{\eta}$ for convenient calculation.
The above ((ref)) moving MLE can be easily directly optimized for $\sigma$ scale parameter of EPD (exponential power distribution) $\rho(x)\sim \exp(-|x|^\kappa)$ adaptive thin tail family containing e.g. Gauss and Laplace distributions, from absolute central moments: $m_p= E[|x-\mu|^p]$, for adaptation evolving with EMA (exponential moving average):
using $p=\kappa$ for EPD, and $\mu_t$ as constant or also adapted using EMA. Here we will take it to Student's t-distribution, this time not through direct MLE due to lack of explicit formula, but through method of moments instead - estimating $\sigma$ scale parameter from absolute central moment for a single power, or $\nu$ degrees of freedom from such moments for two powers.
On example 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 evaluation, here for Student's t-distribution slightly better than for EPD adaptive. Also essentially better than standard methods of $\sigma$ prediction like GARCH(1,1) garch - from one side focused on Gaussian distribution, but also arbitrarily assumed dependencies - here replaced with agnostic philosophy of moving estimator optimizing local parameters.
Such adaptive estimation can be combined with other methods, which might be added in later versions of this article. For example online PCA OPCA or adaptive linear regression regression to combine information from multiple sources like companies here or macroeconomical data - e.g. to improve prediction of the moments, used for parameter estimation here. Finally, as discussed in adaptive, we can use such parametric distributions for normalization $y_t = \textrm{CDF}_t(x_t)$, and then model slight distortion from uniform distribution of $\{y_t\}$ with HCR (hierarchical correlation reconstruction) hcr modelling density as a linear combination, in static or adaptive (evolving in time) way.
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).
The Student's t-distribution was first introduced by Friedrich Helmert in 1875 student, and later in 1908 by William Sealy Gosset signed as "Student" student1, leading to the popular name.
Its basic application is for distribution of sum of $\nu+1$ i.i.d. Gaussian random variables: for the difference between the sample mean and the real mean. For $\nu=1$ it is Cauchy distribution, for large $\nu\to\infty$ it approaches Gaussian distribution.
Its PDF (probability density function), shown in Fig. (ref), is:
for $\mu\in\mathbb{R}$ and $\sigma,\nu\in \mathbb{R}^+$, $\Gamma(z)=\int_0^\infty t^{z-1} e^{-t} dt$ gamma function. Crucially, it has one over polynomial heavy tails $\rho(x)\sim |x|^{-\nu-1}$ for $|x|\to \infty$, hence finite moments $E[x^p]$ only for $p < \nu$.
Its CDF (cumulative distribution function) for $\mu=0, \sigma=1$ is below, for the general case substitute $x \to (x-\mu)/\sigma$:
for $F_{1,2}$ hypergeometric function.
For method of moments we will use absolute central moments: $E[|x-\mu|^p]$ for not necessarily integer power $p\in \mathbb{R}^+$. Using Mathematica there was calculated moment formula as the below integral, finite for $p<\nu$:
Having a $\{x_t \}_{t=1..T}$ data sample, fixing $\nu$ 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 $(0,\nu)$ range, where the possibility to use non-integer $p$ might be crucial for the $p<\nu$ requirement.
Additionally, using various $p$ for such $\sigma$ estimation has various uncertainty depending on $\nu$, as shown in Fig. (ref) - suggesting to optimize $p$ e.g. based on the used $\nu$ range, or even modify $p$ dynamically. For large $\nu$ the optimal $p$ is close to $p=2$ variance estimation, standard for $\nu\to\infty$ Gauss distribution limit. For small $\nu$ the optimal $p$ is $\approx \nu/6$.\\
To estimate $\nu$, 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_{\nu p_1}/M_{\nu p_2}$ is monotonous with $\nu$ (examples in Fig. (ref)), we can e.g. put its behavior into a table and interpolate based on the averages to estimate $\nu$, e.g. done as \verb"find"$\nu$ in the code in Fig. (ref).
However, analogously to $1/(n-1)$ standard adjustment in variance estimator, ((ref)) estimation seems biased - needs adjustment by calculating its expected value, preferably with an explicit formula (yet to be found). In Fig. (ref) such slight adjustment was made by just adding (tuned) 0.9 to found $\nu$.\\
Alternatively, as used for stable distribution (est,stablebook), we can use $\nu>p,-p$ opposite powers, multiplying moments instead: $$E\left[|X|^p\right] E\left[|X|^{-p}\right]=\frac{\Gamma\left(\frac{p+1}{2}\right)\Gamma\left(\frac{\nu-p}{2}\right)} {\sqrt{\pi}\, \Gamma(\nu/2)} \frac{\Gamma\left(\frac{-p+1}{2}\right)\Gamma\left(\frac{\nu+p}{2}\right)} {\sqrt{\pi}\, \Gamma(\nu/2)} $$ However, it lead to inferior log-likelihood for DJIA.
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 (0,2)$ power ($p<\min_t(\nu_t)$), this time with (central absolute) moments evolving in time:
Finally, for $\nu$ degrees of freedom 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.003, \eta_2=0.05, \eta_3=0.005$ for correspondingly $\mu,\sigma,\nu$ (much faster for scale parameter $\sigma$).
Figure (ref) shows evaluation using fixed $\mu=0$ center and various fixed $\nu$ for single MLE $\sigma$ parameter, or $\sigma$ adapted using ((ref)) estimation with $p=1$ power and $\eta_2=0.05$ learning rate - e.g. leading to log-likelihood worse only by $\approx 0.004$ than for optimized evolution of all 3 parameters in Fig. (ref). The $\nu$ estimator needs adjustment - here done by just adding tuned parameter, hopefully to be improved, automatized in future.
The $\nu$ evolution, unavailable in standard ARMA-ARCH approaches, evaluates local tail shapes, probability of potentially destabilizing extreme events - suggesting to call it stability, complementing popular volatility evaluation similar to $\sigma$. Figures (ref), (ref) check that indeed 1967-1983 range suggested in Fig. (ref) has much thinner tails. Figure (ref) shows $\nu$ evolution for all the companies ($\eta_3=0.005$) - such analysis might help to localize and understand stability influencing factors/mechanisms, which hopefully could be applied in future to reduce probability of potentially catastrophic extreme events.
Alternative approach considered e.g. for stable distribution (est,stablebook) is estimation from moments of logarithm of absolute value. We start with transformation of the original variable $X$ to $Y=\ln(|X|)$.
Transforming its moment generating function: $$ E\left[\exp(pY)\right]=\sum_{k=0}^\infty E\left[Y^k\right] \frac{p^k}{k!} =E\left[|X|^p\right]=(M_{\nu p} \sigma)^p $$ which allows to express these moments as
Calculating it for $k=1,2$ using ((ref)) we can get: $$ E[Y^2]-E[Y]^2=E[(Y-E[Y])^2]=\frac{\pi^2 + \textrm{PolyGamma}(1,\nu/2)}{8} $$ Allowing to estimate $\nu$ from moments, also in adaptive way by their EMA update. However, for the discussed data, such moving $\nu$ estimator has led to slightly worse log-likelihood.
As especially the tails e.g. in Fig. (ref) analysis are clearly asymmetric, it might be also valuable to include it. There are two classical approaches to asymmetrize Student-t distribution: noncentral t-distribution noncentral with PDF for $\nu,\delta$ parameters: $$ \frac{e^{-\frac{\delta ^2}{2}} 2^{\nu } \nu ^{\frac{\nu }{2}+1} \Gamma \left(\frac{\nu +1}{2}\right) }{\pi }H_{-\nu -1}\left(\frac{-x \delta }{\sqrt{2} \sqrt{x^2+\nu }}\right) \left(\nu +x^2\right)^{-\frac{\nu +1}{2}}$$ and skewed generalized t-distribution skewed with PDF: $$ \frac{p}{2 v \sigma q^{\frac{1}{p}} B(\frac{1}{p},q) \left(1 + \frac{| x-\mu + m |^p}{q (v \sigma)^p (1 + \lambda \text{sgn}(x-\mu + m))^p}\right)^{\frac{1}{p}+q}} $$
However, they are much more complicated, especially if wanting to search for adaptive estimation, what might be worth to consider in the future.
Instead, for simplicity let us consider asymmetry by just gluing in $\mu$ two standard Student t densities of separate parameters: $\overleftarrow{\sigma},\overleftarrow{\nu}$ describing the part on the left $x<\mu$, and $\overrightarrow{\sigma},\overrightarrow{\nu}$ describing the part on the right $x>\mu$. While it is not necessary, wanting continuous PDF we need to enforce that both parts have the same value in $\mu$, what through linear scaling leads to PDF:
It has advantage that we can use adaptive estimation of $\sigma,\nu$ exactly as previously, just separate for values below/above current $\mu$: we maintain two copies of $\sigma,\nu$, and update one of them based on $\textrm{sgn}(x-\mu)$.
However, tests on DJIA has lead to $\approx 0.01$ worse log-likelihood for separate left/right adaptation of both $\overleftarrow{\nu},\overrightarrow{\nu}$ and $\overleftarrow{\sigma},\overrightarrow{\sigma}$. It has returned to the original log-likelihood for separate $\overleftarrow{\nu},\overrightarrow{\nu}$ adaptation and common $\sigma=\overleftarrow{\sigma}=\overrightarrow{\sigma}$. So while $\overleftarrow{\nu},\overrightarrow{\nu}$ can be used to describe shape of separate left/right tail, for $\sigma$ scale parameter it seems better to use a common one. Also we could search for a more sophisticated adaptive estimation, or maybe use the noncentral or skewed t-distribution.
We can use such evolving especially $\nu$ to evaluate the market: e.g. to estimate probability of extreme events separately toward left and right, also use them as local parameters for various models. Figure. (ref) shows its evolution for 107 years of DJIA.
Later improvement of log-likelihood was reached by estimating $\sigma,\nu$ together, but then using deformed asymmetric Student's t-distribution ((ref)) of tuned parameters, like in Fig. (ref). As usually price drop left tail is heavier: has lower $\nu$, and $\overrightarrow{\nu}-\overleftarrow{\nu} \sim 2$ in Fig. (ref) here. Fixing this difference like in Fig. (ref) has finally allowed to improve log-likelihood from 3.3389 to 3.3406, by using common $\sigma=\overleftarrow{\sigma}=\overrightarrow{\sigma}$ estimated using $\nu$ online updated as previously, but with asymmetric $\overleftarrow{\nu}=\nu+\Delta_\nu,\overrightarrow{\nu}=\nu-\Delta_\nu$ for optimized $\Delta_\nu=0.8$. Its further tuning, also slightly deforming $\sigma$ (suggested by left/right probability difference due to rescaling for ((ref)) continuity), has allowed to increase log-likelihood to 3.3413 using source shown in Fig. (ref), its density and parameter evolution is presented in Fig. (ref).
This article introduces looking novel extensions of method of moments - both to absolute central moments with not necessarily natural powers (crucial to work with low $\nu$), but more importantly as EMA moving estimators - for parameters evolving in time, also asymmetrically for left/right tails. Beside better log-likelihood evaluation, it provides evolution of these crucial parameters like in Fig. (ref), (ref), (ref) - including $\nu$ degrees of freedom evaluating probability of extreme events, which dependency understanding might allow to introduce some market stabilizing mechanisms. For example it suggests search for mechanisms of drastic increase of $\nu$ especially in 1967-1983 period for DJIA, confirmed in Fig. (ref), (ref).
This is a general approach which might be worth taking also to other distributions like alpha-stable, and larger models. Also it is worth combining with other especially adaptive models, like online linear regression and HCR (hierarchical correlation reconstruction) - what is planned to be done in further versions of this article.
Examples of plans for further work: