8 · A Whirlwind Tour — Some Ways Forward
In this final lecture we will take a tour of some ways in which the subject of statistical inference can be taken further. We will prove almost nothing, but rather show that the theory we have developed — in particular maximum likelihood and the analysis of the information content of a sample — is the natural point of departure for a number of active areas of modern statistics, some of them very applied, and some quite theoretical.
Most of what follows is still maximum likelihood, only applied to a richer model, to a different target, or to simulated rather than observed data. The estimation machinery of the previous two lectures serves as the foundation to everything (it is more or less the whole point of the course to realize this). We organize our tour around two themes: going beyond the mean, and adopting a new paradigm for inference.
1. Beyond the mean
In nearly every model we have studied, the modelling effort went into the expectation. We assumed \(X_i \sim N(\mu, \sigma^2)\) and estimated \(\mu\). We wrote \(Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i\) and estimated the conditional mean. The variance was a nuisance parameter, the shape of the distribution was fixed, and the tails were whatever the normal happened to give us. In many applications, and in finance especially, this is the wrong emphasis. The mean there is small and difficult to predict, while the variance, the asymmetry, and the size of the rare events are where the interesting structure lies. The first three stops on the tour each look at ways to deal with such structures.
Volatility models
Consider a long series of daily returns \(r_t\) on a stock index. Their mean is tiny and close to unpredictable, but their variability is not. Turbulent days follow turbulent days and calm follow calm, a phenomenon known as volatility clustering that was noted by Mandelbrot already in the 1960s. The sample autocorrelation of \(r_t\) is essentially zero, but the autocorrelations of \(r_t^2\) and \(|r_t|\) are clearly positive and slowly decaying. The returns are therefore close to uncorrelated but far from independent. The magnitude of today’s return tells us something about the magnitude of tomorrow’s return (but not the sign!).
A model with a single constant variance \(\sigma^2\) cannot represent this. Engle (1982) proposed instead to let the conditional variance depend on the recent past. Write \(r_t = \mu + \varepsilon_t\) with \(\varepsilon_t = \sigma_t z_t\), where the \(z_t\) are iid with mean \(0\) and variance \(1\) (standard normal, or a \(t\)-distribution if we want heavier tails), and let the conditional variance follow \[ \sigma_t^2 = \omega + \alpha \varepsilon_{t-1}^2, \] so that a large shock yesterday raises the variance today. This is the (first-order) ARCH model. Bollerslev (1986) added a feedback term in the variance itself and obtained the GARCH(1,1) model, which has become a standard tool today: \[ \sigma_t^2 = \omega + \alpha \varepsilon_{t-1}^2 + \beta \sigma_{t-1}^2, \qquad \omega > 0,\quad \alpha, \beta \geq 0. \]
Today’s variance is a combination of a baseline level \(\omega\), the most recent squared shock \(\varepsilon_{t-1}^2\) and yesterday’s variance \(\sigma_{t-1}^2\). The process is weakly stationary when \(\alpha + \beta < 1\), in which case the unconditional variance is \(\omega/(1 - \alpha - \beta)\), and values of \(\alpha + \beta\) close to one produce the long, slowly decaying memory in volatility that we see in real data.
The parameters are estimated by maximum likelihood as before, but with one important complication: the observations are no longer independent, so we cannot write the joint density as a product of marginals. Instead we factor it into one-step-ahead conditional densities using the law of total probability: \[ f(r_1, \ldots, r_T) = f(r_1) \prod_{t=2}^{T} f(r_t \mid r_{t-1}, \ldots, r_1), \] and since each conditional distribution, if we use normal innovations, is \(N(\mu, \sigma_t^2)\) — with \(\sigma_t^2\) being a known function of the past and the parameters — the log-likelihood is \[ \ell(\mu, \omega, \alpha, \beta) = -\frac{1}{2}\sum_{t=1}^{T}\left(\log(2\pi\sigma_t^2) + \frac{(r_t - \mu)^2}{\sigma_t^2}\right), \] which we maximise numerically. Everything we proved about the score, the Fisher information and the asymptotic normality of maximum likelihood estimation carries over under suitable conditions on the dependence in the series.
From here one can move on to more generally formulated stochastic volatility models where \(\sigma_t\) is itself an unobserved random process and the likelihood involves an integral that is not available in closed form, or to multivariate GARCH models for whole portfolios, and to the use of all of these in computing Value-at-Risk and in pricing options.
Tsay’s Analysis of Financial Time Series is a good place to continue.
Flexible distributions with GAMLSS
In ordinary regression, and in the generalized linear models that extend it, the covariates are allowed to move the mean of the response and nothing else. In the normal linear model \(E(Y \mid X) = \beta_0 + \beta_1 X\) while the variance is constant \(\sigma^2\) and the shape is fixed at the normal. A generalized linear model relaxes the distribution to other members of the exponential family, such as the binomial, Poisson or gamma, but the same restriction remains. The covariates act on a single parameter, the mean, through a link function, and any spread or skewness is either constant or a fixed function of the mean. In the Poisson model, for instance, the variance is not free at all since it is exactly equal to the mean.
Real data are often less obliging. The spread of household expenditure grows with income, the tails of a financial return are more heavy in some periods than in others, insurance claims are more skewed for some types of claims than for others. In each case it is not only the mean that depends on the covariates, but the scale and the shape as well.
The GAMLSS framework — the Generalized Additive Models for Location, Scale and Shape of Rigby and Stasinopoulos (2005), and a later book-length treatment — takes this seriously. We choose a response distribution with up to four parameters \[ Y \sim D(\mu, \sigma, \nu, \tau), \] interpreted as location, scale, skewness and kurtosis, and give each of them its own predictor: \[ g_1(\mu) = \eta_1, \quad g_2(\sigma) = \eta_2, \quad g_3(\nu) = \eta_3, \quad g_4(\tau) = \eta_4, \] where each \(g_k\) is a link function and each predictor \[ \eta_k = X_k \beta_k + \sum_{j=1}^{J} s_{kj}(x_j) \] combines linear terms with smooth additive functions \(s_{kj}\) of the covariates, for example splines. These smooth terms are the “additive” in the name, and they are estimated nonparametrically.
We are no longer estimating a conditional mean. We are estimating the entire conditional distribution of \(Y\) as a function of \(X\). Estimation is once again via maximum likelihood, but with a “penalty” on the roughness of the smooth terms, so that we maximize a penalized log-likelihood of the form \[ \ell(\beta) - \sum_{jk} \lambda_{jk} \int \big(s_{jk}''(u)\big)^2 du, \] in which the smoothing parameters \(\lambda_{jk}\) trade fidelity to the data against smoothness of the fitted functions. That penalty is an expression of a very common regularization idea.
Extreme value theory
The limit theory in this course is, effectively, a theory of averages. The central limit theorem tells us something about the fluctuations of \(\bar{X}\) around \(\mu\), but it says nothing about the behaviour of the largest observation in the sample. But many questions are really about the maximum, the rare and extreme event: the highest water level that a dam must withstand, the largest claim an insurer must survive, the worst daily loss against which a bank must hold capital. Using a normal model for such questions is not merely inconvenient, but dangerous, since the thin tails of the normal will systematically understate the probability of the very events we are trying to guard against.
The financial crisis of 2008 supplied an expensive illustration. To price mortgage-backed securities, banks and rating agencies relied on the so-called Gaussian copula of David X. Li, which models how a whole pool of mortgages might default together through a single correlation drawn from a normal dependence structure. But that structure has no tail dependence: it assigns a very small probability to the event that many borrowers default at once, which is exactly what happens in a housing crash. This is the multivariate cousin of thin tails; it is not only about marginal thin tails, but also the dependence between tails, that are assumed away in a “normal” world.
Anyway, there is a limit theory for the maximum that runs closely parallel to that of the mean. Let \(M_n = \max(X_1, \ldots, X_n)\) for an iid sample. Just as we centre and scale a sum by \(\sqrt{n}\) before it settles down to a normal limit, we center and scale the maximum by sequences \(a_n > 0\) and \(b_n\). The classical theorem by Fisher, Tippett and Gnedenko says, essentially, that if \[ \frac{M_n - b_n}{a_n} \overset{d}{\to} G \] for some non-degenerate limit \(G\), then \(G\) must belong to a single family: the generalized extreme value distribution \[ G(x) = \exp\left\{-\left[1 + \xi\left(\frac{x - \mu}{\sigma}\right)\right]^{-1/\xi}\right\}, \] defined for those \(x\) with \(1 + \xi(x-\mu)/\sigma > 0\). The shape parameter \(\xi\) alone governs the tail. The value \(\xi = 0\) gives a Gumbel type with a light, exponential-type tail. \(\xi > 0\) gives the Fréchet type with a heavy, polynomially decaying tail, which is the case relevant for financial losses. \(\xi < 0\) gives the Weibull type with a finite upper end-point.
That the limit can only be one of these three types, whatever the underlying distribution of the \(X_i\), is the same kind of universality that makes the CLT so useful, and it is what allows us to extrapolate, in a principled way, beyond the range of the data we have already seen.
In practice one of two strategies is used:
The block maxima approach divides the record into blocks, for instance years, and fits the generalized extreme value distribution to the block maxima.
The peak-over-threshold approach instead keeps every observation above a high threshold \(u\) and uses a companion result, the theorem of Pickands, Balkema and de Haan (1974–75), which states that the distribution of the exceedance \(X - u\), conditional on \(X > u\), converges as \(u\) grows to a generalized Pareto distribution.
The peaks-over-threshold method usually makes better use of the data, since it does not discard a large observation simply because it occurs in a block with an even larger observation. Coles (2001) gives a very readable account, and Embrechts, Klüppelberg and Mikosch (1997) treat the subject with finance and insurance in mind.
The setup is the same iid setup we used for the CLT, with the same flavour of universal limiting family, but with the maximum in place of the sum. Together with the volatility models above, it forms a small unit of knowledge/skills on tail risk, which we have not dealt with properly otherwise in this course.
2. New paradigms
The three topics so far introduced more sophisticated models, but we stuck with the same philosophy of inference: write down a likelihood, maximise it, and describe the uncertainty through the frequentist sampling distribution of the estimator. In the following three excursions, we will keep the models simple and recognizable but rather alter the philosophy of statistical inference.
Bayesian inference
The Bayesian treats the unknown parameter \(\theta\) not as a fixed constant but as a random quantity, and encodes whatever is known about it before the data are seen into a prior distribution \(\pi(\theta)\). The model supplies the likelihood \(L(\theta \mid x)\) exactly as before, and the two are combined by Bayes’ Theorem into the posterior distribution: \[ \pi(\theta \mid x) = \frac{L(\theta \mid x)\,\pi(\theta)}{\int L(\theta \mid x)\,\pi(\theta)\,d\theta} \;\propto\; L(\theta \mid x)\,\pi(\theta). \]
The likelihood — the central object of the latter part of this course — is the same as before. It is multiplied by the prior and renormalised so that the posterior integrates to one. All inference then follows from this single distribution: a point estimate is a summary of the posterior such as its mean or its mode, and the posterior mode coincides with the maximum likelihood estimate when the prior is flat, which is a nice relative position to what we already know. An interval estimate is a credible interval, a set of posterior probability of, say, \(0.95\), which is a direct probability statement about the value of \(\theta\) (which a confidence interval is not).
Example. Suppose \(X \mid \theta \sim \text{Binomial}(n, \theta)\) and we take a \(\text{Beta}(a,b)\) prior for \(\theta\). Then \[ \pi(\theta \mid x) \propto \theta^x (1-\theta)^{n-x} \cdot \theta^{a-1}(1-\theta)^{b-1} = \theta^{x+a-1}(1-\theta)^{n-x+b-1}, \] which we recognize as a \(B(x+a,\, n-x+b)\) distribution. The posterior is in the same family as the prior, a property called conjugacy, and its mean \[ \frac{x+a}{n+a+b} \] is a weighted average of the sample proportion \(x/n\) and the prior mean \(a/(a+b)\). As \(n\) grows, the data come to dominate and the influence of the prior fades. We see the same phenomenon for the normal mean under a normal prior and for the Poisson rate under a gamma prior.
Conjugacy is the exception and not the rule. For most models, the integral in the denominator cannot be done, and the posterior is known only up to its normalising constant. What made Bayesian inference widely usable was the realization that we do not need that constant in order to sample from the posterior. MCMC methods — Metropolis–Hastings, Gibbs sampling and HMC implemented in Stan — generate a dependent series whose stationary distribution is exactly \(\pi(\theta \mid x)\), and we summarize the posterior from the draws.
Prediction is particularly clean under Bayesian inference. Rather than plug a point estimate into the model, the Bayesian averages the model over the posterior, giving the posterior predictive distribution \[ p(x_{\text{new}} \mid x) = \int p(x_{\text{new}} \mid \theta)\,p(\theta \mid x)\,d\theta. \] The uncertainty about \(\theta\) is carried through into the prediction instead of being thrown away.
Want more? See Statistical Rethinking by McElreath (with accompanying YouTube lectures), Bayesian Data Analysis by Gelman et al. (modern classic), or relevant sections in C&B, starting with Section 7.2.3.
Predictive statistics
Every method that we have studied so far takes a model and a parameter as its central objects. We estimate \(\theta\), we summarise the information that a sample carries about \(\theta\), and we test hypotheses about \(\theta\). Clarke & Clarke, in Predictive Statistics (2018), argue that this is the wrong primary emphasis for much of modern data analysis, and that the central object instead should be the prediction of the next observation. Their argument is simple and hard to “unsee” once you have seen it. A prediction can be compared with reality and shown to be right or wrong, whereas a parameter inside a model that we already know to be a simplification cannot be checked in the same direct way. Predictive accuracy is an externally verifiable measuring tape, while closeness to a parameter in a possibly false model is not.
And furthermore: does it even make sense to analyze a model that does not predict well? Some might be provoked by the following argument: drawing conclusions based on a model without checking the predictive performance is nothing more than a circular argument, where the researcher first encodes his/her opinion about some effect of \(X\) on \(Y\) into a model specification (“it is linear, so let’s use linear regression”) and then concludes afterwards that “yes, indeed, it looks like changing \(X\) with one unit is associated with \(\beta\) units of change in \(Y\)!”. Nothing important can be learned by this exercise, according to the predictive view.
Clarke & Clarke organize problems by how far we are willing to trust our list of candidate models. In the \(\mathcal{M}\)-closed case the true data generating process is one of the candidates we have written down, and classical model selection is unproblematic. In the \(\mathcal{M}\)-complete case a true model exists in principle but is too complex for us to write down or to use, so we reason with simpler surrogate models while admitting that none of them is correct. In the \(\mathcal{M}\)-open case, which is the honest description of most messy problems, there is no usable true model at all and the only thing we can reasonably ask of a method is that it predicts well. The further we move from \(\mathcal{M}\)-closed toward \(\mathcal{M}\)-open, the less it makes sense to speak of estimating the “true” parameter, and the more natural it becomes to judge a procedure by its performance on data it has not yet seen.
This predictive emphasis puts much of what we have done so far in a new light: methods should be assessed by how they work out-of-sample, by scoring each observation against the prediction that was made for it before it was revealed, by cross-validation and proper scoring rules. The boundary between classical statistical inference and modern predictive modelling is therefore pretty soft. The book by Clarke & Clarke develops this view with great care.
The bootstrap
In our treatment of sampling and limit theory, we went to some trouble to find the sampling distribution of a statistic. We derived the distribution of \(\bar{X}\) and of \(S^2\), we used the CLT and the delta method to obtain approximate distributions of more complicated quantities (such as the MLE), but each new case requires its own analysis of this kind. For an awkward statistic (the sample median, a trimmed mean, the ratio of two estimates, a correlation coefficient, etc.) that analysis can be very hard or even impossible (for a given skill-set ☺).
Efron’s bootstrap gives a different answer to the same question. We do not know the true distribution \(F\) that generated the data, but we have a very natural estimate for it, namely the empirical distribution \(F_n\), which places probability \(1/n\) on each observed value. Much of the difficulty earlier was that the sampling distribution of a statistic depends on the unknown \(F\). The bootstrap simply substitutes \(F_n\) for \(F\) and reads off the consequences. Since drawing a sample from \(F_n\) is nothing more than drawing from the observed values with replacement, this substitution is something we can always carry out:
Draw a resample \(x_1^{*}, \ldots, x_n^{*}\) of size \(n\) from the data, with replacement, and compute the statistic on it. Repeat this \(B\) times. The spread of the \(B\) resampled values approximates the sampling distribution of the statistic, and from it we read off standard errors, biases, confidence intervals, and whatever we want, really.
A percentile confidence interval is simply the interval running from the 2.5th percentile to the 97.5th percentile of the bootstrap values. The justification is the same plug-in reasoning that underlies the method of moments, supported by the fact that \(F_n\) converges to \(F\): if \(F_n\) is close to \(F\), then the sampling behavior of the statistic computed under \(F_n\) is close to its behavior under \(F\). The argument is asymptotic, and there are statistics for which the bootstrap fails, but across a wide class of problems it works very well and is an alternative to figuring out analytical expressions for the asymptotic properties of a complicated statistic.
The basic idea has many extensions. The parametric bootstrap samples from a fitted parametric model rather than from the raw data. The block bootstrap resamples contiguous blocks of observations instead of single points, so as to preserve dependence in time series. Efron and Tibshirani’s An Introduction to the Bootstrap is the standard text on the subject.
And the name? It comes from the feeling that this method provides of “lifting oneself up from the bootstraps”.