GME (along with others such as AMZ) have made huge moves this past week. These were really the outliers or Black Swans Taleb warned us about… right? Not quite: all of these returns are in line with their pre-“shock” power-law distributions, which means that there is a method which does tell you that this is not a one-in-a-gazillion-years event.

We will need the tidyquant, tidyverse, plotly, and evir packages.

Fetch the data with the tidyquant functions and get daily returns. Also define a breakpoint date, which we will use to analyze how statistics change before this period and after.

GME = tq_get('GME', from = '2002-02-13', to = '2021-01-30') %>%
  tq_transmute(select = adjusted,
               mutate_fun = dailyReturn,
               col_rename = 'returns')

bpd = '2021-01-22'

First, a notional view of the returns since GME’s IPO on February 13th, 2002. We will denote the daily arithmetic return as \(r_t\), with \(t=1,\dots,T\).

plot_ly(GME, x=~date, y=~returns, type = 'scatter', mode = 'lines') %>%
  layout(
    yaxis = list(title = 'Daily Return', tickformat = '%'),
    xaxis = list(title = 'Date'),
    title = 'GME Daily Returns'
  )

The traditional measure of risk in finance is the volatility of these returns.

# before
sd(GME %>% filter(date <= bpd) %>% pull(returns))
## [1] 0.03515004
sd(GME %>% filter(date <= bpd) %>% pull(returns))*sqrt(252) # annualized
## [1] 0.5579896
# after
sd(GME %>% pull(returns))
## [1] 0.04402463
sd(GME %>% pull(returns))*sqrt(252) # annualized
## [1] 0.6988693
# total observations
nrow(GME)
## [1] 4774
# fraction of observations
5/nrow(GME)
## [1] 0.00104734

So, 0.1% of observations (five days, this week) is enough to increase the daily returns’ standard deviation by almost 1 pp, and annualized volatility by about 14 pp. We have more than 4700 observations, which could be more than enough to apply the approximations implied by the Central Limit Theorem, yet so few observations impact the estimate so much…

Let’s have some fun and assume a Gaussian distribution. How frequently can we expect a return larger than the return observed on a given date?

m = mean(GME %>% pull(returns))
sd = sd(GME %>% pull(returns))

GME %>%
  filter(returns > 0) %>%
  mutate(gaussProb = 1/(252*(pnorm(returns, m, sd, lower.tail = F)))) %>%
  plot_ly(x=~date, y=~gaussProb, type = 'scatter', mode = 'lines') %>%
  layout(
    yaxis = list(title = 'Once in ... years',
                 type = 'log'),
    xaxis = list(title = 'Date'),
    title = 'Expected Frequency of GME Returns (Assuming Gaussian)'
  )

The January 27th, 2021 observation should only happen once in some \(10^{202}\) years. This is a huge number: for reference, the lifetime of the universe is approximately 13.8 billion years, or some \(10^{10}\) years, and this number does not get us anywhere close to \(10^{200}\).

Let’s start with what we found previously regarding the volatility. Is the standard deviation of returns a stable estimate? An MS plot is a way of visualizing the finite-sample stability of the moment of a distribution. The value of the MS series at time \(t\) is the ratio

\[ \frac{\max_{\tau \leq t} r_\tau^p}{\sum_{\tau \leq t} r_\tau^p} \]

In this case, we are interested in the second moment (corresponding to the variance/standard derivation), so \(p=2\). If moment \(p\) is finite, then the ratio should tend to zero (and stay at zero).

GME %>%
  mutate(ms = cummax(returns^2)/cumsum(returns^2)) %>%
  plot_ly(x=~date, y=~ms, type = 'scatter', mode = 'lines') %>%
  layout(
    yaxis = list(title = 'M/S'),
    xaxis = list(title = 'Date'),
    title = 'GME 2nd Moment Maximum-Sum Plot '
  )

