The threshold is a dial

Six hundred records of two thousand readings from a normal parent, fitted seven times each at seven thresholds. The mean squared error of the estimated shape is smallest at the 0.925 quantile, at 0.03402, and 80.6% of it is still bias there rather than spread. A flat rule that always takes the 0.90 quantile scores 0.03411 on the same records — worse by a quarter of a per cent, which is to say by nothing.

Then the two diagnostics that exist so that the threshold need not be a flat rule are scored against that flat rule, one record at a time. The mean residual life plot comes out 1.590 times worse on average and is the worse answer on 55.2% of records. The shape-stability plot comes out 11.881 times worse and loses on 61.0%.

A knob whose best setting is indistinguishable from the default, and two ways of turning it that both make things worse: that is the finding, and it is not the interesting half. The interesting half is that the trade the knob is supposed to make is not a property of the estimator at all. Run the same code on a heavy-tailed parent and the error curve is monotone, with no interior optimum anywhere and 0.37% of its error being bias. The dial exists because of what the data are, not because of what is done to them.

One knob, and what is on each side of it #

A block-maxima analysis throws away every reading but the largest in each block. A peaks-over-threshold analysis keeps everything above a level $u$ and fits the excesses $y = x - u$ with a generalised Pareto law,

$$P(Y \le y) = 1 - \left(1 + \frac{\xi y}{\sigma}\right)^{-1/\xi}$$

whose shape $\xi$ is the same $\xi$ as the corresponding law for maxima. That shared shape is the point of the pair rather than a coincidence: a threshold fit and a block-maxima fit are two estimators of one number, and the essay that counts how often that number’s family is named correctly is about the other one of them.

The theorem says the excesses are generalised Pareto in the limit as the threshold rises. So a low threshold leaves many exceedances that are not quite generalised Pareto, and a high one leaves few that nearly are. Bias falls, spread rises, and somewhere between them there is a setting. This is the same shape as choosing a bandwidth for an estimated covariance, where both ends of the dial are wrong for different reasons, and it is the same shape as the block length that minimises a window’s bias not being the one that minimises its error.

The fit is Grimshaw’s exact reduction rather than a two-parameter search. Setting $\theta = \xi/\sigma$, the shape that maximises the likelihood at fixed $\theta$ is available in closed form,

$$\hat\xi(\theta) = \frac{1}{k}\sum_{i=1}^{k} \log\left(1 + \theta y_i\right)$$

so the whole two-parameter problem is a one-dimensional profile on a bracket that respects the only constraint there is. A one-dimensional search on a bracket cannot walk off a boundary that moves with the parameters, which is exactly what a simplex does here.

The trade, counted #

The exceedance count runs 399, 299, 199, 149, 99, 49, 19 across the seven thresholds — the whole usable range of a two-thousand-reading record, and the top of it is already thin. The bias runs −0.1977, −0.1841, −0.1708, −0.1656, −0.1674, −0.1816, −0.3771, and the spread runs 0.0490, 0.0559, 0.0701, 0.0812, 0.1025, 0.1665, 0.5268.

Two things in those rows are worth stopping on. The bias improves by about a sixth over the whole usable range and then reverses: at the 0.99 quantile it is more than twice the worst value anywhere else. That reversal is not a theorem failing; it is nineteen exceedances being too few for a maximum-likelihood shape to be a sensible object, and the fit collapsing toward the boundary of its own support on the records where the largest excess happens to be small. And the spread rises by a factor of eleven across the same range, which is the price the theorem’s limit is bought at.

The sum has an interior minimum, so the dial is real. What the sum does not have is a sharp one: the error at 0.90 is 0.03411, at 0.925 it is 0.03402 and at 0.95 it is 0.03854. Three neighbouring settings within a seventh of each other, with the best of them beating the middle of them by a quarter of a per cent. Any rule that lands anywhere in that stretch is as good as the best rule available, and any rule that leaves it is worse by a lot.

The trade belongs to the parent, not to the estimator #

The obvious suspicion about a bias–variance story is that the bias is the fitting code. It is not, and the control that says so is a second parent run through the identical machinery.

A Pareto’s exceedances are exactly generalised Pareto at every threshold above one. The conditional survival function of the excess over $u$ is $(1 + y/u)^{-\alpha}$, which is a generalised Pareto with shape $1/\alpha$ and scale $u/\alpha$, with nothing asymptotic in it. So that parent has no bias to buy, and its error is spread alone: it rises monotonically from 0.00574 at the 0.80 quantile to 0.21743 at the 0.99, and at its own best threshold — the lowest one on the grid — only 0.37% of the error is bias, against 80.6% for the light tail.

That single comparison is what makes the section’s claim sayable. Identical code, identical thresholds, identical record length, and one parent has a dial while the other has a slope. Whatever is being traded, it is a property of how far the parent’s own tail sits from the law being fitted to it, and it is the same distance the rate at which a maximum reaches its limit measures on the other side of the pair.

Even at the top of the grid, where nineteen exceedances make the estimate wild, the heavy-tailed parent’s error is almost all spread: a bias of −0.1561 against a spread of 0.4394. Compare the light tail at the same threshold — a bias of −0.3771 against a spread of 0.5268 — and the difference is the whole content of the word “trade”.

The plot that is read for straightness #

The first diagnostic is the mean residual life plot: the average excess over a threshold, plotted against that threshold, read for the level above which it becomes straight. The reason it is read that way is exact rather than heuristic — for a generalised Pareto the mean excess is

$$e(u) = \frac{\sigma + \xi u}{1 - \xi}$$

a straight line in $u$. So a stretch where the empirical version is straight is a stretch where the law fits.

The heavy-tailed parent’s plot is straight from the first point, and it is straight in a way that can be checked against arithmetic: for $\alpha = 2$ the formula gives $e(u) = u$ exactly, and the sweep reads 2.2395 at a threshold of 2.2361 and 10.0887 at a threshold of 10. Two routes, one line.

The light-tailed parent’s plot falls, from 0.5593 at the 0.80 quantile to 0.3389 at the 0.99, and it curves throughout. There is no straight stretch in it. A practitioner reading that plot is being asked to find a feature that is not present in the range the data reach, and what they will find instead is wherever the noise happens to look flattest.

Both rules lose, and one loses badly #

A plot cannot be scored. So each diagnostic is written as the decision a careful reader of its plot would make, and then run on every record. The mean residual life rule takes the lowest level from which the remaining mean excesses sit within one standard error of the least-squares line through them — which is what “straight from here up” means when it is made precise. The stability rule takes the lowest level whose 95% interval for the shape contains every shape estimated above it, which is what “flat from here up” means.

The two rules fail in opposite directions and for the same reason. The mean residual life rule piles up at the bottom — 59.0% of records at the 0.80 quantile, and a median pick of 0.80 — because it is looking for the lowest level from which the plot is straight, and on a light tail the plot is straight nowhere, so the criterion is met first by whatever noise it can accommodate. The stability rule piles up at the top: 37.7% at 0.99 and 21.2% at 0.975, because a wide interval at a high threshold contains everything above it trivially. A criterion of the form “this estimate’s interval covers the ones above it” is satisfied best by the estimate with the widest interval.

Scored against the fixed rule on the same records: the mean residual life rule is 1.590 times worse and loses on 55.2%; the stability rule is 11.881 times worse and loses on 61.0%. The twelvefold figure is the 0.99 column doing all the work — a mean squared error of 0.41979 sits at that threshold, twelve times the error anywhere else on the grid, and a rule that lands there on three records in eight inherits it.

And the rules do not improve where the model is exactly right. On the heavy-tailed parent, where every threshold is a correct threshold and the only thing to lose is exceedances, the stability rule is 15.245 times worse than the fixed quantile and loses on 65.0% of records. It is worse there than on the parent it was supposed to be diagnosing. That is what a rule with a systematic pull toward the top of its grid does when the top of the grid is simply thinner: it is not diagnosing anything, it is selecting for wide intervals, and the parent has nothing to do with it.

The constant that won is not a good rule either #

The fixed 0.90 quantile beat both diagnostics, and it would be a mistake to read that as an endorsement of it.

On the light-tailed parent it wins because it lands inside the flat stretch, and the flat stretch is wide: 0.90, 0.925 and 0.95 score 0.03411, 0.03402 and 0.03854, so the constant is right by about a quarter of a per cent and would have been right by the same margin at three of the seven settings. On the heavy-tailed parent it is not right at all. There the error is monotone, the best threshold is the lowest one on the grid at 0.00574, and the fixed rule’s 0.01142 is almost exactly twice it. A constant cannot follow a curve that has no minimum.

