applying Bayesian analysis to Extreme Value Theory (EVT)

[Colin Rust]

The following came up in an EVT (Extreme Value Theory) analysis of data on extreme geomagnetic storms.  These matter because a really extreme storm could fry a lot of telecom, etc.  So it's helpful to have an estimate of how the likelihood of storms of various magnitudes.

Basically you have a data set (in this case with 67 points), and you know/believe the true distribution is drawn from a two parameter family (in this case Generalized Pareto).  You want to estimate the probability of an event above some extreme cutoff C, which is more that any point in your data set.  How do you go about this?  

1. One approach would be to fit the ML distribution in the family and use that. This gives a point estimate of the probability P(X>C).

2. #1 can be extended using bootstrapping: Draw from the data set with replacement a bunch of samples of the same size (67 points).  From this, you can infer a distribution for the probability P(X>C) around the point estimate from #1.

That (#1 and #2) exactly what David Roodman does here:

http://blog.givewell.org/2015/07/13/geomagnetic-storms-using-extreme-value-theory-to-gauge-the-risk/

(I also posted in the comments there.)

That approach generated some weird results.  In particular, for some smaller cutoffs D that actually are reached in the data, it assigns non-zero probability to P(X>D) being 0 (hope that's not too confusing).

So I was wondering if there's a better way than the bootstrap:  Instead assign a prior to the two parameter family of distributions (but what prior?).  Then, using the data, compute a posterior density over the two parameter family.  From there, you can then infer a distribution for P(X>C).

Any thoughts appreciated. This may all be a little confusing with densities on densities so to speak, but hopefully it makes sense.

[Michael Betancourt]

The Bayesian solution is to

a) Construct a likelihood for such events given a set of parameters, theta,

    p(data | theta)

b) Construct priors encoding any information about the parameters.  See,

    for example, https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations.

c) The likelihood, prior, and data specify a posterior distribution for the parameters.

    p(theta | data) \propto p(data | theta) p(theta)

d) The posterior distribution can then be used to compute the posterior predictive

    distribution,

    p(new data | data) = \int dtheta p(new data | theta) p(theta | data)

    and then the expected probability of a tail event,

    E [ P ] = \int_C^infty d(new data) p(new data | data).

In Stan this is accomplished by

i) Specifying the posterior and data in a Stan model

ii) Drawing theta samples from the posterior using HMC

iii) For each theta sample compute p(new data > c | theta)

iv) Average those values together to get an estimate of E [ P ].

You can take the generalized Pareto for your likelihood but you’ll have to

be careful with your priors.  The biggest issue is that these kinds of 

survival models are really hard to fit — lots of different parameters will

be consistent with the data so you will have lots of uncertainty in the

tail probabilities.  This uncertainty is one of the reasons why frequentist

survival models are so brittle.

You’ll have to brush up on Bayesian methods before attempting something like

this.  I recommend going through the examples in the manual and reading

some good texts.

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Aki Vehtari]

Do you have the data in easily accessible format? I was interested to experiment with GPD model, but couldn't find the 67 point data nor 373 point data mentioned in at the page you linked.

Aki

[Aki Vehtari]

As I'm interested in generalized Pareto for other applications, I made a Stan model for the extreme geomagnetic storms discussed at 

http://blog.givewell.org/2015/07/13/geomagnetic-storms-using-extreme-value-theory-to-gauge-the-risk/

David Roodman kindly provided the n=373 observations shown in the previous (non-Bayesian) analysis.

I've attached the Stan code, Matlab code and a figure showing the results. I hope my Matlab code is clear enough to be translated to other environments.

Stan doesn't have yet builtin generalized Pareto, so I wrote the necessary log density and ccdf functions, too.

I used uniform priors, since the model has only two parameters and n=373 (data are informative), and I don't have any expert knowledge about geomagnetic storms..

A more interesting model would take into account the variation in the activity over time...

[David Roodman]

Excellent, Aki.

