Climate of the Past Extracting a common high frequency signal from Northern Quebec black spruce tree-rings with a Bayesian hierarchical model

Abstract. One basic premise of dendroclimatology is that tree rings can be viewed as climate proxies, i.e. rings are assumed to contain some hidden information about past climate. From a statistical perspective, this extraction problem can be understood as the search of a hidden variable which represents the common signal within a collection of tree-ring width series. Classical average-based techniques used in dendrochronology have been applied to estimate the mean behavior of this latent variable. Still, depending on tree species, regional factors and statistical methods, a precise quantification of uncertainties associated to the hidden variable distribution is difficult to assess. To model the error propagation throughout the extraction procedure, we propose and study a Bayesian hierarchical model that focuses on extracting an inter-annual high frequency signal. Our method is applied to black spruce (Picea mariana)tree-rings recorded in Northern Quebec and compared to a classical averagebased techniques used by dendrochronologists ( Cook and Kairiukstis, 1992).


Dendrochronology
In our changing climate, the search for accurate information about the past remains essential to understand and link past, present and future climate variations.Direct measurements are missing beyond the length of instrumental records and proxies are necessary to reconstruct chronologies of Correspondence to: J.-J.Boreux (jj.boreux@ulg.ac.be) past temperatures and precipitation for a given region and/or period for which direct observations are unavailable.One of the most widely known proxy consists in tree-ring widths that possess good skill in representing climate information at the interannual to decadal time scale.An overview on this topic can be found in Cook and Kairiukstis (1992).The fundamental assumption in dendroclimatology is that a climatic signal can be hidden into tree-ring growths.Since the pioneering work of Douglass (1936), dendrochronologists have developed various methods to extract such common signals for different species.A required step in dendrochronology, called standardization, is classically needed to transform ring-width series, that are non-stationary due to tree aging processes, into relative tree-ring indices with unit mean and a constant variance.This can be accomplished by dividing each measured ring width by its expected value, i.e. the growth trend is modelled as a regression function of tree ages.Then a common signal is derived by averaging the ensemble of such tree-ring indices across series for each year.Several methods exist to calculate indices averages (e.g., Melvin et al., 2007).Esper et al. (2002) noticed that the low-frequency climate component can be highly sensitive to the standardization method.Recently, Nicault et al. (2009) proposed a neural network approach to remove the age effect and to estimate regional growth curve via explanatory variables such as tree age and their productivity.They developed a standardization procedure to preserve long-term fluctuations.
In contrast with these past methods, our goal in this paper is neither to reconstruct a series of temperatures or precipitation, nor to propose novel regression schemes based on well-chosen explanatory variables as in Nicault et al. (2009).We prefer to focus on the problem of extracting a common inter-annual high frequency signal from a given tree species and region, without regressing on possible predictants.One Published by Copernicus Publications on behalf of the European Geosciences Union.
reason for such a choice is based on the intrinsic difficulties in linking tree-ring growths to specific explanatory variables and in interpreting these relationships.Depending on the tree species under study, it is not always clear to dendrochronologists, even today, what are the precise contributions of precipitation, temperatures, soil and hydrological characteristics, competition and other factors, to tree-ring growths.This is particularly true for black spruces in Northern Quebec.Long and reliable instrumental records of precipitation and temperature are not available for this region.By bypassing this selection, our strategy is to let the raw data "speak" for themselves.Of course, our extracted common signal could be interpreted with respect to local measurements of temperatures, precipitation and other hydroclimatological variables, whenever such information would be available.Hence, independently of the estimation step, explanatory variables could be employed in a validation scheme.To some extent, our extraction strategy could be thought within a "blind experiments" framework.The latter is classically used in medical studies to remove the experimenter bias.In dendrochronology, the "experimenter" could be viewed as the dendroclimatologist who can select his/her favorite explanatory variables (precipitation, temperatures, soil information, etc.).In our statistical modeling, neither climatological nor environmental data (besides tree rings) are used in the analysis.Consequently, our hidden signal should not be influenced by the choice of the "experimenter".Ideally, the extracted common signal should be then linked to climate variables by independent researchers who did not participate in our data analysis.Such a reasoning about design experiments is very common in many research fields (medical studies, nuclear physics, etc.) but it is rarely a topic of interest in climatology.Here we believe that our blind extraction represents an advantage because it can reduce the subjective link between the extracted common tree signal and the climate.

Bayesian Hierarchical Modeling
Assessing uncertainties in any statistical dendrochronological procedure has to be carefully addressed.To tackle this important statistical issue, we opt to work within a Bayesian Hierarchical Modeling (BHM) framework.The main idea of BHMs is to statistically model a complex process and its relationships to observations in several simple components throughout a hierarchy of layers.BHMs handle efficiently the uncertainty assessment of each layer by clearly identifying prior and posterior distributions of underlining processes.Schematically, the prior corresponds to a probability distribution representing knowledge or belief about an unknown quantity a priori, that is, before any data have been observed.Then, in the light of relevant data, the prior probability is updated via Bayes' theorem and becomes the posterior.For an introduction to such models, see e.g.Gelman et al. (2003).In environmental sciences, BHM has become more and more popular during the last two decades.For example, Berliner et al. (2000) studied long-lead predictions of Pacific Sea Surface Temperatures via Bayesian Dynamic Modeling.Cooley et al. (2005) implemented a BHM to infer glacial retreats in Bolivia using lichen growths as a proxy.Cooley et al. (2007) estimated extreme precipitation return levels by combining BHM and extreme value theory.Concerning dendrochronology, Hooten and Wikle (2007) recently investigated with a BHM shifts in the spatio-temporal growth dynamics of shortleaf pine.
The uncertainty in BHM is spread over different layers, usually three.The base level, called the data layer, characterizes observations, e.g.tree ring areas in our case.The second level in the hierarchy, called the process layer, models the latent process that drives the growth of such rings, tree-to-tree and regional variations.In this second layer, one can start incorporating temporal processes, e.g. the tree memories.The third level, called the parameter layer, consists of the information concerning prior parameters distributions that control the latent process.
What is the interest of BHMs for dendrochronologists?The choice of the Bayesian paradigm allows the use of unobserved variables in a hierarchical structure, while easily modeling uncertainties at each different level of this structure.In particular, expert information can be integrated via probability densities (the priors).In other words, past knowledge, even diffuse or imperfect, from scientists can be taken advantage of.More precisely, each parameter of a Bayesian hierarchical model can be viewed as a random variable and hence, a dialogue with dendrochronologists can be engaged to set the prior distribution of this random variable.If the expert has no prior knowledge then the distribution is set to be very wide (a diffuse prior), otherwise the uncertainty of the parameter can be reduced by using knowledge from past studies.In a following step, the incoming data (tree-ring areas here) are used to update all the parameters of our model.The Bayes' theorem provides the mathematical formula to perform this updating, i.e. to derive the posterior distributions.In summary, one can see the above Bayesian strategy as an assembly of elementary parts.Its modular character makes it possible to replace prior uncertainty knowledge (set by experts) by posterior distributional information, throughout the incoming data.In this sense, it is an evolutionary construction.
The paper is organized as follows.Section 2 describes the data and the regional characteristics of the site under study.The details of our latent model are presented in Sect.3. A short discussion about our application is proposed in Sect. 4. Perspectives are given in the conclusion.

Data and region of interest
To extract a common tree signal, the dendrochronologist has to make a series of important decisions about the tree species, the region of interest and the sampling procedure ities.The black spruce (Picea mariana) was selected because it is a widespread species in northern Quebec.Fifteen trees covering a period of 158 years were sampled.These trees were carefully chosen by an expert who removed singular individuals (sick trees, dominated trees, etc.).Each tree provided a ring width series from which annual growth ring areas were estimated.This transformation from ring width to ring area diminishes the geometrical effect impact, basically older trees have thiner rings.The last ring of all sampled trees, albeit missing rings, should correspond to the calendar year.Hence the youngest tree determines the common period length of all trees.
To illustrate the type of dendrochronological times series under study, Figure 2 shows the temporal behavior of three ring area series, randomly chosen from fifteen trees.The right panels represent those three ring area series.From these three right panels, it is clear that each tree has a different trend and it seems difficult to find a common hidden signal in the low frequency domain.In addition, the variability around the cubic-spline trend in the right panels seems to be stronger after 1880 for trees 1 and 2. This example illustrates the high complexity of separating tree ring areas into their individual growth component and their common hidden component in the low frequency part of these signals.Different techniques (e.g.working with residuals after fitting a reference growth curve) exist to deal with this important issue.In this paper we do not address directly this issue.Instead we apply a simple non-parametric transformation to remove trends and to work with stationary time series.This implies that we only focus on inter-annual high frequencies in tree rings.The simple non-parametric transformation is defined as with t = 2, . . ., T and s = 1, . . ., S and where X ts represent the measured annual ring area produced during year t by tree s and T is the length of the temporal sequence and S the number of trees.Transformation (1) is extensively used in finance (Gencay et al., 2002).Besides its simplicity of implementation, this log-difference has the advantage of removing any smooth (i.e.polynomial) trend, see the right panels of Figure 2. In addition the change of variability aforementioned in trees 1 and 2 is less pronounced in the right panels of Figure 2. The drawbacks of using (1) are that, if present, the low frequency part of a possible common signal has been removed and that the time unit t in Y ts does not correspond to a year anymore but to a one-year increment.The latter has to be kept in mind when interpreting our results.The former implies that our model described below will only focus on the high frequency part of a possible common signal.
Concerning the interpretation of Y ts defined as a logdifference between two consecutive ring area values, the following simple facts need to be recalled.Whenever the relative ratio of two inter-annual consecutive ring areas from the same tree is close to one, then Y ts is close to zero.If this relative ratio is very large (ie the ring area from year t is much larger than the one formed during year t − 1), then Y ts has to strongly positive.Conversely, a negative Y ts represents a large decrease in ring areas between two consecutive years.As exemplified by Figure 3, working with Y ts instead of the raw ring areas X ts allows us to remove long-term trends, to focus on the inter-annual relative variability and to work with time series that can be assumed to be stationary and Gaussian.One drawback is that we have lost the absolute value of X ts , ie working with the couple (X ts , X t−1,s ) is equivalent to analyzing the couple (aX ts , aX t−1,s ) for any a > 0, independently of the value of a. Keeping in mind this drawback and those advantages, the correlation meaning in Y ts and Z t can be viewed as the short term memory in the relative log- (e.g., George et al., 2008).Concerning the region choice, Hydro-Quebec, one of the founding agencies involved in this project, has had a strong interest in Northern Quebec because of its hydro-electrical capacities.With this constraint in mind, a mesic site, i.e. with a moderate supply of moisture, close to lake Hurault (54 • 15 N, 70 • 47 W) was chosen, see the red star called HM-1 in the right panel of Fig. 1.This site has the advantages to belong to a climatic homogenous region and of being far away from most human activities.The black spruce (Picea mariana) was selected because it is a widespread species in Northern Quebec.Fifteen trees covering a period of 158 years were sampled.These trees were carefully chosen by an expert who removed singular individuals (sick trees, dominated trees, etc.).Each tree provided a ring width series from which annual growth ring areas were estimated.This transformation from ring width to ring area diminishes the geometrical effect impact, basically older trees have thiner rings.The last ring of all sampled trees, albeit missing rings, should correspond to the calendar year.Hence the youngest tree determines the common period length of all trees.
To illustrate the type of dendrochronological times series under study, Fig. 2 shows the temporal behavior of three ring area series, randomly chosen from fifteen trees.The right panels represent those three ring area series.From these three right panels, it is clear that each tree has a different trend and it seems difficult to find a common hidden signal in the low frequency domain.In addition, the variability around the cubic-spline trend in the right panels seems to be stronger after 1880 for trees 1 and 2. This example illustrates the high complexity of separating tree ring areas into their individual growth component and their common hidden component in the low frequency part of these signals.Different techniques (e.g.working with residuals after fitting a reference growth curve) exist to deal with this important issue.In this paper we do not address directly this issue.Instead we apply a simple Boreux et al.: Extracting a common signal in tree-rings 3 This site has the advantages to belong to a climatic homogenous region and of being far away from most human activities.The black spruce (Picea mariana) was selected because it is a widespread species in northern Quebec.Fifteen trees covering a period of 158 years were sampled.These trees were carefully chosen by an expert who removed singular individuals (sick trees, dominated trees, etc.).Each tree provided a ring width series from which annual growth ring areas were estimated.This transformation from ring width to ring area diminishes the geometrical effect impact, basically older trees have thiner rings.The last ring of all sampled trees, albeit missing rings, should correspond to the calendar year.Hence the youngest tree determines the common period length of all trees.
To illustrate the type of dendrochronological times series under study, Figure 2 shows the temporal behavior of three ring area series, randomly chosen from fifteen trees.The right panels represent those three ring area series.From these three right panels, it is clear that each tree has a different trend and it seems difficult to find a common hidden signal in the low frequency domain.In addition, the variability around the cubic-spline trend in the right panels seems to be stronger after 1880 for trees 1 and 2. This example illustrates the high complexity of separating tree ring areas into their individual growth component and their common hidden component in the low frequency part of these signals.Different techniques (e.g.working with residuals after fitting a reference growth curve) exist to deal with this important issue.In this paper we do not address directly this issue.Instead we apply a simple non-parametric transformation to remove trends and to work with stationary time series.This implies that we only focus on inter-annual high frequencies in tree rings.The simple non-parametric transformation is defined as with t = 2, . . ., T and s = 1, . . ., S and where X ts represent the measured annual ring area produced during year t by tree s and T is the length of the temporal sequence and S the number of trees.Transformation ( 1) is extensively used in finance (Gencay et al., 2002).Besides its simplicity of implementation, this log-difference has the advantage of removing any smooth (i.e.polynomial) trend, see the right panels of Figure 2. In addition the change of variability aforementioned in trees 1 and 2 is less pronounced in the right panels of Figure 2. The drawbacks of using (1) are that, if present, the low frequency part of a possible common signal has been removed and that the time unit t in Y ts does not correspond to a year anymore but to a one-year increment.The latter has to be kept in mind when interpreting our results.The former implies that our model described below will only focus on the high frequency part of a possible common signal.
Concerning the interpretation of Y ts defined as a logdifference between two consecutive ring area values, the following simple facts need to be recalled.Whenever the relative ratio of two inter-annual consecutive ring areas from the same tree is close to one, then Y ts is close to zero.If this relative ratio is very large (ie the ring area from year t is much larger than the one formed during year t − 1), then Y ts has to strongly positive.Conversely, a negative Y ts represents a large decrease in ring areas between two consecutive years.As exemplified by Figure 3, working with Y ts instead of the raw ring areas X ts allows us to remove long-term trends, to focus on the inter-annual relative variability and to work with time series that can be assumed to be stationary and Gaussian.One drawback is that we have lost the absolute value of X ts , ie working with the couple (X ts , X t−1,s ) is equivalent to analyzing the couple (aX ts , aX t−1,s ) for any a > 0, independently of the value of a. Keeping in mind this drawback and those advantages, the correlation meaning in Y ts and Z t can be viewed as the short term memory in the relative log- non-parametric transformation to remove trends and to work with stationary time series.This implies that we only focus on inter-annual high frequencies in tree rings.The simple non-parametric transformation is defined as with t=2, . .., T and s=1, . .., S and where X ts represent the measured annual ring area produced during year t by tree s and T is the length of the temporal sequence and S the number of trees.Transformation (1) is extensively used in finance (Gencay et al., 2002).Besides its simplicity of implementation, this log-difference has the advantage of removing any smooth (i.e.polynomial) trend, see the right panels of Fig. 2.
In addition the change of variability aforementioned in trees 1 and 2 is less pronounced in the right panels of Fig. 2. The drawbacks of using (1) are that, if present, the low frequency part of a possible common signal has been removed and that the time unit t in Y ts does not correspond to a year anymore but to a one-year increment.The latter has to be kept in mind when interpreting our results.The former implies that our model described below will only focus on the high frequency part of a possible common signal.
Concerning the interpretation of Y ts defined as a logdifference between two consecutive ring area values, the following simple facts need to be recalled.Whenever the relative ratio of two inter-annual consecutive ring areas from the same tree is close to one, then Y ts is close to zero.If this relative ratio is very large (i.e. the ring area from year t is much larger than the one formed during year t−1), then Y ts has to strongly positive.Conversely, a negative Y ts represents a large decrease in ring areas between two consecutive years.As exemplified by Fig. 3, working with Y ts instead of the raw ring areas X ts allows us to remove long-term trends, to focus on the inter-annual relative variability and to work with time series that can be assumed to be stationary and Gaussian.One drawback is that we have lost the absolute value of X ts , i.e. working with the couple (X ts , X t−1,s ) is equivalent to analyzing the couple (aX ts , aX t−1,s ) for any a>0, independently of the value of a. Keeping in mind this drawback and those advantages, the correlation meaning in Y ts and Z t can be viewed as the short term memory in the relative logtransfom rate between two consecutive ring areas.
Before closing this section we would like to emphasize that our detrending choice represented by ( 1) is not unique and others techniques could be used to provide stationary signals.For example, we could have worked with the residuals obtained from the cubic spline fit shown in the left panels of Fig. 2.

An additive latent model
The random variable Y ts defined by Eq. ( 1) is assumed to follow an additive model with a latent variable Z t with t=2, . .., T and s=1, . .., S, and where µ s corresponds to the mean level of tree s, Z t represents the hidden regional signal common to all trees and st describes local fluctuations of tree s during year t.Tree-to-tree variations captured by ε st can be due to reserves accumulated by tree s and other factors that are not directly linked to environmental causes, the latter ones should be represented by Z t .For each calendar year t, the product λ s Z t measures how the hidden factor Z t contributes to the growth of tree s.We assume that Z t and ε st are independent processes.With respect to the BHMs described in Sect.1.2, the random variables Y ts corresponds to the data layer and Z t belongs to the process layer.Before describing the probabilistic structure within Z t and st , it is advantageous to rewrite model (2) with obvious vectorial notations where 1 is the unit vector of length T −1.Each tree s may have a temporal memory that should depend on the hydrological stress or other conditions that are particular to this tree location.Although these tree-to-tree effects can be complex, to keep the inference simple and the risk of over-parametrization low, we opt for a simple zero-mean Gaussian auto-regressive process of order one for , i.e.
s =φ s −s +V s .The notation −s corresponds to s shifted by one year, i.e. −s = ε s1 , . .., ε s(T −1) , φ s represents the auto-regressive coefficient of tree s, and the random vector V s of length T −1 follows a zero-mean multivariate Gaussian distribution with precision η s ×I where I is the identity matrix of size T −1.In other words, all components of vector V s correspond to a standardized normal independent random noise.
To allow the common regional factor Z t to have a short year-to-year memory, we assume that the latent Z t can be modeled as a zero-mean Gaussian auto-regressive process of order one, i.e.Z=ρZ − +U where Z − = Z 1 , . .., Z (T −2) and U represents a zero-mean multivariate normal vector of length T −1 with precision τ ×I.
Our full model counts 2+4 S parameters, namely (ρ, τ ) and θ s = (λ s , µ s , φ s , η s ) with s=1, 2, . .., S. We assume that the priors distributions [ρ, τ ], [θ 1 ] , . .., and [θ S ] are mutually independent.By writing the joint distribution as a product of conditional distributions with a marginal distribution, the prior for (ρ, τ ) can take the following form . In a classical way, we assume that the precision parameter τ follows a gamma distribution with two hyperparameters that must be fixed to reflect prior beliefs.In our application, a diffuse prior is chosen by setting the two gamma parameters to zero.
The choice of the auto-regressive coefficient prior [ρ|τ ] is more delicate.Classically, it is assumed that auto-regressive processes are a priori stationary.This implies that autoregressive coefficients have to belong to the interval [−1, 1].As Bayesian statisticians, we defend idea that the underlying characteristics of the hidden process Z t should not be imposed but arise form the data via the Bayes' rule or via prior knowledge.For this reason, we assume that [ρ|τ ] follows a zero-mean Gaussian distribution with a precision proportional to τ .This multiplicative factor must be fixed between zero and one, mainly to degrade the precision a little.In our application, we work with a diffuse prior by equaling the multiplicative factor to zero.
Concerning the prior of the random vector θ s = (λ s , µ s , φ s , η s ), we assume conditional independence, i.e. [θ s ] = [λ s |η s ] [µ s |η s ] [φ s |η s ] [η s ] where the variable η s follows a gamma distribution with two hyperparameters (set to zero in our application).The distributions [λ s |η s ], [µ s |η s ] and [φ s |η s ] are assumed to be diffuse Gaussian priors in this paper.As for the auto-regressive coefficient of Z t , this means that the auto-regressive coefficient of st are not a priori assumed to be in the interval [−1, 1].
To compute the posteriors of the latent vector Z t and of the 2+4 S parameters, we implement the Gibbs sampler described in the Appendix.The Bayesian inference was carried out with the open source R statistical software (R Development Core Team, 2009).Our programs are available upon request.

Results and discussion
The solid line in Fig. 3 shows the estimated posterior median value of the common factor Z t over the period 1846-2003.The shaded area corresponds to the 90% credible regions (CR).Note that the value of Z t and λ s are estimated up to a constant because it is always possible in (2) to multiply Z t by a constant and divide the λ s by the same constant tracted signals makes us believe that our BHM approach is capable of providing meaningful outputs for dendrochronologists because they do not contradict past results and offer another statistical approach to this community of scientists.Concerning the memory within Z t , the posterior distribution The dashed line represents the so-called tree-growth index which is an arithmetic mean of ratios over all trees.Each ratio is derived by dividing ring thickness over a temporally smoothed tree signal for each tree (e.g., Cook and Kairiukstis, 1992).
of the autoregressive coefficient ρ indicates a negative correlation because its 25%, 50% and 75% posterior quantiles are equal to −0.40, −0.36 and −0.32, respectively.In addition to CRs, our methods allow the practitioner to derive a finer analysis of her/his tree ring data.For example, an analysis tree-by-tree can be undertaken.For each of the fifteen trees, Figure 4 displays the posterior mean and 90% CRs of the parameters µ s , λ s and φ s , respectively.The mean posterior value of µ s mostly oscillates around zero for all trees.Overall, each tree but tree 2 appears to have a mild negative inter-annual memory, all autoregressive coefficients (but tree 2) shown in the bottom panel of Figure 4 have a φ s posterior median around -0.4.The central panel clearly points out tree 1 which seem to contribute the most to Z t .To check the quality of our estimation, Figure 5 displays for trees 1, 2 and 3 (shown in Figure 2), the observed Y ts versus the naive estimate Ŷts obtained by plugging our median posterior parameter values in (2) without noise.As expected, the relationships appear to be linear.The same result holds for the other trees.Fig. 4. For each of the fifteen trees, the posterior mean of µ s , λ s and ρ s are represented by circles in the top, middle and bottom panels, respectively.Vertical bars correspond to the 90% credible regions.
Fig. 5.For each of the three randomly chosen trees described in Figure 2, the observed Y ts versus its estimator from model ( 2) is plotted.The dashed line represents the so-called tree-growth index which is an arithmetic mean of ratios over all trees.Each ratio is derived by dividing ring thickness over a temporally smoothed tree signal for each tree (e.g., Cook and Kairiukstis, 1992).
without being able to identify this multiplicative factor.In Fig. 3, we compare our BHM results with a classical technique employed by dendrochronologists.The output of this procedure is represented by the dashed line, a so-called treegrowth index which is an arithmetic mean of ratios over all trees.Each ratio is derived by dividing ring thickness over a temporally smoothed tree signal for each tree (e.g., Cook and Kairiukstis, 1992).Up to a constant (this explains the two different scales for the y-axis), the classical tree-growth index behaves similarly to Z t by staying in the CR over a long time period.From about 1875 to 1900, there is a discrepancy between Z t and the classical tree-growth index, the latter producing higher values during this period.Although fairly localized in time, this difference indicates that this classical technique by not providing confidence intervals shows its limitations.Still, this comparison between the two extracted signals makes us believe that our BHM approach is capable of providing meaningful outputs for dendrochronologists because they do not contradict past results and offer another statistical approach to this community of scientists.
Concerning the memory within Z t , the posterior distribution of the autoregressive coefficient ρ indicates a negative correlation because its 25%, 50% and 75% posterior quantiles are equal to −0.40, −0.36 and −0.32, respectively.In addition to CRs, our methods allow the practitioner to derive a finer analysis of her/his tree ring data.For example, an analysis tree-by-tree can be undertaken.For each of the fifteen trees, Fig. 4 displays the posterior mean and 90% CRs of the parameters µ s , λ s and φ s , respectively.The mean posterior value of µ s mostly oscillates around zero for all trees.Overall, each tree but tree 2 appears to have a mild negative inter-annual memory, all autoregressive coefficients (but tree 2) shown in the bottom panel of Fig. 4  and Kairiukstis, 1992).Up to a constant (this explains the two different scales for the y-axis), the classical tree-growth index behaves similarly to Z t by staying in the CR over a long time period.From about 1875 to 1900, there is a discrepancy between Z t and the classical tree-growth index, the latter producing higher values during this period.Although fairly localized in time, this difference indicates that this classical technique by not providing confidence intervals shows its limitations.Still, this comparison between the two extracted signals makes us believe that our BHM approach is capable of providing meaningful outputs for dendrochronologists because they do not contradict past results and offer another statistical approach to this community of scientists.
Concerning the memory within Z t , the posterior distribution The dashed line represents the so-called tree-growth index which is an arithmetic mean of ratios over all trees.Each ratio is derived by dividing ring thickness over a temporally smoothed tree signal for each tree (e.g., Cook and Kairiukstis, 1992).
of the autoregressive coefficient ρ indicates a negative correlation because its 25%, 50% and 75% posterior quantiles are equal to −0.40, −0.36 and −0.32, respectively.In addition to CRs, our methods allow the practitioner to derive a finer analysis of her/his tree ring data.For example, an analysis tree-by-tree can be undertaken.For each of the fifteen trees, Figure 4 displays the posterior mean and 90% CRs of the parameters µ s , λ s and φ s , respectively.The mean posterior value of µ s mostly oscillates around zero for all trees.Overall, each tree but tree 2 appears to have a mild negative inter-annual memory, all autoregressive coefficients (but tree 2) shown in the bottom panel of Figure 4 have a φ s posterior median around -0.4.The central panel clearly points out tree 1 which seem to contribute the most to Z t .To check the quality of our estimation, Figure 5 displays for trees 1, 2 and 3 (shown in Figure 2), the observed Y ts versus the naive estimate Ŷts obtained by plugging our median posterior parameter values in (2) without noise.As expected, the relationships appear to be linear.The same result holds for the other trees.Fig. 4. For each of the fifteen trees, the posterior mean of µ s , λ s and ρ s are represented by circles in the top, middle and bottom panels, respectively.Vertical bars correspond to the 90% credible regions.To check the quality of our estimation, Fig. 5 displays for trees 1, 2 and 3 (shown in Fig. 2), the observed Y ts versus the naive estimate Ŷts obtained by plugging our median posterior parameter values in (2) without noise.As expected, the relationships appear to be linear.The same result holds for the other trees.

Conclusions
To summarize our findings, we have implemented a hierarchical Bayesian model to estimate a common hidden signal in high frequency component of trees.This latent signal should be viewed as a representation of the regional pressure affecting black spruce trees over our studied area in Northern Quebec.The hierarchical structure provides another way to model the temporal structure associated to tree memories at the regional and tree-to-tree levels.This model attempts to quantify the contribution of a high frequency common hidden signal to each tree growth.This could help selecting trees with regard to a possible climatological interpretation in a reconstruction context.Compared with past approaches, our hidden signal was strongly correlated to the estimate obtained with the most traditional procedure.This confirms a past method derived by dendrochronologists, while bringing the benefits of a BHM approach.As a further step in this analysis, it would be of interest to integrate low frequency respectively.Vertical bars correspond to the 90% credible regions.   in Eq. ( 2).One possibility is to bypass transformation (1) by making the term µ s in (2) varying in time.For example, µ ts could be modeled by Bayesian splines.Besides the complexity of such an approach, the main difficulty is our limited sample size (fifteen trees).Another aspect is the handling of missing values and consequently avoid the limitation brought by the age of the youngest tree.In addition, ongoing field trips should provide a much larger sample of tree rings and allows us to extend our BHM procedure in future research.In this context, our present work should rather be viewed as an addition of a simple statistical procedure to the mathematical toolbox of dendroclimatologists rather than a comprehensive study of black spruce trees in Northern Quebec.
The combination of the additive model described by (2) and the Bayesian paradigm allows the practitioner to easily generate the full posterior distribution of the hidden signal, and consequently it is possible to simulate realizations of the relative log-transform ratio between two consecutive ring areas.Such simulations could help simulating the way tree growth in response to climatic forces that drives the common inter-annual variations.Another interesting perspective of the Bayesian approach resides in the possibility to compute predictive posterior densities for future years.Since we can derive the posterior density of the extracted signal, the predictive posterior [z t+1 |y t , θ] for an unobserved year t+1 can obtained by computing the hidden state posterior at time t+1.where k ψ et m ψ are prior parameters which are invariant from tree to tree (e.g.k ψ =m ψ =0).The vectors f s and g s depend on handling parameters.

-Step
3.2: Draw precision η s |µ s , λ s , ϕ s , y s , z, z 0 , y 0s from a gamma distribution with parameters c+ T +3 2 and [d+0.5vv] −1 where c and d are prior parameters (e.g.c=d=0) -Step 4.: Draw vector U|τ, η s , L s , R s from a multivariate normal distribution with mean ω and covariance −1 and set z=z 0 B 0 +Bu.The mean ω and matrix −1 relate vector L and matrix R which depend on previous parameters.
Note that when the Gamma hyper parameters are theoretically equal to zero, this means that they are set to a very small value in real computations, e.g.a=b=.001 in our case.

Fig. 1 .
Fig. 1.The upper panel corresponds to northern Quebec.The lower panel is a zoom near the Caniapiscau region and the red star called HM-1 represents the site from which fifteen trees have been sampled.

Fig. 2 .
Fig. 2. Temporal behavior of three ring area time series (randomly chosen from a set of fifteen trees) over the period 1846-2003.The left panels correspond to the measured tree ring areas with a fitted cubic spline trend.The right panels indicate the log difference of the same ring areas, see Equation (1).

Fig. 1 .
Fig. 1.The upper panel corresponds to Northern Quebec.The lower panel is a zoom near the Caniapiscau region and the red star called HM-1 represents the site from which fifteen trees have been sampled.

