Abstract: This study presents a new characterization on optimal harvesting problem. By relying on both Gompertz shape stochastic growth model and the Shanon, Tsallis, Kullback, Fisher information measures, the solution of biologically optimal rotation problem is presented and exemplified. As an present experience a real data set is used from the repeated measurements on permanent sample plots of pine stands in Lithuania. All results are implemented in symbolic algebra system MAPLE.
INTRODUCTION
Forest ecosystems contain 91% of the earth`s estimated biomass (Begon et al., 1996) and experience discontinuities dynamics in their growth. Harvesting strategies are affected by many factors such as temporal development of forest stands, site quality, desired product, stocking and intensity of forest management activity. What is the best timber-harvesting solution? There is no single answer to this question, because forest owners and managers must hold two opposite paradigms; the forest protection and felling. Formulating basic principles of timber-harvesting theories necessarily involves keeping a balance between conserving forest biodiversity, which is socially useful and timber production, which is economically useful. It is evident that a good answer must reflect a growth process of forest stands. Appropriate application of principles derived from knowledge of the dynamics of stands is an important element in determination of an optimal harvesting procedure. A tree diameter distribution is one of the most descriptive and important stand characteristics. Therefore, knowledge the evolution of age-dependent trajectory of tree diameter size at breast height is fundamental to developing a biological optimal rotation period.
Over the last few decades the problem of rational harvesting planning has significantly increased interest in the forest scientific literature. For centuries sustainable forestry has been interpreted in the forestry literature as a situation while the volume of the timber harvest does not exceed the natural growth of the forest stand. In this context, foresters usually examine two paradigms: biologically optimal rotation age and economically optimal rotation age. Optimal rotation age in biological sense can be defined by the concept of maximum sustainable yield. Maximum sustainable yield is the largest yield that can be harvested and which leaves the forest stand in good shape for future uses. Nowadays forest management strives to compound both optimal rotation periods. Koskela et al. (2007) introduced a socially optimal period which is based on Hartman`s (1976) approach and found that biodiversity conservation benefits prolong the economically optimal rotation period.
Historically, in forestry practice the optimal rotation period was usually based on economic principles. When the forest growth process is deterministic, the optimal rotation period was primarily defined by the German forester (Faustmann, 1849). A weakness of Faustmann`s model is its deterministic nature.
The complexity of forest ecosystems and the related issue of variability in growth of forest stand present difficulties to mathematical modeling of stand dynamics. The randomness is usually caused by a limitation of our knowledge of analyzed growth process. The consideration of randomness phenomena is central in modern ecology and management (Spagnolo et al., 2004). It is widely recognized that the forest stands operate in a highly uncertain environment (Garcia, 1983, 1994, 2005; Clarke and Read, 1989, 1990; Willassen, 1998; Alvarez and Koskela, 2003, 2007; Rupšys, 2003, 2007; Zhou and Boungiorno, 2004, 2006; Lohmander and Mohammadi, 2007; Rupšys et al., 2007 and references therein). When the forest growth is continuous stochastic process driven by the geometric Brownian motion the periodicity of the harvests for the case of even-aged forest were determined by Nordstrom (1975). Further investigations of the continuous stochastic version of Faustmann`s optimal rotation problem have been made by Clarke and Read (1989, 1990), Willassen (1998), Alvarez and Koskela (2003, 2007), Penttinen (2006), Chladna (2007), Gong and Lofgren (2008) and references therein.
In discrete space the studies of stochastic Faustmann`s optimal rotation problem have been made by Lembersky and Johnson (1975), Boungiorno (2001), Francois et al. (2005), Mohammadi and Lohmander (2007) and references therein.
In present analysis, the forest manager is a non-industrial forest owner in the sense that main earnings do not come from forest harvestings. It can be generalized to managers of publicly owned forests who value benefits such as various forms of recreation, urban sequestration, biodiversity preservation, reducing global CO2 emission, or others. The aim of present study is to establish the biologically optimal timber harvesting period which focused on the stochastic growth and yield dynamics of forest stands. Present objective is to derive a qualitative description of the optimal rotation problem of even-aged stands from only time series data and a specified growth law of tree diameter at breast height. For examining diameter development is used age-dependent transition probability density function of tree diameter at breast height. It is received that the Gompertz growth law (Gompertz, 1825) provides an adequate description of stochastic tree diameter growth (Rupšys, 2003, 2007; Rupšys et al., 2007). In this study a drift term of stochastic tree diameter growth process is modeled by a logistic Gomperz shape growth function and a diffusion term is modeled by a multiplicative noise. Other logistic growth models such as Verhulst, von Bertalanffy, Richards and much more, are equally plausible, but they lead to a numerical solution of transition probability density (Rupšys, 2007). The mathematical reasons for the usage of stochastic model we motivate by the distinction between the stochastic logistic growth model and its deterministic counterpart. For the stochastic logistic shape models of diameter growth we can derive the first two moments, namely, the mean mt and variance st of diameter size (Rupšys, 2007). These equations showed that the mean tree diameter is extremely sensitive with respect to the form and size of the coefficient of volatility (amplitude of noise) and mostly deterministic logistic shape models overestimate the true tree diameter size in the presence of stochastic perturbations (Rupšys, 2007).
Over the past 20 years, the study of the concept of information has become increasingly significant for understanding the fundamental laws of science (Goswami et al., 2005; Frieden, 2006). The information theoretical measures play a crucial theoretical role in physics of macroscopic equilibrium systems. The Shannon`s (1948), Tsallis` (1988), Kullback`s (1997) entropy and Fisher`s (1922) information represent promising tools to show the behavior of multidimensional systems in biology, ecology, biomedicine (Moniz et al., 2007; Al-Khafari et al., 2007; Aoki, 2008). Many researches have emphasized that information and entropy are interrelated and that biological systems necessarily maintain a stable entropy state (Gatenby and Frieden, 2004). A stable steady state in biology is important for an understanding of biological systems as dynamical complex processes. A departure from the steady state indicates a negative unhealthy situation of biological systems. The importance of stable steady state as a criteria for biological well-being, emphasize many researches (Salthe, 2003; Karunanithi et al., 2008). So, it is reasonable to define a biologically optimal rotation period by the age T when the age-dependent transition probability density function p(t,x) of tree diameter size at breast height of forest stand depolarizes with the steady-state transition probability density function p(x). We will compare both non-steady state and steady state densities using the information entropy, defined by Shannon (1948) and Tsallis (1988), the generalized entropy, defined by Kullback (1997) and Fisher (1922) information. In present study four information measures are used as a measure of dynamics order of tree diameter growth process.
The main aim of this contribution is to present a tree diameter density function method of biologically optimal rotation problem. This approach is rather different from other economically optimal rotation procedures defined by Clarke and Read (1989, 1990), Willassen (1998), Boungiorno (2001), Alvarez and Koskela (2003, 2007), Penttinen (2006), Chladna (2007), Koskela et al. (2007), Gong and Lofgren (2008) and Lohmander and Mohammadi (2007), since these procedures are based on economic principles of Faustmann discounting dynamics.
MATERIALS AND METHODS
Biology literature uses allometric laws to describe the relationship between various parameters of growth processes in living organisms. In this study we will focus on dynamics of tree diameter distribution of forest stands. We suppose that dynamics of tree diameter growth is expressed in terms of the Gompertz stochastic ordinary differential equation with multiplicative noise (Rupšys, 2007; Rupšys et al., 2007)
(1) |
Where:
| t | = | Age of a forest stand |
| r | = | Diameter intrinsic growth rate |
| K | = | Diameter carrying capacity and forms a numerical upper bound on the diameter size |
| X(t) | = | Breast height diameter at the age t |
| σ | = | Intensity of noise |
| W(t) | = | Standard Brownian motion (white noise), which is a random process whose increments are independent and normally distributed with zero mean and with variance equal to the length of the time interval over which the increment take place |
This diameter growth model, expressed as the ordinary stochastic differential Eq. 1, holds the transition probability density function p(t, x), which project the distribution of tree diameter size subject to the age t. This transition probability density function of tree diameter distribution can be obtained from the Fokker-Plank or forward Kolmogorov equation (Gikhman and Skorokhod, 1965), which relates the variation of diameter with the age t. Each solution of the ordinary stochastic differential Eq. 1 describes one path of diameter development. The ensemble of realizations, described by the transitional probability density function p(t, x), satisfies the following Fokker-Plank equation.
(2) |
In the sequel we assume, that initial diameter size is not random and known exactly x(t0) = x0. The transformation Y(t) = ln(X(t)) and the Ito`s formula, the diffusion term σX(t) of the stochastic growth model (1) converts to constant volatility σ and the nonlinear process (1) transforms into the Ornstein-Uhlenbeck process. So, we can write an exact age dependent solution of the Fokker-Planck Eq. 2. The exact non-steady state solution of Eq. 2 has the following form:
(3) |
The steady-state solution p(x) of Eq. 2 has the form:
(4) |
As was mentioned above, we will define the biologically optimal rotation period from the point of view of infodynamics. The important information measures of a noise-driven dynamical system are Shannon`s (1948), Tsallis` (1988), Kullback`s (1997) generalized entropy and Fisher`s (1922) information and much more.
The spreading of the transition probability density function p(t, x) and the steady state probability density function p(x) is best measured by the Shannon information entropy, also named differential entropy:
(5) |
The Shannon entropy is a well-known method for estimating the degree of disorder in dynamical system. The Shannon entropy conveys uncertainty and information measure. The less is the uncertainty of the system the larger is the information that we acquire.
The Fisher information measure of p(t, x) and p(x) is given by:
(6) |
which measures sharpness or concentration of dynamical system. The Fisher information measure is always positive and reflects the localization characteristics of probability density function more sensitively than the Shannon measure.
The Tsallis entropy, defined by:
(7) |
is an extension of Shannon entropy with one-real-parameter of q. In the limit of q→1, the Tsallis entropy (Eq. 7) reduces to the Shannon entropy (Eq. 5), since (p(t, x))q-1 = e(q-1)lnp(t, x) ≈ 1+(q-1)(t, x).
For the study of difference between the non-steady state and steady state structure of probability density functions of tree diameter size, the Shanon`s information entropy (Eq. 5) we replace by the generalized entropy or the Kullback (1997) information entropy in the following form:
(8) |
The Kullback (1997) information entropy allows us to focus on the non-steady state probability density function p(t, x) of tree diameter size, defined by Eq. 3 and then test the scenario that non-steady state probability density function p(t, x) described the process of diameter growth as well as the steady state probability density function p(x), defined by Eq. 4. Larger Kullback entropy Eq. 8 reflects poorer performance p(x) as an approximation to p(t, x).
The distances S-S(t), I-I(t), S2-S2(t), -H(t) describe the loss of information when the non-steady state probability density p(t, x) is used to approximate the steady state probability density p(x) and measure the difference between the tree diameter maximum entropy and its entropy at given age t. While smaller these distances clearly indicate less information is lost by the non-steady state probability density function, there is no absolute scale against which to judge the significance of distances S-S(t), I-I(t), S2-S2(t), -H(t). So, it is worth noting that significance of these distances might be redefined through one-sided confidence intervals, at first providing lower (upper) bound for the steady state Shannon, Tsallis, Kullback, (Fisher) information measures (Eq. 5-8) and then by comparing values of the steady state lower (upper) bound and the non-steady state curve.
Next we discuss a comparison between the non-steady state transition probability density function Eq. 3 of tree diameter size and the steady state probability density Eq. 4, namely, how far from each other these two probability densities are? To analyze dissimilarity of probability densities (Eq. 3, 4), we will use a Monte Carlo simulation approach. The Monte Carlo method (Rubinstein and Kroese, 2008) allows us to find out a lower (upper) bound of the all used steady state information measures (Eq. 5-8).
The simulation procedure for calculating the steady state lower (upper) bound of the Shannon, Tsallis, Kullack, (Fisher) information measures can be described in the following steps:
| • | Estimate parameters of Eq. 1 using a real data set of n observations |
| • | Generate large number of samples (e.g., s = 1000) from the uniformly distributed on (0,1) of size n equal to the number of observed data, uij, j = 1,2,...,s, i = 1,2,...,n |
| • | Calculate percentiles qij, j = 1,2,…,s,
i = 1,2,…,n by the following form qij =
F-1 (uij), where |
| • | Estimate parameters of Eq. 1 using the simulated data sets { q1j, q2j,… qnj}, j = 1,2,…,s; |
| • | Calculate 95, 99 and 99.5% lower (upper) bounds for the Shannon, Tsallis, Kullback, (Fisher) steady state information measures |
| • | Find 95, 99 and 99.5% hitting ages |
RESULTS AND DISCUSSION
The numerical analysis of tree diameter growth is based on experiments in Pine stands at Dubrava district in Lithuania. The data were provided by the Lithuanian National Forest Inventory, which consisted of a systematic sample plots distributed on a square grid of 5 km, with a 5 year measurement interval. The sample method used circular plots of fixed radius with the area 500 m2 (Kuliešis et al., 2003). The data analyzed were collected during 1976. During this time these stands. have been remeasured 5 times on stand variables: age, number of trees per hectare, breast height diameter, trees position co-ordinates, age and height. The measurements have been conducted in 34 permanent treatment plots. The observed data of study plots are shown in Fig. 1.
| Fig. 1: | Plot of the diameter including data from pine forests in Lithuania |
The estimation method that we followed is first estimate by the least squares estimate method the parameters of the deterministic model (σ = 0), which represents the drift term of the stochastic model (1) and then use the maximum likelihood procedure to estimate the parameter σ, which represents the amplitude of diffusion (Rupšys et al., 2007). Practically, we could estimate simultaneously all parameters appearing in Eq. 1 by the maximum likelihood procedure. Unfortunately, the precise of such estimate suffers from the adjusted logarithmic transformation of the original process X(t) and the linearization of the drift term of a transformed process.
In this study the analysis of models` precision is based on the data used to fit them. We will examine two statistics: coefficient of determination for nonlinear regression (R2) and relative error (RE%) as follows:
where, p is the number of model parameters, the mean mt and variance st of tree diameter size are defined by the following system of ordinary differential equations:
(9) |
The parameter estimates and their corresponding goodness of fit statistics are shown in Table 1. The age-dependent evolution of mean and variance of tree diameter size showed that they increase to the steady state mean, 37.9902 and variance, 455.3135 (Fig. 2). As we can see from system Eq. 9, the mean tree diameter size mt of a forest stand is extremely sensitive with respect to the size of the coefficient of volatility σ. An analysis of goodness of fit statistics shown in Table 1 contrasts the performance of the stochastic Gompertz shape diameter growth law against the deterministic Gompertz shape diameter growth law.
| Table 1: | Estimations of parameters |
| Fig. 2: | (a) Plot of the mean trajectory (solid line-stochastic, dash line-deterministic) and (b) Variance trajectory |
Now we define an optimal rotation period as the age T when the age-dependent non-steady state differential information measures Eq. (5-8) fall the steady state lower (upper) bounds. Under this assumption, the optimal rotation period can be formulated as the minimal value (age) T which satisfies the inequality:
(10) |
where, b(S), b(T2), b(H), b(F) are lower (upper) bound of the Shannon, Tsallis, Kullback, (Fisher) steady state information measure (Eq. 5-8). An exact solution of inequality (Eq. 10) does not exist. This solution is simulated numerically and displayed in the graphic form.
As we already pointed out the defined test of optimal rotation period is managed on four information measures. A simulation study was conducted to investigate the order of convergence of the differential information measures (Eq. 5-8) to their steady state lower (upper) bound values. In order to evaluate the steady state lower (upper) bound of the Shannon, Tsallis, Kullback, (Fisher) information measures, a Monte Carlo study was performed. Based on the considerations of the above, were simulated 1000 samples (n = 1581) from the steady state probability density function Eq. 4. Using these simulations were computed the steady state lower (upper) bounds of the Shannon, Tsallis, Kullback, (Fisher) information measures (Table 2) and plotted jointly with the curve of the differential information measures (Eq. 5-8). We considered three cases (99.5, 99 and 95% confidence levels) of the lower (upper) bounds. The non-steady state information measures (Eq. 5-8) and their steady state lower (upper) bounds are shown in Fig. 3. As we can show in Fig. 3 the curves of non-steady state information measures (Eq. 5-8) fall their steady state lower (upper) bounds. These intersection points are shown in Table 3.
| Table 2: | Simulated values of information measures |
| ***99.5% confidence level, **99% confidence level, *95% confidence level | |
| Table 3: | Optimal rotation periods |
| ***99.5% confidence level, **99% confidence level, *95% confidence level | |
| Fig. 3: | Plot of information measures vs. age t and corresponding steady state lower (upper) bound values (99.5% confidence level) |
All aforementional information measures are closely related. As we show in Table 3, the optimal harvesting defined by Shannon`s entropy has minimal periods and optimal harvesting defined by Fisher`s information has maximal periods.
CONCLUSION
In this study is suggested a new approach in investigating the dynamics of the forest stands on the basis of the Shannon, Tsallis, Kullback, Fisher information measures (Eq. 5-8).
We focused on the forest stands in which human activity is based by uneconomical motives.
The smaller Shannon`s, Tsallis` differential entropy is, the more concentrated is the probability density function of tree diameter size of forest stands; the larger the Fisher`s information is, the smaller is uncertainty and the higher is the accuracy in predicting the localization of the tree diameter size of forest stands. Moreover, these three functionals of tree diameter size are closely connected with other diameter density functionals, which characterize the first and second moments of diameter size.
Most forest stands reach their economically optimal harvesting period prior to the defined biologically optimal harvesting period.
The optimal harvesting procedure defined by the Shannon entropy has periods of minimal size and optimal harvesting defined by the Fisher information has periods of maximal size.
One additional useful observation can be made from the non-steady state probability density (Eq. 3). Site quality affects the biologically optimal rotation periods. Hence, on better sites, where the diameter intrinsic growth rates r are higher, less biologically optimal periods as larger (lesser) the value of the Shannon, Tsallis, Kullback, Fisher information measures (Eq. 5-8).
The results presented in this study encourage the use of these techniques for analyzing other renewable resources such as fish. Exploitation of fish and forest resources has common parallels in maintaining equilibrium of growth and harvest.