Reconstructing past temperatures from natural proxies and estimated climate forcings using short and longmemory models
Abstract
We produce new reconstructions of Northern Hemisphere annually averaged temperature anomalies back to 1000 AD, and explore the effects of including external climate forcings within the reconstruction and of accounting for shortmemory and longmemory features. Our reconstructions are based on two linear models, with the first linking the latent temperature series to three main external forcings (solar irradiance, greenhouse gas concentration and volcanism), and the second linking the observed temperature proxy data (tree rings, sediment record, ice cores, etc.) to the unobserved temperature series. Uncertainty is captured with additive noise, and a rigorous statistical investigation of the correlation structure in the regression errors is conducted through systematic comparisons between reconstructions that assume no memory, shortmemory autoregressive models, and longmemory fractional Gaussian noise models.
We use Bayesian estimation to fit the model parameters and to perform separate reconstructions of landonly and combined landandmarine temperature anomalies. For model formulations that include forcings, both exploratory and Bayesian data analysis provide evidence against models with no memory. Model assessments indicate that models with no memory underestimate uncertainty. However, no single line of evidence is sufficient to favor shortmemory models over longmemory ones, or to favor the opposite choice. When forcings are not included, the longmemory models appear to be necessary. While including external climate forcings substantially improves the reconstruction, accurate reconstructions that exclude these forcings are vital for testing the fidelity of climate models used for future projections.
Finally, we use posterior samples of model parameters to arrive at an estimate of the transient climate response to greenhouse gas forcings of 2.5C (95% credible interval of [, ]C), which is on the high end of, but consistent with, the expertassessmentbased uncertainties given in the recent Fifth Assessment Report of the IPCC.
10.1214/14AOAS785 \volume8 \issue4 2014 \firstpage1966 \lastpage2001 \docsubtyFLA
Paleoclimate reconstruction
A]\fnmsLuis \snmBarboza\corref\thanksrefM1,T1label=e1]luisalberto.barboza@ucr.ac.cr,
B]\fnmsBo \snmLi\thanksrefM2,T2label=e2]libo@illinois.edu,
C]\fnmsMartin P. \snmTingley\thanksrefM3,M31,T3label=e3]mpt14@psu.edu
and D]\fnmsFrederi G. \snmViens\thanksrefM4,T4label=e4]viens@purdue.edu
\thankstextT1Supported in part by NSF Grant DMS0907321.
\thankstextT2Supported in part by NSF Grants DMS1007686 and DPP1418339.
\thankstextT3Supported in part by NSF Grant 1304309.
\thankstextT4Supported in part by NSF Grant DMS1407762 and by Government of Chile Conicyt MEC competition Grant 0011232.
External forcings \kwdlongmemory \kwdproxies \kwdtemperature reconstruction
1 Introduction
An understanding of recently observed and projected future climate changes [Stocker et al. (2013)] within the context of the natural variability and dynamics of the climate system requires accurate and precise reconstructions of past climate. As spatially widespread instrumental temperature observations extend back to only about 1850, it is necessary to turn to the noisy and sparsely distributed paleoclimate record to characterize natural climate variability on longer time scales. In addition, reconstructions of past climate allow for important outofsample assessments of the Atmosphere–Ocean General Circulation Models (GCM) that are used to project future climate under various emissions scenarios [MassonDelmotte et al. (2013), Flato et al. (2013)]. While there is now a rich tradition of inferring past climate from natural proxies, such as tree rings, corals, ice cores, lake floor sediment cores and measurement on speleothems [for recent reviews, see NRC (2006); Jones et al. (2009), Tingley et al. (2012)], many scientific and statistical challenges remain.
1.1 Paleoclimatology context
Reconstructions of past surface temperatures from networks of multiple proxy types are prevalent in the climate science literature of the last 15 years—notable examples include Overpeck et al. (1997), Mann, Bradley and Hughes (1998, 1999), Luterbacher et al. (2004), Moberg et al. (2005), Juckes et al. (2006), Mann et al. (2008a, 2009), Kaufman et al. (2009), Tingley and Huybers (2013) and PAGES 2k Consortium (2013). While these studies have substantially increased our understanding of past climate, limitations remain in terms of the statistical treatment and uncertainty quantification. As described in Tingley et al. (2012), the most commonly used approaches to paleoclimate reconstruction are all variants of multiple linear regression [see, e.g., Table 1 of Christiansen, Schmith and Thejll (2009)], regularized in some fashion to account for the “” problem in the estimation procedure. Examples of particularly popular estimation approaches include regularized variants of the Expectation–Maximization algorithm [Dempster, Laird and Rubin (1977), Mann et al. (2007, 2005), Rutherford et al. (2003, 2005), Schneider (2001), Zhang, Mann and Cook (2004), Steig et al. (2009)] and principal component regression [Cook, Briffa and Jones (1994), mannbradley (1998), Luterbacher et al. (2004), Wahl and Smerdon (2012)], which is sometimes combined with canonical correlation analysis [Smerdon et al. (2010)]. A common shortcoming of these studies lies in the limited propagation of parameter uncertainty through the model, including uncertainty in the estimation of regularization parameters; for further discussion see Schneider (2001), Smerdon et al. (2010) and the supplement to Wahl and Smerdon (2012).
Recently, hierarchical modeling and Bayesian inference techniques have been proposed and employed to reconstruct past climate from proxies[Haslett et al. (2006), Li, Nychka and Ammann (2010), Tingley and Huybers (2010a, 2010b, 2013), Werner, Luterbacher and Smerdon (2013)]. Hierarchical modeling is a natural framework for including the available scientific understanding of both the target climate process (e.g., annual surface temperature anomalies) and how the various natural proxies are causally affected by variations in the climate system. Bayesian inference, in turn, provides a cohesive framework for propagating uncertainty, while the posterior draws of the target climate quantity are a more statistically precise and scientifically useful result than a point estimate and associated uncertainty interval [Tingley et al. (2012)].
In this paper, we reconstruct Northern Hemisphere (NH) temperature anomalies over the past millennium using a hierarchical Bayesian model that describes temperature as linearly dependent on three important climate forcings: green house gas concentrations, volcanic aerosol concentrations and variations in solar irradiance. The proxies, in turn, are modeled as linear in the latent temperature process. Motivated by existing evidence of longrange correlation in temperature series [e.g., Brody, Syroka and Zervos (2002), Benth and ŠaltytėBenth (2005), Huybers and Curry (2006), Imbers et al. (2014)], we explore the effects of specifying white noise (no memory), autoregressive (short memory) and longmemory correlation structures for the two error processes. To our knowledge, this is the first ensemblebased paleoclimate reconstruction that includes the effects of climate forcings, and the first systematic investigation of error structure in the temperature reconstruction. As our method involves first reducing the proxy data set to a single time series, and then inferring hemispheric average temperature anomalies, rather than the spatial pattern, our analysis is a form of compositeplusscaling [Tingley et al. (2012)].
The external forcings used in the analysis are closely related to global temperature evolution. The Intergovernmental Panel on Climate Change (IPCC) has steadily increased its certainty level on stating the causal relationship between increasing atmospheric concentrations of anthropogenic greenhouse gases and increasing average global temperatures, reaching the “extremely likely” level of 95% confidence in 2013 [Bindoff et al. (2013)]. The relationship between solar irradiance and surface temperatures is studied in Crowley and Kim (1996), Lean, Beer and Bradley (1995), while Briffa et al. (1998), Crowley and Kim (1993), Crowley, Criste and Smith (1993) and Landrum et al. (2013) analyzed the effect of volcanic activity on global temperatures.
The conceptual study of Li, Nychka and Ammann (2010) demonstrated that temperature reconstructions are improved when information about the climate forcing is included in the reconstruction. We therefore explore the effects of including these three major external forcings in our reconstructions, reporting results for both cases. While the forcings are expected to improve the reconstructions, we note that reconstructions that exclude the forcings are necessary for the evaluations of GCMs [MassonDelmotte et al. (2013), Flato et al. (2013)] to avoid the circularity of using the same forcings in the simulation of past climate and the reconstruction used to assess the simulation.
1.2 Longmemory modeling and estimation challenges
To our knowledge, the error terms in all previous models for multiproxy climate reconstructions are assumed to be white or autoregressive [AR; see, e.g., Tingley et al. (2012)]. For instance, Li, Nychka and Ammann (2010), Tingley and Huybers (2010a) and McShane and Wyner (2011) use AR(1) or AR(2) errors, while reconstructions based on the Expectation–Maximization algorithm or principal component regression have generally not explicitly modeled temporal autocorrelation [see Section 8.7.4 of Tingley et al. (2012)].
The assessment of longmemory behavior in hierarchical models is complicated by the fact that graphs of the autocorrelation and partial autocorrelation functions (a.c.f. and p.a.c.f.) are generally not adequate diagnostic tools. In addition, the short data streams we are faced with disallow reliance on known asymptotic properties, while lack of selfsimilarity means that inference on one range of frequencies cannot apply to another. These issues are well known for widely used longmemory time series models, such as fractional autoregressive integrated moving average (FARIMA) models [Beran (1994)]. Misspecification of a longmemory process with a shortmemory model can lead to erroneously attributing longmemory effects to deterministic trends or external forcings, and thus will affect uncertainty quantification. Specifically, since longmemory models can exhibit larger asymptotic variances than their relatively shortmemory model analogues [see Chronopoulou, Viens and Tudor (2009) and references therein], reported uncertainty levels under memory misspecification can be lower than the nominal values.
Motivated by the limitations of the data, and our goal of using a robust model, we focus on a simple longmemory model: linear regression with fractional Gaussian noise (fGn) errors. The theoretical question of estimating memory length for nonselfsimilar models, such as our hierarchical linear model, is notoriously difficult. Asymptotic theory is still under development, and current work on highfrequency or increasinghorizon versions of our model cannot yet be considered definitive. Online Supplement A.1 in Barboza et al. (2014) provides brief background information on longmemory estimation, while further details can be found in references therein; see, in particular, Gneiting and Schlather (2004).
In the context of annual paleoclimate observations, time intervals cannot be assumed small, and the calibration period is short. On account of the long time intervals, we cannot use the local path behavior of the data (e.g., Hölder continuity) as a proxy for long memory—an approach that is possible for fGndriven models where high frequency data exists. Such models are asymptotically Höldercontinuous in the limit of ultrahigh frequency, with a single parameter that also governs long memory. On account of the short calibration period, methodologically sound results from low frequency increasinghorizon asymptotics [see Tudor and Viens (2007)] cannot be used to measure longrange dependence in our case, as there is simply not enough data. Instead we resort to a fully Bayesian framework to estimate all parameters, including those responsible for memory length, with the added benefit of a complete evaluation and propagation of uncertainty.
This article is structured as follows. Section 2 describes the data sets used in the reconstruction, and Section 3 gives the details of the hierarchical Bayesian models. Section 4 presents the results of our Bayesian reconstructions, including parameter posterior distributions and model validation metrics; it compares models with different error structures and which include or exclude the climate forcings. We also compare our results with previous reconstructions and discuss the estimation of transient climate response in Section 5 before summarizing our quantitative conclusions and discussing remaining challenges in Section 6. Two online supplements provide further details on the modeling framework and additional quantitative results [see Barboza et al. (2014)].
2 Data sets
The analysis makes use of three distinct data sources: instrumentally observed temperature anomalies (in Celsius) over the period 1900–1998; a suite of temperaturesensitive proxies over the period 1000–1998 taken from the database originally described in Mann et al. (2008a) and used additionally in Mann et al. (2009); and estimates of external climate forcings from 1000–1998 AD.
We make use of two different instrumental estimates of NH temperature anomalies, both developed by the Climate Research Unit of the University of East Anglia [Brohan et al. (2006)]. The CRUTEM3v data set (abbreviated hereafter as CRU) is an estimate of air surface temperature anomalies over land, while HadCRUT3v (hereafter abbreviated as HAD) is an estimate of combined land air and marine seasurface temperatures. These data sets are widely used for the calibration of proxybased climate reconstructions [e.g., Mann et al. (2008a), Luterbacher et al. (2004), Rutherford et al. (2005), Kaufman et al. (2009), McShane and Wyner (2011), Tingley and Huybers (2013)]. We make use of the varianceadjusted version of each data set to facilitate comparisons with results from Mann et al. (2008a). While both instrumental data sets extend back to 1850, we choose 1900–1998 as our calibration period, as the sparsity of instrumental observations results in less trustworthy hemispheric estimates prior to about 1900 [Smith (2010)].
The proxies used in our analysis are selected from the 1209 climatesensitive proxies originally compiled in Mann et al. (2008a). ^{4}^{4}4For more details on the data set, see the NOAAPaleoclimatology/World Data Center at \surlhttp://www.ncdc.noaa.gov/paleo/pubs/pcn/pcnproxy.html. This compilation brings together a wide array of proxy types, including tree ring widths and densities, marine sediment cores, speleothems (cave deposits), lacustrine sediment cores, ice cores, coral records and historical documentary information [see NRC (2006) and Jones et al. (2009) for further descriptions of each of these data types]. The proxy data are not raw observations, but are rather processed to remove nonclimatic variability, such as age effects associated with tree ring data. This type of processing results in a data product which may be more directly interpreted as “climate sensitive,” according to the paleoclimatology community. While it is common to base climate reconstructions on the postprocessed data, as is done here, we acknowledge that doing so does neglect the uncertainty inherent in the processing steps. We set aside for future research the challenge of including the processing of raw climate proxy observations into climatesensitive series within the hierarchical framework developed here. For further details concerning the processing of raw proxy observations, see, for example, NRC (2006), Jones et al. (2009).
Estimates of the external climate forcings—atmospheric greenhouse gas concentrations (C), solar irradiance (S) and volcanism (V)—are described and plotted in Li, Nychka and Ammann (2010) and described more fully in Ammann et al. (2007). The original greenhouse gas concentration time series is in units of CO equivalent in parts per million; the solar irradiance series is in Wattm and the volcanic series is an estimate of the climate forcing, in , derived sulphate measurements on ice cores [see Ammann et al. (2007) for further details].
3 Model specification
Hierarchical Bayesian models typically consist of three levels. The data level describes the likelihood of the observations conditional on a latent stochastic process. In our context, the latent process is the time series of NH mean temperature anomalies, and the observations are the proxies. The process level describes the parametric structure of the latent process—often with recourse to prior scientific information, such as knowledge of the underlying physical dynamics [e.g., Berliner, Wikle and Cressie (2000)]. Finally, the prior level provides closure and allows for Bayesian inference by providing prior distributions for all unknown parameters in the data and processlevels. For a general description of hierarchical modeling and Bayesian inference in the paleoclimate context, see Tingley et al. (2012). Following Li, Nychka and Ammann (2010), the datalevel models the proxies as a normal distribution with mean equal to a linear function of the latent, unobserved true temperatures, while the processlevel models the latent temperature process as normal with mean given by a linear function of the external forcings [Li, Nychka and Ammann (2010)]. We add to previous work by applying the model to actual proxy data, as opposed to using pseudo proxy experiments derived from climate model output [Li, Nychka and Ammann (2010)], as well as identifying appropriate memory lengths in the error structures of the residuals at both levels.
The Bayesian modeling framework is closely related to stochastic filtering methods. An interesting application of classical Kalman filtering [see Kalman and Bucy (1961)] to climatic reconstruction is in Lee, Zwiers and Tsao (2008), where the authors use forcings and a smaller proxy data set to reconstruct temperatures on a decadal basis. However, there are, to our knowledge, no practical tools for filtering with fGn errors and, in addition, stochastic filters, which are adapted to tracking moving signals dynamically in time, are notoriously poor at estimating fixed parameters; see Yang et al. (2008) and Chronopoulou and Viens (2012). Thus, they are not an optimal choice for our exploration of short versus longmemory models in paleoclimate reconstructions. In contrast, the Bayesian approach adopted here allows for all parameters to be estimated simultaneously while avoiding the known estimation difficulties inherent to filtering. Moreover, since the proxy observations are not being updated over time, the sequential updating property of filtering is not advantageous.
3.1 Proxy data reduction
It is desirable for several reasons to reduce the dimensionality of the proxy data set, which consists of 1209 time series. First, as there are only a limited number of years in the calibration interval, dimension reduction can lead to a more parsimonious model, avoid overfitting, and lead to more robust temperature reconstructions. Second, our interest in inferring global mean temperatures rather than spatial fields motivates a reduction, prior to fitting a hierarchical model, to a single time series that reflects the shared variability between the proxies that is likely attributable to a common, climatic origin. Third, the proxy reduction is important in limiting the computational burden of estimating parameters describing long memory; for a comparison between computational and asymptotic efficiency for various longmemory parameter estimators, see Chronopoulou and Viens (2010). We therefore apply a sequence of steps to reduce the number of proxies while attempting to retain as much climatically useful information as possible.
Following Mann et al. (2008a), we first select only those proxies that are recorded at least as far back as 1000 AD and, in addition, have a significant correlation with their closest instrumental time series (marine or land) over their period of mutual overlap. We use local temperature information in the screening procedure, as any proxy that might correlate to hemispheric temperature with some degree of accuracy should relate to its local temperature with higher precision [Mann et al. (2008a)]. Such a criterion does not take into account the possibility of exploiting physical teleconnections that exist in the actual climate system [mannbradley (1998), Tingley et al. (2012), Werner, Luterbacher and Smerdon (2013)]. This screening procedure yields 38 proxies whose distribution by type and location is given in Table 1. Tree rings represent the majority of proxies that pass the screening criteria, consistent with the ubiquitous use of tree ring information in annual resolution temperature reconstructions [NRC (2006), Jones et al. (2009), PAGES 2k Consortium (2013) and references therein].
Type  #  Locations 

