Up-to-date probabilistic temperature climatologies

We considered homogeneity-adjusted, quality-controlled time series of monthly mean surface air temperature from the 7279 stations in the Global Historical Climatology Network (GHCN) [11, 12]. To evaluate climatology construction methods, we constructed hindcasts for each of the last 30 years (1984–2013) for each station that had at least 50 years of valid previous data. The 3222 stations meeting this criterion were distributed worldwide, although concentrated in densely populated and industrialized areas and in the United States (figure 1), for a total of 674 964 station-months with valid hindcasts. Each hindcast used only observations from the years preceding the hindcast year.

Zoom In Zoom Out Reset image size

We assumed that the yearly time series of temperature for a given calendar month T(t) can be represented as the sum of a smooth trend component $\bar{T}(t)$ and a zero-mean high-frequency component $\epsilon (t)$ [13]:

$$\begin{eqnarray}&&T(t)=\bar{T}(t)+\epsilon (t).\end{eqnarray} \tag{ 1 }$$

Using this framework, our goal was to estimate a probability distribution for T at a given year t_f, given observations from previous years. We took the expectation of $T({{t}_{f}})$, or T_f for short, to be equal to $\bar{T}({{t}_{f}}),$ thus neglecting any information that might be available, particularly for short lead times, about the high-frequency component $\epsilon ({{t}_{f}}),$ which various methods of seasonal forecasting attempt to estimate [10, 14–16]. The probability distribution of $T({{t}_{f}})$ therefore includes uncertainty due to $\epsilon ({{t}_{f}}),$ which will be the same for all the methods considered here, as well as uncertainty in estimating $\bar{T}({{t}_{f}})$ from available previous observations.

A common class of methods for estimating T_f is the moving average $MA(n),$ where the estimate is the average value of T over the previous n years for which observations are available. Previous comparisons focusing on the United States [6, 7] found that $n\approx 15$ minimizes mean square error for seasonal temperature series over recent years. The 30-year standard WMO averaging corresponds to $MA$(30) (although here we update the averaging period annually rather than only every 10–30 years as is typically done for climatological standard normals), whereas with a value of n at least equal to the observational record length $MA\;$ reduces to the average of all previous observations, which would be the optimal estimate of the expected value in a stationary climate. Here we considered values of n equal to 5, 10, 15, 30, 60, and 90 years.

A second method considered for estimating T_f is the exponentially weighted moving average $EW(\tau )$. Unlike the moving average, this uses all the past observations T but gives more weighting to the more recent observations, which offers theoretically optimal predictions if the trend follows a random walk [17]. The behavior of the exponentially weighted moving average depends on its e-folding weighting timescale $\tau .$ Here we considered timescales τ equal to 5, 10, 15, 30, 60, and 90 years, allowing direct comparison with the simple moving average results. Large τ correspond to weighting all past observations equally, whereas small τ correspond to giving significant weights only to the most recent observations.

A third method considered is the hinge fit [2]. This uses only observations since 1940 and estimates T_f by piecewise linear regression with a break point at 1975. At the break point, the fitted temperature is continuous, but its first derivative is discontinuous. One version of the hinge fit (${\rm H}{{{\rm F}}_{1}}$) assumes that $\bar{T}$ was constant between 1940 and 1975, whereas a second version (${\rm H}{{{\rm F}}_{2}}$) fits a slope for $\bar{T}$ between 1940 and 1975 as well as after 1975 [7].

We also tested several more complex methods, including polynomial fits, linear regression using atmospheric CO₂ concentration as a predictor [18], and a smoothing cubic spline fit to all previous observations with smoothing parameter p [19, chapter XIV]. Since in our preliminary tests these methods did not outperform the simpler methods just described for monthly temperature prediction at individual stations, we do not consider them further here.

Each of the above methods, for a given choice of parameter (e.g., n or τ), estimates $\bar{T}({{t}_{f}})$ as a linear combination of the observed $T.$ If $\bar{T}$ actually follows the assumed linear model and the terms are independent identically distributed Gaussian variables, the predictive distribution of T_f given the observations T can be calculated as a Student t distribution centered at the estimated $\bar{T}({{t}_{f}}),$ giving the required probability distribution for the expected temperature at t_f [20].

The t distribution with mean $\mu ,$ scale parameter $\sigma ,$ and ν degrees of freedom is given by

