Autoregressive Modeling of Forest Dynamics

11/20/2019 ∙ by Olga Rumyantseva, et al. ∙ Washington State University University of Nevada, Reno 0

In this work, we employ autoregressive models developed in financial engineering for modeling of forest dynamics. Autoregressive models have some theoretical advantage over currently employed forest modeling approaches such as Markov chains and individual-based models, as autoregressive models are both analytically tractable and operate with continuous state space. We perform time series statistical analysis of forest biomass and basal area recorded in Quebec provincial forest inventories in 1970-2007. The geometric random walk model adequately describes the yearly average dynamics. For individual patches, we fit an AR(1) process capable to model negative feedback (mean-reversion). Overall, the best fit also turns out to be geometric random walk, however, the normality tests for residuals fail. In contrast, yearly means are adequately described by normal fluctuations, with annual growth, on average, 2.3 with standard deviation of order 40 uneven number of observations per year. This work demonstrates that autoregressive models represent a valuable tool for modeling of forest dynamics. In particular, they quantify stochastic effects of environmental disturbances and develop predictive empirical models on short and intermediate temporal scales.



There are no comments yet.


page 10

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

1.1. Background

Understanding the dynamics and self-organization of ecosystems is one of the most challenging problems in modern ecology [19]. Self-organization occurs simultaneously on several levels of hierarchical ecosystem organization and involves dynamic processes operating on different temporal and spatial scales [21]

. Forest dynamics refers to temporal and spatial changes that occur simultaneously at different levels of ecosystem organization. Various modeling approaches employed to understand and predict these changes include a number of discrete and continuous stochastic and deterministic models, such as Markov chains, individual-based models, ordinary and partial differential equations

[41]. A large number of forest models have been developed over the last decades [2, 37, 33, 34, 31, 40]. Still, there are fundamental questions that existing quantitative approaches fail to fully address. One of the major challenges is the understanding of forest succession and biomass dynamics under the non-stationary disturbance regimes related to climatic factors and anthropogenic activities. An incomplete list of disturbances that substantially affect tree survival and lead to forest biomass decrease include wind, frost, hurricanes, harvesting and forest fires. Markov chains are able to capture effects of all these disturbances on forest biomass dynamics [28]. However, their application is based on the discretization of the state space [39]. Spatially-explicit individual-based models are able to simulate effects of disturbances in continuous time and state space [32, 40, 41]. However, these models are typically analytically intractable, i.e. the model predictions are produced by computer simulations only, and it limits our ability to understand what model prediction in general [40, 26].

Another challenge to our understanding of forest dynamics is the multidimensional nature of this complex adaptive system [21, 27]. Forest disturbances are traditionally associated with a loss of biomass. However, Markov chain models based only on biomass do not capture forest succession comprehensively [39, 27]. This limitation motivates the need for alternative formulations that are able to consider several forest dimensions instead of only one. In our previous study, we combined multivariate statistical analysis with Markov chain approach to develop a multidimensional Markov Chain [27]. However, simultaneous discretization of several independent forest characteristics of a different nature substantially reduced our ability to understand the ecological meaning of model predictions, which was one of the major advantages of the Markov Chain approach [34, 5, 39, 28].

1.2. Forests and stock markets

In this work, we employ time series (autoregressive models and random walk) to quantify disturbance regimes and to build a predictive stochastic model of multiple disturbance classes. This type of models can overcome both major shortcomings of previous models outlined above. Autoregressive models operate with continuous space, they are analytically tractable, and they can operate with multidimensional characteristics of complex adaptive systems. Similar approaches have been successfully applied, for example, in modeling stock market fluctuations. We develop a stochastic theory of forest dynamics using an analogy to stock market theory in financial mathematics. A stock market is another complex system with random fluctuations due to multiple difficult-to-predict factors. Each individual stock has fluctuations with heavy tails. But the total stock market, measured by an index (such as Standard & Poor 500) has long-run fluctuations (3–4 years) which follow Gaussian distribution. These fluctuations depend on various factors, such a price-earnings ratio (this measures whether the stock market is underpriced or overpriced; one can informally think of this as temperature of the stock market) and Treasury rates (long-term and short-term). These factors, in turn, can be modeled as various autoregressive models. Our main idea is to imagine that individual patches behave like stocks, and an average over a particular region is a stock market index.

1.3. Understanding and modeling of forest patch dynamics

In this work we propose to employ autoregressive models to understand and predict forest dynamics at the patch level. The patch-mosaic concept [46] was actively developed in the second half of the twentieth century after suggestion in [45] that forested ecological systems can be considered a collection of patches at different successional stages. Dynamic equilibrium arises at the level of the patch mosaic rather than at the level of individual patches. The classic patch-mosaic methodology assumes that patch dynamics can be represented by changes in macroscopic variables characterizing the state of the patch as a function of time [20]. Conservation law modeling in the case of continuous time and patch state results in the reaction-advection-diffusion model [20].

Patch dynamics concept can be applied for understanding and predicting of forest dynamics at different levels of forest organisation within the hierarchical patch dynamics paradigm [46, 41]. At the level of individual trees patch dynamics concept is traditionally called the forest gap dynamics [37, 36, 30]. Individual-based forest models capture gap dynamics by simulating growth, competition and mortality of individual trees [33, 32, 3, 8]. Individual-based models and analytically-tractable models approximating gap dynamics [18, 40] provide scaling from individual-level dynamics to the next level of forest hierarchical organisation, the stand-level.