If one looked at this plot sometime in 2016 or 2017, and ignored the December 2002 spike, one could have (naively) concluded that the second moment is finite. However, the December 2002 spike points towards a possible infinite second moment, and the latest spike reconfirms this. However, more importantly, an MS plot does not serve to confirm the finiteness of a moment, but rather to confirm its infiniteness. Epistemologically speaking, one should never conclude moment finiteness from an MS graph. In any case, this result suggests that \(\alpha < 2\).

We look at the Hill plot for estimates of the tail index \(\alpha\) before and after the breakpoint date. The Hill estimator is the closed-form maximum likelihood estimator for the Pareto distribution with CDF \(F(x) = 1-\left(\frac{x_m}{x}\right)^\alpha\). This estimator is \[ H(x; k) = \left( \frac{1}{k} \sum_{\kappa \leq k} \left(\log(r_{n-\kappa:n} - \log(r_{\kappa:n})) \right)\right)^{-1} \] where \(r_{i:n}\) is the order statistic \(i\), that is, \(r_{1:n}\) is the minimum and \(r_{n:n}\) is the maximum. We subset the daily returns by selecting only those which are positive, and so \(n\) is the total number of positive returns.

hill.ests = function(x){
  n=length(x)
  k = 1:n
  loggs = log(sort(x, decreasing = T))
  avesumlog = cumsum(loggs)/k
  alphahat = 1/c(NA, (avesumlog - loggs)[2:n])
  return(alphahat)
}

hill.ci = function(hill.ests, p = 0.05, type = 'lower'){
  if (type == 'lower'){
    ci = hill.ests + qnorm(p/2)*hill.ests/sqrt(1:length(hill.ests))
  } else if (type == 'upper'){
    ci = hill.ests + qnorm(1-p/2)*hill.ests/sqrt(1:length(hill.ests))
  }
  return(ci)
}

tibble(
  n = 1:nrow(GME %>% filter(returns > 0, date <= bpd)),
  hill = hill.ests(GME %>% filter(returns > 0, date <= bpd) %>% pull(returns)),
  ci.up = hill.ci(hill, type = 'upper'),
  ci.lo = hill.ci(hill, type = 'lower')
) %>%
  plot_ly(x=~n, y=~hill, type = 'scatter', mode = 'lines', name = 'Hill Estimate') %>%
  add_trace(y=~ci.up, name = 'Upper 5% CI', line = list(color = 'red', dash = 'dash')) %>%
  add_trace(y=~ci.lo, name = 'Lower 5% CI', line = list(color = 'red', dash = 'dash')) %>%
  add_trace(y=1.8, name = 'alpha = 1.8', line = list(color = 'black', dash = 'dot')) %>%
  layout(
    yaxis = list(title = 'alpha'),
    xaxis = list(title = 'k'),
    title = 'Hill Plot - Before "Shock"'
  )
tibble(
  n = 1:nrow(GME %>% filter(returns > 0)),
  hill = hill.ests(GME %>% filter(returns > 0) %>% pull(returns)),
  ci.up = hill.ci(hill, type = 'upper'),
  ci.lo = hill.ci(hill, type = 'lower')
) %>%
  plot_ly(x=~n, y=~hill, type = 'scatter', mode = 'lines', name = 'Hill Estimate') %>%
  add_trace(y=~ci.up, name = 'Upper 5% CI', line = list(color = 'red', dash = 'dash')) %>%
  add_trace(y=~ci.lo, name = 'Lower 5% CI', line = list(color = 'red', dash = 'dash')) %>%
  add_trace(y=1.8, name = 'alpha = 1.8', line = list(color = 'black', dash = 'dot')) %>%
  layout(
    yaxis = list(title = 'alpha'),
    xaxis = list(title = 'k'),
    title = 'Hill Plot - After "Shock"'
  )