$$\begin{eqnarray}&&\begin{array}{l} p({{T}_{f}})=\frac{\Gamma (\frac{\nu +1}{2})}{\Gamma (\frac{\nu }{2})\sqrt{\pi \nu }\sigma } \\ \quad \times {{\left( 1+\frac{1}{\nu }{{\left( \frac{{{T}_{f}}-\mu }{\sigma } \right)}^{2}} \right)}^{-\frac{\nu +1}{2}}}. \\ \end{array}\end{eqnarray} \tag{ 2 }$$

For MA(n), μ is simply the mean of the past n observations, σ is the sample standard deviation of the observations, and ν is equal to $n-1.$

For EW(τ), define $\gamma =\frac{1}{{{\tau }^{2}}+1},$ quantifying the relative magnitude of the trend compared with interannual variability [21]. Then we get

$$\begin{eqnarray}&&\mu ={{({{{\bf X}}^{T}}{\bf PX})}^{-1}}{\bf XPT},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\sigma =\sqrt{\frac{1-\gamma +{{({{{\bf X}}^{T}}{\bf PX})}^{-1}}}{n-1}{{({\bf T}-\mu )}^{T}}{\bf P}({\bf T}-\mu )},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\nu =n-1,\end{eqnarray} \tag{ 5 }$$

where ${\bf T}$ is the $n\times 1$ vector of observed past temperatures, ${\bf X}$ is an $n\times 1$ vector of ones, and ${\bf P}$ is the inverse of the model covariance matrix ${\bf Q},$ defined as

$$\begin{eqnarray}&&{\bf Q}=(1-\gamma ){\bf I}+\gamma {\bf R},\end{eqnarray} \tag{ 6 }$$

with ${\bf I}$ the $n\times n$ identity matrix and ${\bf R}$ the $n\times n$ matrix with elements

$$\begin{eqnarray}&&{{R}_{ij}}={\rm min} (\mid {{t}_{f}}+1-{{t}_{i}}\mid ,\mid {{t}_{f}}+1-{{t}_{j}}\mid ),\end{eqnarray} \tag{ 7 }$$

where t_f is the forecast time and t_i are the n past observation times [21].

For a general linear regression such as the hinge fit, define ${\bf X}$ to be the $n\times k$ matrix of predictor values at the observation times and ${{X}^{*}}$ to be the $1\times k$ vector of predictor values at the forecast time ${{t}_{f}}.$ Then

$$\begin{eqnarray}&&\mu ={{X}^{*}}{{({{{\bf X}}^{T}}{\bf X})}^{-1}}{\bf XT},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\sigma =\sqrt{\frac{1+{{X}^{*}}{{({{{\bf X}}^{T}}{\bf X})}^{-1}}{{X}^{*T}}}{n-k}{{({\bf T}-\mu )}^{T}}({\bf T}-\mu )},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\nu =n-k,\end{eqnarray} \tag{ 10 }$$

where k is the number of regression coefficients (2 for ${\rm H}{{{\rm F}}_{1}}$ and 3 for ${\rm H}{{{\rm F}}_{2}}$).

Bias–variance trade-off is characteristic of a wide variety of inference problems and methods [22, 23]. The temperature trend at t_f is best estimated by using recent observations, but these observations are not numerous enough to confidently separate the trend from interannual variability subsumed under $\epsilon .$ On the other hand, averaging observations over a long period will reduce the influence of $\epsilon$ but also obscure recent climate changes. This trade-off can be studied quantitatively with a bias–variance decomposition of mean square hindcast error (MSE):

$$\begin{eqnarray}&&\begin{array}{l} {\rm MSE}={\rm bia}{{{\rm s}}^{2}} \\ \qquad \quad +{\rm variance} \\ \end{array}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\begin{array}{l} \langle {{({{y}^{*}}-y)}^{2}}\rangle ={{\langle {{y}^{*}}-y\rangle }^{2}} \\ \quad +\langle {{({{y}^{*}}-y-\langle {{y}^{*}}-y\rangle )}^{2}}\rangle , \\ \end{array}\end{eqnarray} \tag{ 12 }$$

where $y*$ denotes the predicted expectation, y is the actual value, and $\langle \cdot \rangle$ denotes averaging over hindcasts. Figure 2 (see [24] for a somewhat similar concept) shows these bias and variance components with the GHCN data for $MA(n)$ and $EW(\tau )$ with averaging timescales $n,\tau$ ranging from 5 to 90 years. Small $n,\tau$ corresponds to little bias but high variance, whereas high $n,\tau$ leads to reduced variance but mean bias of $\gt 0.6$ K since current conditions are substantially warmer than the average of a long series of past observations. MSE is thus minimized at intermediate $n,\tau$ of $\approx 15$ years. The $EW$ curve is slightly below and to the left of the $MA$ curve on this plot, indicating that $EW$ offers improved bias–variance characteristics compared with what is possible with $MA.$ On the other hand, the hinge fits offer low bias but mean variance higher than any of the $MA$ or $EW$ cases (figure 2).