One thing I think you need to consider though: Theory says that the tails of most distribution converge to a Generalized Pareto distribution as one moves along their tails toward infinity. Of course our data do not extend to infinity, so we have to stop moving to the right before we run out of data. This then raises the question of whether the GP well approximates the chunk of the tail we examine. So one issue in GP fitting is the choice of the threshold defining the tail. Too high and you have too little data; too low and you are too far from the theoretical GP distribution.

One popular diagnostic in choosing the threshold is the mean residual life plot, which is the average exceedance of all points above a threshold, as a function of possible thresholds. (E.g., see https://books.google.com/books?id=SonbBwAAQBAJ&lpg=PR4&dq=stuart%20coles%20extreme%20value%20theory&pg=PA80#v=onepage&q&f=false.) If a distribution is GP, the MRL will be linear, as shown in the Coles book linked to. So in practice, we look at the MRL plot and ask when it starts to become linear as we move to the right. As you can imagine, one hopes one's conclusions are robust to reasonable changes in the threshold selected; and there is a literature (maybe some of it Bayesian) on more formal methods for choosing the threshold. In effect, it is a third parameter you need to take on board in some way.

Attached is the MRL I used, made with my "extreme" module for Stata. It includes 95% confidence intervals. From this, I chose a threshold of 150, as noted in a footnote in the blog post. Below 150, it seems like there is a statistically significant bend.

So at a minimum it would be interesting to know what happens if you use that threshold, i.e., only fit to the data above 150. Certainly that choice is debatable, and I did not systematically explore alternatives. It appears that you used 100, right?

This matters because, as seems clear, at lower Dst levels the distribution is more power-law like, which if extrapolated produces a fatter tail. So the more of the distribution you fit to, the higher the risk of extremes you will project. To some extent the uncertainty here may be irreducible, since we don't know when GP behavior sets in. But I think some attention to the threshold issue would probably be worthwhile.

--David

[David Roodman]

By the way, attached is what I get if I change my threshold from 150 as in my posted analysis to the 100 that I think is used in yours. It brings my risk estimates much more into line with yours. I think this shows that the threshold matters more than the method here...though of course the Bayesian approach may well improve on my thinking about the choice of threshold, and the influence thereof.

Question: a subset of GP distributions have absolute maxima (those with \xi<0). So the posterior distribution could assign non-vanishing probability to the hypothesis that a Carrington-sized event is impossible. And yet it appears not to. Does the prior used above exclude the possibility of a finite maximum on geomagnetic storm size?

--David

[Aki Vehtari]

On Tuesday, November 10, 2015 at 6:19:54 PM UTC+2, David Roodman wrote:

One thing I think you need to consider though: Theory says that the tails of most distribution converge to a Generalized Pareto distribution as one moves along their tails toward infinity. Of course our data do not extend to infinity, so we have to stop moving to the right before we run out of data. This then raises the question of whether the GP well approximates the chunk of the tail we examine. So one issue in GP fitting is the choice of the threshold defining the tail. Too high and you have too little data; too low and you are too far from the theoretical GP distribution.

Yes, I'm aware of this bias/variance problem and I know there are heuristic methods for choosing the threshold.

 

So at a minimum it would be interesting to know what happens if you use that threshold, i.e., only fit to the data above 150. Certainly that choice is debatable, and I did not systematically explore alternatives. It appears that you used 100, right?

I used 100 as the plots on the web page showed observations from 100, and I didn't notice that some other threshold was used. I've attached the result with threshold 150. I don't think the posterior mean changes a lot, but the posterior variance increases a lot.

Question: a subset of GP distributions have absolute maxima (those with \xi<0). So the posterior distribution could assign non-vanishing probability to the hypothesis that a Carrington-sized event is impossible. And yet it appears not to. Does the prior used above exclude the possibility of a finite maximum on geomagnetic storm size?

 The model I used with threshold 100 had the possibility that maxima is the max(y). I changed the model for threshold 150, to use Carrington-size maxima. It would be easy also to restrict k>0 for infinite maxima.

Aki

[David Roodman]

That does look like a good fit.

If I understand right, the new graph embodies two changes from the old:

* Changed threshold from 100 to 150.

* Moved from pure uniform prior to imposition of assumption that the absolute maximum of the distribution exceeds 850 (or really, 850-150).

Is that right?

If so, would you be able to do it with just the first change? I think this would provide a useful comparator to what I did, and help me understand the influence of various choices on the conclusions.

--David

[David Roodman]

Aki and Colin,

This exchange has been very stimulating for me. It prompted me this afternoon to develop more fully the reason I hesitated to incorporate the information. But that also led me to want to revise my published estimate upward a bit.

My concern is this. I have an estimator for a parameter. The estimator is imprecise and so has some distribution. My (historically) preferred way of estimating the distribution can put non-zero probability on the belief that the Carrington event was impossible, contradicting reality. Colin seems to propose that I chop off the point mass off the low end of my distribution, which is at 0. Then, evidently, I take the mean or median of the remaining distribution. But I worry that I will systematically bias my central estimate upward by first chopping off one end of the distribution. I will be incorporating a highly nonrepresentative selection of pre-1957 information. Thus while the point mass may look impossible in itself, it may still be better to leave it there in order to keep my central estimate properly centered.

It would be different if we had complete Dst data back to 1859. Or really, it should go back even earlier, since if we keep the data set as short as possible while still embracing that historic extreme, the selection would still be biased. At any rate, if we had all that data, then we could incorporate not only the peak of 1859, but also all the information about the lows that occurred in 1859-57. Right now we are doing the first but not the second.

My theory that imposing Pr[|Dst|>850]>0 induces bias generates a testable prediction wIthin the estimation framework I have used and understand well. If I draw smaller and smaller samples from the data set we are working with, fit the GP to them, and throw out those that predict Pr[|Dst|>850]=0, the median Pr[|Dst|>850] for those that remain should rise. This is because as the samples shrink the estimator's variance grows, producing more extreme probability on the high and low sides. Those on the low side will all be at 0. So as this "tail" grows, chopping it off should induce more bias. On the other hand, if I don't throw away the implausible Pr[|Dst|>850]=0 estimates, the median shouldn't become more biased as the samples shrink.

I tested this with simulations. I found it to be true--after introducing one tweak. Maximum likelihood is biased in small samples; more precisely it has O(n^-1) bias. So shrinking the samples per se tends to change results. To compensate, I applied the Cox-Snell correction for the ML Generalized Pareto estimator. (See eqs 53 & 54 here.) Early on I had thought about using this in my published analysis as well. The paper just linked to shows through Monte Carlo simulations that the correction does work. But the paper does not check whether it reduces bias in estimates of return levels and return probabilities, nor whether it improves inference about them. And some preliminary Monte Carlo simulations of my own were discouraging. So I thought I should be careful about grabbing the shiny new thing.

Attached is a graph of the results. Recall that the full event data set has 134 observations (episodes of |Dst|>=100). For each sample size N from 134 down to 134/4 I drew 10,000 datasets of size N, with replacement. For each, I fit the GP, made the Cox-Snell bias correction, and computed the predicted probability Pr[|Dst|>=850 | |Dst|>=100]. I computed the mean and median probability with and without filtering for Pr[|Dst|>=850 | |Dst|>=100]>0.

The results suggest to me that only the unfiltered median is unbiased, or at least that it is less-biased. If I understand right, this would correspond in Aki's work to removing the requirement that Pr[Carrington]>0.

The good experience with the Cox-Snell bias-correction made me more inclined to include it my posted analysis. Here is a graph that incorporates it. If you compare it carefully to the one at http://blog.givewell.org/2015/07/13/geomagnetic-storms-using-extreme-value-theory-to-gauge-the-risk/, you'll see the GP line is slightly higher now. The 95% CI per decade for -850 is now [0%, 5.8%] instead of [0%,4.0%] and the median is 0.56%/decade.

Dinner calls...

--David

[Colin Rust]

David,

Do you see an advantage to your analysis based on finding the ML fit over the Bayesian analysis as Aki has done?   To me the latter seems much more natural.

I suspect you're right that in the context of the analysis you've done, rejecting resamplings for which the ML fit satisfies P(Quebec)=0 or even just the weaker condition that P(Carrrington)=0 doesn't make sense. 

In the context of a Bayesian analysis of the data set that includes Quebec of course you automatically have P(Quebec)>0.  Also, it's natural to impose on the prior that P(Carrington)>0.

Colin.

[David Roodman]

Hi Colin,

                No, I don’t see any strong advantage of the frequentist ML fit over the Bayesian approach, though my understanding of the latter is incomplete. My prior, as it were, is that the two should give similar answers if given similar inputs. Possibly Aki’s next answer will prove me wrong on that.

                I think there’s a difference between effectively imposing P(Quebec)>0 and P(Carrington)>0 because the first bit of information comes embedded in a complete data series for several decades, chosen with minimal arbitrariness, whereas the latter is selected precisely for its extremity and so isn’t representative of the larger data range it is part of.

                (In the bootstrapping method I use, it is of course possible for P(Quebec)=0 in some replications.)

--David

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/C6IxqRU1Cuw/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Aki Vehtari]

