NOTE: see the UPDATE at the end of the post.

Jeff Masters at Wunderblog (part of Weather Underground) reported that for the lower-48 states of the USA, every one of the last 13 months was in the top third of its historical distribution. He calculated the odds of that happening by random chance, in an unchanging climate, being only 1/3 to the 13th power, or a mere 1 chance out of about 1.6 million. Pretty small odds.



However, the calculation isn’t correct because month-to-month temperatures for the USA48 are not independent, they show autocorrelation. If last month was hotter-than-average, this month is more likely to be as well. This increases the probability for a run of 13 top-third months in a row. Michael Tobis pointed this out, among others.

Enter Lucia, who tries to get a better estimate by modelling the monthly temperature as an AR(1) (red noise) process and says the chance that the most recent 13 months will all be in the top third of their distribution is about 10% — way more than 1 in 1.6 million! But Lucia estimates the lag-1 autocorrelation not by using USA48 temperature, but by using global temperature. As a result she uses the rather ridiculous figure of 0.936 for the autocorrelation parameter.

Enter Anthony Watts, who states that he had “started on an essay to describe this meteostatistical failure last night, but Lucia beat me to it.” He pronounces that “Lucia reduces this probability estimate calculation to rubble” and quotes her estimate that the chance is about 10%.

Jeff Masters isn’t the only one to have miscalculated the odds, so did Lucia. At least she recognized this — but as of this writing, Watts still hasn’t mentioned Lucia’s update:



Update Wow! I didn’t realize the US temperatures had such low serial auto-correlation! I obtained data for the lower 48 states here: http://www7.ncdc.noaa.gov/CDO/CDODivisionalSelect.jsp Based on this, the lag 1 autocorrelation is R=.150, which is much lower than R=0.936. So ‘white noise’ isn’t such a bad model. I am getting a probability less than 1 in 100,000. I have to run the script longer to get the correct value!



She is now saying that the odds are around 1 in a million. Pretty small. Jeff Masters also noted that his calculation wasn’t correct, but that the odds of 13 in a row in the top third is still pretty darn small.

Last but not and in fact least, enter Willis Eschenbach, contributing what must be the dumbest “analysis” of this situation yet. He suggests that the number of months out of 13 which will be among the top third will follow a Poisson distribution. He then counts how many 13-month periods show each number of top-third months and uses that to fit a Poisson distribution, apparently by least-squares regression. According to Eschenbach, using just 116 June-to-June periods he estimates the Poisson parameter at 5.206, using all 1374 13-month periods he estimates 5.213. The first estimate gives a probability of 13 out of 13 being in the top third as 0.001817, the second gives a probability of 0.001836. Both estimates are much larger than Jeff Masters’ and Lucia’s revised estimates. In either case, we would expect about 2.5 such occurences in the period of record, but only one has been observed (the last 13 months).

First let’s point out that it’s a mistake for Eschenbach to estimate the Poisson parameter by fitting that distribution to observations using least-squares. Deviations from the Poisson distribution aren’t even close to following the normal distribution — they’re deviations from the Poisson distribution. There’s a much easier and better way to estimate the Poisson parameter, and Eschenbach’s method overestimates it by a sizeable amount.

Second, let’s point out that if we accept Eschenbach’s hypothesis, then during the period of record we would expect to have seen about 1 instance of 14 out of 13 months with temperature in the top third of the distribution. What would you guess is the actual probability of that?

This much is clear: in an unchanging climate the probability of the observed event — the last 13 months all in the top third of their probability distribution — is larger than Jeff Masters’ original estimate, but still pretty small.

How small?

We can’t be sure because we don’t know for sure what process governs the noise in temperature data for USA48. But we can come up with a better estimate than 1 out of 1.6 million. We can come up with a better estimate than Lucia’s original 10%, in fact she has already done so. We can certainly come up with a better model than one which predicts we should already have seen a 13-month period with 14 months in the top third.

First let’s take note of a relevant fact. The only reason we’re talking about a 13-month run, rather than a 12-month or 14-month or some other length run, is that the run has been 13 months; the number 13 was chosen because that’s what was observed. If we really want to know how unusual this is (in an unchanging climate) I think it’s more realistic to put the question this way: given that 13 months ago was in the top third distribution-wise, what’s the probability that all of the 12 which followed would be also?

I’ll use data from the National Climate Data Center. USA48 temperature differs from global temperature in many respects. For one thing, signal-to-noise ratio is smaller (i.e., the noise-to-signal ratio is bigger). For another thing, the autocorrelation is smaller — much smaller. In addition, the noise level is month-dependent. There’s much more inherent variation in temperature for winter months than for summer months:

We can compensate to some degree by transforming temperature, not to anomaly — the departure from average for the given month — but to normalized anomaly: departure from average for the given month scaled by the standard deviation for the given month.

The standardized anomalies follow the normal distribution, at least approximately (the Shapiro-Wilk test does not reject that hypothesis), so we can use it to estimate probabilities. But we must still account for autocorrelation. The lag-1 autocorrelation of the standardized anomalies is 0.2185, so we’ll use that figure in an AR(1) model for the noise. This isn’t correct, but it’s a reasonable approximation.

For a standardized anomaly to be in the top third of the normal distribution it must be at least as large as 0.4307. The chance of that happening is just 1/3! But — what if the previous month was also in the top third? For an AR(1) process, a given value is a multiple (the autocorrelation parameter) times the previous value plus white noise

.

Furthermore, for the variance of the red noise to be equal to (which equals 1 for our standardized anomalies) the variance of the white-noise must be

.

The observed standardized anomaly from June 2011 was 0.9724. Let’s treat each of the following 12 months as a normal variable standard deviation given by , and with mean given by its expected value given that June 2010 was 0.9724. That expected value will be 0.9724 times to the power of the number of months after June 2011.

The we can compute the probability of being in the top third for each of the following 12 months, and approximate the probability they will all be so. The result: about 1 out of 458000. That’s pretty small.

The calculation is far from perfect. Mainly, the individual monthly results aren’t independent. A better approach would be to model the following 12 months using a multivariate normal distribution. But then we have to integrate that distribution over a nontrivial region of 12-dimensional space! Maybe there’s some clever way to do so — but I’m not sufficiently motivated to spend effort on that.

Another good approach is, of course, Monte Carlo simulation (Lucia’s approach). Her latest (that I’ve seen) indicates a probability somewhere in the 1-in-a-million range.

This much is clear: the odds of what we’ve seen having happened in an unchanging climate are pretty small. Jeff Masters’ original estimate wasn’t right, but it does appear to be within an order of magnitude. Willis Eschenbach’s is not, quite apart from the fact that it doesn’t make sense.

P.S. If you really want a laugh, check out this comment on Eschenbach’s post.

UPDATE

Lucia’s latest estimate is that the probability (given no climate change) is 1 in 134381. That’s 0.0000074415.

I think the Monte Carlo estimation is a better method than my own calculation. The actual probability remains uncertain, but is surely a heckuva lot less than 10%, or Willis Eschenbach’s Poisson-distribution nonsense.