The graph after the shock does produce lower estimates of the tail index for the relevant range of order statistics. Taleb suggests that \(\alpha \approx 1.8\). Both Hill plots support this, and we can also get a threshold: \(k \approx 40\). This is a better fit for the second plot, but still definitely plausible for the first. This helps us determine the scale/location parameter of a Pareto distribution, \(x_m\).

GME %>% filter(returns > 0) %>% arrange(desc(returns)) %>% slice(c(39,40,41))
## # A tibble: 3 x 2
##   date       returns
##   <date>       <dbl>
## 1 2007-03-27   0.114
## 2 2020-03-19   0.111
## 3 2008-12-16   0.111

The threshold \(k \approx 30\) corresponds to \(x_m \approx 0.11\).

Basically, what we have just done is taken all returns higher than 13% and estimated a Pareto distribution using these data. Now, to make sure this makes sense, we use a Zipf plot. This is a plot of the log of the empirical survival function against the log of the observation. Manipulating the Pareto CDF a bit, we get: \[ \begin{aligned} F(x) &= 1 - \left(\frac{x_m}{x} \right)^\alpha \\ 1 - F(x) &= \left(\frac{x_m}{x} \right)^\alpha \\ \log \left(1 - F(x) \right) &= \alpha \log\left(\frac{x_m}{x} \right) \\ \log \left(1 - F(x) \right) &= \alpha \log\left(x_m\right) - \alpha \log\left(x\right) \\ \end{aligned} \]

So the plot should show a straight line with slope \(-\alpha\), and intercept \(\alpha \log(x_m)\). We do not use R’s built-in ecdf function, since it assigns a 1 to the maximum and a non-zero to the minimum. We want a zero for the minimum and a non-one to the maximum.

alpha = 1.8
x_m = 0.11

GME %>%
  filter(returns > x_m) %>%
  mutate(logsurv = log(1-(rank(returns)-1)/length(returns)),
         logret = log(returns),
         bef_aft = if_else(date > bpd, 'After', 'Before'),
         fit = alpha*log(x_m) - alpha*logret,
         textlab = paste0('Return: ', round(100*returns, 1), '%\n',
                          'Date: ', as.character(date), '\n')) %>%
  plot_ly(x=~logret, y=~logsurv, color = ~bef_aft, type = 'scatter', mode = 'markers',
          text = ~textlab, hoverinfo = 'text', colors = c('darkred', 'darkgreen')) %>%
  add_trace(y=~fit, name = 'Pareto', mode = 'lines', line = list(color = 'darkblue'),
            hoverinfo = 'skip', showlegend = F) %>%
  layout(xaxis = list(title = 'log(x)', zeroline = F),
         yaxis = list(title = 'log(1-F(x))', zeroline = F),
         title = 'GME Zipf Plot')

The tail decay exponent has remained the same pretty much since IPO (the earliest observation on this graph is February 18th, 2003). Currently, though, we just have more data for a better tail index estimate, but the tail was always there.

Let’s revisit the expected frequency plot, but this time with the Pareto distribution.

mu = x_m
beta = x_m/alpha

GME %>%
  filter(returns > 0) %>%
  mutate(gaussProb = 1/(252*(1-pgpd(returns, xi = 1/alpha, mu=mu, beta=beta)))) %>%
  plot_ly(x=~date, y=~gaussProb, type = 'scatter', mode = 'lines') %>%
  layout(
    yaxis = list(title = 'Once in ... years'),
    xaxis = list(title = 'Date'),
    title = 'Expected Frequency of GME Returns (Assuming Pareto)'
  )

This plot is still inaccurate: the January 27th, 2021 observation should occur once every 0.36 years, or about \(10^{-1}\). However, we have around 18 years of data (Feb 2003 - Jan 2021), so the correct expected frequency is probably between \(10^{0}\) and \(10^{2}\). The error of the expected frequency provided by the Pareto distribution is one to three orders of magnitude. The error of the Gaussian is 202 to 200 orders of magnitude.