Tree ring  USA, Argentina, Norway, New Zealand, Poland, Sweden  
Lacustrine  Mexico, Ecuador, Finland  
Speleothem  China, Scotland, Yemen, Costa Rica, South Africa  
Ice cores  Peru, Greenland, Canada  
Other\tabnotereftt1  China, Mongolia, Tasmania, New Zealand 
tt1The category named “Other” contains data from composite temperature reconstructions and historical documentary series.
A number of the 38 proxy series in Table 1 show undesirable properties given our assumption of a stationary relationship between the proxies and temperatures. In particular, several of the lacustrine and speleothem records feature much greater variability in the early portion of the time interval than in the calibration period. On such bases, we exclude 13 proxies, leaving a total of 25; see Figure B.1 and Table B.1 in the Online Supplement B for details [see Barboza et al. (2014)]. The single lacustrine proxy included in the reconstructions is the tiljander_2003_darksum series from Finland [Tiljander et al. (2003)]. We apply a logtransformation on this series in order to dampen the few years that feature very thick varves [Loso (2009)], and to produce a series that is inline with the assumption of normal errors in our statistical models. Figure 1 shows the spatial locations of the 25 proxies.
To increase computational tractability, and to ensure that the heterogenous spatial distribution of the proxies does not bias estimates of the spatial average, we further reduce the 25 proxies into a single series, termed the “reduced proxy,” via a weighted averaging procedure. Intuitively, we seek a reduced proxy series that captures the common signal of globally averaged climate reflected in the shared variability between the proxies. We estimate the averaging weights used to form the reduced proxy using least squares regression, first centering and scaling each of the 25 proxy series over the period 1000–1998. Denoting these scaled proxies as and and the HAD or CRU series as (mean temperature anomalies), we estimate the weights via an ordinary least squares fit to , where is white noise. Since most of the proxies end after 1982, here we fit the model using only the data from 1900 to 1982. The least squares parameter estimates provide a weighted average of proxies that maximizes the explained variance. Denote the reduced proxy as , then
(1) 
The percentage of variation in temperatures that can be explained by the reduced proxy is for the HAD data set and for the CRU data set; note that the is higher for the HAD data set despite all proxies being terrestrial. The proxies are selected on the basis of local correlations, and the higher percentage of explained variation with the HAD data set is indicative of the fact that temperature observations at the locations of the proxies (many of which are coastal) are better at predicting global land and sea temperatures than global landonly temperatures. Note that colinearity is not an issue, as the do not feature strong correlations with one another, and, in addition, our interest lies in the linear combination of rather than the coefficients .
The geophysical distribution of the weights (in percentage of absolute value) is displayed in Tables B.3 and B.4 in Online Supplement B of Barboza et al. (2014). For both HAD and CRU data sets, proxies in the United States are most heavily weighted, followed by the Mongolian composite. The remaining countries have a fairly uniform distribution, with no single country exceeding the 8% level (HAD) or 7% level (CRU). Our selected proxies therefore have broad spatial coverage, inasmuch as possible with the available data. The weights heavily concentrate on the “Tree rings” and “Other” categories, consistent once more with the prevalence of tree ring series in climate reconstructions [e.g., Overpeck et al. (1997), mannbradley (1998), Luterbacher et al. (2004), Moberg et al. (2005), Tingley and Huybers (2013), PAGES 2k Consortium (2013)]. The weight for the single lacustine series, from Tiljander et al. (2003), is less than 8% for both HAD and CRU data sets, indicating that it exerts a limited control on the overall reconstructions. The limited influence of this lacustrine series is of particular importance given the known difficulties in calibrating it, due to the potential of anthropogenic impact on the lake catchment [Tiljander et al. (2003), Mann et al. (2008b)]; we return to this point in Section 4.3.
The modeling approach taken here, based on a weighted average of proxies that pass a local screening condition, does not explicitly consider longrange spatial dependencies, or teleconnections, within the climate system. Another option would be to set the reduced proxy to the leading principal component of the 25 proxies that pass the screening test. Such an approach would extract the dominant common signal shared by the proxies, whereas for the purposes of this analysis we are more interested in retaining the common temperature signal they share. While methods based on principal component or canonical correlation analysis are prevalent in paleoclimatology, both for the reconstruction of spatial patterns and (as here) spatial averages, there is ongoing debate as to the merits of such methods; see Cook, Briffa and Jones (1994), NRC (2006), Wahl and Smerdon (2012), Tingley et al. (2012), Werner, Luterbacher and Smerdon (2013), PAGES 2k Consortium (2013) for discussion.
3.2 Examination of longmemory correlation in the proxy data
While the temperature–proxy relationship is almost universally assumed to be linear [e.g., Luterbacher et al. (2004), Rutherford et al. (2005), Li, Nychka and Ammann (2010), Tingley and Huybers (2010b), Kaufman et al. (2009), McShane and Wyner (2011), Christiansen (2011), Smerdon et al. (2010), and each of the methods in Table 1 of Christiansen, Schmith and Thejll (2009) and discussed in Section 8 of Tingley et al. (2012)], the correlation structure in the error term has not been thoroughly studied. The choice of model for the correlation structure is of particular importance, as its adequacy directly affects the accuracy and precision of the uncertainty quantification associated with the reconstruction. Here we consider models of the form
(2) 
where is a zeromean, unitvariance stationary stochastic process and a constant variance parameter. We fit model (2) using leastsquares over the 1900–1982 interval, using either the HAD or CRU as , and examine the correlation structure of the resulting residuals.
We first explore the correlation structure of using estimates of the spectral density, , of the empirical residuals. If the residuals have longmemory behavior, then the logarithm of the spectrum will feature a negative slope with respect to logfrequency. More specifically, a stationary stochastic process is generally said to have long memory when its autocovariance function decays at the rate for large time lag , where is the longmemory parameter. This behavior is essentially equivalent to requiring that have a singular behavior for small frequencies [see Beran (1994)]. Since for longmemory models, the plot of against for a longmemory model will be approximately a straight line with negative slope . While spectral methods are not generally accepted as a formal way to estimate , save for very simple models, they do offer a useful diagnostic tool to evaluate the longmemory structure in the data [see Beran (1994)].
Based on the regression residuals from equation (2), we compute two widely used estimators of the spectral density: the periodogram and the adaptive multitaper estimator [see Online Supplement A in Barboza et al. (2014) for a brief description for each estimator]. Figure 2 shows both estimators on a log–log scale for the HAD and CRU data sets, respectively. In both cases, the multitaper spectral estimator features a clear negative slope on the log–log scale, indicating possible longmemory behaviors. Results for the periodogram are less striking than the multitaper estimate, but still show a negative slope in log–log space.
To examine more formally the long memory behavior of the residuals, we employ the test developed by Robinson (1995) (Section 3.4 presents results of alternative tests). To introduce the idea of this method briefly, consider a stationary process with spectral density . The may satisfy the power law as for a positive value and some . The socalled Hurst parameter measures the length of the correlation as illustrated by the negative slope of the spectrum in Figure 2. Typical examples that follow this power law include FARIMA and fGn. This fGn is the discretetime stationary Gaussian process that is the firstorder difference process of the socalled fractional Brownian motion (fBm) process evaluated at integer times. The spectrum of the distributional derivative of the fBm process is proportional to . The spectrum of fGn has the same behavior asymptotically for small . Historically, the parameter first made its appearance when fBm was introduced by Kolmogoroff (1940); the name Hurst arose after Mandelbrot proposed that fBm might be a good model to explain the power behavior of a statistic introduced by the hydrologist H. E. Hurst to study yearly levels of the Nile river; see Mandelbrot (1965), Mandelbrot and Van Ness (1968) and the account in Taqqu (2013). More information on fGn can be found in Online Supplement A in Barboza et al. (2014). The FARIMA model depends on a parameter usually denoted by and features a spectral density with the same lowfrequency and longmemory asymptotics as fGn.
(a) HAD data set  (b) CRU data set 
The null hypothesis for the Robinson (1995) test is (no memory), while the alternative hypothesis is (longmemory). The test is based on the asymptotic normality of the semiparametric Gaussian estimate of . Other tests for the memory length are reviewed in Murphy and Izzeldin (2009), who recommend Robinson’s test due to its power properties and its good performance for relatively small samples when combined with bootstrap resampling.
We perform Robinson’s test on the regression residuals in (2), resulting in values of for HAD and for CRU. Both data sets therefore show strong evidence, according to Robinson’s test, in favor of rejecting the null hypothesis of . Note that, the test, while consistent with long memory, does not provide evidence in favor of longmemory correlations over shorter nonzero ones; in the modelcomparison exercises below (Section 4.2), we also consider models which contain short memory, AR(1) errors.
3.3 Examination of longmemory behavior in the temperature anomalies
In the specification of the process level of the hierarchical model, we follow Li, Nychka and Ammann (2010) and model the latent temperature as linear in the external forcings. We apply the following transformations to the forcings, where , and are, respectively, the time series of solar irradiance, volcanism and greenhouse gases:

. Exploratory data analysis indicated that this transformation increases the explanatory power of volcanism. From a physical standpoint, it dampens the effects of very large events, and thus provides a form of regularization given the larger uncertainties associated with the larger values [Li, Nychka and Ammann (2010)].

. Following Hegerl et al. (2007), we use a logtransformation to approximate the radiative forcing due to changes in the equivalent CO concentration.
The resulting processlevel model is
(3) 
where denotes a stationary stochastic process with zero mean and unit variance, and is a constant variance parameter. Li, Nychka and Ammann (2010) employ an AR(2) for the error term, based on an examination of auto and partial autocorrelation functions. However, in a similar situation, Beran (1994) shows that the residuals are appropriately modeled as FARIMA (0, , 0), with Hurst parameter . Benth and ŠaltytėBenth (2005) and Brody, Syroka and Zervos (2002) also provide examples of estimation of longmemory parameters over regression residuals on temperature series for specific locations in Norway and England, respectively, while Huybers and Curry (2006) provides statistical evidence of a powerlaw behavior in the spectrum of surface temperatures. Finally, Imbers et al. (2014) use a longmemory fractionaldifferencing process that is very similar to fGn in terms of its asymptotic longmemory behavior, in order to test the presence of an anthropogenic impact on presentday temperatures.
We repeat the same diagnostic procedure and hypothesis testing as in Section 3.2 to assess the long memory behavior of . We first fit model (3) using the ordinary leastsquares criterion, and find values of 73% for HAD and 66% for CRU, indicating the strong explanatory power of the forcings. Figure 3 plots spectral density estimates in log–log space, for both HAD and CRU, and shows that HAD, but not CRU, exhibits a negative slope.
(a) HAD data set  (b) CRU data set 
The value associated with Robinson’s test is for HAD and for CRU, indicating strong evidence against no memory for the HAD data set, but not for the CRU data set.
As there is value in reconstructions that exclude the forcings (e.g., for the purpose of General Circulation Model assessment), we also consider a reduced form of the processlevel model that exclude the forcings, and models climate variability as a purely stochastic process. Applying Robinson’s test to the CRU and HAD data sets results in values of for both cases, where we can note that the amount of evidence against no memory increases when we exclude the forcings.
3.4 Other tests
We briefly discuss results for several alternatives to Robinson’s test. Beran’s test [see Beran (1992)] evaluates the goodness of fit of a particular stochastic process model (e.g., fGn) to a realization of a time series. Let be a stationary Gaussian process with spectral density , whose realization one observes. When testing for fGn, for instance, if is the spectral density of an fGn process with Hurst parameter , then the null hypothesis for Beran’s test is and the alternative is . Both the Robinson and Beran tests base their test statistics on the Whittle estimator of , which enjoys the desirable property of insensitivity to certain changes of scale [see Online Supplement A in Barboza et al. (2014) for additional technical details].
We performed Beran’s test on six data sets: the four residuals from the HAD and CRU data sets, for both the proxy [equation (2)] and instrumental [equation (3)] equations, and the two HAD and CRU temperature data series themselves with no forcings. To test the presence of memory, we use three distinct memory structures: fractional Gaussian noise, AR(1) and AR(2). The null hypothesis in each test is that the data comes from a spectral density equal to that of the given memory structure. Thus, a nonrejection of the null is not inconsistent with the tested memory structure. For our eighteen Beran’s tests, the corresponding values are shown in Table 2. The results indicate that Beran’s test cannot reject the null in any of the eighteen cases; this is consistent with the presence of memory, but the tests do not point to a preferred memory structure.
Model  fGn  AR(1)  AR(2) 