The parameters had the following constraints

  real<lower=0> sigma; 

  real<lower=-sigma/ymax> k; 

in the first model ymax was the largest observation 589-100 and in the

second model ymax was 850-150 as I interpreted your previous messages

that you want Carrington level event to be possible (the Carrington

level event was not includes as an observation). I also tested with

threshold 150 and ymax=589-150 as you wished and there is not much

difference (the figure attached).

Aki

[David Roodman]

Thank you, Aki.

I did not intend to imply that I thought we ought to require Carrington/850 to be possible. I just asked about that to try to understand better what you did. I do think however that this is pretty much what Colin was looking for.

 

In my bootstrapping, the only constraint is sigma>0. Except in the “filtered” regressions in the graph in my last message, I require 850 to be possible.

 

I tend to think now that just introducing this requirement introduces a kind of bias. Maybe there is a name for this in Bayesian analysis, right? If I am interested in inference about a parameter k and I introduce the prior k > kmin based on the historical record, but choose not to impose a max based on the historical record, then this in itself tends to raise my central estimate of k, right? Is there a name for this phenomenon?

 

--David

 

From: stan-...@googlegroups.com [mailto:stan-...@googlegroups.com] On Behalf Of Aki Vehtari
Sent: Thursday, November 12, 2015 3:00 PM
To: Stan users mailing list <stan-...@googlegroups.com>
Subject: [stan-users] Re: applying Bayesian analysis to Extreme Value Theory (EVT)

 