In the present work we apply autoregressive models to the stand-level forest patch dynamics. At this level of forest organisation we operate with forest patches (forest stands) consisting of a large number of individual trees [39]. Forest stand is affected not by individual-level tree dynamics, as well as by large-scale disturbances such as forest fires, drought and hurricanes [15, 30], which affect many trees in the stand at the same time. The stand-level dynamics scales up to the next hierarchical level of particular forest type or regional patch mosaic (level 3 in the hierarchical patch mosaic Matreshka model [41]).

The interplay of individual-level and stand-level changes and disturbances creates complex dynamical patterns at the stand-level. One particular source of complexity is related to a large number of intermediate level disturbances affecting only a fraction of trees in the patch [39]. As the consequence of this system complexity a classical linear patch dynamics model does not capture patch dynamics of the US and Canadian forests [39, 23, 28]. This classical patch-dynamics model can be represented in continuous case by advection-reaction patch-dynamics conservation law model [20], and in discrete case by birth and death process that can be written as a Markov Chain [41], or as a simple forest fire model [43]. As the result, in order to capture patch dynamics of the US and Canadian forests, we need to consider more complicated models. In particular, if we discretize forest dynamics with respect to both time and state variable (biomass) we can achieve an adequate representation of forest patch dynamics within the framework of Markov chains [39, 27, 28]. Markov chains provide an analytically tractable representation of forest stand dynamics, while they have a discretization error that is challenging to quantify. This work introduces an autoregressive modeling approach in application to the forest patch dynamics in Quebec. Theoretically, our modeling approach will deliver stochastic and analytically tractable models operating with continuous state-space and -time, without discretization errors.

1.4. Our contributions

Here, we model dynamics of Quebec forests using a traditional AR(1) process borrowed from quantitative finance without modifications. We select the Quebec forest inventory for this proof of concept work as it is a long term dataset collected over more than 3 decades using the same field survey protocol [27]. We operate with the same biomass and basal area data derived from Quebec forest inventories in our previous publication on Markov chain modeling (data-mining protocol is available in [27, 25]

). For both individual patches and the Quebec region, we model logarithms of biomass or basal area as autoregressive process. The best fit, in a certain sense, turns out to be a random walk, with independent increments, which allows to quantify forest disturbance regime overall at the regional level. Regional averages are well-described by normal distribution, while individual patches are not: Fluctuations have heavy tails. This is similar to financial markets, with individual stocks having non-normal fluctuations, and stock indices (in 3–4 years or more) having normal fluctuations. To account for differing number of observations each year, we use Bayesian analysis for annual averages.

2. Materials and Methods

2.1. Data mining of Quebec provincial forest inventories

We base our research on Quebec forest inventory data 1970–2007 Each permanent forest inventory plot has a circular form of 400 . The database consists of 32552 plot re-measurements at 11660 different locations. The Quebec forest inventory is designed to comprehensive describe patch mosaic of Quebec forests and plots cover the Quebec territory practically uniformly. The GIS-based map of forest inventory plots is published as Figure 1 in [28]. Forest inventory plots cover hardwood and mixed forests in the northern temperate zone (9621 and 7663 measurements, respectively) and continuous boreal forests in the boreal zone (11969 measurements). These forest patches (forest inventory plots) are remeasured every few years often with irregular time intervals between measurements. The inventory plots were affected by natural and anthropogenic disturbances including fire and harvesting. The statistical analysis of the measurement dynamics and re-measurement intervals are published in [39] (see Figure 2 in Appendix to [39]). In this inventory, each patch observation includes diameter of each tree, its species, soil moisture, and other characteristics. This is the raw data which is then converted to a more tractable data series. In particular, we are interested in biomass and basal area. Calculations of biomass and basal area were previously done according to [17]. The computations of biomass and basal area (as well as other characteristics, such as shade tolerance index, and biodiversity measured by Shannon entropy) is done in articles [27, 25, 23, 39]. The biomass is this article refers to the plot biomass, which is the sum of biomasses of all trees computed using formulas from [17] (see section 3.1 in [27] for the details). The code used for this article is available on GitHub repository asarantsev/Quebec.

2.2. Autoregressive model for individual forest patches

We propose a new method of modeling the biomass of an individual patch: autoregressive model, when each next year’s logarithm of biomass is a weighted sum of the previous year’s logarithm of biomass and a random Gaussian noise. See the primer on autoregressive models in Appendix B. We measure biomass on a logarithmic scale since it is always positive. That is,


where are constants, and are i.i.d. (independent identically distributed) random variables. If , this sequence exhibits mean reversion to its long-term average . That is, if , then is likely to be smaller than , and vice versa. Examples of earlier use of such models for forest modeling include [11, 22, 12]. They include also spatial models (incorporating distance between patches). We shall not attempt it here, instead treating every patch as effectively isolated. Building a spatial model for Quebec forest is left for future research.

Since data is collected on irregular time intervals, we apply (1) multiple times to itself to get the expression of from if and are consecutive years for which this patch was observed. Then we try various and obtain for each

the maximum likelihood estimate via linear regression. We compare these likelihoods and find the best fit for

. It turns out to be . That is, this sequence actually does not exhibit any mean reversion, but behaves like a random walk, when each next increment is independent from the past:

The biomass itself is a geometric random walk: a process whose logarithm is random walk.

However, the residuals for

do not pass the normality test. This model does not actually fit well, and we cannot find confidence intervals for

using standard statistical techniques. This is due to noise in measurements of individual patches. Later in the article, we find that the average biomass over all patches exhibits more regular behavior, with normal increments.

We perform two versions of this computation: for biomass and for basal area. For each version, we do it in two ways: (a) original logarithms of biomass/basal area; (b) with logarithm of mean biomass/basal area for this year subtracted. In both cases, the maximum likelihood estimate gives us random walk .

The biomass and the basal area are highly dependent, and one can plausibly use one of these metrics instead of the other.

We inherit these techniques from quantitative finance. In particular, the geometric random walk model is a classic model for the stock market movements, going back to classic research by [10]. Mean reversion is commonly observed in financial ratios such as earnings yield or dividend yield, [1, 16]. See also [4, 7, 13] for this research and influence of financial ratios on stock market performance. However, these techniques are less known in mathematical biology.

2.3. Annual averages

We are also interested in each year’s average over all patches. We have only 38 observations for the mean value. Let be the logarithm of this mean. We find that the random walk adequately describes this:

In particular, we find that the increments indeed have normal distribution, and not heavy tails. However, we cannot simply take yearly averages for every year , since they would have different precision. Reason: each year

has a different number of observations. To account for this, we use Bayesian statistics. We put a prior distribution on the values of

and (which corresponds to the lack of any existing information), and then compute the posterior distribution from the likelihood. Bayesian techniques are increasingly used in ecology, [9] as well as in medical statistics [14], and quantitative finance [35].

3. Results and Discussion

We have 32552 observations of 11660 patches of Quebec forests. Each patch is observed at most once a year, during 1970–2007. On average, each patch is observed around 3 times: . Out of these 11660 patches, 10215 have more than one observation, which allows us to model time dynamics. Each observation consists of 4 numbers: Patch ID; year; biomass; and basal area. Various patches are observed in different years: Patch 7000406701 is observed only in 1970; patch 7000406901 is observed in 1970 and 1978; patch 8509702201 is observed in 1985, 1993, 2003; and patch 7000406902 is observed in 1970, 1978, 1985, 1997. Let be the set of observed patches. Each has observations at years , where . Let . Yearly means are defined as:

And we define . The main difficulty is that for almost all patches, a gap between consecutive observations is more than one year, and it differs from patch to patch. In particular, there are 3334 pairs of patch-year observations with the same patch and the gap 8 years, 1923 such pairs with the gap 9 years, but only 66 pairs with gap 22 years. More detailed data is in Appendix D.

(a) Annual Means of Biomass
(b) Annual Means of Basal Area
Figure 1. Annual means of biomass (measured in ) and basal area (measured in )

3.1. Autoregressive model for individual patches

Consider the following time series:


where is an AR(1) parameter, and are i.i.d. . As mentioned earlier, our main difficulty is that we do not have observations for all and , only for . Assume are subsequent time points in . Then


We do this both for (raw data), (centered data), and centered data without and . As discussed above, these years have only very few observations, and we do not have much confidence in these values. We could write the log-likelihood of (3) and apply maximum likelihood estimate using gradient descent. However, since we do not have many data points, we can use an equivalent method which is computationally inefficient but easy to implement: Fix an and run regression with respect to . Then choose an

such that the standard error of this regression is smallest. To properly apply this linear regression, divide (

3) by a constant to make the standard error in error terms in (3) the same for all :


For every which is a multiple of , fit a simple linear regression (without intercept) to find and the standard error

. Appendix A explains why minimizing standard error given normalized residual variance is equivalent to maximizing likelihood. We do this analysis three times: for original values, centered values, and centered values with years 1982 and 2004 removed. We repeat this for biomass and basal area. For all six cases, the parameter

minimizes the standard error (thus maximizing the likelihood). Thus, the dynamics in (2) is given by , which is a simple random walk with Gaussian increments.

(a) Standard error vs
(b) QQ plot for residuals
Figure 2. Non-centered biomass data for individual patches, with measured in

Assuming that the model is, in fact, a random walk, let us make quantile-quantile (QQ) plots for residuals in (

4) (we have simple linear regression without intercept):

These QQ plots of residuals in this regression (4) are not normal. Thus more analysis is needed.

For non-centered biomass, the minimal standard error is achieved when , see Figure 2 (A). This corresponds to a random walk model, and we get and . However, the residuals are not normal: See them on the QQ plot in Figure 2 (B). Then we repeat the computation above for centered values: instead of . Similarly, the standard error is minimized for . For this value, , . The QQ plot of residuals is still not normal. For centered values with years 1982 and 2004 removed, the standard error once again is minimized for , with and . The QQ plot of residuals is still not normal. We omit the standard error graph and the QQ plot for the last two cases: centered data with all years, and centered data without 1982 and 2004, since these plots are very similar to their counterparts for original (non-centered) data.

Modeling basal area data as in (2), and computing standard error of the regression (4), we again get , , and . Again, the best-fitting model among AR(1) according to the maximum likelihood estimation is the random walk. If we center the basal area data, and consider all years, then again, the standard error is minimized for , with and . Centering the basal area data and removing years 1980 and 2004, we get: , , . In all these cases, similarly to the case of biomass, the QQ plots of residuals for basal area are not normal, with both tails fat. We do not provide pictures of QQ plots, since they are very similar to that for biomass. Correlation between biomass and basal area: Take all patches and corresponding years in Denote by the logarithm of biomass, and by the logarithm of basal area for patch and year . If has more than one year, order them in increasing order: Compute correlation coefficient between

It is equal to 0.983. With years 1982 and 2004 removed, this number does not change (in the first four decimal digits). Thus biomass and basal area for individual patches are very correlated. For practical purposes, this means we can use either measure as a size of patch.

(a) Biomass 2007
(b) Basal area 2007
Figure 3. Histograms for biomass (measured in ) and basal area (measured in ) in 2007
(a) Biomass 2019
(b) Basal area 2019
Figure 4. Histograms for biomass (measured in ) and basal area (measured in ) in 2019 superimposed upon 2007

3.2. Yearly averages, frequentist analysis

We model the logarithm of yearly average using random walk:


However, years 1982 and 2004 have only a few observations, see the table from Appendix D. Thus we do not have much confidence in these results. Thus we do the second test: Remove these years, and consider increments (5) for . Then we test normality of increments using the QQ plot and Shapiro-Wilk test. We confirm that, indeed, increments are i.i.d. normal. Rewrite the model (5) as the random walk:


Or, equivalently, we can rewrite (6) in the original scale, instead of the logarithmic scale:


From (7), we can compute the mean and variance of :

For the biomass, we get: and . These increments passes Shapiro-Wilk normality test with . With removed two years 1982 and 2004, the estimates will not change much: and . This still passes Shapiro-Wilk normality test with . Repeating this analysis for basal area instead of biomass, we get: for all years, and for years without 1982 and 2004. Shapiro-Wilk test gives us for all years, and for years without 1982 and 2004. The QQ plots in Figure 5 show that, indeed, the residuals are close to normal.

In Figure 3, we plot histograms of the biomass and basal area for an individual patch in 2007. In Figure 4, we simulate biomass and basal area as in (6) until 2019, starting from a patch randomly selected among observed patches in 2007. We superimpose this histogram upon the one for 2007.

Taking increments of logarithms of yearly means for biomass and basal area (37 data points), we get correlation 0.977. For years 1982 and 2004 removed, we get correlation 0.980. Previously, we got very high correlation between increments of logarithms for individual patches, thus we conclude that biomass and basal area are the same for practical purposes, as measures of size. Now we see that the same is true for yearly averages.

(a) All Years
(b) Without 1982 and 2004
Figure 5. QQ plot for biomass yearly average logarithm increments

3.3. Yearly averages, Bayesian analysis

As mentioned earlier, the analysis in the previous section is deficient: The means have different precision for different years, because the quantity of patches observed differs from year to year. Thus we apply Bayesian statistics in this section.

A primer on Bayesian statistics for normal distribution can be found in Appendix C. Let for each fixed year the logarithms of observations are . For these values we performed analysis above: The posterior mean and the posterior variance were generated. The simulation of and was performed times. In Figure 6, we have histograms of 1000 simulations for and for (year 1970), for both biomass and basal area. Hence we obtain 38 sequences of numbers: , . The average growth rate based on simulated results is


The mean increments are: . Assuming these are i.i.d. we estimate and as in (8) and :

We compute the point estimates of and for biomass and basal area:

Thus the growth rate of forest, measured by the biomass or basal area (on the logarithmic scale), is per year. These numbers are close to the ones from frequentist analysis from the previous subsection: with years 1982 and 2004. However, with removed years 1982 and 2004, this estimate changes to . For basal area, we have a similar comparison with the previous subsection: with 1982 and 2004, without these years. As discussed earlier, we view Bayesian analysis as the more statistically sound. Thus growth per year seems more reasonable.

From year to year, however, there are a lot of fluctuations, or, to use a stock market term, volatility: The standard deviation for yearly fluctuations is estimated as , that is, per year for the biomass, and , that is, per year for the basal area. Similarly to the mean estimates, we view these as more scientifically sound that the ones from frequentist analysis from the previous subsection.

Figure 6. Histograms of simulations for means and variances of the logarithm for biomass and basal area, 1970. The biomass is measured in , and the basal area is measured in .

4. General discussion

4.1. Towards autoregressive theory of forest dynamics

We have applied AR(1) autoregressive model to patch/stand dynamics of Quebec forests. Overall, we followed major steps of application of autoregressive model in financial engineering considering stand biomass and basal area as variables similar to a stock market index. However, forest and stock market are different complex systems despite the observed similarities. We consider this work as the first step, and the further discussion is required: it is hard to expect that the forest dynamics modeling will simply mirror stock market theory.

One interesting result reported in the Section 2.2 is that and the residuals do not pass the normality test. In simple terms, it means that each forest patch is highly volatile, and its dynamics cannot be well described by the normal distribution. Still, among all AR(1), random walk has the maximum likelihood. This suggests for possible use random-walk models with tails heavier than normal. This dichotomy of average patch vs individual patch reminds us of a stock market dynamics, where each individual stock is highly volatile, with non-Gaussian fluctuations [6], but a stock market overall (measured by Standard & Poor 500 or another stock market index) in the long run follows geometric random walk with normal increments (if the time step is large enough, say 3-4 years or more).

From the perspective of random processes means that we are dealing with the divergent AR(1) process, traditionally called the random walk, which does not have a stationary distribution. The AR(1) converges to its stationary distribution when . This means that biomass and basal area at every step will drift away like in both positive and negative directions. We can search for mechanistic meaning of this result in both fundamental patterns of forested ecosystem dynamics and in general modeling assumptions related to the autoregressive modeling approach and the patch dynamics modeling framework.

Patch dynamics modeling framework assumes that each patch, a forest stand in our case, has the same underlying transition probabilities related to internal changes such as tree growth as well as with respect to natural disturbances

[39]. These are typical background assumptions of practically all forest models, however the validity of these assumptions should be evaluated. In particular, in our application we consider the combined data set consisting of hardwood and mixed temperate forest stands as well as continuous boreal forests. In addition, some of the stands even in the same forest types are located in different climatic conditions and might be affected by various disturbance regimes. The overall robustness of modeling predictions with respect to this type of data variability is typically known for traditional models. We have previously investigated effects of these assumptions in the Markov chain modeling framework [39, 27, 28]. Similarly, individual-based models are often applied for forest simulations in different regions with the same growth/mortality characteristics of particular tree species. Autoregressive modeling is a new approach in forest modeling, and it is not yet known how robust our predictions with respect to variability of forest patch mosaic. We can hypothesise that outcomes reported in the Section 2.2 ( and the normality test for residuals) might be related to the variability in forest inventories.

Another common underlying assumption in forest modeling is that the random process is time homogeneous or stationary [39]. It is often assumed that climatic and disturbance regime changes occur at the slower scale and can be ignore for the short- and intermediate-term modeling [27]. Time inhomogeneous Markov chains allow to relax this assumption and take into account both climate driven changes in tree growth and frequency of environmental disturbances [28]. Similarly to inhomogeneous Markov chain theory, autoregressive models can be generalized to include directional changes in disturbance dynamics. This leads us to autoregressive integrated moving average (ARIMA) models, which can capture non-stationary (time inhomogenous) dynamics. In particular, we previously detected some non-stationary effects related to the disturbance regimes recorded in the Quebec forest inventory, such as small decrease of the total disturbance rate in boreal and hardwood forests (see Table 1 in [28]). However, the relatively small number of the measurements in the data set did not allow us to report some definite trend. Nevertheless, we anticipate that the pameterization and validation of ARIMA models can be an important step towards time series theory of forest patch dynamics.

Time series analysis can also be applied to forest dynamics at different hierarchical levels. In a certain way autoregressive models are continuous time and state counterpart of discrete Markov chain models. Markov chains were applied at the individual tree and species levels [44, 38], stand-level [42], and landscape-level [29]. The autoregressive and ARIMA models can be also employed to capture similar dynamics. In addition, time series models might be used in approximation of the forest gap dynamics and individual-based models.

4.2. Future research

That the random walk with Gaussian increments performed better than any other AR(1) suggests that there is not much dependence of annual change of the biomass on the current biomass. However, testing this statement requires further research.

Time-series analysis is a novel modeling approach for forest modeling at the large scale. We wish to compare this approach with state-of-art mathematical models developed on other principles, in particular with Markov chains [39, 27, 25, 23] and individual-based models [41, 26], developed using the same forest inventory data sets. We also need to understand the merits of discrete vs continuous space. This is critical for practical applications of forest models for natural resource management, and risk assessment of forest vulnerability to changes in environmental disturbance regimes.

We can generalize time series approach to model forest dynamics and climate at the regional scale. Such generalization in financial engineering resulted in the development of vector autoregressive integrated moving average (VARIMA) models, which are particularly suited for multidimensional forest modeling under the dynamic disturbance regimes. This can result in an enhanced tolerance-disturbance model build upon our past analyses of tolerance patterns in North American forests, including observed shade tolerance patterns

[23], associations between shade tolerance and soil moisture [25], and relative trade-offs between shade, drought, and water-logging tolerance at the continental scale [24].

We can extend our research to the USA Forest Inventories: Analysis of current disturbance regimes across US ecoregions, and how disturbance regimes are connected with other forest macroscopic characteristics. In particular, forest tolerance is a key ecological driver that controls the response of forested ecosystems to climate change-related disturbances such as drought- and fires [24]. We can link climatic models and forest tolerances [23, 25, 24]. Spatially-explicit information from the USDA Forest Service Forest Inventory and Analysis Program (FIA), Canadian Forest Inventory, and climate data sets and models (WorldClim and PRISM), can be integrated to quantify the effect climate variables in North America have on forest tolerances and disturbances via recently proposed methodology [24].

5. Conclusion

Individual patches have biomass and basal area whose dynamics (on logarithmic scale) are not well described by classic autoregressive model with normal noise terms , because the tails of these noise terms turn out to be heavier than normal. However, among these AR(1) models, the best fit (according to the maximum likelihood) is random walk, with , and increments independent of the past . Thus one can try a heavy-tailed random walk, in which increments have tails heavier than Gaussian. This topic is left for future research.

