In this final lecture we take a tour of some ways in which statistical inference can be taken further.
Some of them are very applied, some quite theoretical.
Most of what follows is still maximum likelihood, only applied to
The estimation machinery of the previous two lectures is the foundation for everything — it is more or less the whole point of the course to realise this.
We organise the tour around two themes: going beyond the mean, and adopting a new paradigm for inference.
In nearly every model we have studied, the modelling effort went into the expectation:
The variance was a nuisance parameter, the shape was fixed, and the tails were whatever the normal happened to give us.
In many applications — in finance especially — this is the wrong emphasis. The mean is small and hard to predict, while the variance, the asymmetry and the size of rare events are where the interesting structure lies.
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 follows calm: volatility clustering, noted by Mandelbrot already in the 1960s.
So returns are close to uncorrelated but far from independent: the magnitude of today’s return says something about tomorrow’s magnitude — but not its sign.
A model with a single constant variance \(\sigma^2\) cannot represent this. Engle (1982) 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 \(t\) if we want heavier tails), and let \[\sigma_t^2 = \omega + \alpha \varepsilon_{t-1}^2.\]
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:
The GARCH(1,1) model. \[\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 combines a baseline \(\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\), with unconditional variance \(\omega/(1 - \alpha - \beta)\). 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 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).\]
Each conditional distribution is \(N(\mu, \sigma_t^2)\) under normal innovations, with \(\sigma_t^2\) a known function of the past and the parameters.
\[\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 the MLE carries over, under suitable conditions on the dependence in the series.
From here one can move on to
Tsay’s Analysis of Financial Time Series is a good place to continue.
In ordinary regression — and in the generalized linear models that extend it — the covariates move the mean of the response and nothing else.
Covariates act on a single parameter through a link function. Any spread or skewness is constant, or a fixed function of the mean — in the Poisson model the variance is not free at all, being exactly equal to the mean.
In each case it is not only the mean that depends on the covariates, but the scale and the shape as well.
The Generalized Additive Models for Location, Scale and Shape of Rigby and Stasinopoulos (2005) take this seriously. 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 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 \[\eta_k = X_k \beta_k + \sum_{j=1}^{J} s_{kj}(x_j)\] combines linear terms with smooth additive functions \(s_{kj}\) — for example splines.
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 by maximum likelihood, but with a penalty on the roughness of the smooth terms — we maximise \[\ell(\beta) - \sum_{jk} \lambda_{jk} \int \bigl(s_{jk}''(u)\bigr)^2 du,\]
where 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.
The limit theory in this course is, effectively, a theory of averages. The CLT tells us about the fluctuations of \(\bar{X}\) around \(\mu\) — but says nothing about the largest observation in the sample.
Many questions are really about the maximum:
Using a normal model here is not merely inconvenient but dangerous: its thin tails systematically understate the probability of the very events we are guarding against.
The financial crisis of 2008. To price mortgage-backed securities, banks and rating agencies relied on the Gaussian copula of David X. Li, modelling 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 — exactly what happens in a housing crash.
This is the multivariate cousin of thin tails. It is not only about marginal thin tails, but the dependence between tails, assumed away in a “normal” world.
There is a limit theory for the maximum running 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 centre and scale the maximum by sequences \(a_n > 0\) and \(b_n\).
Fisher–Tippett–Gnedenko. 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.
\[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:
That the limit can only be one of three types, whatever the distribution of the \(X_i\), is the same universality that makes the CLT so useful — and it lets us extrapolate, in a principled way, beyond the range of the data we have seen.
Block maxima. Divide the record into blocks — say years — and fit the generalized extreme value distribution to the block maxima.
Peaks over threshold. Keep every observation above a high threshold \(u\), and use the companion theorem of Pickands, Balkema and de Haan (1974–75): the distribution of the exceedance \(X - u\), conditional on \(X > u\), converges as \(u\) grows to a generalized Pareto distribution.
Peaks-over-threshold usually makes better use of the data — it does not discard a large observation simply because it sits in a block with an even larger one.
Coles (2001) is a very readable account; Embrechts, Klüppelberg and Mikosch (1997) treat the subject with finance and insurance in mind.
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 next three excursions we keep the models simple and recognisable, and instead alter the philosophy of statistical inference.
The Bayesian treats \(\theta\) not as a fixed constant but as a random quantity, encoding what is known before the data are seen into a prior \(\pi(\theta)\).
The model supplies the likelihood \(L(\theta \mid x)\) exactly as before, and the two combine by Bayes’ theorem into the posterior:
\[\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 this course — is unchanged. It is multiplied by the prior and renormalised.
All inference follows from this single distribution:
That is a direct probability statement about the value of \(\theta\) — which a confidence interval is not.
Suppose \(X \mid \theta \sim \text{Binomial}(n, \theta)\) with a \(\text{Beta}(a, b)\) prior. 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 recognise as a \(\text{Beta}(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.
Conjugacy is the exception, 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 realisation that we do not need that constant in order to sample from the posterior.
MCMC methods — Metropolis–Hastings, Gibbs sampling, HMC as implemented in Stan — generate a dependent series whose stationary distribution is exactly \(\pi(\theta \mid x)\). We then summarise the posterior from the draws.
Rather than plug a point estimate into the model, the Bayesian averages the model over the posterior:
\[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? Statistical Rethinking by McElreath (with accompanying lectures), Bayesian Data Analysis by Gelman et al., or C&B from Section 7.2.3.
Every method we have studied takes a model and a parameter as its central objects: we estimate \(\theta\), summarise the information a sample carries about \(\theta\), and test hypotheses about \(\theta\).
Clarke & Clarke (Predictive Statistics, 2018) argue this is the wrong primary emphasis for much of modern data analysis — the central object should instead be the prediction of the next observation.
The argument is simple and hard to “unsee”: a prediction can be compared with reality and shown right or wrong, whereas a parameter inside a model we already know to be a simplification cannot be checked so directly.
Does it even make sense to analyse a model that does not predict well?
Drawing conclusions from a model without checking predictive performance is arguably a circular argument: the researcher first encodes an opinion about the 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\) by one unit is associated with \(\beta\) units of change in \(Y\)!”
Nothing important has been learned, according to the predictive view.
Clarke & Clarke organise problems by how far we trust our list of candidate models.
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.
Methods should be assessed by how they work out-of-sample:
The boundary between classical statistical inference and modern predictive modelling is therefore pretty soft.
In our treatment of sampling and limit theory we went to some trouble to find the sampling distribution of a statistic.
But each new case requires its own analysis. For an awkward statistic — the sample median, a trimmed mean, the ratio of two estimates, a correlation coefficient — that analysis can be very hard, or impossible.
We do not know the true distribution \(F\) that generated the data — but we have a very natural estimate: the empirical distribution \(F_n\), placing probability \(1/n\) on each observed value.
Much of the difficulty 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.
And 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 else we want.
A percentile confidence interval is simply the interval running from the 2.5th 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 \(F_n \to F\): if \(F_n\) is close to \(F\), the behaviour of the statistic computed under \(F_n\) is close to its behaviour 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 it is an alternative to deriving analytical asymptotics for a complicated statistic.
Efron and Tibshirani’s An Introduction to the Bootstrap is the standard text.
It comes from the feeling this method provides of “lifting oneself up from the bootstraps”.
Beyond the mean:
New paradigms: