1. Introduction
1.1. Background
Understanding the dynamics and selforganization of ecosystems is one of the most challenging problems in modern ecology [19]. Selforganization 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, individualbased 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 nonstationary 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]. Spatiallyexplicit individualbased 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 difficulttopredict factors. Each individual stock has fluctuations with heavy tails. But the total stock market, measured by an index (such as Standard & Poor 500) has longrun fluctuations (3–4 years) which follow Gaussian distribution. These fluctuations depend on various factors, such a priceearnings 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 (longterm and shortterm). 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 patchmosaic 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 patchmosaic 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 reactionadvectiondiffusion 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]. Individualbased forest models capture gap dynamics by simulating growth, competition and mortality of individual trees [33, 32, 3, 8]. Individualbased models and analyticallytractable models approximating gap dynamics [18, 40] provide scaling from individuallevel dynamics to the next level of forest hierarchical organisation, the standlevel.
In the present work we apply autoregressive models to the standlevel 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 individuallevel tree dynamics, as well as by largescale disturbances such as forest fires, drought and hurricanes [15, 30], which affect many trees in the stand at the same time. The standlevel 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 individuallevel and standlevel changes and disturbances creates complex dynamical patterns at the standlevel. 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 patchdynamics model can be represented in continuous case by advectionreaction patchdynamics 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 statespace 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 (datamining 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 welldescribed by normal distribution, while individual patches are not: Fluctuations have heavy tails. This is similar to financial markets, with individual stocks having nonnormal 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 www.mffp.gouv.qc.ca. Each permanent forest inventory plot has a circular form of 400 . The database consists of 32552 plot remeasurements 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 GISbased 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 remeasurement 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,
(1) 
where are constants, and are i.i.d. (independent identically distributed) random variables. If , this sequence exhibits mean reversion to its longterm 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 patchyear 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.
3.1. Autoregressive model for individual patches
Consider the following time series:
(2) 
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
(3) 
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 loglikelihood 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 :(4) 
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.Assuming that the model is, in fact, a random walk, let us make quantilequantile (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 noncentered 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 (noncentered) data.
Modeling basal area data as in (2), and computing standard error of the regression (4), we again get , , and . Again, the bestfitting 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.
3.2. Yearly averages, frequentist analysis
We model the logarithm of yearly average using random walk:
(5) 
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 ShapiroWilk test. We confirm that, indeed, increments are i.i.d. normal. Rewrite the model (5) as the random walk:
(6) 
Or, equivalently, we can rewrite (6) in the original scale, instead of the logarithmic scale:
(7) 
From (7), we can compute the mean and variance of :
For the biomass, we get: and . These increments passes ShapiroWilk normality test with . With removed two years 1982 and 2004, the estimates will not change much: and . This still passes ShapiroWilk normality test with . Repeating this analysis for basal area instead of biomass, we get: for all years, and for years without 1982 and 2004. ShapiroWilk 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.
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
(8) 
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.
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 randomwalk 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 nonGaussian 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 34 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, individualbased 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 intermediateterm 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 nonstationary (time inhomogenous) dynamics. In particular, we previously detected some nonstationary 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], standlevel [42], and landscapelevel [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 individualbased 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.
Timeseries analysis is a novel modeling approach for forest modeling at the large scale. We wish to compare this approach with stateofart mathematical models developed on other principles, in particular with Markov chains [39, 27, 25, 23] and individualbased 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 tolerancedisturbance 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 tradeoffs between shade, drought, and waterlogging 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 changerelated disturbances such as drought and fires [24]. We can link climatic models and forest tolerances [23, 25, 24]. Spatiallyexplicit 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 heavytailed 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, 1year 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 noninformative and the model is Gaussian, we do not need to do any Monte Carlo or MetropolisHastings computations: There are explicit formulas for the posterior distribution of parameters.
Acknowledgments
We thank Dr. Adam Erickson for help with proofreading. Initial datamining 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:
(9) 
Solve for by maximum likelihood estimation. The Gaussian density for is given by
(10) 
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
(11) 
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 onestep 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 noninformative 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(12) 
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:
(13) 
Recall that Gamma distribution with shape
and scale has density and expectation:(14) 
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 patchyear 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 
References
 [1] Alicia Barrett, Peter Rappoport (2011). PriceEarnings 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 LongRun 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: NonGaussian Fluctuations and Financial Factors. Review of Quantitative Finance and Accounting 30 (4), 355–370.
 [7] Joseph Davis, Roger AliagaDí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 SpatiallyExplicit 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, JeanNo’́el Candau (2002). Scalingup an Autoregressive TimeSeries 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). NationalScale 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 PatchMosaic 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 ClimateRelated 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 IndividualBased 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). DataIntensive 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 TimeInhomogeneous 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 SelfOrganization in an OldGrowth MixedConifer 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). IndividualBased 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). AgeClass 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.
Comments
There are no comments yet.