Zoom In Zoom Out Reset image size

It is possible to sharply reduce the bias by consideration of the global mean warming trend. This is because the amplitude of the interannual variability is much smaller for global mean temperatures than for individual stations, so it is possible to estimate the global trend accurately. We estimated the global trend by spline smoothing of a global mean temperature series, with the smoothing parameter value p chosen to minimize the corrected Akaike information criterion [25, 26]. Accurate estimates of the expected bias, for example, for an $EW(\tau )$ or $MA(n)$ estimate, can be obtained using the less-variable global series, and this bias can then be subtracted from the $EW(\tau )$ or $MA(n)$ expectation μ for each station. As expected, this bias is greater for longer averaging timescales and shows some variability and a general increase in magnitude over time, particularly for longer τ (figure 3). The MA case is similar, with an averaging timescale even as short as 15 years [7] having a consistent, albeit modest $(\sim 0.2\;{\rm K})$, negative bias because it does not fully capture recent warming.

Zoom In Zoom Out Reset image size

Adjusting the forecast probability distributions by an additive factor as shown in figure 3 does greatly reduce the bias term even for large τ ('adjusted' curves in figure 2). After bias adjustment, the minimum hindcast root mean square error decreases ∼2% from ∼1.92 K to ∼1.89 K, and this minimum is reached at $n\approx 30$ years (figure 3).

An example of the hindcast expected values for $EW(15)$ with and without bias adjustment is shown in figure 4(a). The unadjusted $EW(15)$ method has a clear low bias, which adjustment largely removed through translating the probability distribution (i.e., changing the mean μ) by 0.16 K for 1984, increasing to 0.39 K by 2013 (see figure 3). The adjusted $EW(15)$ probabilistic hindcasts (figure 4(b)) show the large effect of the warming trend particularly on the expected frequencies of temperatures at the hot and cold extremes of the distributions.

Zoom In Zoom Out Reset image size

MSE (or its square root, RMSE) and its bias and variance components are useful and widely employed for evaluating the skill of point estimates of climatology. Additional metrics are needed, however, to evaluate the skill of probabilistic climatologies. The skill of a probabilistic forecast may be evaluated by how probable the actual outcomes are in the forecast. This can be quantified with the mean negative log likelihood (NLL):

$$\begin{eqnarray}&&{\rm NLL}=\langle -{\rm log} p(y)\rangle ,\end{eqnarray} \tag{ 13 }$$

where p is the hindcast probability distribution and y is the observed temperature. The NLL values for different forecast methods have units of information (bits or nats, depending on the base of the logarithm taken—here we take natural logarithms) and can be related to the methods' ability to reduce uncertainty in a decision-making framework [10, 27–29]. NLL is a strictly proper forecast-evaluation skill score which is sensitive to both the mean of the forecast distribution and its variance [30, 31].

Particularly for assessing the likelihood of extreme events, it is also important that the method produce probability distributions that are consistent with the observed values. This can be assessed graphically by plotting the histogram of the hindcast cumulative distribution function at the observed temperature (known as the rank verification histogram [32]), which should approximate the uniform distribution, and numerically by calculating the Kolmogorov-Smirnov statistic D_n for the maximum deviation of the empirical cumulative distribution from the standard uniform distribution. D_n ranges between 0 and 1, and for n independent trials and a well-calibrated forecast method should tend to zero as ${{n}^{-1/2}}$ [33]. In practice we cannot compute exactly how small we expect this statistic to be under the assumptions of each method because this depends on the correlation of temperature data across different stations and months [32], which is difficult to estimate; still, lower values of this statistic generally indicate better-performing probabilistic forecast methods.

Here we averaged metrics such as MSE, NLL, and D_n across years and stations to produce global skill measures. The same metrics could also be averaged for spatiotemporal subsets in order to study, for example, whether the ranking of methods is consistent across regions or seasons.