Fig. 1 .
Fig. 1.The upper panel corresponds to northern Quebec.The lower panel is a zoom near the Caniapiscau region and the red star called HM-1 represents the site from which fifteen trees have been sampled.

Fig. 2 .
Fig. 2. Temporal behavior of three ring area time series (randomly chosen from a set of fifteen trees) over the period 1846-2003.The left panels correspond to the measured tree ring areas with a fitted cubic spline trend.The right panels indicate the log difference of the same ring areas, see Equation (1).

Fig. 2 .
Fig. 2. Temporal behavior of three ring area time series (randomly chosen from a set of fifteen trees) over the period 1846-2003.The left panels correspond to the measured tree ring areas with a fitted cubic spline trend.The right panels indicate the log difference of the same ring areas, see Eq. (1).

Fig. 3 .
Fig. 3.The solid line corresponds to the estimated posterior median value of the common signal Z t from (2) over the period 1846-2003.The shaded area corresponds to the 90% credible regions.The dashed line represents the so-called tree-growth index which is an arithmetic mean of ratios over all trees.Each ratio is derived by dividing ring thickness over a temporally smoothed tree signal for each tree (e.g.,Cook and Kairiukstis, 1992).

Fig. 3 .
Fig. 3.The solid line corresponds to the estimated posterior median value of the common signal Z t from (2) over the period 1846-2003.The shaded area corresponds to the 90% credible regions.The dashed line represents the so-called tree-growth index which is an arithmetic mean of ratios over all trees.Each ratio is derived by dividing ring thickness over a temporally smoothed tree signal for each tree (e.g.,Cook and Kairiukstis, 1992).
have a φ s posterior Boreux et al.: Extracting a common signal in tree-rings 5

Fig. 3 .
Fig. 3.The solid line corresponds to the estimated posterior median value of the common signal Z t from (2) over the period 1846-2003.The shaded area corresponds to the 90% credible regions.The dashed line represents the so-called tree-growth index which is an arithmetic mean of ratios over all trees.Each ratio is derived by dividing ring thickness over a temporally smoothed tree signal for each tree (e.g.,Cook and Kairiukstis, 1992).

Fig. 5 .
Fig. 5.For each of the three randomly chosen trees described in Figure2, the observed Y ts versus its estimator from model (2) is plotted.

Fig. 4 .
Fig. 4.For each of the fifteen trees, the posterior mean of µ s , λ s and ρ s are represented by circles in the top, middle and bottom panels, respectively.Vertical bars correspond to the 90% credible regions.

Fig. 5 .
Fig. 5.For each of the three randomly chosen trees described in Figure2, the observed Y ts versus its estimator from model (2) is plotted.

Fig. 5 .
Fig. 5.For each of the three randomly chosen trees in Fig.2, the observed Y ts versus its estimator from model (2) is plotted.