And on that parent the mean residual life rule produces the one genuinely awkward reading in the sweep: it is 1.177 times worse than the constant on average and is the worse answer on only 35.7% of records. It wins nearly two records in three and loses the aggregate, which means the records it loses on it loses badly. A comparison reported only as a mean would have called it worse; one reported only as a win rate would have called it better; both were counted for exactly that reason, and neither alone would have been an honest summary.

The standard error stops working before the grid does #

There is a quieter number in the sweep, and it is the one a threshold rule actually reads. The closed-form asymptotic standard error of a generalised Pareto shape is $(1 + \xi)/\sqrt{k}$, and it is what every interval in a stability plot is drawn from.

Compared with the counted spread across records, it tracks well and then stops. At 199 exceedances it reports 0.0588 against a counted 0.0701. At nineteen it reports 0.1429 against a counted 0.5268 — a factor of 3.7 too small, at exactly the threshold the stability rule chooses most often. So the rule that lands at the 0.99 quantile is not merely landing somewhere noisy; it is landing somewhere its own error bars are wrong by a factor of nearly four, and it is landing there because they are.

This is the same failure as an interval computed from a formula that assumes independent observations, arriving from the other direction: there the count is right and the information in it is less than the count says, and here the count is simply too small for the asymptotics that turn it into an error. Both produce a reported precision that is not there, and the check that decides whether every interval in an analysis is trustworthy has no counterpart here at all.

The half of the fit that was not there #

While this was being built, the generalised Pareto fit had no negative half, and nothing said so.

The profile was written as $k \log \theta - k \log \hat\xi$, which is the correct arithmetic when the shape is positive. When the shape is negative, $\hat\xi$ and $\theta$ are both negative and both logarithms are of negative numbers, so the profile evaluated to a not-a-number, which compares false against every candidate the search considers, so the search never improved on its $\xi = 0$ fallback and returned it. Written as $k \log(\hat\xi/\theta)$ — the same arithmetic for a positive shape, and the only form that exists for a negative one, since the ratio is the scale and the scale is positive — it works.

Before the fix, every bounded parent was fitted at a shape of exactly zero. Nothing threw. The numbers that came out were plausible, the figures drew, and the entire Weibull half of the subject — the half in which a fitted shape asserts that the quantity has a ceiling — simply did not exist. The symptom was that a whole domain of the parameter space was absent, and the absence is what nothing looks for: every check asks whether a number that is there is right. A walk that reaches exactly half its admissible set and reports the two-sided p-value correctly for ever is the same defect in a different subject, and it survived for the same reason.

Where this does not hold, and what was not tried #

Two limits, both real.

The parent is known here and is not known anywhere else. Every bias above is a difference from a number that is arithmetic. A record whose parent is unknown is exactly the case the diagnostics exist for, and nothing measured here can say what they cost there — only what they cost when the answer is known. The finding that they lose to a fixed quantile is therefore a finding about a controlled comparison, and the honest reading of it is that a diagnostic which cannot beat a constant when the truth is available has not earned the benefit of the doubt when it is not.

And the whole section prices the threshold rather than the estimator. Eighty per cent of the error at the best threshold is bias, so a corrected shape estimator would move more than any threshold rule can — a maximum-penalised likelihood, or the Zhang–Stephens estimator, either of which changes the number being minimised rather than the point it is minimised at. Neither was run. Adding one would have made the comparison two-dimensional and would not have answered the question that was asked, which was whether a diagnostic beats a constant. It does not, and the reason it does not is that on a light tail there is nothing in the plot to read.

The last thing worth carrying is that the two halves of the subject agree. A block-maxima fit on a light-tailed parent carries a bias that falls with the block size and not with the number of blocks; a threshold fit on the same parent carries a bias that falls with the threshold and not with the record length. Both are the distance between the parent’s own tail and the law being fitted to it, and in both cases the dial that reduces it is the dial that spends the data. The essay in this collection that measured a rule against the length nobody has named the same trouble in a different setting: the setting that would be best is not a setting anybody can aim at, and a rule that claims to aim at it should be scored rather than believed.