HADproxy  0.77  0.40  0.76 
CRUproxy  0.91  0.72  0.92 
HADtemp.  0.56  0.58  0.59 
CRUtemp.  0.73  0.33  0.40 
HADtemp. (no forcings)  0.61  0.63  0.67 
CRUtemp. (no forcings)  0.46  0.19  0.47 
Finally, we apply the test proposed by Davies and Harte (1987); see Section A.1.2 for technical details. The fGn is used as the underlying parametric model for this test, and the null and alternative hypotheses are identical to Robinson’s test: (no memory) versus (longmemory). Thus, in contrast with Beran’s test, rejection of the null is evidence against no memory. As in Beran’s test, we use the four residuals from the HAD and CRU data sets, for both the proxy [equation (2)] and instrumental [equation (3)] equations, and the two HAD and CRU series with no forcings. values in Table 3 show that the null can be rejected in three out of four cases when we include forcings within the models, and in the two cases without forcings. In fact, the evidence against no memory increases when we exclude forcings.
No single method employed here is a perfect indicator for the presence or absence of memory in our error processes. Taken together, however, the spectral density estimates and applications of the tests of Robinson (1995), Beran (1992) and Davies and Harte (1987) indicate to us that the possibility of memory, long or short, cannot be ignored in developing models for the residuals or for the HAD and CRU series themselves. In Section 4, we further investigate the memory properties of the residual processes, via Bayesian parameter estimates and reconstruction validation measures.
3.5 Hierarchical Bayesian model with long or shortmemory errors
Given the statistical evidence for long or shortmemory correlation in the empirical residuals from equations (2) and (3), and the implication for fGn or AR model by Beran’s test, we explore the results of modeling the errors using either fGn or AR processes. As the strategy for fitting the hierarchical Bayesian reconstruction is similar in each case, we present details for the more computationally involved fGn error assumption. Comparisons between various modeling choices (long memory vs. short memory vs. no memory; with or without forcings) are given in Section 4.3. A summary of the data and process levels of the hierarchical model is as follows:
(4)  
where and are independent fGn processes with respective parameters and which control the longmemory behavior. We assume these models hold throughout the entire prediction period (1000–1899) and calibration period (1900–1998). Independence between and is a reasonable assumption, as represents the stochastic aspect of the proxies that is not explained by the climate, while is the longmemory aspect of the climate not attributable to the forcings.
The modeling framework [equation (4)] is based on the assumption that the relationship between the proxies and temperatures is invariant through time. While stationarity may be an idealized assumption, we note that our data selection procedure ensures that stationarity is at the very least not an unreasonable assumption, while the short calibration period precludes a more indepth study of possible nonstationarity in the temperature–proxy relationship. Moreover, we note that the modeling framework could be made more realistic by specifying a (possibly independent) error structure for each individual proxy series. We do not pursue these specifics here, but rather focus on exploring the effects of long memory and forcings on the reconstruction.
Following Li, Nychka and Ammann (2010), we define the following prior distributions for the parameters , , , , and :

; ;

, ;

; ;
where is the identity matrix of dimension .
Let denote the vector of unknown temperatures and the vector of instrumental temperatures. Our goal is to infer based on , , , and . The full conditional posterior distributions of and all unknown parameters save and can be derived explicitly, thus allowing for standard Gibbs sampling in the Markov chain Monte Carlo (MCMC) method. We resort to Metropolis–Hasting steps to sample and . The derivation of full conditional distributions can be found in Online Supplement A in Barboza et al. (2014). We implement the MCMC using a number of R packages: MCMCpack [Martin, Quinn and Park (2011)], mvtnorm [Genz and Bretz (2009)], ltsa [McLeod, Yu and Krougly (2007)] and msm [Jackson (2011)].
4 Numerical results
The diagnostic tests in Section 3, while providing no conclusive evidence for the presence of long or short memory, indicate the possibility of certain correlations. In order to further investigate appropriate models for error structures and to assess the benefit of incorporating external forcings in the reconstruction, we compare eight model variants on the basis of their parameter estimates and reconstruction validation metrics: {longlist}[B]
Possible long memory ( and not fixed), with external forcings.
AR(1) errors in (4), with external forcings.
No memory (), with external forcings.
Possible long memory ( and not fixed), no external forcings ().
AR(1) errors in (4), no external forcings ().
No memory (), no external forcings (). The AR(1) model is included, as it features short memory—an intermediate model between assuming fGn and assuming uncorrelated white noise. We refer to Scenarios E and H as having no memory, as they are based on Gaussian white noise errors that are independent and thus have no memory. Scenario B allows for a longmemory model for the proxies while assuming short memory in the temperature residuals, while Scenario C reverses the assumptions of Scenario B.
For reconstructions using both the HAD and CRU instrumental records, we sample 5000 times from the posterior distribution and discard the first 1000 replicates to account for the burnin period. The details of posterior samples are shown in Online Supplement B [see Barboza et al. (2014)]. Here we summarize the results and show a selection of representative plots and focus on reconstructions using the HAD data set.
4.1 Bayesian parameter estimates
We first examine parameter estimates using the HAD data set and including the forcings. Figure 4 shows trace plots and histograms of the and parameters that are responsible for long memory in Scenario A. Visually, the posterior draws quickly stabilize; see Section 4.3 for a formal assessment of convergence for these and other parameters. The histograms of and for the HAD reconstruction clearly indicate that both parameters are significantly greater than 0.5, suggesting that the data are consistent with a longrange correlation model. Figure 5 shows the posterior distribution of and for the CRU reconstruction. The distribution of (memory structure of the proxy residuals) is similar to that arising from the HAD analysis, whereas the posterior distribution of for the CRU analysis is centered on smaller values than for HAD, but still remains significantly greater than 0.5. The larger value of for the HAD data set, which includes the oceans, is in line with intuition, on account of the larger heat capacity of the oceans resulting in a longer timescale response to changes in the forcings.
(a) Traceplots  (b) Histograms 
(a) Traceplots  (b) Histograms 
For Scenarios B, C and D, both HAD and CRU show that all AR(1) parameter estimates are significantly greater than zero, and all longmemory parameter estimates are significantly greater than 0.5 [see Figures B.2, B.3 and B.4 in the Online Supplement B of Barboza et al. (2014)]. For Scenarios F and G, which exclude the forcings, and have, respectively, long and short memory, Bayesian posteriors for the memory parameters provide evidence against models with no memory at higher levels of certainty than for models that include forcings, especially in the CRU case; see Figure B.5 in Online Supplement B [see Barboza et al. (2014)]. These results indicate that while there is a certain amount of memory in the error structures, there is insufficient evidence to select between short or longmemory assumptions. In the subsequent section, we resort to reconstruction validation metrics to compare different models.
Posterior samples for the processlevel regression coefficients (the ) for Scenario A show that the transformed volcanic and greenhouse gas forcing series are meaningful predictors of the temperature evolution for both HAD and CRU, while solar irradiance is less influential (Figures B.8 and B.12). While the forcings are useful predictors of past temperatures, we stress that the reconstructions that exclude the forcings are also of scientific interest. Such reconstructions may not provide the most accurate estimates of past climate fluctuations, but provide necessary test beds for assessing the GCMs used to project future climate, since comparisons between forcingsbased reconstructions and GCM simulations, which are based on the same forcings, would pose circularity issues.
4.2 Validation measures
We provide quantitative assessments of the eight reconstructions using a number of statistical measures: squared bias (squared sample mean of differences between the posterior mean and the observed anomalies); variance (sample variance of the differences used in bias calculation); root mean squared error (RMSE); empirical coverage probabilities (ECP) of the credible intervals at the and levels; Interval Scores (IS) at the and levels; and, since we obtain MCMC samples from the predictive distribution, the Continuous Ranked Probability Score (CRPS). The ECP measures the accuracy of the uncertainty quantification and values closer to nominal level are more desirable, while the IS and CRPS provide more nuanced assessments of the posterior predictive distributions, rewarding both calibration and sharpness simultaneously; details of these scoring rules are available in Gneiting and Raftery (2007), Gneiting, Balabdaoui and Raftery (2007), Gschlößl and Czado (2007) and Online Supplement A.3 in Barboza et al. (2014). For convenience, we report the negative IS and CRPS so that smaller values indicate higher quality predictions.
Scenarios  Sq. bias  Variance  RMSE  ECP  ECP  I  IS  CRPS  

HAD  A  0.016  0.012  0.168  92.9  74.7  0.062  0.178  0.208 
B  0.017  0.013  0.171  92.9  74.7  0.062  0.179  0.205  
C  0.015  0.011  0.160  90.9  72.7  0.064  0.179  0.212  
D  0.015  0.011  0.162  90.9  74.7  0.063  0.176  0.209  
E  0.014  0.010  0.154  90.9  69.7  0.060  0.171  0.195  
F  0.055  0.072  0.356  99.0  84.8  0.110  0.323  0.229  
G  0.081  0.071  0.390  94.9  75.8  0.118  0.389  0.259  
H  0.113  0.059  0.415  82.8  59.6  0.168  0.511  0.304  
CRU  A  0.032  0.025  0.238  91.9  73.7  0.084  0.251  0.245 
B  0.031  0.025  0.235  91.9  75.8  0.081  0.245  0.237  
C  0.033  0.024  0.238  91.9  71.7  0.090  0.258  0.252  
D  0.032  0.023  0.234  90.9  70.7  0.087  0.255  0.245  
E  0.031  0.024  0.235  91.9  73.7  0.085  0.250  0.242  
F  0.089  0.097  0.432  97.0  78.8  0.131  0.416  0.274  
G  0.120  0.095  0.464  90.9  75.8  0.150  0.482  0.303  
H  0.148  0.080  0.477  84.8  62.6  0.206  0.570  0.335 
[]tt4RMSE: Root Mean Square Error; : Empirical Coverage Probability at confidence level; : Interval Score at confidence level; CPRS: Continuous Ranked Probability Score. \tabnotetext[]tt4 HAD and CRU refer to the two instrumental data sets, with HAD including the oceans.
Table 4 summarizes the quantitative assessments of the reconstructions for both the HAD and CRU data sets. The benefit of the external forcings are readily apparent [cf. Li, Nychka and Ammann (2010)], as their inclusion substantially reduces the squared bias, variance and, consequently, the RMSE, as well as the IS and CRPS (compare Scenario A to F, Scenario D to G, and Scenario E to H). This corroborates the fact that the posterior distributions of the coefficients for both the volcanic and green house gas forcing series are significant. Moreover, the widths of the 95% credible intervals are likewise narrower when the external forcings are included (see Figure 8, below, and Figure B.16 in the Online Supplement).
When external forcings are included in the reconstruction, the squared biases, variances and RMSEs are generally similar across the different error models for each of the two data sets. For the HAD data set, and among reconstructions that include forcings, Scenarios A and B are optimal in terms of ECP; Scenario E in terms of CRPS; and there are no appreciable differences in IS. Note that Scenario E exhibits the worst ECP, indicating an underestimation of uncertainty compared to Scenarios A and B. This is consistent with the rejection of no memory models in our tests in Section 3. For the CRU data set, Scenario B is optimal in terms of ECP and CRPS, and again there is no appreciable difference in terms of IS. Based on these validation measures, while there continues to be support for memory models, there is no clear indication of a single, best model for the error structures among the reconstructions that include forcings, with Scenarios A and B featuring comparable performance metrics. Indeed, tests for selecting between long and shortmemory models for climate time series are often inconclusive [e.g., Percival, Overland and Mofjeld (2001)].
When forcings are not included, the greater variability of validation metrics across the scenarios allows for more meaningful ranking of the error correlation assumptions. For both data sets, Scenario F is optimal in terms of squared bias, RMSE, IS and CRPS. For the HAD data set, Scenario G is optimal in terms of ECP at the 95% level but is equally distant from the nominal level as Scenario F at the 80% level, while for the CRU data set, ECP favors Scenario F. In general, the results indicate that when forcings are not included the longmemory models play an important role in capturing the correlation structure in proxies and temperature and should be employed in the hierarchical model. As discussed in Li, Nychka and Ammann (2010), reconstructions are improved when information is included at a broad range of frequency scales. In the absence of forcings, which feature longrange correlations and lowfrequency behavior, the inclusion of more highly structured noise processes leads to marked improvements in the reconstructions.
4.3 Temperature reconstruction results
According to validation measures in Table 4, the reconstruction scenarios that include forcings are similar to one another. Here we focus on Scenario B due to slightly better validation measures for both HAD and CRU data sets. Figure 6 shows the Scenario B temperature reconstruction together with 95% pointwise credible intervals, using the HAD data set. The reconstruction shows a slight downward trend during the period 1000–1899 [cf. Kaufman et al. (2009)], and no maxima in the posterior distributions exceed the levels observed after approximately 1950. The reconstruction for the CRU data set (see Figure B.14) is qualitatively similar, but features higher variance due to the more variable CRU temperatures.
In order to evaluate our reconstruction, we use 1900–1998 as an insample validation period. Due to the limited number of available observations and the necessity of inferring the memory parameters, outofsample validation was not feasible. Figure 7 shows the posterior mean and 95% pointwise credible intervals for predictions using the HAD data in Scenarios A, B, E, F and H, as well as the actual HAD observations. The scenarios that include forcings (A, B, E) result in reconstructions that are qualitatively similar to one another and feature good qualitative agreement with the observations, with Scenario E exhibiting slightly narrower credible intervals. Reconstructions resulting from scenarios that exclude the forcings (F and H) feature greater divergence from the observations—particularly for Scenario H, which models the error structure as white noise. Results are similar for the CRU data set (see Figure B.15). Note that the reduced variability of the posterior mean as compared with the observations is akin to the predictions from a linear regression being less variable than the observations. A key advantage of a Bayesian analysis, such as that used here, is that, provided the processlevel model assumptions are reasonable, the temporal variability of individual posterior draws will be similar to that of the actual climate, even while variability of the mean across them is attenuated [see Figure 2 of Tingley and Huybers (2010b) for further discussion]. Repeating the reconstructions with the single lacustrine record excluded from the reduced proxy leads to similar results; see Figure B.17.
4.4 MCMC diagnostics
To establish convergence of the MCMC samples, we examine trace plots (Figures B.6–B.13) and calculate the potential scale reduction factor [PSRF; Gelman and Rubin (1992)] and its multivariate version [Brooks and Gelman (1998)]; see Brooks and Roberts (1997) and Cowles and Carlin (1996) for further details. We present diagnostic results for Scenario A, as it represents the most complex model for estimation. If the PSRF is close to unity for all parameters, then the Markov chain simulation is close to its stationary distribution, while a large PSRF indicates that the chain has not converged [Gelman and Rubin (1992)]. Brooks and Gelman (1998) provide a generalization that allows for the computation of a single PSRF for all model parameters.
Mul.  

HAD PSRF  1.01  1.01  1.01  1.01  1.01  1.00  1.01  1.00  1.01  1.01  1.02 
UB  1.02  1.04  1.01  1.03  1.02  1.00  1.02  1.01  1.01  1.03  – 
CRU PSRF  1.00  1.01  1.03  1.04  1.01  1.01  1.00  1.01  1.00  1.00  1.05 
UB  1.00  1.03  1.09  1.11  1.01  1.02  1.01  1.01  1.00  1.01  – 
For both the HAD and CRU data sets, we run five MCMC simulations, each of length 5000, and discard the first 1000 samples to allow the chain to burn in. We compute PSRFs for the scalar parameters of the model (, , , , , , , ) and the multivariate PSRF, along with their upper 95% confidence bounds, using the coda Rpackage [Plummer et al. (2006)]. Results in Table 5 show that all the individual PSRFs are relatively close to unity, indicating their successful convergence to the stationary distribution. The multivariate PSRF likewise indicates convergence.
5 Comparison with other works
5.1 Comparison with previous reconstructions
We compare our reconstructions to those reported in Mann et al. (2008a), as both use similar proxy and temperature data sets. Mann et al. (2008a) assume no memory in the error processes, do not include the external forcings, and present reconstructions, along with uncertainty bands, based on two regression approaches: composite plus scale (CPS) and errors in variables (EIV). The CPS approach computes a weighted average of the proxy data, and then calibrates this weighted average by matching its mean and variance to those of the instrumental temperature data during their overlap period. The EIV regression approach allows for errors in both the dependent and independent variables, and we refer to Mann et al. (2008a, 2008b) for details. The EIV and CPS reconstructions, and their associated uncertainty estimates, are available online^{5}^{5}5http://www.ncdc.noaa.gov/paleo/pubs/mann2008/mann2008.html. as decadally smoothed time series, as Mann et al. (2008a) focuses on lowfrequency climate variability. In contrast, the reconstructions we present here are available at annual temporal resolution, with no smoothing. In comparisons, we show the posterior mean and uncertainty of our reconstructions at annual resolution, and additionally include the posterior mean that results from first smoothing each posterior draw with a Butterworth filter^{6}^{6}6Our calculations are based on the Matlab code associated with Mann et al. (2008a), posted online at http://www.ncdc.noaa.gov/paleo/pubs/mann2008/mann2008.html. We smooth using the filtfilt command in the R package “signal.” with cutoff frequency equal to 0.1 cyclesyear.
Figure 8 compares our reconstructions using the HAD data, and under Scenarios A, B, E, F and H, to those from Mann et al. (2008a). In all cases, and especially when including the forcings, our reconstructions are generally cooler than both the EIV and CPS reconstructions from Mann et al. (2008a), particularly during the 1000–1400 interval, and feature a smaller amplitude of preinstrumental temperature variability. We are not the first to report a lower variability than Mann et al. (2008a)—for example, PAGES 2k Consortium (2013) report a change in 30 year average temperatures between 1000 AD and the 1800s of about C, compared with about C for Mann et al. (2008a); see Figure 4 of PAGES 2k Consortium (2013).
The model settings of Mann et al. (2008a) are most similar to our Scenario H, which includes neither the forcings nor the longmemory processes. Indeed, the EIV predictions from Mann et al. (2008a) are visually most similar to smoothed Scenario H results, and 88.4% of the EIV predictions from Mann et al. (2008a) fall within the 95% pointwise credible intervals for the smoothed Scenario H results. Results are similar when using the CRU data set (Figure B.16).
To facilitate numerical comparisons with the Mann et al. (2008a) reconstruction, we recalculate validation metrics for Scenario H after first smoothing each posterior draw; results are shown in Table 6 for the 20th century validation interval. The main difference between our smoothed Scenario H results and the Mann et al. (2008a) results is in terms of squared bias, with the Mann et al. (2008a) reconstruction featuring biases that are about an order of magnitude smaller and variances that are about 1.5–2 times larger. The net result is that the Mann et al. (2008a) reconstructions feature smaller RMSE than our smoothed Scenario H, on par with results from our annually resolved Scenarios A, B, C, D and E. As measured by the ECP, the uncertainties for the Mann et al. (2008a) reconstructions are too wide, in the sense that the empirical coverage rate is greater than the nominal rate. The uncertainties for our smoothed Scenario H are smaller than that in Mann et al. (2008a), but due to the relatively large bias, the ECPs appear to be too low compared to their nominal value. The Interval Scores for the smoothed Mann et al. (2008a) reconstructions are much better than those for our smoothed Scenario H and, like the RMSE, are similar to those for our annually resolved Scenarios A, B, C, D and E which carry small squared bias by including the forcings (see Table 4).
Scenarios  Sq. bias  Variance  RMSE  ECP  ECP  IS  IS  