The parameters had the following constraints

--

[Michael Betancourt]

You've got to stop thinking about bias as the end all be all.  Aki's model takes the lower bound into account by assuming compatibility with the assumed generalized Pareto distribution, which is uses the data in a self-consistent way.  Imposing an explicit lower bound lilies by the data need not be self consistent and hence can lose information in the posterior.  

The only worry in Aki's model is not that there is bias, because again it is self consistent and hence optimal conditioned on the modeling assumptions, but rather that the modeling assumptions themselves may be wrong.  Now looking at the predictive check it's not hard to see an elbow that the generalized Pareto is too rigid to fit, but again this an issue with the assumed tail and not the lower bound.

You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Aki Vehtari]

On Friday, November 13, 2015 at 10:20:10 AM UTC+2, Michael Betancourt wrote:

Now looking at the predictive check it's not hard to see an elbow that the generalized Pareto is too rigid to fit

Discrepancy is not bad, it fits inside the uncertainty intervals. My guess is that this kind of elbow is due to changes in the overall activation level in time and that should be modeled for the best accuracy (instead of using a more complex univariate distribution). The attached figure shows the estimated event intensity for abs dst > 1000 using a log Gaussian Cox process model (with GPstuff). There's clearly a change in the amount activation (which I guess is related to the sun spot cycle).

Aki

[Michael Betancourt]

> Discrepancy is not bad, it fits inside the uncertainty intervals.

Except that you’re ignoring correlations — the predictive
distributions are largely random around the median but
the data are systematically correlated away.

> My guess is that this kind of elbow is due to changes in the overall activation level in time and that should be modeled for the best accuracy (instead of using a more complex univariate distribution).

Sure, but then you’re just creating a generative process
which is implicitly creating a more complex univariate
distribution when you marginalize the nuisance parameters
out. But I agree that it is the right way to move forward.