Without adjustment, the MA and EW methods all tend to give temperature hindcasts that are biased cold, meaning that they underestimate the warming seen (table 1). Reducing n in the MA method or τ in the EW method reduces the magnitude of this bias, but at the cost of higher variance. RMSE, which combines bias and variance, is lowest for $EW(10)$, with $EW(15)$ close behind. The information measure NLL gives the same ranking. The negative bias means that EW and MA hindcasts present cold conditions as more probable, and hot conditions as less probable, than was actually the case (figure 5(a)), a mismatch that gives values of the Kolmogorov-Smirnov statistic around 0.07–0.10 for the most skillful EW and MA variants. The hinge fit methods have little bias ($\lt$ 0.1 K) and good Kolmogorov-Smirnov statistics of under 0.01, but higher variance, so their overall performance as measured by RMSE and NLL is worse than that of EW and MA methods (table 1).

Zoom In Zoom Out Reset image size

Table 1.  Temperature hindcast skill metrics for methods without additional bias adjustment. For all metrics shown, values closer to zero indicate a more skillful method. See the text for details of the methods compared. RMSE = root mean square error; NLL = mean negative log likelihood; D_n = Kolmogorov-Smirnov statistic.

Method	Bias (K)	$\sqrt{{\rm variance}}$ (K)	RMSE (K)	NLL (nats)	Dn (%)
MA(5)	−0.089	2.023	2.025	2.273	6.22
MA(10)	−0.174	1.934	1.942	1.974	5.74
MA(15)	−0.256	1.908	1.925	1.928	7.42
MA(30)	−0.433	1.895	1.944	1.926	11.53
MA(60)	−0.509	1.901	1.968	1.947	13.42
MA(90)	−0.571	1.910	1.993	1.966	14.80
EW(5)	−0.159	1.937	1.944	1.894	4.23
$EW(10)$	−0.274	1.896	1.916	1.885	7.28
EW(15)	−0.349	1.886	1.918	1.893	9.26
EW(30)	−0.465	1.886	1.943	1.921	12.22
EW(60)	−0.547	1.900	1.977	1.951	14.15
EW(90)	−0.573	1.908	1.992	1.962	14.76
${\rm H}{{{\rm F}}_{{}}}_{1}$	−0.014	2.042	2.042	1.940	0.44
${\rm H}{{{\rm F}}_{{}}}_{2}$	+0.033	2.044	2.044	1.940	0.69

After bias adjustment, the MA and EW methods have greatly reduced mean bias, generally under 0.1 K regardless of the exact parameter values used, whereas the hindcast variance is little affected (table 2). $EW(15)$ is now the best performer as measured by NLL, with a decrease of some 2% in RMSE and 0.04 bits in NLL due to bias adjustment, followed by $EW(10)$ and $EW(30)$. $EW(30)$ ranks best in terms of RMSE. (It is expected that bias adjustment will tend to increase the optimal values of the averaging timescale n or $\tau$ since now bias is a smaller contributor to the error, allowing us to take advantage of the reduced variance afforded by a longer averaging time.) The Kolmogorov-Smirnov statistic for these methods is around 0.01, indicating a much better calibrated probabilistic hindcast that is comparable in quality to the hinge fit results. This improvement can be clearly seen in rank verification histograms (figure 5(b)).

Table 2.  Temperature hindcast skill metrics. Same as table 1, but with global bias adjustment.

Method	Bias (K)	$\sqrt{{\rm variance}}$ (K)	RMSE (K)	NLL (nats)	Dn (%)
MA(5)	−0.025	2.021	2.022	2.269	5.16
MA(10)	−0.047	1.935	1.936	1.970	2.67
MA(15)	−0.063	1.911	1.912	1.916	2.23
MA(30)	−0.079	1.897	1.898	1.886	1.69
MA(60)	−0.055	1.899	1.899	1.883	0.92
MA(90)	−0.044	1.907	1.908	1.890	0.47
EW(5)	−0.032	1.939	1.939	1.889	0.91
$EW(10)$	−0.047	1.899	1.900	1.870	1.04
EW(15)	−0.052	1.889	1.889	1.867	1.04
EW(30)	−0.038	1.887	1.888	1.871	0.79
EW(60)	+0.044	1.899	1.900	1.886	2.20
EW(90)	+0.114	1.906	1.909	1.899	4.14
${\rm H}{{{\rm F}}_{{}}}_{1}$	−0.019	2.043	2.043	1.942	0.52
${\rm H}{{{\rm F}}_{{}}}_{2}$	−0.007	2.045	2.045	1.940	0.47