In contrast, yearly means (on the log scale) are well described by random walk with normal increments. Bayesian analysis accounts for the fact that different years have different number of observations. We get growth rate (measured on the log scale) per year, with standard deviation for biomass, and for basal area. Thus the forest grows on average in the long run, but from year to year there is a lot of volatility.

An important part of forest ecological modeling is to quantify disturbances from fires, droughts, etc. These events quickly destroy significant parts of the forest. An analogy for the stock market would be a crash, as in 2001 or 2008. Indeed, for the stock market, 1-year fluctuations are not adequately described by the normal distribution (only 3–4 years or more are). But for the forest as a whole (as opposed to particular patches), annual fluctuations are normal. We do not need to introduce separate distributions for modeling disturbances. Since we assume the prior is Jeffrey’s non-informative and the model is Gaussian, we do not need to do any Monte Carlo or Metropolis-Hastings computations: There are explicit formulas for the posterior distribution of parameters.


We thank Dr. Adam Erickson for help with proofreading. Initial data-mining of Quebec forest inventories was done by Dr. Jean Lienard for the Markov chain modeling and published in [27]. We also acknowledge welcoming environment for applied quantitative research in our departments. A.S. is grateful to Dr. Anna Panorska and Dr. Tomasz Kozubowski for mentorship and support. This work was partially supported by a grant from the Simons Foundation (# 283770 to N.S.)

6. Appendix

6.1. Maximal likelihood and minimal standard error

Take a family of linear regressions depending on the parameter , with factors and the intercept:


Solve for by maximum likelihood estimation. The Gaussian density for is given by


The log likelihood of (9) is derived from the Gaussian density (10) and is given by

But the standard error of the regression (9) is estimated as

Thus we can express


We used crucially here that the residuals in (9) are normalized so that they have the same variance . To maximize from (11), we need to first minimize the standard error by computing it for each and then choosing an appropriate ; then to choose the which maximizes (11) for given , which turns out to be .

6.2. Background on autoregressive models and random walk

Consider a time series of random variables :

This is called an autoregressive process of order 1 or AR(1): Regression of this sequence onto itself with a one-step time lag. It has the following properties. For , this sequence converges weakly to the stationary distribution: , with

This means that for every interval ,

In addition, mean and variance of converge to that of the limiting distribution: , . Thus this time series exhibits mean reversion: If then is more likely to decrease compared to than to increase. For , this is random walk: Increments are independent for different . This sequence does not have a limit as . The expectation is constant. But the variance .

6.3. Background on Bayesian inference

Assume we have a sample . Denote the variance by . Random variables are independent, and has density

We set a non-informative prior on , which means we do not have any existing information about these parameters: . This is an infinite measure:

Thus we cannot normalize it (divide it by a constant) to make it a probability measure. But we can still apply Bayesian statistics to this measure. We choose this form of measure because we can get explicit posterior. Bayesian inference works as follows. The

likelihood, that is, density of given is


Unlike the prior, the posterior is a finite measure. We can normalize it by computing its integral and dividing it by this integral. After computation, we get:


Recall that Gamma distribution with shape

and scale has density and expectation:


Using (14), we can rewrite (13) as follows:

That is, has marginal Gamma distribution with shape and expectation ; and has inverse Gamma distribution with shape . The conditional distribution of given is normal. The unconditional (marginal) distribution of is Student (-distribution). A Student distribution has heavier tails than a normal distribution, which implies more uncertainty, resulting from our Bayesian estimation framework.

6.4. Empirical data

Empirical means and variances of biomass and basal area logarithms, for each year

Year Number Yearly Mean Yearly Variance Yearly Mean Yearly Variance
of of Biomass of Biomass of Basal Area of Basal Area
Observations Logarithm Logarithm Logarithm Logarithm
1970 522 1.81 0.77 1.67 0.48
1971 1216 2.05 0.8 1.9 0.5
1972 1286 1.97 0.67 1.87 0.43
1973 335 1.66 0.59 1.64 0.39
1974 304 1.68 0.55 1.66 0.36
1975 902 1.97 0.66 1.96 0.54
1976 1883 2.17 0.61 2.02 0.4
1977 422 1.82 0.37 1.81 0.25
1978 1319 2.77 0.7 2.53 0.48
1979 1339 2.68 0.77 2.52 0.54
1980 1047 2.2 0.7 2.13 0.52
1981 396 1.86 0.7 1.8 0.51
1982 8 1.98 0.12 1.97 0.12
1983 98 2.23 0.66 2.06 0.43
1984 358 2.65 0.57 2.32 0.36
1985 629 2.38 0.75 2.15 0.47
1986 665 2.24 0.62 2.07 0.4
1987 732 2.22 0.62 2.05 0.42
1988 604 1.98 0.56 1.92 0.34
1989 1597 2.49 0.95 2.38 0.58
1990 723 1.46 0.47 1.54 0.34
1991 581 2.02 0.61 1.92 0.38
1992 1782 2.67 0.68 2.54 0.46
1993 865 2.88 0.77 2.65 0.55
1994 647 3.21 0.67 2.93 0.48
1995 625 2.85 0.77 2.66 0.52
1996 858 2.57 0.86 2.43 0.68
1997 2247 3.08 0.81 2.82 0.57
1998 977 2.6 0.68 2.49 0.49
1999 905 2.11 0.67 2.13 0.54
2000 101 2.62 1.0 2.46 0.74
2001 756 2.47 0.38 2.45 0.3
2002 309 2.32 0.47 2.32 0.39
2003 3414 3.08 0.85 2.83 0.61
2004 19 2.55 0.29 2.51 0.23
2005 641 2.71 0.73 2.57 0.57
2006 599 3.42 0.81 3.08 0.53
2007 841 2.67 0.81 2.55 0.62

Quantity of observation patch-year pairs with given time gap

Year Gap 3 4 5 6 7 8 9 10 11 12 13 14
Quantity 1 2 65 915 1381 3334 1923 2214 2543 2677 1972 694
Year Gap 15 16 17 18 19 20 21 22 23 24 25 26
Quantity 922 533 391 569 390 123 8 66 8 17 22 22
Year Gap 27 28 29 30 31 32 33 34 35 36 37 38
Quantity 4 17 14 9 6 5 12 1 0 5 2 0


  • [1] Alicia Barrett, Peter Rappoport (2011). Price-Earnings Investing. Reality in Returns JPMorgan Asset Management, 1, 1–12.
  • [2] Daniel B. Botkin (1993). Forest Dynamics: An Ecological Model. Oxford University Press.
  • [3] Harald Bugmann (2001). A Review of Forest Gap Models. Climatic Change 51 (3–4), 259–305.
  • [4] John Y. Campbell, Robert J. Shiller (2001). Valuation Ratios and the Long-Run Stock Market Outlook: An Update. National Bureau of Economic Research.
  • [5] Hal Caswell (2001). Matrix Population Models: Construction, Analysis, and Interpretation. Sinauer Associates.
  • [6] Bruce D. Craven, Sardar M. N. Islam (2008). A Model for Stock Market Returns: Non-Gaussian Fluctuations and Financial Factors. Review of Quantitative Finance and Accounting 30 (4), 355–370.
  • [7] Joseph Davis, Roger Aliaga-Díaz, Charles J. Thomas (2012). Forecasting Stock Returns: What Signals Matter, and What Do They Say Now. Vanguard Research.
  • [8] Patrick Dube, Marie J. Fortin, Charles D. Canham, Danielle J. Marceau (2001). Quantifying Gap Dynamics at the Patch Mosaic Level Using a Spatially-Explicit Model of a Northern Hardwood Forest Ecosystem. Ecological Modelling 142 (1), 39–60.
  • [9] Aaron M. Ellison (2004). Bayesian Inference in Ecology. Ecology Letters 7 (6), 509–520.
  • [10] Eugene F. Fama (1995). Random Walks in Stock Market Prices. Financial Analysts Journal 51 (1), 75–80.
  • [11] Richard A. Fleming, Hugh J. Barclay, Jean-No’́el Candau (2002). Scaling-up an Autoregressive Time-Series Model (of Spruce Budworm Population Dynamics) Changes its Qualitative Behavior. Ecological Modelling 149 (1–2), 127–142.
  • [12] Julian C. Fox, Huiquan Bi, Peter K. Ades (2008). Modelling Spatial Dependence in an Irregular Natural Forest. Silva Fennica 42 (1), 1–35.
  • [13] Amit Goyal, Ivo Welch (2003). Predicting the Equity Premium with Dividend Ratios. Management Science 49 (5), 639–654.
  • [14] Lyle C. Gurrin, Jennifer J. Kurinczuk, Paul R. Burton (2000). Bayesian Statistics in Medical Research: an Intuitive Alternative to Conventional Data Analysis. Journal of Evaluation in Clinical Practice 6 (2), 193–204.
  • [15] Paul J. Hanson, Jake F. Weltzin (2000). Drought Disturbance from Climate Change: Response of United States Forests. Science of the Total Environment 262 (3), 205–220.
  • [16] Christos Ioannidis, David A. Peel, Michael J. Peel (2003). The Time Series Properties of Financial Ratios: Lev Revisited. Journal of Business Finance & Accounting 30 (5–6), 699–714.
  • [17] Jennifer C. Jenkins, David C. Chojnacky, Linda S. Heath, Richard A. Birdsey (2003). National-Scale Biomass Estimators for United States Tree Species. Forest Science 49 (1), 12–35.
  • [18] Takashi Kohyama, Eizi Suzuki, Tukirin Partomihardjo, Toshihiro Yamada (2001). Dynamic Steady State of Patch-Mosaic Tree Size Structure of a Mixed Dipterocarp Forest Regulated by Local Crowding. Ecological Research 16 (1), 85–98.
  • [19] Simon A. Levin (1999). Fragile Dominion: Complexity and the Commons. Perseus Publishing.
  • [20] Simon A. Levin, Robert T. Paine (1974). Disturbance, Patch Formation, and Community Structure. Proceedings of the National Academy of Sciences 71 (7), 2744–2747.
  • [21] Simon A. Levin (2003). Complex Adaptive Systems: Exploring the Known, the Unknown and the Unknowable. Bulletin of the American Mathematical Society 40, 3–19.
  • [22] Jeremy W. Lichstein, Theodore R. Simons, Susan A. Shriner, Kathleen E. Franzreb (2002). Spatial Autocorrelation and Autoregressive Models in Ecology. Ecological Monographs 72 (3), 445–463.
  • [23] Jean Liénard, Ionut Florescu, Nikolay Strigul (2015). An Appraisal of the Classic Forest Succession Paradigm with the Shade Tolerance Index. PLoS ONE 10 (2), 117–138.
  • [24] Jean Liénard, John Harrison, Nikolay Strigul (2016). US Forest Response to Projected Climate-Related Stress: a Tolerance Perspective. Global Change Biology 22 (8), 2875–2886.
  • [25] Jean Liénard, Nikolay Strigul (2015). Linking Forest Shade Tolerance and Soil Moisture in North America. Ecological Indicators 58, 332–334.
  • [26] Jean Liénard, Nikolay Strigul (2016). An Individual-Based Forest Model Links Canopy Dynamics and Shade Tolerances Along a Soil Moisture Gradient. Royal Society Open Science 3 (2), 13 pages.
  • [27] Jean Liénard, Dominique Gravel, Nikolay Strigul (2015). Data-Intensive Modeling of Forest Dynamics. Environmental Modelling & Software 67, 138–148.
  • [28] Jean Liénard, Nikolay Strigul (2016). Modelling of Hardwood Forest in Quebec Under Dynamic Disturbance Regimes: a Time-Inhomogeneous Markov Chain Approach. Journal of Ecology 104 (3), 806–816.
  • [29] Dmitrii O. Logofet, Ekaterina V. Lesnaya

    (2000). The Mathematics of Markov Models: What Markov Chains Can Really Predict in Forest Successions.

    Ecological Modelling 126 (2), 285–298.
  • [30] John McCarthy (2001). Gap Dynamics of Forest Trees: a Review with Particular Attention to Boreal Forests. Environmental Reviews 9 (1), 1–59.
  • [31] Paul R. Moorcroft, George C. Hurtt, Stephen W. Pacala (2001). A Method for Scaling Vegetation Dynamics: the Ecosystem Demography Model. Ecological Monographs 71 (4), 557–586.
  • [32] Stephen W. Pacala, Charles D. Canham, John Saponara, John A. Silander Jr., Richard K. Kobe, Eric Ribbens (1996). Forest Models Defined by Field Measurements: Estimation, Error Analysis and Dynamics. Ecological Monographs 66 (1), 1–43.
  • [33] Stephen W. Pacala, Charles D. Canham, John A. Silander Jr. (1993). Forest Models Defined by Field Measurements: I. The Design of a Northeastern Forest Simulator. Canadian Journal of Forest Research 23 (10), 1980–1988.
  • [34] John Pastor, Angela Sharp, Peter Wolter (2005). An Application of Markov Models to The Dynamics of Minnesota’s Forests. Canadian Journal of Forest Research 35 (12), 3011–3019.
  • [35] Svetlozar T. Rachev, John S.J. Hsu, Biliana S. Bagasheva, Frank J. Fabozzi (2008). Bayesian Methods in Finance, Wiley.
  • [36] Andrew E. Scholl, Alan H. Taylor (2010). Fire Regimes, Forest Change, and Self-Organization in an Old-Growth Mixed-Conifer Forest in Yosemite National Park. Ecological Applications 20 (2), 362–380.
  • [37] Herman H. Shugart (2003). A Theory of Forest Dynamics: The Ecological Implications of Forest Succession Models. Blackburn Press.
  • [38] George R. Stephens, Paul E. Waggoner (1980). A Half Century of Natural Transitions in Mixed Hardwood Forests. Connecticut Agricultural Experiment Station Bulletin 783.
  • [39] Nikolay Strigul, Ionut Florescu, Alicia R. Welden, Fabian Michalczewski (2012). Modelling of Forest Stand Dynamics Using Markov Chains. Environmental Modelling and Software 31, 64–75.
  • [40] Nikolay Strigul, Denis Pristinski, Drew Purves, Jonathan Dushoff, Stephen Pacala (2008). Scaling from Trees to Forests: Tractable Macroscopic Equations for Forest Dynamics. Ecological Monographs 78 (4), 523–545.
  • [41] Nikolay Strigul (2012). Individual-Based Models and Scaling Methods for Ecological Forestry: Implications of Tree Phenotypic Plasticity. Sustainable Forest Management 359–384, InTech.
  • [42] Michael B. Usher (1979). Markovian Approaches to Ecological Succession. The Journal of Animal Ecology 48 (2), 413–426.
  • [43] C. E. Van Wagner (1978). Age-Class Distribution and the Forest Fire Cycle. Canadian Journal of Forest Research 8 (2), 220–227.
  • [44] Paul E. Waggoner, George R. Stephens (1970). Transition Probabilities for a Forest. Nature 225, 1160–1161.
  • [45] Alex S. Watt (1947). Pattern and Process in the Plant Community. The Journal of Ecology 35 (1–2), 1–22.
  • [46] Jianguo Wu, Orie L. Loucks (1995). From Balance of Nature to Hierarchical Patch Dynamics: A Paradigm Shift in Ecology. Quarterly Review of Biology 70 (4), 439–466.