[Aki Vehtari]

On Saturday, November 14, 2015 at 11:54:03 PM UTC+2, Michael Betancourt wrote:

> My guess is that this kind of elbow is due to changes in the overall activation level in time and that should be modeled for the best accuracy (instead of using a more complex univariate distribution).

Sure, but then you’re just creating a generative process
which is implicitly creating a more complex univariate
distribution when you marginalize the nuisance parameters
out.

I guessed you would answer like this, but I was too lazy to write more accurately :)

Aki 

[Aki Vehtari]

On Saturday, November 14, 2015 at 11:54:03 PM UTC+2, Michael Betancourt wrote:

> Discrepancy is not bad, it fits inside the uncertainty intervals.

Except that you’re ignoring correlations — the predictive
distributions are largely random around the median but
the data are systematically correlated away.

These correlations are due to how the data is plotted in sorted order (in Matlab with loglog(y,(n:-1:1)/n,'o')). I had now time to check how these plot look like if I generate simulated data from a generalised Pareto distribution with known parameters. Most of the time the plots have similar looking elbows or bends (one ore more) in different directions. I've attached one example with k=0.25 and sigma=40.

[Colin Rust]

Thanks Aki, David, Michael, this has been a really interesting discussion.

A key issue for policy is the probability of a really extreme storm, like the Carrington event.  Under the ML distribution, David found a probability of ~0.005% per solar storm, which works out to 0.33% per decade or 6.4% probability of at least one over a 200 year period.  As David noted in his original post, those numbers seem low.  Indeed, eyeballing Aki's charts, we see numbers about a factor of 20x larger on a Bayesian analysis:  a probability of >0.1% per storm which works out to ~6% per decade or a probability of ~80% (computed as 1-(1-6%)^20 = 81%) over a 200 year period. 

That's a sufficiently large difference that I would describe it as qualitative:  Under the ML distribution, one could describe Carrington as a freakishly extreme event, that humanity was unlucky to have encountered in the period since the Industrial Revolution and thus should have limited weight in policy planning.  Under a Bayesian analysis that incorporates contributions from the non-ML distributions, we instead see Carrington as a magnitude of storm that might plausibly recur in the coming decades and so is very important for planning.

[David Roodman]

Hi Colin,

 

I think that summary is correct, but I would add that the simulations I did showed that if the last ~200 years are statistically identical to the observed 1957-2014 Dst history aside from Carrington, then estimates that factor in the Carrington occurrence and nothing else from pre-1957 will be upward-biased. So somehow I think we need to build in some prior about the pre-1957 Dst history.

 