HAD  H (smoothed)  0.100  0.012  0.335  33.3  0.46  0.71  
CPS  0.009  0.024  0.183  96.9  0.06  0.16  
EIV  0.003  0.022  0.157  99.0  0.06  0.23  
CRU  H (smoothed)  0.121  0.016  0.371  36.4  0.45  0.73  
CPS  0.017  0.025  0.207  99.0  0.07  0.25  
EIV  0.006  0.021  0.163  98.0  0.07  0.17 
tt6The statistics for EIV and CPS reconstructions are calculated using the estimated standard deviations associated with Mann et al. (2008a). They are posted as “2sigma uncertainties” (), hence, the formula for their 95% confidence bands is , where is their predicted temperature mean.
We caution against drawing substantive conclusions from the comparison of the validation and scoring metrics between the Mann et al. (2008a) results and the smoothed Scenario H, as numerous lines of evidence indicate that Scenario H is the least appropriate of the eight scenarios explored here: validation metrics and scores (Table 4) are generally worst for Scenario H; the inclusion of the forcings is motivated by the scientific understanding of their connection with temperatures; and the inclusion of the longmemory processes in the absence of forcings is driven by the structure of the data. Indeed, we view Scenario H as a misspecified model, and the high squared bias and associated inadequacies of the ECPs are therefore to be expected. Perhaps the most telling conclusion to be drawn from the numerical comparisons is that our annually resolved Scenarios A, B, C and D, which include the forcings as well as short and/or longmemory processes, are comparable in terms of RMSE and Interval Scores to the decadally resolved Mann et al. (2008a) results while featuring ECPs which are closer to their nominal values.
Finally, we note that the proxy selection and modeling treatments do differ between our Scenario H and the reconstructions in Mann et al. (2008a) so that the comparison remains imperfect. In particular, we note that the Mann et al. (2008a) reconstructions include proxies with decadal resolution, whereas here we focus on proxies with annual resolution. Indeed, the CPS reconstruction is performed after smoothing all proxies to a common decadal resolution, while the EIV reconstruction is based on a “hybrid” frequency approach that involves separate calibrations to infer climate on interdecadal (periods longer than 20 years) and interannual (periods shorter than 20 years) timescales [Mann et al. (2007, 2008b)]. Due to the differing methods and the focus on lower frequency variability in Mann et al. (2008a), the differing validation metrics between our smoothed Scenario H and those for the Mann et al. (2008a) reconstructions are not surprising.
5.2 Transient climate response
The Fourth Assessment Report of the IPCC [see page 723 in Hegerl et al. (2007)] refers to the “transient climate response” (TCR) as the “global mean temperature change that is realized at the time of CO doubling TCR is therefore indicative of the temperature trend associated with external forcing, and can be constrained by an observable quantity, the observed warming trend that is attributable to greenhouse gas (GHG) forcing.” In our model, the transient response to a doubling of GHG is functionally related to the parameter , and the resulting estimates of TCR are based on the instrumental temperature record since 1900 and proxy and forcing information over the past millennium. We believe that our Bayesian approach to computing the transient response to GHG forcing from both instrumental and proxy observations, without recourse to global climate models, is new to the field.
Taking into account the transformations applied to the CO series, we define TCR in terms of as
where is the standard deviation of the logarithm of the GHG series and is computed over the entire period 1000–1998. An important advantage of Bayesian estimation is the possibility of obtaining a sample estimate of the marginal posterior distribution of given the data, from which we can compute a nonparametric estimator of the probability density function for TCR that accounts for the uncertainties in all other parameters in the model.
We present results of TCR estimates using the global land and marine HAD data set, for the five scenarios that include the forcings: Scenarios A, B, C, D and E (Figure 9). There is substantial variability between the TCR estimates from the five scenarios. TCR estimates are the lowest and most sharply peaked for the memoryfree Scenario E, with a median around and an approximate 95% credible interval of C. The TCR distributions become progressively broader as more memory is included, in Scenarios D, then in B and C, and finally in the fully longmemory Scenario A which features the broadest 95% credible interval of C. A quantitative explanation for this increasing uncertainty behavior can be found by inspecting the formula for the covariance matrix of the posterior distribution of the vector given : from formula (A.4) in the Online Supplement A [see Barboza et al. (2014)], one sees that is the inverse of a matrix which is affine in , that is, affine in the inverse of the covariance matrix for the noise model being used in each scenario. It is known [see, e.g., Palma and Bondon (2003)^{7}^{7}7In this paper, the authors provide the estimate for the top eigenvalue of the covariance matrix of a vector of contiguous terms of a stationary sequence whose covariance matrix satisfies , which is the case for our fGn sequence. Thus, indeed, is roughly increasing in for all . Palma and Bondon (2003) state this result in the case of the ARFIMA process, in Example 2 on pages 99–100, but an inspection of their proof shows that the result holds for all satisfying the above asymptotics.] that the magnitude (e.g., the operator norm) of increases with memory length; this and the formula for can explain the increasing behavior we observe in Figure 9.
On the other hand, the progression of posterior medians for the TCRs is not monotone. Scenario D [AR(1) errors] features the largest median TCR value: C, followed by Scenario C with short memory in the proxy model: C. The two scenarios with long memory in the proxy model, A and B, have lower median TCR values, respectively, C and C, with Scenario E (no memory) having the lowest median of C, as reported above. The formula for the posterior distribution of the vector given is again helpful: the posterior mean is the product of the increasing , as discussed above, and of a matrix which is affine in , thus with presumably decreasing magnitude with respect to memory length; the competition between these two effects could induce nonmonotonicity with respect to memory length.
To arrive at a best estimate of TCR, we mix with equal weights the posterior estimates from Scenarios A and B, yielding a median TCR of about C with a combined credible interval of about C. Scenarios A and B feature superior validation metrics as compared with Scenarios C and D, while Scenario E presumably underreports TCR uncertainty because it is based on a model that lacks memory. We therefore focus on TCR estimates derived from an equally weighted mixture of estimates from Scenarios A and B. It is not possible to select between these two scenarios using model diagnostics and, as both include memory, our estimates are conservative with respect to uncertainty.
It is instructive to compare our five TCR distributions reported here with the consensus (expert assessment) recently released in the IPCC’s Fifth Assessment Report [Bindoff et al. (2013), Collins et al. (2013)], where TRC is reported as “likely” within the interval C and “extremely unlikely” to exceed C. Using the IPCC’s definitions/guidance on uncertainty language, these expert assessments can be interpreted as meaning that the probability of the estimated TCR from one of our scenarios falling within the C interval should exceed 0.66, while the probability that the estimated TCR exceeds C should not exceed . For all scenarios reported in Figure 9, the posterior probability that the TCR exceeds C is in all cases lower than . This exceedance probability is essentially zero for Scenario E, which features the narrowest TCR distribution: Scenario E presumably underestimates uncertainty by using no memory for modeling errors. As for falling in the interval C with probability around 0.66, Scenario E does satisfy this condition; Scenario B nearly does; for Scenario A the probability is closer to 0.5; but Scenarios C and D, with their significantly higher median values, fail the condition by some margin. Our best estimate, derived from mixing Scenarios A and B, meets the TCR upper bound condition: the probability that it exceeds C is about . It falls slightly short of meeting the confidence interval condition: the probability that it falls within the interval C is about .
All of our reported TCRs are on the high side compared to the latest IPCC consensus, and as compared with several specific recent studies which have arrived at TCR estimates by combining information from models and the instrumental temperature record. Gillett et al. (2012) produce a TCR estimate of 1.3–1.8C using the global HAD data set and a single global climate model, but note that this TCR estimate may be unrealistically narrow as it results from a single climate model. A more recent study [Gillett et al. (2013)] that combines information from an ensemble of models and the instrumental record results in a wider range of TCR estimates, 0.9–2.3C, featuring greater overlap with our results. Otto et al. (2013) use global, decadal averages of the HAD data set over the 1970–2009 to arrive at a databased TCR estimate in the range of 0.7–2.5C, but caution against strong conclusions based on a such a short time interval. A particularly high estimate of TCR, of at least 2.5–3.6C, is reported by Tung, Zhou and Camp (2008), based on an analysis of the 11year solar cycle.
Hence, both the specific modeldata fusion studies discussed in the previous paragraph and in Bindoff et al. (2013), as well as the synthesis provided by the IPCC Fifth Assessment Report, generally feature broader uncertainties and are peaked at lower values as compared to our posterior estimates of TCR. Indeed, only one of the estimated TCR distributions shown in Figure 10.20 of Bindoff et al. (2013) is peaked at a value greater than C, while the high estimate of Tung, Zhou and Camp (2008) is explained as resulting from solar forcing having a different mechanistic effect on climate [Bindoff et al. (2013)]. Interestingly, the single plotted TCR distribution peaked at greater than C is that of Harris et al. (2013), which estimates TCR using a Bayesian approach that combine information from GCMs and recently observed temperature changes. A possible cause for the narrower uncertainties and higher TCR values estimated here is the more extensive use of data, in terms of both variety (instrumental temperatures, proxies, and estimates of CO, volcanic and solar forcings) and duration (observations over the last millennium).
6 Conclusions and discussion
We use a comprehensive multiproxy data set to produce new reconstructions of NH temperature anomaly time series back to 1000 AD and systematically evaluate the effects of including or excluding external drivers of climate variability, and of assuming the error processes feature long, short or no memory, by considering eight modeling scenarios. Hierarchical Bayesian models are used throughout as they provide a natural framework for integrating the different information sources—proxy and instrumental temperatures observations, and time series of solar, greenhouse gas and volcanic forcings. Bayesian inference additionally permits for estimation of all unknown quantities, including past temperatures, and facilitates uncertainty propagation.
While the possibility of long memory was suggested by exploratory data analysis, and the significance of longmemory parameters verified by Bayesian estimation, model diagnostics indicated that short and longmemory models yield comparable results provided the climate forcings are incorporated into the reconstructions. The inclusion of the external forcings is motivated from physical principles and the conclusions of Li, Nychka and Ammann (2010), and additionally allows for estimation of the transient climate response. While our TCR estimates are near the upper bound of the expertderived “extremely likely” interval provided in the IPCC Fifth Assessment Report [Bindoff et al. (2013)], they do not violate this uncertainty consensus, and we note that our estimate is based on both the instrumental and paleoclimate records, and does not rely on GCMs.
If the forcings are excluded from the reconstruction, as is necessary for reconstructions to be suitable for GCM assessment exercises, the longmemory processes substantially improve the quality of the reconstructions. The scenario with neither forcings nor memory is similar to the benchmark reconstruction of Mann et al. (2008a), though we note that there remain differences in both method and data usage. Our reconstructions generally indicate cooler temperatures than those of Mann et al. (2008a), particularly before the year 1400.
The basic framework presented in this paper can be extended in several directions, and we anticipate that doing so will produce further insights into the climate of the late Holocene. An obvious extension is to incorporate a spatial element, by combining the model used here with the space–time model in Tingley and Huybers (2010a). Doing so would require generalizing the reducedproxy framework and instead specifying a separate longmemory error model for each proxy time series, or perhaps a common model for each proxy type [cf. Tingley and Huybers (2010a)]. Such an implementation would pose technical challenges, as the estimation of the longmemory parameters is the most numerically demanding component of the analysis. Prior scientific understanding of the mechanisms by which the proxies record variations in the climate may be helpful in selecting appropriate temporal correlation models for the residuals, and can potentially be used to simplify calculations. Such a computationally demanding generalization may be a more scientifically defensible use of the proxies and may allow for further insights into the proxy–climate relationship.
Acknowledgments
The authors thank the Editor, three anonymous referees, and Michael Stein for their helpful comments and suggestions.
Supplement to: “Reconstructing past temperatures from natural proxies and estimated climate forcings using short and longmemory models” \slink[doi]10.1214/14AOAS785SUPP \sdatatype.pdf \sfilename10.1214/14AOAS785_supp.pdf \sdescriptionWe provide a background on longmemory models, the multitaper estimator and scoring rules together with some calculations of our model’s posterior distributions. Finally, we include additional plots and tables.
References
 Ammann et al. (2007) {barticle}[author] \bauthor\bsnmAmmann, \bfnmCaspar M.\binitsC. M., \bauthor\bsnmJoos, \bfnmFortunat\binitsF., \bauthor\bsnmSchimel, \bfnmDavid S.\binitsD. S., \bauthor\bsnmOttoBliesner, \bfnmBette L.\binitsB. L. \AND\bauthor\bsnmTomas, \bfnmRobert A.\binitsR. A. (\byear2007). \btitleSolar influence on climate during the past millennium: Results from transient simulations with the NCAR climate system model. \bjournalProc. Natl. Acad. Sci. USA \bvolume104 \bpages3713–3718. \bptokimsref\endbibitem
 Barboza et al. (2014) {bmisc}[author] \bauthor\bsnmBarboza, \binitsL., \bauthor\bsnmLi, \binitsB., \bauthor\bsnmTingley, \binitsM. P. and \bauthor\bsnmViens, \binitsF. G. (\byear2014). \bhowpublishedSupplement to “Reconstructing past temperatures from natural proxies and estimated climate forcings using short and longmemory models.” DOI:\doiurl10.1214/14AOAS785SUPP. \bptokimsref \endbibitem\bptokimsref\endbibitem
 Benth and ŠaltytėBenth (2005) {barticle}[author] \bauthor\bsnmBenth, \bfnmFred Espen\binitsF. E. \AND\bauthor\bsnmŠaltytėBenth, \bfnmJūratė\binitsJ. (\byear2005). \btitleStochastic modelling of temperature variations with a view towards weather derivatives. \bjournalAppl. Math. Finance \bvolume12 \bpages53–85. \bptokimsref\endbibitem
 Beran (1992) {barticle}[mr] \bauthor\bsnmBeran, \bfnmJan\binitsJ. (\byear1992). \btitleA goodnessoffit test for time series with long range dependence. \bjournalJ. Roy. Statist. Soc. Ser. B \bvolume54 \bpages749–760. \bidissn=00359246, mr=1185220 \bptokimsref\endbibitem
 Beran (1994) {bbook}[mr] \bauthor\bsnmBeran, \bfnmJan\binitsJ. (\byear1994). \btitleStatistics for LongMemory Processes. \bseriesMonographs on Statistics and Applied Probability \bvolume61. \bpublisherChapman & Hall, \blocationNew York. \bidmr=1304490 \bptokimsref\endbibitem
 Berliner, Wikle and Cressie (2000) {barticle}[author] \bauthor\bsnmBerliner, \bfnmL Mark\binitsL. M., \bauthor\bsnmWikle, \bfnmChristopher K.\binitsC. K. \AND\bauthor\bsnmCressie, \bfnmNoel\binitsN. (\byear2000). \btitleLonglead prediction of pacific SSTs via Bayesian dynamic modeling. \bjournalJ. Climate \bvolume13 \bpages3953–3968. \bptokimsref\endbibitem
 Bindoff et al. (2013) {bincollection}[author] \bauthor\bsnmBindoff, \bfnmN. L.\binitsN. L., \bauthor\bsnmStott, \bfnmP. A.\binitsP. A., \bauthor\bsnmAchutaRao, \bfnmK. M.\binitsK. M., \bauthor\bsnmAllen, \bfnmM. R.\binitsM. R., \bauthor\bsnmGillett, \bfnmN.\binitsN., \bauthor\bsnmGutzler, \bfnmD.\binitsD., \bauthor\bsnmHansingo, \bfnmK.\binitsK., \bauthor\bsnmHegerl, \bfnmG.\binitsG., \bauthor\bsnmHu, \bfnmY.\binitsY., \bauthor\bsnmJain, \bfnmS.\binitsS., \bauthor\bsnmMokhov, \bfnmI. I.\binitsI. I., \bauthor\bsnmOverland, \bfnmJ.\binitsJ., \bauthor\bsnmPerlwitz, \bfnmJ.\binitsJ., \bauthor\bsnmSebbari, \bfnmR.\binitsR. \AND\bauthor\bsnmZhang, \bfnmX.\binitsX. (\byear2013). \btitleDetection and attribution of climate change: From global to regional. In \bbooktitleClimate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (\beditor\bfnmT. F.\binitsT. F. \bsnmStocker, \beditor\bfnmD.\binitsD. \bsnmQin, \beditor\bfnmG. K.\binitsG. K. \bsnmPlattner, \beditor\bfnmM.\binitsM. \bsnmTignor, \beditor\bfnmS. K.\binitsS. K. \bsnmAllen, \beditor\bfnmJ.\binitsJ. \bsnmBoschung, \beditor\bfnmA.\binitsA. \bsnmNauels, \beditor\bfnmY.\binitsY. \bsnmXia, \beditor\bfnmV.\binitsV. \bsnmBex \AND\beditor\bfnmP. M.\binitsP. M. \bsnmMidgley, eds.). \bpublisherCambridge Univ. Press, \blocationCambridge. \bptokimsref\endbibitem
 Briffa et al. (1998) {barticle}[author] \bauthor\bsnmBriffa, \bfnmK R.\binitsK. R., \bauthor\bsnmJones, \bfnmP D.\binitsP. D., \bauthor\bsnmSchweingruber, \bfnmF H.\binitsF. H. \AND\bauthor\bsnmOsborn, \bfnmT J.\binitsT. J. (\byear1998). \btitleInfluence of volcanic eruptions on Northern Hemisphere summer temperature over the past 600 years. \bjournalNature \bvolume393 \bpages450–455. \biddoi=10.1038/30943 \bptokimsref\endbibitem
 Brody, Syroka and Zervos (2002) {barticle}[mr] \bauthor\bsnmBrody, \bfnmDorje C.\binitsD. C., \bauthor\bsnmSyroka, \bfnmJoanna\binitsJ. \AND\bauthor\bsnmZervos, \bfnmMihail\binitsM. (\byear2002). \btitleDynamical pricing of weather derivatives. \bjournalQuant. Finance \bvolume2 \bpages189–198. \biddoi=10.1088/14697688/2/3/302, issn=14697688, mr=1913206 \bptokimsref\endbibitem
 Brohan et al. (2006) {barticle}[author] \bauthor\bsnmBrohan, \bfnmP.\binitsP., \bauthor\bsnmKennedy, \bfnmJ. J.\binitsJ. J., \bauthor\bsnmHarris, \bfnmI.\binitsI., \bauthor\bsnmTett, \bfnmS. F. B.\binitsS. F. B. \AND\bauthor\bsnmJones, \bfnmP. D.\binitsP. D. (\byear2006). \btitleUncertainty estimates in regional and global observed temperature changes: A new data set from 1850. \bjournalJ. Geophys. Res. \bvolume111 \bpagesD12106+. \biddoi=10.1029/2005JD006548 \bptokimsref\endbibitem
 Brooks and Gelman (1998) {barticle}[mr] \bauthor\bsnmBrooks, \bfnmStephen P.\binitsS. P. \AND\bauthor\bsnmGelman, \bfnmAndrew\binitsA. (\byear1998). \btitleGeneral methods for monitoring convergence of iterative simulations. \bjournalJ. Comput. Graph. Statist. \bvolume7 \bpages434–455. \biddoi=10.2307/1390675, issn=10618600, mr=1665662 \bptokimsref\endbibitem
 Brooks and Roberts (1997) {barticle}[author] \bauthor\bsnmBrooks, \bfnmStephen P.\binitsS. P. \AND\bauthor\bsnmRoberts, \bfnmGareth O.\binitsG. O. (\byear1997). \btitleAssessing convergence of Markov chain Monte Carlo algorithms. \bjournalStatist. Comput. \bvolume8 \bpages319–335. \bptokimsref\endbibitem
 Christiansen (2011) {barticle}[author] \bauthor\bsnmChristiansen, \bfnmBo\binitsB. (\byear2011). \btitleReconstructing the NH mean temperature: Can underestimation of trends and variability be avoided? \bjournalJ. Climate \bvolume24 \bpages674–692. \bptokimsref\endbibitem
 Christiansen, Schmith and Thejll (2009) {barticle}[author] \bauthor\bsnmChristiansen, \bfnmB.\binitsB., \bauthor\bsnmSchmith, \bfnmT.\binitsT. \AND\bauthor\bsnmThejll, \bfnmP.\binitsP. (\byear2009). \btitleA surrogate ensemble study of climate reconstruction methods: Stochasticity and robustness. \bjournalJ. Climate \bvolume22 \bpages951–976. \bptokimsref\endbibitem
 Chronopoulou and Viens (2010) {bincollection}[mr] \bauthor\bsnmChronopoulou, \bfnmAlexandra\binitsA. \AND\bauthor\bsnmViens, \bfnmFrederi G.\binitsF. G. (\byear2010). \btitleHurst index estimation for selfsimilar processes with longmemory. In \bbooktitleRecent Development in Stochastic Dynamics and Stochastic Analysis (\beditor\bfnmJinqiao\binitsJ. \bsnmDuan, \beditor\bfnmShunlong\binitsS. \bsnmLuo \AND\beditor\bfnmCaishi\binitsC. \bsnmWang, eds.). \bseriesInterdiscip. Math. Sci. \bvolume8 \bpages91–117. \bpublisherWorld Sci. Publ., \blocationHackensack, NJ. \biddoi=10.1142/9789814277266_0006, mr=2807815 \bptnotecheck year \bptokimsref\endbibitem
 Chronopoulou and Viens (2012) {barticle}[mr] \bauthor\bsnmChronopoulou, \bfnmAlexandra\binitsA. \AND\bauthor\bsnmViens, \bfnmFrederi G.\binitsF. G. (\byear2012). \btitleEstimation and pricing under longmemory stochastic volatility. \bjournalAnn. Finance \bvolume8 \bpages379–403. \biddoi=10.1007/s1043601001564, issn=16142446, mr=2922802 \bptokimsref\endbibitem
 Chronopoulou, Viens and Tudor (2009) {barticle}[mr] \bauthor\bsnmChronopoulou, \bfnmAlexandra\binitsA., \bauthor\bsnmViens, \bfnmFrederi G.\binitsF. G. \AND\bauthor\bsnmTudor, \bfnmCiprian A.\binitsC. A. (\byear2009). \btitleVariations and Hurst index estimation for a Rosenblatt process using longer filters. \bjournalElectron. J. Stat. \bvolume3 \bpages1393–1435. \biddoi=10.1214/09EJS423, issn=19357524, mr=2578831 \bptokimsref\endbibitem
 Collins et al. (2013) {bincollection}[author] \bauthor\bsnmCollins, \bfnmM. R.\binitsM. R., \bauthor\bsnmKnutti, \bfnmR.\binitsR., \bauthor\bsnmArblaster, \bfnmJ.\binitsJ., \bauthor\bsnmDufresne, \bfnmJ.–L.\binitsJ.L., \bauthor\bsnmFichefet, \bfnmT.\binitsT., \bauthor\bsnmFriedlingstein, \bfnmP.\binitsP., \bauthor\bsnmGao, \bfnmX.\binitsX., \bauthor\bsnmGutowski, \bfnmW. J.\binitsW. J., \bauthor\bsnmJohns, \bfnmT.\binitsT., \bauthor\bsnmKrinner, \bfnmG.\binitsG., \bauthor\bsnmShongwe, \bfnmM.\binitsM., \bauthor\bsnmTebaldi, \bfnmC.\binitsC., \bauthor\bsnmWeaver, \bfnmA. J.\binitsA. J. \AND\bauthor\bsnmWehner, \bfnmM.\binitsM. (\byear2013). \btitleLongterm climate change: Projections, commitments and irreversibility. In \bbooktitleClimate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (\beditor\bfnmT. F.\binitsT. F. \bsnmStocker, \beditor\bfnmD.\binitsD. \bsnmQin, \beditor\bfnmG. K.\binitsG. K. \bsnmPlattner, \beditor\bfnmM.\binitsM. \bsnmTignor, \beditor\bfnmS. K.\binitsS. K. \bsnmAllen, \beditor\bfnmJ.\binitsJ. \bsnmBoschung, \beditor\bfnmA.\binitsA. \bsnmNauels, \beditor\bfnmY.\binitsY. \bsnmXia, \beditor\bfnmV.\binitsV. \bsnmBex \AND\beditor\bfnmP. M.\binitsP. M. \bsnmMidgley, eds.). \bpublisherCambridge Univ. Press, \blocationCambridge. \bptokimsref\endbibitem
 Cook, Briffa and Jones (1994) {barticle}[author] \bauthor\bsnmCook, \bfnmEdward R.\binitsE. R., \bauthor\bsnmBriffa, \bfnmKeith R.\binitsK. R. \AND\bauthor\bsnmJones, \bfnmPhilip D.\binitsP. D. (\byear1994). \btitleSpatial regression methods in dendroclimatology: A review and comparison of two techniques. \bjournalInt. J. Climatol. \bvolume14 \bpages379–402. \bptokimsref\endbibitem
 Cowles and Carlin (1996) {barticle}[mr] \bauthor\bsnmCowles, \bfnmMary Kathryn\binitsM. K. \AND\bauthor\bsnmCarlin, \bfnmBradley P.\binitsB. P. (\byear1996). \btitleMarkov chain Monte Carlo convergence diagnostics: A comparative review. \bjournalJ. Amer. Statist. Assoc. \bvolume91 \bpages883–904. \biddoi=10.2307/2291683, issn=01621459, mr=1395755 \bptokimsref\endbibitem
 Crowley, Criste and Smith (1993) {barticle}[author] \bauthor\bsnmCrowley, \bfnmThomas J.\binitsT. J., \bauthor\bsnmCriste, \bfnmTamara A.\binitsT. A. \AND\bauthor\bsnmSmith, \bfnmNeil R.\binitsN. R. (\byear1993). \btitleReassessment of Crete (Greenland) ice core acidity/volcanism link to climate change. \bjournalGeophys. Res. Lett. \bvolume20 \bpages209–212. \biddoi=10.1029/93GL00207 \bptokimsref\endbibitem
 Crowley and Kim (1993) {barticle}[author] \bauthor\bsnmCrowley, \bfnmT.\binitsT. \AND\bauthor\bsnmKim, \bfnmK.\binitsK. (\byear1993). \btitleTowards development of a strategy for determining the origin of decadal–centennial scale climate variability. \bjournalQuat. Sci. Rev. \bvolume12 \bpages375–385. \biddoi=10.1016/S02773791(05)800034 \bptokimsref\endbibitem
 Crowley and Kim (1996) {barticle}[author] \bauthor\bsnmCrowley, \bfnmThomas J.\binitsT. J. \AND\bauthor\bsnmKim, \bfnmKwangYul\binitsK.Y. (\byear1996). \btitleComparison of proxy records of climate change and solar forcing. \bjournalGeophys. Res. Lett. \bvolume23 \bpages359–362. \biddoi=10.1029/96GL00243 \bptokimsref\endbibitem
 Davies and Harte (1987) {barticle}[mr] \bauthor\bsnmDavies, \bfnmR. B.\binitsR. B. \AND\bauthor\bsnmHarte, \bfnmD. S.\binitsD. S. (\byear1987). \btitleTests for Hurst effect. \bjournalBiometrika \bvolume74 \bpages95–101. \biddoi=10.1093/biomet/74.1.95, issn=00063444, mr=0885922 \bptokimsref\endbibitem
 Dempster, Laird and Rubin (1977) {barticle}[mr] \bauthor\bsnmDempster, \bfnmA. P.\binitsA. P., \bauthor\bsnmLaird, \bfnmN. M.\binitsN. M. \AND\bauthor\bsnmRubin, \bfnmD. B.\binitsD. B. (\byear1977). \btitleMaximum likelihood from incomplete data via the EM algorithm. \bjournalJ. Roy. Statist. Soc. Ser. B \bvolume39 \bpages1–38. \bidissn=00359246, mr=0501537 \bptnotecheck related \bptokimsref\endbibitem
 Flato et al. (2013) {bincollection}[author] \bauthor\bsnmFlato, \bfnmG.\binitsG., \bauthor\bsnmMarotzke, \bfnmJ.\binitsJ., \bauthor\bsnmAbiodun, \bfnmB.\binitsB., \bauthor\bsnmBraconnot, \bfnmP.\binitsP., \bauthor\bsnmChou, \bfnmS. C.\binitsS. C., \bauthor\bsnmCollins, \bfnmW.\binitsW., \bauthor\bsnmCox, \bfnmP.\binitsP., \bauthor\bsnmDriouech, \bfnmF.\binitsF., \bauthor\bsnmEmori, \bfnmS.\binitsS., \bauthor\bsnmEyring, \bfnmV.\binitsV., \bauthor\bsnmForest, \bfnmC.\binitsC., \bauthor\bsnmGleckler, \bfnmP.\binitsP., \bauthor\bsnmGuilyardi, \bfnmE.\binitsE., \bauthor\bsnmJakob, \bfnmC.\binitsC., \bauthor\bsnmKattsov, \bfnmV.\binitsV., \bauthor\bsnmReason, \bfnmC.\binitsC. \AND\bauthor\bsnmRummukainen, \bfnmM.\binitsM. (\byear2013). \btitleEvaluation of climate models. In \bbooktitleClimate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (\beditor\bfnmT. F.\binitsT. F. \bsnmStocker, \beditor\bfnmD.\binitsD. \bsnmQin, \beditor\bfnmG. K.\binitsG. K. \bsnmPlattner, \beditor\bfnmM.\binitsM. \bsnmTignor, \beditor\bfnmS. K.\binitsS. K. \bsnmAllen, \beditor\bfnmJ.\binitsJ. \bsnmBoschung, \beditor\bfnmA.\binitsA. \bsnmNauels, \beditor\bfnmY.\binitsY. \bsnmXia, \beditor\bfnmV.\binitsV. \bsnmBex \AND\beditor\bfnmP. M.\binitsP. M. \bsnmMidgley, eds.). \bpublisherCambridge Univ. Press, \blocationCambridge. \bptokimsref\endbibitem
 Gelman and Rubin (1992) {barticle}[author] \bauthor\bsnmGelman, \bfnmA.\binitsA. \AND\bauthor\bsnmRubin, \bfnmD B.\binitsD. B. (\byear1992). \btitleInference from iterative simulation using multiple sequences. \bjournalStatist. Sci. \bvolume7 \bpages457–472. \bptokimsref\endbibitem
 Genz and Bretz (2009) {bbook}[mr] \bauthor\bsnmGenz, \bfnmAlan\binitsA. \AND\bauthor\bsnmBretz, \bfnmFrank\binitsF. (\byear2009). \btitleComputation of Multivariate Normal and Probabilities. \bseriesLecture Notes in Statistics \bvolume195. \bpublisherSpringer, \blocationDordrecht. \biddoi=10.1007/9783642016899, mr=2840595 \bptokimsref\endbibitem
 Gillett et al. (2012) {barticle}[author] \bauthor\bsnmGillett, \bfnmN P.\binitsN. P., \bauthor\bsnmArora, \bfnmV K.\binitsV. K., \bauthor\bsnmFlato, \bfnmG M.\binitsG. M., \bauthor\bsnmScinocca, \bfnmJ F.\binitsJ. F. \AND\bauthor\bparticlevon \bsnmSalzen, \bfnmK.\binitsK. (\byear2012). \btitleImproved constraints on 21stcentury warming derived using 160 years of temperature observations. \bjournalGeophys. Res. Lett. \bvolume39 \bpagesL01704. \biddoi=10.1029/2011GL050226 \bptokimsref\endbibitem
 Gillett et al. (2013) {barticle}[author] \bauthor\bsnmGillett, \bfnmNathan P.\binitsN. P., \bauthor\bsnmArora, \bfnmVivek K.\binitsV. K., \bauthor\bsnmMatthews, \bfnmDamon\binitsD. \AND\bauthor\bsnmAllen, \bfnmMyles R.\binitsM. R. (\byear2013). \btitleConstraining the ratio of global warming to cumulative CO emissions using CMIP5 simulations. \bjournalJ. Climate \bvolume26 \bpages6844–6858. \bptokimsref\endbibitem
 Gneiting, Balabdaoui and Raftery (2007) {barticle}[mr] \bauthor\bsnmGneiting, \bfnmTilmann\binitsT., \bauthor\bsnmBalabdaoui, \bfnmFadoua\binitsF. \AND\bauthor\bsnmRaftery, \bfnmAdrian E.\binitsA. E. (\byear2007). \btitleProbabilistic forecasts, calibration and sharpness. \bjournalJ. R. Stat. Soc. Ser. B Stat. Methodol. \bvolume69 \bpages243–268. \biddoi=10.1111/j.14679868.2007.00587.x, issn=13697412, mr=2325275 \bptokimsref\endbibitem
 Gneiting and Raftery (2007) {barticle}[mr] \bauthor\bsnmGneiting, \bfnmTilmann\binitsT. \AND\bauthor\bsnmRaftery, \bfnmAdrian E.\binitsA. E. (\byear2007). \btitleStrictly proper scoring rules, prediction, and estimation. \bjournalJ. Amer. Statist. Assoc. \bvolume102 \bpages359–378. \biddoi=10.1198/016214506000001437, issn=01621459, mr=2345548 \bptokimsref\endbibitem
 Gneiting and Schlather (2004) {barticle}[mr] \bauthor\bsnmGneiting, \bfnmTilmann\binitsT. \AND\bauthor\bsnmSchlather, \bfnmMartin\binitsM. (\byear2004). \btitleStochastic models that separate fractal dimension and the Hurst effect. \bjournalSIAM Rev. \bvolume46 \bpages269–282 (electronic). \biddoi=10.1137/S0036144501394387, issn=00361445, mr=2114455 \bptokimsref\endbibitem
 Gschlößl and Czado (2007) {barticle}[mr] \bauthor\bsnmGschlößl, \bfnmSusanne\binitsS. \AND\bauthor\bsnmCzado, \bfnmClaudia\binitsC. (\byear2007). \btitleSpatial modelling of claim frequency and claim size in nonlife insurance. \bjournalScand. Actuar. J. \bvolume3 \bpages202–225. \biddoi=10.1080/03461230701414764, issn=03461238, mr=2361126 \bptokimsref\endbibitem
 Harris et al. (2013) {barticle}[author] \bauthor\bsnmHarris, \bfnmG. R.\binitsG. R., \bauthor\bsnmSexton, \bfnmD .M. H.\binitsD. M. H., \bauthor\bsnmBooth, \bfnmB. B. B.\binitsB. B. B., \bauthor\bsnmCollins, \bfnmM.\binitsM. \AND\bauthor\bsnmMurphy, \bfnmJ. M.\binitsJ. M. (\byear2013). \btitleProbabilistic projections of transient climate change. \bjournalClim. Dyn. \bvolume40 \bpages2937–2972. \bptokimsref\endbibitem
 Haslett et al. (2006) {barticle}[mr] \bauthor\bsnmHaslett, \bfnmJ.\binitsJ., \bauthor\bsnmWhiley, \bfnmM.\binitsM., \bauthor\bsnmBhattacharya, \bfnmS.\binitsS., \bauthor\bsnmSalterTownshend, \bfnmM.\binitsM., \bauthor\bsnmWilson, \bfnmSimon P.\binitsS. P., \bauthor\bsnmAllen, \bfnmJ. R. M.\binitsJ. R. M., \bauthor\bsnmHuntley, \bfnmB.\binitsB. \AND\bauthor\bsnmMitchell, \bfnmF. J. G.\binitsF. J. G. (\byear2006). \btitleBayesian palaeoclimate reconstruction. \bjournalJ. Roy. Statist. Soc. Ser. A \bvolume169 \bpages395–438. \biddoi=10.1111/j.1467985X.2006.00429.x, issn=09641998, mr=2236914 \bptokimsref\endbibitem
 Hegerl et al. (2007) {bincollection}[author] \bauthor\bsnmHegerl, \bfnmG. C.\binitsG. C., \bauthor\bsnmZwiers, \bfnmF. W.\binitsF. W., \bauthor\bsnmBraconnot, \bfnmP.\binitsP., \bauthor\bsnmGillett, \bfnmN. P.\binitsN. P., \bauthor\bsnmLuo, \bfnmY.\binitsY., \bauthor\bsnmMarengo Orsini, \bfnmJ. A.\binitsJ. A., \bauthor\bsnmNicholls, \bfnmN.\binitsN., \bauthor\bsnmPenner, \bfnmJ. E.\binitsJ. E. \AND\bauthor\bsnmStott, \bfnmP. A.\binitsP. A. (\byear2007). \btitleUnderstanding and attributing climate change. In \bbooktitleClimate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change (\beditor\bfnmS.\binitsS. \bsnmSolomon, \beditor\bfnmD.\binitsD. \bsnmQin, \beditor\bfnmM.\binitsM. \bsnmManning, \beditor\bfnmZ.\binitsZ. \bsnmChen, \beditor\bfnmM.\binitsM. \bsnmMarquis, \beditor\bfnmK B.\binitsK. B. \bsnmAveryt, \beditor\bfnmM.\binitsM. \bsnmTignor \AND\beditor\bfnmH L.\binitsH. L. \bsnmMiller, eds.). \bpublisherCambridge Univ. Press, \blocationCambridge. \bptokimsref\endbibitem
 Huybers and Curry (2006) {barticle}[pbm] \bauthor\bsnmHuybers, \bfnmPeter\binitsP. \AND\bauthor\bsnmCurry, \bfnmWilliam\binitsW. (\byear2006). \btitleLinks between annual, Milankovitch and continuum temperature variability. \bjournalNature \bvolume441 \bpages329–332. \biddoi=10.1038/nature04745, issn=14764687, pii=nature04745, pmid=16710417 \bptokimsref\endbibitem
 Imbers et al. (2014) {barticle}[author] \bauthor\bsnmImbers, \bfnmJara\binitsJ., \bauthor\bsnmLopez, \bfnmAna\binitsA., \bauthor\bsnmHuntingford, \bfnmChris\binitsC. \AND\bauthor\bsnmAllen, \bfnmMyles\binitsM. (\byear2014). \btitleSensitivity of climate change detection and attribution to the characterization of internal variability. \bjournalJ. Climate \bvolume27 \bpages3477–3491. \bptokimsref\endbibitem
 Jackson (2011) {barticle}[author] \bauthor\bsnmJackson, \bfnmChristopher\binitsC. (\byear2011). \btitleMultistate models for panel data: The msm package for R. \bjournalJ. Stat. Softw. \bvolume38 \bpages1–28. \bptokimsref\endbibitem
 Jones et al. (2009) {barticle}[author] \bauthor\bsnmJones, \bfnmP D.\binitsP. D., \bauthor\bsnmBriffa, \bfnmK R.\binitsK. R., \bauthor\bsnmOsborn, \bfnmT J.\binitsT. J., \bauthor\bsnmLough, \bfnmJ M.\binitsJ. M., \bauthor\bparticlevan \bsnmOmmen, \bfnmT D.\binitsT. D., \bauthor\bsnmVinther, \bfnmB M.\binitsB. M., \bauthor\bsnmLuterbacher, \bfnmJ.\binitsJ., \bauthor\bsnmWahl, \bfnmE R.\binitsE. R., \bauthor\bsnmZwiers, \bfnmF W.\binitsF. W., \bauthor\bsnmMann, \bfnmM E.\binitsM. E., \bauthor\bsnmSchmidt, \bfnmG A.\binitsG. A., \bauthor\bsnmAmmann, \bfnmC M.\binitsC. M., \bauthor\bsnmBuckley, \bfnmB M.\binitsB. M., \bauthor\bsnmCobb, \bfnmK M.\binitsK. M., \bauthor\bsnmEsper, \bfnmJ.\binitsJ., \bauthor\bsnmGoosse, \bfnmH.\binitsH., \bauthor\bsnmGraham, \bfnmN.\binitsN., \bauthor\bsnmJansen, \bfnmE.\binitsE., \bauthor\bsnmKiefer, \bfnmT.\binitsT., \bauthor\bsnmKull, \bfnmC.\binitsC., \bauthor\bsnmKüttel, \bfnmM.\binitsM., \bauthor\bsnmMosleyThompson, \bfnmE.\binitsE., \bauthor\bsnmOverpeck, \bfnmJ T.\binitsJ. T., \bauthor\bsnmRiedwyl, \bfnmN.\binitsN., \bauthor\bsnmSchulz, \bfnmM.\binitsM., \bauthor\bsnmTudhope, \bfnmA W.\binitsA. W., \bauthor\bsnmVillalba, \bfnmR.\binitsR., \bauthor\bsnmWanner, \bfnmH.\binitsH., \bauthor\bsnmWolff, \bfnmE.\binitsE. \AND\bauthor\bsnmXoplaki, \bfnmE.\binitsE. (\byear2009). \btitleHighresolution palaeoclimatology of the last millennium: A review of current status and future prospects. \bjournalHolocene \bvolume19 \bpages3–49. \biddoi=10.1177/0959683608098952 \bptokimsref\endbibitem
 Juckes et al. (2006) {barticle}[author] \bauthor\bsnmJuckes, \bfnmM N.\binitsM. N., \bauthor\bsnmAllen, \bfnmM R.\binitsM. R., \bauthor\bsnmBriffa, \bfnmK R.\binitsK. R., \bauthor\bsnmEsper, \bfnmJ.\binitsJ., \bauthor\bsnmHegerl, \bfnmG C.\binitsG. C., \bauthor\bsnmMoberg, \bfnmA.\binitsA., \bauthor\bsnmOsborn, \bfnmT J.\binitsT. J., \bauthor\bsnmWeber, \bfnmS L.\binitsS. L. \AND\bauthor\bsnmZorita, \bfnmE.\binitsE. (\byear2006). \btitleMillennial temperature reconstruction intercomparison and evaluation. \bjournalClim. Past Discuss. \bvolume2 \bpages1001–1049. \bptokimsref\endbibitem
 Kalman and Bucy (1961) {barticle}[mr] \bauthor\bsnmKalman, \bfnmR. E.\binitsR. E. \AND\bauthor\bsnmBucy, \bfnmR. S.\binitsR. S. (\byear1961). \btitleNew results in linear filtering and prediction theory. \bjournalTrans. ASME Ser. D, J. Basic Engrg. \bvolume83 \bpages95–108. \bidmr=0234760 \bptokimsref\endbibitem
 Kaufman et al. (2009) {barticle}[author] \bauthor\bsnmKaufman, \bfnmD S.\binitsD. S., \bauthor\bsnmSchneider, \bfnmD P.\binitsD. P., \bauthor\bsnmMcKay, \bfnmN P.\binitsN. P., \bauthor\bsnmAmmann, \bfnmC M.\binitsC. M., \bauthor\bsnmBradley, \bfnmR S.\binitsR. S., \bauthor\bsnmBriffa, \bfnmK R.\binitsK. R., \bauthor\bsnmMiller, \bfnmG H.\binitsG. H., \bauthor\bsnmOttoBliesner, \bfnmB L.\binitsB. L., \bauthor\bsnmOverpeck, \bfnmJ T.\binitsJ. T., \bauthor\bsnmVinther, \bfnmB M.\binitsB. M. \AND\bauthor\bsnmOthers (\byear2009). \btitleRecent warming reverses longterm arctic cooling. \bjournalScience \bvolume325 \bpages1236–1239. \bptokimsref\endbibitem
 Kolmogoroff (1940) {barticle}[mr] \bauthor\bsnmKolmogoroff, \bfnmA. N.\binitsA. N. (\byear1940). \btitleWienersche Spiralen und einige andere interessante Kurven im Hilbertschen Raum. \bjournalC. R. (Doklady) Acad. Sci. USSR (N.S.) \bvolume26 \bpages115–118. \bidmr=0003441 \bptokimsref\endbibitem
 Landrum et al. (2013) {barticle}[author] \bauthor\bsnmLandrum, \bfnmLaura\binitsL., \bauthor\bsnmOttoBliesner, \bfnmBette L.\binitsB. L., \bauthor\bsnmWahl, \bfnmEugene R.\binitsE. R., \bauthor\bsnmConley, \bfnmAndrew\binitsA., \bauthor\bsnmLawrence, \bfnmPeter J.\binitsP. J., \bauthor\bsnmRosenbloom, \bfnmNan\binitsN. \AND\bauthor\bsnmTeng, \bfnmHaiyan\binitsH. (\byear2013). \btitleLast millennium climate and its variability in CCSM4. \bjournalJ. Climate \bvolume26 \bpages1085–1111. \biddoi=10.1175/JCLID1100326.1 \bptokimsref\endbibitem
 Lean, Beer and Bradley (1995) {barticle}[author] \bauthor\bsnmLean, \bfnmJudith\binitsJ., \bauthor\bsnmBeer, \bfnmJuerg\binitsJ. \AND\bauthor\bsnmBradley, \bfnmRaymond\binitsR. (\byear1995). \btitleReconstruction of solar irradiance since 1610: Implications for climate change. \bjournalGeophys. Res. Lett. \bvolume22 \bpages3195–3198. \biddoi=10.1029/95GL03093 \bptokimsref\endbibitem
 Lee, Zwiers and Tsao (2008) {barticle}[author] \bauthor\bsnmLee, \bfnmTerry C K.\binitsT. C. K., \bauthor\bsnmZwiers, \bfnmFrancis W.\binitsF. W. \AND\bauthor\bsnmTsao, \bfnmMin\binitsM. (\byear2008). \btitleEvaluation of proxybased millennial reconstruction methods. \bjournalClim. Dyn. \bvolume31 \bpages263–281. \biddoi=10.1007/s0038200703519 \bptokimsref\endbibitem
 Li, Nychka and Ammann (2010) {barticle}[mr] \bauthor\bsnmLi, \bfnmBo\binitsB., \bauthor\bsnmNychka, \bfnmDouglas W.\binitsD. W. \AND\bauthor\bsnmAmmann, \bfnmCaspar M.\binitsC. M. (\byear2010). \btitleThe value of multiproxy reconstruction of past climate. \bjournalJ. Amer. Statist. Assoc. \bvolume105 \bpages883–895. \biddoi=10.1198/jasa.2010.ap09379, issn=01621459, mr=2752583 \bptokimsref\endbibitem
 Loso (2009) {barticle}[author] \bauthor\bsnmLoso, \bfnmM G.\binitsM. G. (\byear2009). \btitleSummer temperatures during the Medieval Warm Period and Little Ice Age inferred from varved proglacial lake sediments in southern Alaska. \bjournalJ. Paleolimnol. \bvolume41 \bpages117–128. \bptokimsref\endbibitem
 Luterbacher et al. (2004) {barticle}[pbm] \bauthor\bsnmLuterbacher, \bfnmJürg\binitsJ., \bauthor\bsnmDietrich, \bfnmDaniel\binitsD., \bauthor\bsnmXoplaki, \bfnmElena\binitsE., \bauthor\bsnmGrosjean, \bfnmMartin\binitsM. \AND\bauthor\bsnmWanner, \bfnmHeinz\binitsH. (\byear2004). \btitleEuropean seasonal and annual temperature variability, trends, and extremes since 1500. \bjournalScience \bvolume303 \bpages1499–1503. \biddoi=10.1126/science.1093877, issn=10959203, pii=303/5663/1499, pmid=15001774 \bptokimsref\endbibitem
 Mandelbrot (1965) {barticle}[mr] \bauthor\bsnmMandelbrot, \bfnmBenoît\binitsB. (\byear1965). \btitleUne classe processus stochastiques homothétiques à soi; application à la loi climatologique H. E. Hurst. \bjournalC. R. Math. Acad. Sci. Paris \bvolume260 \bpages3274–3277. \bidmr=0176521 \bptokimsref\endbibitem
 Mandelbrot and Van Ness (1968) {barticle}[mr] \bauthor\bsnmMandelbrot, \bfnmBenoit B.\binitsB. B. \AND\bauthor\bsnmVan Ness, \bfnmJohn W.\binitsJ. W. (\byear1968). \btitleFractional Brownian motions, fractional noises and applications. \bjournalSIAM Rev. \bvolume10 \bpages422–437. \bidissn=00361445, mr=0242239 \bptokimsref\endbibitem
 (54) {barticle}[author] \bauthor\bsnmMann, \bfnmMichael E.\binitsM. E., \bauthor\bsnmBradley, \bfnmRaymond S.\binitsR. S. \AND\bauthor\bsnmHughes, \bfnmMalcolm K.\binitsM. K. (\byear1998). \btitleGlobalscale temperature patterns and climate forcing over the past six centuries. \bjournalNature \bvolume392 \bpages779–787. \biddoi=10.1038/33859 \bptokimsref\endbibitem
 Mann, Bradley and Hughes (1999) {barticle}[author] \bauthor\bsnmMann, \bfnmMichael E.\binitsM. E., \bauthor\bsnmBradley, \bfnmRaymond S.\binitsR. S. \AND\bauthor\bsnmHughes, \bfnmMalcolm K.\binitsM. K. (\byear1999). \btitleNorthern Hemisphere temperatures during the past millennium: Inferences, uncertainties, and limitations. \bjournalGeophys. Res. Lett. \bvolume26 \bpages759–762. \biddoi=10.1029/1999GL900070 \bptokimsref\endbibitem
 Mann et al. (2005) {barticle}[author] \bauthor\bsnmMann, \bfnmM E.\binitsM. E., \bauthor\bsnmRutherford, \bfnmS.\binitsS., \bauthor\bsnmWahl, \bfnmE.\binitsE. \AND\bauthor\bsnmAmmann, \bfnmC.\binitsC. (\byear2005). \btitleTesting the fidelity of methods used in proxybased reconstructions of past climate. \bjournalJ. Climate \bvolume18 \bpages4097–4107. \bptokimsref\endbibitem
 Mann et al. (2007) {barticle}[author] \bauthor\bsnmMann, \bfnmMichael E.\binitsM. E., \bauthor\bsnmRutherford, \bfnmScott\binitsS., \bauthor\bsnmWahl, \bfnmEugene\binitsE. \AND\bauthor\bsnmAmmann, \bfnmCaspar\binitsC. (\byear2007). \btitleRobustness of proxybased climate field reconstruction methods. \bjournalJ. Geophys. Res. \bvolume112 \bpagesD12109+. \biddoi=10.1029/2006JD008272 \bptokimsref\endbibitem
 Mann et al. (2008a) {barticle}[author] \bauthor\bsnmMann, \bfnmMichael E.\binitsM. E., \bauthor\bsnmZhang, \bfnmZhihua\binitsZ., \bauthor\bsnmHughes, \bfnmMalcolm K.\binitsM. K., \bauthor\bsnmBradley, \bfnmRaymond S.\binitsR. S., \bauthor\bsnmMiller, \bfnmSonya K.\binitsS. K., \bauthor\bsnmRutherford, \bfnmScott\binitsS. \AND\bauthor\bsnmNi, \bfnmFenbiao\binitsF. (\byear2008a). \btitleProxybased reconstructions of hemispheric and global surface temperature variations over the past two millennia. \bjournalProc. Natl. Acad. Sci. USA \bvolume105 \bpages13252–13257. \biddoi=10.1073/pnas.0805721105 \bptokimsref\endbibitem
 Mann et al. (2008b) {barticle}[author] \bauthor\bsnmMann, \bfnmMichael E.\binitsM. E., \bauthor\bsnmZhang, \bfnmZhihua\binitsZ., \bauthor\bsnmHughes, \bfnmMalcolm K.\binitsM. K., \bauthor\bsnmBradley, \bfnmRaymond S.\binitsR. S., \bauthor\bsnmMiller, \bfnmSonya K.\binitsS. K., \bauthor\bsnmRutherford, \bfnmScott\binitsS. \AND\bauthor\bsnmNi, \bfnmFenbiao\binitsF. (\byear2008b). \btitleSupporting information for “Proxybased reconstructions of hemispheric and global surface temperature variations over the past two millennia.” \bjournalProc. Natl. Acad. Sci. USA \bvolume105 \bpages13252–13257. \bptokimsref\endbibitem
 Mann et al. (2009) {barticle}[author] \bauthor\bsnmMann, \bfnmM E.\binitsM. E., \bauthor\bsnmZhang, \bfnmZ.\binitsZ., \bauthor\bsnmRutherford, \bfnmS.\binitsS., \bauthor\bsnmBradley, \bfnmR S.\binitsR. S., \bauthor\bsnmHughes, \bfnmM K.\binitsM. K., \bauthor\bsnmShindell, \bfnmD.\binitsD., \bauthor\bsnmAmmann, \bfnmC M.\binitsC. M., \bauthor\bsnmFaluvegi, \bfnmG.\binitsG. \AND\bauthor\bsnmNi, \bfnmF.\binitsF. (\byear2009). \btitleGlobal signatures and dynamical origins of the Little Ice Age and Medieval Climate Anomaly. \bjournalScience \bvolume326 \bpages1256–1260. \bptokimsref\endbibitem
 Martin, Quinn and Park (2011) {barticle}[author] \bauthor\bsnmMartin, \bfnmAndrew D.\binitsA. D., \bauthor\bsnmQuinn, \bfnmKevin M.\binitsK. M. \AND\bauthor\bsnmPark, \bfnmJong H.\binitsJ. H. (\byear2011). \btitleMCMCpack: Markov chain Monte Carlo in R. \bjournalJ. Stat. Softw. \bvolume42 \bpages1–21. \bptokimsref\endbibitem
 MassonDelmotte et al. (2013) {bincollection}[author] \bauthor\bsnmMassonDelmotte, \bfnmV.\binitsV., \bauthor\bsnmSchulz, \bfnmM.\binitsM., \bauthor\bsnmAbeOuchi, \bfnmA.\binitsA., \bauthor\bsnmBeer, \bfnmJ.\binitsJ., \bauthor\bsnmGanopolski, \bfnmA.\binitsA., \bauthor\bsnmRouco, \bfnmJ. F. González\binitsJ. F. G., \bauthor\bsnmJansen, \bfnmE.\binitsE., \bauthor\bsnmLambeck, \bfnmK.\binitsK., \bauthor\bsnmLuterbacher, \bfnmJ.\binitsJ., \bauthor\bsnmNaish, \bfnmT.\binitsT., \bauthor\bsnmOsborn, \bfnmT.\binitsT., \bauthor\bsnmOttoBliesner, \bfnmB.\binitsB., \bauthor\bsnmQuinn, \bfnmT.\binitsT., \bauthor\bsnmRamesh, \bfnmR.\binitsR., \bauthor\bsnmRojas, \bfnmM.\binitsM., \bauthor\bsnmShao, \bfnmX.\binitsX. \AND\bauthor\bsnmTimmermann, \bfnmA.\binitsA. (\byear2013). \btitleInformation from paleoclimate archives. In \bbooktitleClimate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (\beditor\bfnmT. F.\binitsT. F. \bsnmStocker, \beditor\bfnmD.\binitsD. \bsnmQin, \beditor\bfnmG. K.\binitsG. K. \bsnmPlattner, \beditor\bfnmM.\binitsM. \bsnmTignor, \beditor\bfnmS. K.\binitsS. K. \bsnmAllen, \beditor\bfnmJ.\binitsJ. \bsnmBoschung, \beditor\bfnmA.\binitsA. \bsnmNauels, \beditor\bfnmY.\binitsY. \bsnmXia, \beditor\bfnmV.\binitsV. \bsnmBex \AND\beditor\bfnmP. M.\binitsP. M. \bsnmMidgley, eds.). \bpublisherCambridge Univ. Press, \blocationCambridge. \bptokimsref\endbibitem
 McLeod, Yu and Krougly (2007) {barticle}[author] \bauthor\bsnmMcLeod, \bfnmA Ian\binitsA. I., \bauthor\bsnmYu, \bfnmHao\binitsH. \AND\bauthor\bsnmKrougly, \bfnmZinovi L.\binitsZ. L. (\byear2007). \btitleAlgorithms for linear time series analysis: With R package. \bjournalJournal of Statistical Software \bvolume23 \bpages1–26. \bptokimsref\endbibitem
 McShane and Wyner (2011) {barticle}[author] \bauthor\bsnmMcShane, \bfnmBlakeley B.\binitsB. B. \AND\bauthor\bsnmWyner, \bfnmAbraham J.\binitsA. J. (\byear2011). \btitleA statistical analysis of multiple temperature proxies: Are reconstructions of surface temperatures over the last 1000 years reliable? \bjournalAnn. Appl. Stat. \bvolume5 \bpages5–44. \bptokimsref\endbibitem
 Moberg et al. (2005) {barticle}[pbm] \bauthor\bsnmMoberg, \bfnmAnders\binitsA., \bauthor\bsnmSonechkin, \bfnmDmitry M.\binitsD. M., \bauthor\bsnmHolmgren, \bfnmKarin\binitsK., \bauthor\bsnmDatsenko, \bfnmNina M.\binitsN. M., \bauthor\bsnmKarlén, \bfnmWibjörn\binitsW. \AND\bauthor\bsnmLauritzen, \bfnmSteinErik\binitsS.E. (\byear2005). \btitleHighly variable Northern Hemisphere temperatures reconstructed from low and highresolution proxy data. \bjournalNature \bvolume433 \bpages613–617. \biddoi=10.1038/nature03265, issn=14764687, pii=nature03265, pmid=15703742 \bptokimsref\endbibitem
 Murphy and Izzeldin (2009) {barticle}[mr] \bauthor\bsnmMurphy, \bfnmA.\binitsA. \AND\bauthor\bsnmIzzeldin, \bfnmM.\binitsM. (\byear2009). \btitleBootstrapping long memory tests: Some Monte Carlo results. \bjournalComput. Statist. Data Anal. \bvolume53 \bpages2325–2334. \biddoi=10.1016/j.csda.2008.07.040, issn=01679473, mr=2665831 \bptokimsref\endbibitem
 NRC (2006) {bbook}[author] \bauthor\bsnmNRC (\byear2006). \btitleSurface Temperature Reconstructions for the Last 2000 Years. \bpublisherThe National Academies Press, \blocationWashington, DC. \bptokimsref\endbibitem
 Otto et al. (2013) {barticle}[author] \bauthor\bsnmOtto, \bfnmAlexander\binitsA., \bauthor\bsnmOtto, \bfnmFriederike E L.\binitsF. E. L., \bauthor\bsnmBoucher, \bfnmOlivier\binitsO., \bauthor\bsnmChurch, \bfnmJohn\binitsJ., \bauthor\bsnmHegerl, \bfnmGabi\binitsG., \bauthor\bsnmForster, \bfnmPiers M.\binitsP. M., \bauthor\bsnmGillett, \bfnmNathan P.\binitsN. P., \bauthor\bsnmGregory, \bfnmJonathan\binitsJ., \bauthor\bsnmJohnson, \bfnmGregory C.\binitsG. C., \bauthor\bsnmKnutti, \bfnmReto\binitsR., \bauthor\bsnmLewis, \bfnmNicholas\binitsN., \bauthor\bsnmLohmann, \bfnmUlrike\binitsU., \bauthor\bsnmMarotzke, \bfnmJochem\binitsJ., \bauthor\bsnmMyhre, \bfnmGunnar\binitsG., \bauthor\bsnmShindell, \bfnmDrew\binitsD., \bauthor\bsnmStevens, \bfnmBjorn\binitsB. \AND\bauthor\bsnmAllen, \bfnmMyles R.\binitsM. R. (\byear2013). \btitleEnergy budget constraints on climate response. \bjournalNat. Geosci. \bvolume6 \bpages415–416. \bptokimsref\endbibitem
 Overpeck et al. (1997) {barticle}[author] \bauthor\bsnmOverpeck, \bfnmJ.\binitsJ., \bauthor\bsnmHughen, \bfnmK.\binitsK., \bauthor\bsnmHardy, \bfnmD.\binitsD., \bauthor\bsnmBradley, \bfnmR.\binitsR., \bauthor\bsnmCase, \bfnmR.\binitsR., \bauthor\bsnmDouglas, \bfnmM.\binitsM., \bauthor\bsnmFinney, \bfnmB.\binitsB., \bauthor\bsnmGajewski, \bfnmK.\binitsK., \bauthor\bsnmJacoby, \bfnmG.\binitsG., \bauthor\bsnmJennings, \bfnmA.\binitsA., \bauthor\bsnmLamoureux, \bfnmS.\binitsS., \bauthor\bsnmLasca, \bfnmA.\binitsA., \bauthor\bsnmMacDonald, \bfnmG.\binitsG., \bauthor\bsnmMoore, \bfnmJ.\binitsJ., \bauthor\bsnmRetelle, \bfnmM.\binitsM., \bauthor\bsnmSmith, \bfnmS.\binitsS., \bauthor\bsnmWolfe, \bfnmA.\binitsA. \AND\bauthor\bsnmZielinski, \bfnmG.\binitsG. (\byear1997). \btitleArctic environmental change of the last four centuries. \bjournalScience \bvolume278 \bpages1251–1256. \bptokimsref\endbibitem
 PAGES 2k Consortium (2013) {barticle}[author] \bauthor\bsnmPAGES 2k Consortium (\byear2013). \btitleContinentalscale temperature variability during the past two millennia. \bjournalNature Geoscience \bvolume6 \bpages339–346. \biddoi=10.1038/ngeo1797 \bptokimsref\endbibitem
 Palma and Bondon (2003) {barticle}[mr] \bauthor\bsnmPalma, \bfnmWilfredo\binitsW. \AND\bauthor\bsnmBondon, \bfnmPascal\binitsP. (\byear2003). \btitleOn the eigenstructure of generalized fractional processes. \bjournalStatist. Probab. Lett. \bvolume65 \bpages93–101. \biddoi=10.1016/j.spl.2003.07.008, issn=01677152, mr=2017253 \bptokimsref\endbibitem
 Percival, Overland and Mofjeld (2001) {barticle}[author] \bauthor\bsnmPercival, \bfnmD. B.\binitsD. B., \bauthor\bsnmOverland, \bfnmJ. E.\binitsJ. E. \AND\bauthor\bsnmMofjeld, \bfnmH. O.\binitsH. O. (\byear2001). \btitleInterpretation of North Pacific variability as a short and longmemory process. \bjournalJ. Climate \bvolume14 \bpages7–11. \bptokimsref\endbibitem
 Plummer et al. (2006) {barticle}[author] \bauthor\bsnmPlummer, \bfnmMartyn\binitsM., \bauthor\bsnmBest, \bfnmNicky\binitsN., \bauthor\bsnmCowles, \bfnmKate\binitsK. \AND\bauthor\bsnmVines, \bfnmKaren\binitsK. (\byear2006). \btitleCODA: Convergence diagnosis and output analysis for MCMC. \bjournalR News \bvolume6 \bpages7–11. \bptokimsref\endbibitem
 Robinson (1995) {barticle}[mr] \bauthor\bsnmRobinson, \bfnmP. M.\binitsP. M. (\byear1995). \btitleGaussian semiparametric estimation of long range dependence. \bjournalAnn. Statist. \bvolume23 \bpages1630–1661. \biddoi=10.1214/aos/1176324317, issn=00905364, mr=1370301 \bptokimsref\endbibitem
 Rutherford et al. (2003) {barticle}[author] \bauthor\bsnmRutherford, \bfnmS.\binitsS., \bauthor\bsnmMann, \bfnmM E.\binitsM. E., \bauthor\bsnmDelworth, \bfnmT L.\binitsT. L. \AND\bauthor\bsnmStouffer, \bfnmR J.\binitsR. J. (\byear2003). \btitleClimate field reconstruction under stationary and nonstationary forcing. \bjournalJ. Climate \bvolume16 \bpages462–479. \bptokimsref\endbibitem
 Rutherford et al. (2005) {barticle}[author] \bauthor\bsnmRutherford, \bfnmS.\binitsS., \bauthor\bsnmMann, \bfnmM E.\binitsM. E., \bauthor\bsnmOsborn, \bfnmT J.\binitsT. J., \bauthor\bsnmBriffa, \bfnmK R.\binitsK. R., \bauthor\bsnmJones, \bfnmP D.\binitsP. D., \bauthor\bsnmBradley, \bfnmR S.\binitsR. S. \AND\bauthor\bsnmHughes, \bfnmM K.\binitsM. K. (\byear2005). \btitleProxybased Northern Hemisphere surface temperature reconstructions: Sensitivity to method, predictor network, target season, and target domain. \bjournalJ. Climate \bvolume18 \bpages2308–2329. \bptokimsref\endbibitem
 Schneider (2001) {barticle}[author] \bauthor\bsnmSchneider, \bfnmT.\binitsT. (\byear2001). \btitleAnalysis of incomplete climate data: Estimation of mean values and covariance matrices and imputation of missing values. \bjournalJ. Climate \bvolume14 \bpages853–871. \bptokimsref\endbibitem
 Smerdon et al. (2010) {barticle}[author] \bauthor\bsnmSmerdon, \bfnmJ E.\binitsJ. E., \bauthor\bsnmKaplan, \bfnmA.\binitsA., \bauthor\bsnmChang, \bfnmD.\binitsD. \AND\bauthor\bsnmEvans, \bfnmM N.\binitsM. N. (\byear2010). \btitleA pseudoproxy evaluation of the CCA and RegEM methods for reconstructing climate fields of the last millennium. \bjournalJ. Climate \bvolume23 \bpages4856–4880. \bptokimsref\endbibitem
 Smith (2010) {bmisc}[author] \bauthor\bsnmSmith, \bfnmRichard L.\binitsR. L. (\byear2010). \bhowpublishedUnderstanding sensitivities in paleoclimatic reconstructions. Preprint. \bptokimsref\endbibitem
 Steig et al. (2009) {barticle}[pbm] \bauthor\bsnmSteig, \bfnmEric J.\binitsE. J., \bauthor\bsnmSchneider, \bfnmDavid P.\binitsD. P., \bauthor\bsnmRutherford, \bfnmScott D.\binitsS. D., \bauthor\bsnmMann, \bfnmMichael E.\binitsM. E., \bauthor\bsnmComiso, \bfnmJosefino C.\binitsJ. C. \AND\bauthor\bsnmShindell, \bfnmDrew T.\binitsD. T. (\byear2009). \btitleWarming of the Antarctic icesheet surface since the 1957 international geophysical year. \bjournalNature \bvolume457 \bpages459–462. \biddoi=10.1038/nature07669, issn=14764687, pii=nature07669, pmid=19158794 \bptokimsref\endbibitem
 Stocker et al. (2013) {bincollection}[author] \bauthor\bsnmStocker, \bfnmThomas F.\binitsT. F., \bauthor\bsnmDahe, \bfnmQin\binitsQ., \bauthor\bsnmPlattner, \bfnmGianKasper\binitsG.K., \bauthor\bsnmAlexander, \bfnmLisa\binitsL., \bauthor\bsnmAllen, \bfnmSimon\binitsS., \bauthor\bsnmBindoff, \bfnmNathaniel\binitsN., \bauthor\bsnmBreon, \bfnmFrancoisMarie\binitsF.M., \bauthor\bsnmChurch, \bfnmJohn\binitsJ., \bauthor\bsnmCubasch, \bfnmUlrich\binitsU., \bauthor\bsnmEmori, \bfnmSeita\binitsS., \bauthor\bsnmForster, \bfnmPiers\binitsP., \bauthor\bsnmFriedlingstein, \bfnmPierre\binitsP., \bauthor\bsnmGillett, \bfnmNathan\binitsN., \bauthor\bsnmGregory, \bfnmJonathan\binitsJ., \bauthor\bsnmHartmann, \bfnmDennis\binitsD., \bauthor\bsnmJansen, \bfnmEystein\binitsE., \bauthor\bsnmKirtman, \bfnmBen\binitsB., \bauthor\bsnmKnutti, \bfnmReto\binitsR., \bauthor\bsnmKumar Kanikicharla, \bfnmKrishna\binitsK., \bauthor\bsnmLemke, \bfnmPeter\binitsP., \bauthor\bsnmMarotzke, \bfnmJochem\binitsJ., \bauthor\bsnmMassonDelmotte, \bfnmValerie\binitsV., \bauthor\bsnmMeehl, \bfnmGerald\binitsG., \bauthor\bsnmMokhov, \bfnmIgor\binitsI., \bauthor\bsnmPiao, \bfnmShilong\binitsS., \bauthor\bsnmRamaswamy, \bfnmVenkatachalam\binitsV., \bauthor\bsnmRandall, \bfnmDavid\binitsD., \bauthor\bsnmRhein, \bfnmMonika\binitsM., \bauthor\bsnmRojas, \bfnmMaisa\binitsM., \bauthor\bsnmSabine, \bfnmChristopher\binitsC., \bauthor\bsnmShindell, \bfnmDrew\binitsD., \bauthor\bsnmTalley, \bfnmLynne\binitsL., \bauthor\bsnmVaughan, \bfnmDavid\binitsD. \AND\bauthor\bsnmXie, \bfnmShangPing\binitsS.P. (\byear2013). \btitleTechnical summary. In \bbooktitleClimate Change 2013: The Physical Science Basis (\beditor\bfnmThomas F.\binitsT. F. \bsnmStocker, \beditor\bfnmQin\binitsQ. \bsnmDahe, \beditor\bfnmGianKasper\binitsG.K. \bsnmPlattner, \beditor\bfnmM.\binitsM. \bsnmTignor, \beditor\bfnmSimon\binitsS. \bsnmAllen, \beditor\bfnmJ.\binitsJ. \bsnmBoschung, \beditor\bfnmA.\binitsA. \bsnmNauels, \beditor\bfnmY.\binitsY. \bsnmXia, \beditor\bfnmV.\binitsV. \bsnmBex \AND\beditor\bfnmP. M.\binitsP. M. \bsnmMidgley, eds.). \bpublisherCambridge Univ. Press, \blocationCambridge. \bptokimsref\endbibitem
 Taqqu (2013) {barticle}[mr] \bauthor\bsnmTaqqu, \bfnmMurad S.\binitsM. S. (\byear2013). \btitleBenoît Mandelbrot and fractional Brownian motion. \bjournalStatist. Sci. \bvolume28 \bpages131–134. \biddoi=10.1214/12STS389, issn=08834237, mr=3075342 \bptokimsref\endbibitem
 Tiljander et al. (2003) {barticle}[author] \bauthor\bsnmTiljander, \bfnmM.\binitsM., \bauthor\bsnmSaarnisto, \bfnmM.\binitsM., \bauthor\bsnmOjala, \bfnmA E K.\binitsA. E. K. \AND\bauthor\bsnmSaarinen, \bfnmT.\binitsT. (\byear2003). \btitleA 3000year palaeoenvironmental record from annually laminated sediment of lake Korttajärvi, central Finland. \bjournalBoreas \bvolume32 \bpages566–577. \bptokimsref\endbibitem
 Tingley and Huybers (2010a) {barticle}[author] \bauthor\bsnmTingley, \bfnmMartin P.\binitsM. P. \AND\bauthor\bsnmHuybers, \bfnmPeter\binitsP. (\byear2010a). \btitleA Bayesian algorithm for reconstructing climate anomalies in space and time. Part I: Development and applications to paleoclimate reconstruction problems. \bjournalJ. Climate \bvolume23 \bpages2759–2781. \biddoi=10.1175/2009JCLI3015.1 \bptokimsref\endbibitem
 Tingley and Huybers (2010b) {barticle}[author] \bauthor\bsnmTingley, \bfnmMartin P.\binitsM. P. \AND\bauthor\bsnmHuybers, \bfnmPeter\binitsP. (\byear2010b). \btitleA Bayesian algorithm for reconstructing climate anomalies in space and time. Part II: Comparison with the regularized expectationmaximization algorithm. \bjournalJ. Climate \bvolume23 \bpages2782–2800. \biddoi=10.1175/2009JCLI3016.1 \bptokimsref\endbibitem
 Tingley and Huybers (2013) {barticle}[pbm] \bauthor\bsnmTingley, \bfnmMartin P.\binitsM. P. \AND\bauthor\bsnmHuybers, \bfnmPeter\binitsP. (\byear2013). \btitleRecent temperature extremes at high northern latitudes unprecedented in the past 600 years. \bjournalNature \bvolume496 \bpages201–205. \biddoi=10.1038/nature11969, issn=14764687, pii=nature11969, pmid=23579678 \bptokimsref\endbibitem
 Tingley et al. (2012) {barticle}[author] \bauthor\bsnmTingley, \bfnmMartin P.\binitsM. P., \bauthor\bsnmCraigmile, \bfnmPeter F.\binitsP. F., \bauthor\bsnmHaran, \bfnmMurali\binitsM., \bauthor\bsnmLi, \bfnmBo\binitsB., \bauthor\bsnmMannshardt, \bfnmElizabeth\binitsE. \AND\bauthor\bsnmRajaratnam, \bfnmBala\binitsB. (\byear2012). \btitlePiecing together the past: Statistical insights into paleoclimatic reconstructions. \bjournalQuat. Sci. Rev. \bvolume35 \bpages1–22. \biddoi=10.1016/j.quascirev.2012.01.012 \bptokimsref\endbibitem
 Tudor and Viens (2007) {barticle}[mr] \bauthor\bsnmTudor, \bfnmCiprian A.\binitsC. A. \AND\bauthor\bsnmViens, \bfnmFrederi G.\binitsF. G. (\byear2007). \btitleStatistical aspects of the fractional stochastic calculus. \bjournalAnn. Statist. \bvolume35 \bpages1183–1212. \biddoi=10.1214/009053606000001541, issn=00905364, mr=2341703 \bptokimsref\endbibitem
 Tung, Zhou and Camp (2008) {barticle}[author] \bauthor\bsnmTung, \bfnmK. K.\binitsK. K., \bauthor\bsnmZhou, \bfnmJ. S.\binitsJ. S. \AND\bauthor\bsnmCamp, \bfnmC. D.\binitsC. D. (\byear2008). \btitleConstraining model transient climate response using independent observations of solarcycle forcing and response. \bjournalGeophys. Res. Lett. \bvolume35 \bpagesL17707. \bptokimsref\endbibitem
 Wahl and Smerdon (2012) {barticle}[author] \bauthor\bsnmWahl, \bfnmEugene R.\binitsE. R. \AND\bauthor\bsnmSmerdon, \bfnmJason E.\binitsJ. E. (\byear2012). \btitleComparative performance of paleoclimate field and index reconstructions derived from climate proxies and noiseonly predictors. \bjournalGeophysical Research Letters \bvolume39 \bpagesL06703. \biddoi=10.1029/2012GL051086 \bptokimsref\endbibitem
 Werner, Luterbacher and Smerdon (2013) {barticle}[author] \bauthor\bsnmWerner, \bfnmJohannes P.\binitsJ. P., \bauthor\bsnmLuterbacher, \bfnmJuerg\binitsJ. \AND\bauthor\bsnmSmerdon, \bfnmJason E.\binitsJ. E. (\byear2013). \btitleA pseudoproxy evaluation of Bayesian hierarchical modeling and canonical correlation analysis for climate field reconstructions over Europe. \bjournalJ. Climate \bvolume26 \bpages851–867. \biddoi=10.1175/JCLID1200016.1 \bptokimsref\endbibitem
 Yang et al. (2008) {barticle}[mr] \bauthor\bsnmYang, \bfnmXiaojun\binitsX., \bauthor\bsnmXing, \bfnmKeyi\binitsK., \bauthor\bsnmShi, \bfnmKunlin\binitsK. \AND\bauthor\bsnmPan, \bfnmQuan\binitsQ. (\byear2008). \btitleJoint state and parameter estimation in particle filtering and stochastic optimization. \bjournalJ. Control Theory Appl. \bvolume6 \bpages215–220. \biddoi=10.1007/s1176800860782, issn=16726340, mr=2418951 \bptokimsref\endbibitem
 Zhang, Mann and Cook (2004) {barticle}[author] \bauthor\bsnmZhang, \bfnmZhihua\binitsZ., \bauthor\bsnmMann, \bfnmMichael E.\binitsM. E. \AND\bauthor\bsnmCook, \bfnmEdward R.\binitsE. R. (\byear2004). \btitleAlternative methods of proxybased climate field reconstruction: Application to summer drought over the conterminous United States back to AD 1700 from treering data. \bjournalHolocene \bvolume14 \bpages502–516. \biddoi=10.1191/0959683604hl727rp \bptokimsref\endbibitem