Autoregressive Modeling of Forest Dynamics

Abstract

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 performed a time series statistical analysis of forest biomass and basal areas recorded in Quebec provincial forest inventories from 1970 to 2007. The geometric random walk model adequately describes the yearly average dynamics. For individual patches, we fit an autoregressive process (AR) of order 1 capable to model negative feedback (mean-reversion). Overall, the best fit also turned out to be geometric random walk; however, the normality tests for residuals failed. In contrast, yearly means were adequately described by normal fluctuations, with annual growth on average of 2.3%, but with a standard deviation of order of 40%. We used a Bayesian analysis to account for the uneven number of observations per year. This work demonstrates that autoregressive models represent a valuable tool for the modeling of forest dynamics. In particular, they quantify the stochastic effects of environmental disturbances and develop predictive empirical models on short and intermediate temporal scales.

1. Introduction

1.1. Background

Understanding the dynamics and self-organization of ecosystems is one of the most challenging problems in modern ecology [1]. Self-organization occurs simultaneously on several levels of hierarchical ecosystem organization and involves dynamic processes operating on different temporal and spatial scales [2]. 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 [3]. A large number of forest models have been developed over the last decades [4,5,6,7,8,9]. 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 the effects of all these disturbances on forest biomass dynamics [10]. However, their application is based on the discretization of the state space [11]. Spatially-explicit individual-based models are able to simulate the effects of disturbances in continuous time and state space [3,9,12]. 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 [9,13].

Another challenge to our understanding of forest dynamics is the multidimensional nature of this complex adaptive system [2,14]. Forest disturbances are traditionally associated with a loss of biomass. However, Markov chain models based only on biomass do not capture forest succession comprehensively [11,14]. 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 a Markov chain approach to develop a multidimensional Markov chain [14]. 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 [7,10,11,15].

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 the major shortcomings of the 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 the S& P 500 (Standard & Poor’s 500) index that measures stock performance of 500 large companies on stock exchanges in the United States) 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 the 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 [16] was actively developed in the second half of the twentieth century after suggestion in [17] 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 [18]. Conservation law modeling in the case of continuous time and patch state results in the reaction–advection–diffusion model [18].

The patch dynamics concept can be applied for understanding and predicting of forest dynamics and different levels of forest organisation within the hierarchical patch dynamics paradigm [3,16]. At the level of individual trees, the patch dynamics concept is traditionally called forest gap dynamics [5,19,20]. Individual-based forest models capture gap dynamics by simulating growth, competition, and the mortality of individual trees [6,12,21,22]. Individual-based models and analytically-tractable models approximating gap dynamics [9,23] 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 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 [11]. Forest stand is affected by individual-level tree dynamics, as well as by large-scale disturbances such as forest fires, drought, and hurricanes [20,24], which affect many trees in the stand at the same time. The stand-level dynamics scale up to the next hierarchical level of particular forest type or regional patch mosaic (level 3 in the hierarchical patch mosaic Matreshka [3]).

The interplay of individual-level and stand-level changes and disturbances creates complex 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 [11]. As the consequence of this system complexity, a classical linear patch dynamics model does not capture patch dynamics of the US and Canadian forests [10,11,25]. This classical patch-dynamics model can be represented in the continuous case by the advection-reaction patch-dynamics conservation law model [18], and in the discrete case by a birth and death process that can be written as a Markov chain [3], or as a simple forest fire [26]. As a result, in order to capture the 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 variables (biomass), we can achieve an adequate representation of forest patch dynamics within the framework of Markov chains [10,11,14]. Markov chains provide an analytically tractable 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 the 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 three decades using the same field survey protocol [14]. We operate with the same biomass and basal area data derived from Quebec forest inventories in our previous publication on Markov chain modeling (the data-mining protocol is available in [14,27]). 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 the quantifying of the 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 normal fluctuations, and stock indices (in 3–4 years or more) having normal fluctuations. To account for the 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 forest inventory plot has a circular form of 400 m² [28]. The database consists of 32,552 plot re-measurements at 11,660 different locations. The Quebec forest inventory is designed to comprehensive describe the patch mosaic of Quebec forests and plots cover the Quebec territory practically uniformly [28]. The GIS-based map of forest inventory plots is published as Figure 1 in [10]. Forest plots cover hardwood and mixed forests in the northern temperate zone (9621 and 7663, respectively) and continuous boreal forests in the boreal zone (11,969 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 [11] (see Figure 2 in Appendix to [11]). In this inventory, each patch observation includes the 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 [29,30]. The computations of biomass and basal area (as well as other characteristics, such as shade tolerance index, and biodiversity measured by Shannon entropy) are done in [11,14,25,27]. The biomass is this article refers to the plot biomass, which is the sum of biomasses of all trees computed using formulas from [29] (see Section 3.1 in [14] 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: An autoregressive model, where each next year’s logarithm of biomass y(t + 1) is a weighted sum of the previous year’s logarithm of biomass y(t) 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,

y(t + 1) = r + ay(t) + ε(t)

(1)

where r, a are constants, and ε(t) are i.i.d. (independent identically distributed) N (0, σ²) random variables. If 0 < a < 1, this sequence y(0), y(1),…exhibits mean reversion to its long-term average m = r/(1 − a). That is, if y(t) > m, then y(t + 1) − m is likely to be smaller than y(t) − m, and vice versa. Examples of earlier use of such models for forest modeling include [31,32,33]. They also include 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 forests is left for future research.

Since data are collected on irregular time intervals, we apply (1) multiple times to itself to get the expression of y(t₁) from y(t₀) if t₀ and t₁ are consecutive years for which this patch was observed. Then we try various a and obtain for each a the maximum likelihood estimate via linear regression. We compare these likelihoods and find the best fit for a. It turns out to be a = 1. That is, this sequence does not actually exhibit any mean reversion, but behaves like a random walk, when each next increment is independent from the past:

y(t + 1) = y(t) + r + ε(t), ε∼N(0, σ²) i.i.d.

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

However, the residuals for a = 1 do not pass the normality test. This model does not actually fit well, and we cannot find confidence intervals for a 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 basal area. For each version, we do it in two ways: (a) Original logarithms of biomass/basal area; (b) with a logarithm of mean biomass/basal area for this year subtracted. In both cases, the maximum likelihood estimate gives us random walk a = 1.

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 is a classic model for stock market movements, going back to classic research by [34]. Mean reversion is commonly observed in financial ratios such as earnings yield or dividend yield [35,36]. See also [37,38,39] for this research and the 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, see Figure 1. Let c(t) be the logarithm of this mean. We find that the random walk adequately describes this:

c(t) = c(t − 1) + µ + ε(t), ε(t) ∼ N (0, ρ²) i.i.d.

In particular, we find that the increments ε(t) indeed have a normal distribution and not heavy tails. However, we cannot simply take yearly averages for every year t, since they would have different precision. Reason: Each year, t, 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 [40], as well as in medical statistics [41] and quantitative finance [42].

3. Results and Discussion

We have 32,552 observations of 11,660 patches of Quebec forests. Each patch is observed at most once a year, during 1970–2007. On average, each patch was observed around 3 times: 32552/11660 ≈ 3. Out of these 11,660 patches, 10,215 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 were observed in different years: patch 7000406701 was observed only in 1970; patch 7000406901 was observed in 1970 and 1978; patch 8509702201 was observed in 1985, 1993, and 2003; and patch 7000406902 was observed in 1970, 1978, 1985, and 1997. Let P be the set of observed patches. Each p $ϵ$ P has observations x(t, p) at years t $ϵ$ T(p), where T(p) ⊆ {1970,…, 2007}. Let yp (t): = ln x(t, p). Yearly means are defined as

$$m \left(t\right)=\frac{1}{\#\{p \in P | t \in T\left(p\right)\}}{\displaystyle \sum }_{p: t \in T\left(p\right)}x\left(t,p\right), t=1970, \dots , 2007$$

Additionally, we define c(t) = ln m(t). 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 of 8 years, 1923 such pairs with a gap of 9 years, but only 66 pairs with a gap of 22 years. More detailed data are in Appendix D.

3.1. Autoregressive Model for Individual Patches

Consider the following time series:

$$yp\left(t\right)=r+ayp(t-1)+\mathsf{\epsilon }p\left(t\right),t=1970,\dots ,2007,p\in P$$

(2)

where a is an AR(1) parameter, and εp(t) are i.i.d. N(0, σ²). As mentioned earlier, our main difficulty is that we do not have observations for all t and p, only for t $\in$ T(p). Assume t, t + u are subsequent time points in T(p). Then

$$y\left(\mathrm{t}+\mathrm{u}\right)=a^u y\left(t\right)+\left(1+a+\dots +a^{u-1}\right)\mathrm{r}+\sigma {{\displaystyle \sum }}_{t=1}^s{a}^{v-1}\mathsf{\epsilon }p(t+\mathrm{v})$$

(3)

We do this both for yp(t) (raw data), $\tilde{y}p(t)$ = yp(t) − c(t) (centered data), and centered data without t = 1982 and t = 2004. 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 Equation (3) and apply a 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 a and run regression with respect to r. Then choose an a such that the standard error of this regression is smallest. To properly apply this linear regression, divide Equation (3) by a constant to make the standard error in error terms in Equation (3) the same for all u:

C(a, u)(y(t + u) − a^uy(t)) = D(a, u)r + δp (t + u), δp(t + u) ∼ N(0, σ²) i.i.d.
C(a, u): = [1 + a² +…+ a^2(u−1)]^−1/2
D(a, u): = C(a, u)[1 + a +…+ a^u−1]

(4)

For every a $\in$ (−2, 2) which is a multiple of 0.01, fit a simple linear regression (without intercept) to find r 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 a = 1 minimizes the standard error (thus maximizing the likelihood). Thus, the dynamics in Equation (2) are given by yp(t) = r + yp(t − 1) + ε(t), which is a simple random walk with Gaussian increments.

Assuming that the model is, in fact, a random walk, let us make quantile–quantile (QQ) plots for residuals in Equation (4) (we have simple linear regression without intercept):

$$u^{-\frac{1}{2}}\left(y\left(t+u\right)-y\left(t\right)\right)=r\sqrt{u}+\mathsf{\sigma }{\epsilon }_p\left(t\right)$$

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

For non-centered biomass, the minimal standard error is achieved when a = 1, see Figure 2a. This corresponds to a random walk model, and we get σ = 0.223 and r = 0.0378. However, the residuals are not normal: see them on the QQ plot in Figure 2b. Then we repeated the computation above for centered values: y˜p(t) = yp(t) − z(t) instead of yp(t). Similarly, the standard error is minimized for a = 1. For this value, r = 0.0569, σ = 0.224. 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 a = 1, with σ = 0.224 and r = 0.00570. The QQ plot of residuals is still not normal. We omited 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 Equation (2), and computing standard error of the regression Equation (4), we again get a = 1, r = 0.0337, and σ = 0.207. Again, the best-fitting model among AR(1), according to 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 a = 1, with r = 0.00340 and σ = 0.197. Centering the basal area data and removing years 1980 and 2004, we get: a = 1, r = 0.00341, σ = 0.197. 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 p and corresponding years t in T(p). Denote by y’_p(t) the logarithm of biomass, and by y’’_p(t) the logarithm of basal area for patch p and year t. If T(p) = {t₀, t₁, t₂, …} has more than one year, order them in increasing order: t₀ < t₁ < t₂ <… Compute correlation coefficient between

y′_p(t_k) − y′_p(t_k−1) and y″_p(t_k) − y″_p(t_k−1)

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:

Δc(t): = c(t + 1) − c(t) ∼ N (µ, ρ²) i.d.d.

(5)

However, the 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 Equation (5) for t = 1970,…, 1981, 1984,…, 2003, 2006, 2007. Then we test the normality of increments using the QQ plot and the Shapiro–Wilk test. We confirm that, indeed, increments are i.i.d. normal. Rewrite Equation (5) as the random walk:

c(t) = c(0) + ξ₁ +…+ ξ_t, ξ_i ∼ N (µ, ρ²) i.i.d.

(6)

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

m(t) = m(0) exp(ξ₁ +…+ ξ_t), t = 0, 1, 2,…

(7)

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

E[m(t)] = m(0)· exp(t(µ + σ²/2)); Var[m(t)] = m²(0) exp(2µt + σ²t) (exp(σ²t) − 1)

For the biomass, we get: µ = 0.021 and ρ = 0.512. These increments pass the Shapiro–Wilk test with p = 0.80. With the two years 1982 and 2004 removed, the estimates will not change much: µ = −0.011 and ρ = 0.507. This still passes the Shapiro–Wilk normality test with p = 0.60. With the analysis for basal area instead of biomass, we get: µ = 0.0240, ρ = 0.367 for all years, and µ = 0.00455, ρ = 0.370 for years without 1982 and 2004. The Shapiro–Wilk test gives us p = 0.89 for all years, and p = 0.58 for years without 1982 and 2004. The QQ plots in Figure 3 show that, indeed, the residuals are close to normal.

In Figure 4, we plot histograms of the biomass and basal area for an individual patch in 2007. In Figure 5, we simulate biomass and basal area as in Equation (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 a correlation of 0.977. For years 1982 and 2004 removed, we get a correlation of 0.980. Previously, we had very high correlations 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. For each fixed year t, the logarithms of observations are x₁(t), …, x_n(t). For these values, we performed the analysis above: the posterior mean µ(t) and the posterior variance v(t) were generated. The simulation of µ(t) and v(t) was performed N = 1000 times. In Figure 6, we have histograms of 1000 simulations for µ(t) and σ(t) for t = 0 (year 1970), for both biomass and basal area. Hence we obtained 38 sequences of N numbers: µ₁(t), …, µ_N (t), t = 1970, …, 2007. The average growth rate based on simulated results is

$$\widehat{g}=\frac{1}{NT}{\displaystyle \sum }_{i=1}^N[u_i\left(T\right)-u_i\left(0\right)]$$

(8)

The mean increments are: Dµ_i(t) = µ_i(t) −µ_i(t − 1), t = 1, …, T, i = 1, …, N. Assuming these are N (g, σ²) i.i.d. we estimate g and σ² as $\widehat{g}$ in Equation (8) and ${\widehat{\sigma }}^2$:

$${\widehat{\mathsf{\sigma }}}^2=\frac{1}{NT}{\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{t=1}^T{(∆u_i\left(t\right)-\widehat{g})}^2$$

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

$${\widehat{g}}_{bio}\approx 0.023, {\widehat{\sigma }}_{bio}{}^2\approx 0.214, {\widehat{g}}_{BA}\approx 0.023, {\widehat{\sigma }}_{BA}{}^2\approx 0.134$$

Thus the growth rate of forest, measured by the biomass or basal area (on the logarithmic scale), is 2.3% per year. These numbers are close to the ones from frequentist analysis from the previous subsection: 2.1% with years 1982 and 2004. However, with the years 1982 and 2004 are removed, this estimate changes to −1.1%. For basal area, we have a similar comparison with the previous subsection: 2.4% with years 1982 and 2004, 0.46% without those years. As discussed earlier, we view Bayesian analysis as the more statistically sound. Thus 2.1% 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 $\sqrt{0.214}=0.46$, that is, 46% per year for the biomass, and $\sqrt{0.314}=0.37$, that is, 37% per year for the basal area. Similar 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 an AR(1) autoregressive model to patch/stand dynamics of Quebec forests. Overall, we followed major steps of application of the 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 as it is hard to expect that the forest dynamics modeling will simply mirror stock market theory.

One interesting result reported in Section 2.2 is that a = 1 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 the possible use of random-walk models with tails heavier than normal. This dichotomy of average patch vs. individual patch reminds us of a stock market dynamic, where each individual stock is highly volatile, with non-Gaussian fluctuations [43], but a stock market overall (measured by Standard & Poor 500 or another stock market index) in the long run follows a 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, a = 1 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 −1 < a < 1. This means that biomass and basal area at every step will drift away in both positive and negative directions. We can search for mechanistic meaning of this result in both fundamental patterns of forested ecosystem dynamics and general modeling assumptions related to the autoregressive modeling approach and the patch dynamics modeling framework.

A 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 [11]. 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 the effects of these assumptions in the Markov chain modeling framework [10,11,14]. Similarly, individual-based models are often applied for forest simulations in different regions with the same growth/mortality characteristics of a particular tree species. Autoregressive modeling is a new approach in forest modeling, and it is not yet known how robust our predictions are with respect to variability of forest patch mosaic. We can hypothesise that outcomes reported in Section 2.2 (a = 0 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 [11]. 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 [14]. Time inhomogeneous Markov chains allow us to relax this assumption and take into account both climate driven changes in tree growth and frequency of environmental disturbances [10]. Similar 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 a small decrease of the total disturbance rate in boreal and hardwood forests (see Table 1 in [10]). However, the relatively small number of 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 a 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 a continuous-time and state counterpart of discrete Markov chain models. Markov chains have been applied at the individual tree and species levels [44,45], stand-level [46], and landscape-level [47]. The autoregressive and ARIMA models can also be employed to capture similar dynamics. In addition, time series models might be used in the approximation of 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 the 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 on a large scale. We wish to compare this approach with state-of-art mathematical models developed on other principles, in particular with Markov chains [11,14,25,27] and individual-based models [3,13], 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 assessments of forest vulnerability to changes in environmental disturbance regimes.

We can generalize a 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. built upon our past analyses of tolerance patterns in North American forests, including observed shade tolerance patterns [25], associations between shade tolerance and soil moisture [27], and relative trade-offs between shade, drought, and water-logging tolerance at the continental scale [48].

We can extend our research to the USA Forest Inventories in an 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 [48]. We can link climatic models and forest tolerances [25,27,48]. Spatially-explicit information from the USDA Forest Service Forest Inventory and Analysis Program (FIA), the Canadian Forest Inventory, and climate data sets and models (WorldClim and PRISM), can be integrated to quantify the effect that climate variables in North America have on forest tolerances and disturbances via the recently proposed methodology [48].

5. Conclusions

Individual patches have biomass and basal areas whose dynamics (on logarithmic scale) are not well described by the classic autoregressive model y(t + 1) = ay(t) + r + ε(t) with normal noise terms ε(t), 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 a = 1, and increments y(t + 1) − y(t) independent of the past y(s), s ≤ t. 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 a random walk with normal increments. A Bayesian analysis accounts for the fact that different years have different numbers of observations. We calculated a growth rate (measured on the log scale) of 2.3% per year, with a standard deviation of 46% for biomass, and 37% 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.