Something that might help is the aa index (https://www.ngdc.noaa.gov/stp/geomag/aastar.html), which I believe is based on two antipodal, non-equatorial observatories and goes back to 1868. Maybe it can be used to estimate pre-1957 Dst values.

 

Aki, I’m not sure if this speaks to what you were thinking about exploring time dependency, but two cycles I would include are the semi-annual cycle, with peak Dst activity at the equinoxes, and the sunspot cycle, which I treat as taking 11 years and reaching a low on 1 Jan 2008.

 

--David

 

From: stan-...@googlegroups.com [mailto:stan-...@googlegroups.com] On Behalf Of Colin Rust
Sent: Sunday, November 22, 2015 10:19 AM
To: Stan users mailing list <stan-...@googlegroups.com>
Subject: [stan-users] Re: applying Bayesian analysis to Extreme Value Theory (EVT)

 

Thanks Aki, David, Michael, this has been a really interesting discussion.

A key issue for policy is the probability of a really extreme storm, like the Carrington event.  Under the ML distribution, David found a probability of ~0.005% per solar storm, which works out to 0.33% per decade or 6.4% probability of at least one over a 200 year period.  As David noted in his original post, those numbers seem low.  Indeed, eyeballing Aki's charts, we see numbers about a factor of 20x larger on a Bayesian analysis:  a probability of >0.1% per storm which works out to ~6% per decade or a probability of ~80% (computed as 1-(1-6%)^20 = 81%) over a 200 year period. 

That's a sufficiently large difference that I would describe it as qualitative:  Under the ML distribution, one could describe Carrington as a freakishly extreme event, that humanity was unlucky to have encountered in the period since the Industrial Revolution and thus should have limited weight in policy planning.  Under a Bayesian analysis that incorporates contributions from the non-ML distributions, we instead see Carrington as a magnitude of storm that might plausibly recur in the coming decades and so is very important for planning.

[Aki Vehtari]

I think that summary is correct, but I would add that the simulations I did showed that if the last ~200 years are statistically identical to the observed 1957-2014 Dst history aside from Carrington, then estimates that factor in the Carrington occurrence and nothing else from pre-1957 will be upward-biased. So somehow I think we need to build in some prior about the pre-1957 Dst history.

 

Something that might help is the aa index (https://www.ngdc.noaa.gov/stp/geomag/aastar.html), which I believe is based on two antipodal, non-equatorial observatories and goes back to 1868. Maybe it can be used to estimate pre-1957 Dst values.

Bayesian hierarchical modeling allows easily to combine information from two datasets, and it is possible to take into account that the data generating processes are different from each other. Taking into account additional historical data even without time dependent model would improve the accuracy.

 

Aki, I’m not sure if this speaks to what you were thinking about exploring time dependency, but two cycles I would include are the semi-annual cycle, with peak Dst activity at the equinoxes, and the sunspot cycle, which I treat as taking 11 years and reaching a low on 1 Jan 2008.

These two would be two important parts of such model.

Aki

 

[Aki Vehtari]

On Sunday, November 22, 2015 at 5:18:35 PM UTC+2, Colin Rust wrote:

Thanks Aki, David, Michael, this has been a really interesting discussion.

A key issue for policy is the probability of a really extreme storm, like the Carrington event.  Under the ML distribution, David found a probability of ~0.005% per solar storm, which works out to 0.33% per decade or 6.4% probability of at least one over a 200 year period.  As David noted in his original post, those numbers seem low.  Indeed, eyeballing Aki's charts, we see numbers about a factor of 20x larger on a Bayesian analysis:  a probability of >0.1% per storm which works out to ~6% per decade or a probability of ~80% (computed as 1-(1-6%)^20 = 81%) over a 200 year period. 

My Stan code is available, so you can in addition of eyballing compute more accurate values, too :)

Aki 

[Colin Rust]

On Sunday, November 22, 2015 at 10:40:43 AM UTC-5, David Roodman wrote:

 

I think that summary is correct, but I would add that the simulations I did showed that if the last ~200 years are statistically identical to the observed 1957-2014 Dst history aside from Carrington, then estimates that factor in the Carrington occurrence and nothing else from pre-1957 will be upward-biased. So somehow I think we need to build in some prior about the pre-1957 Dst history.

It's of course good to include more history if you can. But to be clear, I don't think there's a bias in the Bayesian analysis from including Carrington in the minimal way Aki did.  As I understand it, he only required that the process be consistent with the possibility of a Carrington-scale event. 

Fundamentally, what's going on here I think is the ML distribution may not give you a reasonable approximation of the probability of a really extreme event.  I think most of the probability density for a Carrington-scale event comes from lower likelihood distributions with much greater density there.  (I'm not trying to make a mathematically precise statement.  But if I were to try, the flavor would be something like:  Define a "really extreme event" as a function of sample size n as an event that occurs O(1) times.  So the definition of really extreme varies with the sample size.  Then you might see even in the limit n->infinity that the ratio of the probability from the ML distribution over the true probability might converge to something less than 1.  On the other hand, if you had a fixed definition of "extreme" rather than one that depended on n, then it would converge to 1.)

[David Roodman]

Aki, I downloaded the aa index data from ftp://ftp.ngdc.noaa.gov/STP/GEOMAGNETIC_DATA/AASTAR/aaindex. A version that may be easier to work with is at http://1drv.ms/1T8Xfqi.

--David

 

 

From: stan-...@googlegroups.com [mailto:stan-...@googlegroups.com] On Behalf Of Colin Rust

Sent: Sunday, November 22, 2015 5:56 PM
To: Stan users mailing list <stan-...@googlegroups.com>

--