To test the sensitivity of our results to the geographic distribution of stations used, we also calculated RMSE and NLL global skill metrics with an alternative averaging procedure where stations that were close together were weighted less, thus reducing the influence of dense clusters of observations, such as those in the United States. The weight assigned to each station was proportional to the reciprocal of the local station density estimated by an exponential density kernel with a length scale of 3${}^\circ$ (see figure 1). The rankings of the methods did not change greatly under this weighting, but the optimal averaging timescale tended to be somewhat shorter than when stations were weighted equally (tables 3 and 4): for example, the method with lowest RMSE after bias adjustment was $EW(15)$, and the method with lowest NLL was $EW(10)$.

Table 3.  Temperature hindcast skill metrics for methods without additional bias adjustment. For all metrics shown, values closer to zero indicate a more skillful method. Same as table 1 except that here the averaging of metrics across stations downweights stations that are close together. RMSE = root mean square error; NLL = mean negative log likelihood.

Method	Bias (K)	$\sqrt{{\rm variance}}$ (K)	RMSE (K)	NLL (nats)
MA(5)	−0.078	1.856	1.858	2.023
MA(10)	−0.159	1.775	1.783	1.706
MA(15)	−0.244	1.765	1.782	1.670
MA(30)	−0.430	1.756	1.808	1.696
MA(60)	−0.585	1.762	1.856	1.776
MA(90)	−0.643	1.766	1.880	1.806
EW(5)	−0.147	1.780	1.786	1.612
$EW(10)$	−0.269	1.747	1.768	1.622
EW(15)	−0.357	1.741	1.777	1.648
EW(30)	−0.504	1.745	1.816	1.716
EW(60)	−0.606	1.758	1.859	1.776
EW(90)	−0.638	1.764	1.876	1.796
${\rm H}{{{\rm F}}_{{}}}_{1}$	+0.050	1.884	1.885	1.671
${\rm H}{{{\rm F}}_{{}}}_{2}$	+0.054	1.884	1.884	1.661

Table 4.  Temperature hindcast skill metrics for methods. Same as table 3, but with global bias adjustment.

Method	Bias (K)	$\sqrt{{\rm variance}}$ (K)	RMSE (K)	NLL (nats)
MA(5)	−0.009	1.856	1.856	2.024
MA(10)	−0.032	1.776	1.776	1.698
MA(15)	−0.048	1.765	1.766	1.654
MA(30)	−0.067	1.755	1.756	1.633
MA(60)	−0.104	1.758	1.761	1.647
MA(90)	−0.089	1.763	1.765	1.655
EW(5)	−0.016	1.780	1.780	1.606
$EW(10)$	−0.034	1.748	1.748	1.598
EW(15)	−0.049	1.741	1.742	1.603
EW(30)	−0.059	1.744	1.744	1.622
EW(60)	+0.010	1.756	1.756	1.653
EW(90)	+0.078	1.762	1.764	1.677
${\rm H}{{{\rm F}}_{{}}}_{1}$	+0.048	1.884	1.884	1.672
${\rm H}{{{\rm F}}_{{}}}_{2}$	+0.017	1.883	1.884	1.660

We conducted similar analyses using mean-monthly daily minimum and maximum temperatures from GHCN rather than mean temperatures. These gave qualitatively similar results to those seen for mean temperatures, including a bias–variance trade-off for the EW and MA methods and high error variance for the hinge fit methods. In these two cases bias-adjusted $EW(15)$ had the best overall performance as measured by both RMSE and NLL.

The bias-adjusted probabilistic EW(15) hindcasts can be used to estimate how monthly temperature quantiles have shifted over time. For example, over the GHCN stations, the 1st percentile monthly temperature value in 1984 (a one-in-100-year cold event) became on average the 0.32 percentile by 2013, i.e., the recurrence interval lengthened to over 300 years. The 99th percentile monthly temperature value in 1984 (a one-in-100-year hot event) became on average the 95.3 percentile by 2013, i.e., the recurrence interval dropped to under 25 years. Thus, extreme hot episodes have become much more likely in just a few decades, whereas extreme cold episodes have become less likely. For the less extreme 10 and 90 percentiles, the corresponding shifts are to 4.2 and 76.0, so the frequency of cold months has dropped whereas the frequency of hot months has increased by more than a factor of 2. These global average changes were fairly robust to the specific climatology estimation method used. The spatial distribution of the estimated 2013 frequencies of the onetime 100-year events (figure 6) suggests that the changes may tend to be greater in the tropics and in areas with maritime climates, where interannual variability is smaller compared with the warming trend [34].

Zoom In Zoom Out Reset image size