A Multivariate Hybrid Stochastic Differential Equation Model for Whole-Stand Dynamics

Abstract

The growth and yield modeling of a forest stand has progressed rapidly, starting from the generalized nonlinear regression models of uneven/even-aged stands, and continuing to stochastic differential equation (SDE) models. We focus on the adaptation of the SDEs for the modeling of forest stand dynamics, and relate the tree and stand size variables to the age dimension (time). Two different types of diffusion processes are incorporated into a hybrid model in which the shortcomings of each variable types can be overcome to some extent. This paper presents the hybrid multivariate SDE regarding stand basal area and volume models in a forest stand. We estimate the fixed- and mixed-effect parameters for the multivariate hybrid stochastic differential equation using a maximum likelihood procedure. The results are illustrated using a dataset of measurements from Mountain pine tree (Pinus mugo Turra).

1. Introduction

Modeling growth and yield in a forest stand has progressed rapidly, starting from generalized nonlinear regression models of uneven/even-aged stands, and continuing to stochastic differential equations (SDEs) models and artificial neural network (ANN) models [1,2,3]. Nowadays, forest management needs information on the outcomes of different decisions to predict changes in forest structure and yield. All regression models continue to be generalized and refined [4]. Most of the empirical models developed using regression were designed to properly and precisely reproduce observed data sets. In contrast, the adaptation of the SDE for stand growth as an analogy has a shorter history, and aims to relate stand attributes structure to time (age) [5]. Regression and artificial neural network models are able to address questions using the linear or nonlinear relationships between model variables, and they are restricted by a single problem. An individual model is a mathematical equation describing the dynamics of specific components, such as diameter, height, crown width, tree density (number of trees per hectare), basal area, volume, and others. Regression and artificial neural network models suffer from the fixed forms of the cause and effect relationship. The theoretical methodology for both techniques involves trying a variety of equations and choosing the best fitting equation based on particular statistical measures. The limitations of these techniques are that they are laborious and the empirical choices of candidate equations make the results subjective [6].

The multivariate SDE model enables the prediction of the future values of certain outputs, such as diameter, height, tree density, stand basal area, stand volume per hectare, their mean and current annual increases, and more, at a particular age. The dynamics of the transition probability density function governed by the SDE determines the degree of stand variable predictability. Simple quantitative measures of the predictive capability, such as the mean and central moments, of stand variables can be effectively calculated using a derived multivariate probability density function and its marginal univariate and bivariate. In this study, we focused on multivariate continuous-time Markov processes that are of the Ornstein–Uhlenbeck family [7]. Typically, the SDE is considered an ordinary differential equation with a white noise variable, which incorporates an influence that seems random. The possibility of combining SDEs, sophisticated mathematical techniques of parameter estimates, and increased computing power have produced an advanced research methodology for understanding how the stochastic phenomenon affects the predictions in practical applications [8,9,10].

Forestry systems, and particularly those regarding tree density and stem profile, are very noisy [11,12]. The deterministic mechanisms involved in stand variables cannot be fully determined, so they function as a multivariate stochastic process. The discrete time Markov process has been most frequently used in stochastic description of stand variables [13,14]. A few years later, the continuous time univariate Markov process was used to predict the future diameter distributions [15,16]. Continuous time stochastic growth models are useful in the modeling of a stand growth, as the dynamics we are interested in is some interval of the real line. The multivariate stochastic process has been used to study more complex sets of data than can be handled by univariate processes. Understanding the conditional multivariate process possibilities provides the ability to choose response and predictor variables and incorporates the factor of tree size variable dependence using the variance–covariance matrix of tree or stand variables. Historically, the systematic application of the SDE models in forestry uses many different symmetric and asymmetric types that concern sigmoid growth consisting of three key stages: exponential, transitional, and plateau [16,17]. The variables involved in the SDE are therefore treated as diffusion processes and the solutions are described in terms of normal or lognormal probability density functions. Therefore, the central problem is determining how to frame a multivariate SDE whose solution has, e.g., normally-, as well as, lognormally-distributed marginals. In this study, we simultaneously assimilated two sets of the Gompertz and Vasicek SDEs. This is possible through defining a hybrid system of SDEs. Statistical estimation parameters of the observed data sets can be based upon the maximum likelihood procedure by maximizing the maximum log-likelihood (fixed effect scenario) or the approximated log-likelihood function (mixed effect scenario) that consists of the simultaneous assimilation of normal and lognormal probability density functions.

The general methodology of a stationary distribution formalizes an equilibrium structure or steady-state structure of variables, which is commonly used in biology, chemistry, and other branches of science [18]. This concept can be recognized as a constant form probability density function corresponding to the solution of the multivariate SDE, the structure of which remains invariant over time.

The main goal of this paper is to develop a multivariate hybrid mixed-effect parameters SDE model that links different types of the drift and diffusion functions, and to describe the maximum likelihood procedure for fixed-and mixed-effect parameter estimators. First, we find exact form solutions of the multivariate hybrid SDE. Second, we use the newly developed multivariate and conditional hybrid probability density functions for describing the mathematical equations of the dynamics of the stand attributes such as the mean, median, mode, standard deviation, quantiles of the tree density, diameter, height, basal area, and stand volume per hectare.

In the Results and Discussion sections, we consider a possible application of the multivariate hybrid SDE to measurements from Mountain pine (Pinus mugo Turra) stands in Lithuania.

2. Materials and Methods

2.1. Hybrid Model and Its Characteristics

In producing a trivariate hybrid SDE model for predictions of stand attributes over time, it would be instructive to start from fundamental variables related to all other aspects of forest stand development. Traditionally the tree density N, diameter at breast height D, and height H are variables substantially affecting the development of whole stand growth and yield models [19]. Each variable can be analyzed as a separate diffusion process. The multivariate probability density function of the tree density, diameter, and height provides full-scale information about the cause and effect variables. The height–diameter equation or the diameter–height equation can be derived using conditional stationary SDEs. Here, we focus on a multivariate mixed effect parameters hybrid SDE model (Gompertz, Gompertz, and Vasicek, see [17,20,21]), $X^i\left(t\right)={\left(X_1^i\left(t\right),X_2^i\left(t\right),X_3^i\left(t\right)\right)}^T={\left(N^i\left(t\right),D^i\left(t\right),H^i\left(t\right)\right)}^T$ of the tree density, $N^i\left(t\right)$, diameter, $D^i\left(t\right)$, and height, $H^i\left(t\right)$ dynamics as follows:

$$d{X}^i\left(t\right)=A^i\left(X^i\left(t\right)\right)dt+D\left(X^i\left(t\right)\right)B^{\frac{1}{2}}\cdot d{W}^i\left(t\right), P\left(X^i\left(t_0\right)=x_0\right)=1, i=1, \dots , M,$$

(1)

$$A^i\left(x\right)={\left(\left(\left({\alpha }_1+{\phi }_1^i\right)-{\beta }_{1}ln(x_1)\right)x_1,\left(\left({\alpha }_2+{\phi }_2^i\right)-{\beta }_{2}ln(x_2)\right)x_2,{\beta }_3\left({\alpha }_3+{\phi }_3^i-x_3\right)\right)}^T,$$

(2)

$$B=\left(\begin{array}{ccc}{\sigma }_{11}& {\sigma }_{12}& {\sigma }_{13}\\ {\sigma }_{12}& {\sigma }_{22}& {\sigma }_{23}\\ {\sigma }_{13}& {\sigma }_{23}& {\sigma }_{33}\end{array}\right),$$

(3)

$$D\left(x\right)=\left(\begin{array}{ccc}{x}_1& 0& 0\\ 0& x_2& 0\\ 0& 0& 1\end{array}\right),$$

(4)

$$G\left(x\right)=\left(D\left(x\right)B^{\frac{1}{2}}\right){\left(D\left(x\right)B^{\frac{1}{2}}\right)}^T=\left(\begin{array}{ccc}{x}_1& 0& 0\\ 0& x_2& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{ccc}{\sigma }_{11}& {\sigma }_{12}& {\sigma }_{13}\\ {\sigma }_{12}& {\sigma }_{22}& {\sigma }_{23}\\ {\sigma }_{13}& {\sigma }_{23}& {\sigma }_{33}\end{array}\right)\left(\begin{array}{ccc}{x}_1& 0& 0\\ 0& x_2& 0\\ 0& 0& 1\end{array}\right),$$

(5)

where $W^i\left(t\right)={\left(W_1^i\left(t\right),W_2^i\left(t\right),W_3^i\left(t\right)\right)}^T$ are independent multivariate Brownian motions, $t\in \left[t_0;T\right]$ is a finite horizon, T < ∞, $B^{\frac{1}{2}}$ is the Cholesky factorization for a positive definite symmetric matrix B into the product of a lower triangular matrix and its conjugate transpose, and ${\phi }_j^i$, i = 1, …, M, j = 1, …, 3, are independent and normally distributed random variables with zero mean and constant variances (${\phi }_j^i~N\left(0;{\sigma }_j^2\right)$), an initial condition takes the form; if t = t₀, then $X^i\left(t_0\right)=x_0=\left(x_{10},x_{20},x_{30}\right)={\left(n_0,d_0,h_0\right)}^T$, ${\sigma }_{ij}=\sqrt{{\sigma }_{ii}{\sigma }_{jj}}{\rho }_{ij}$, and $\left\{{\alpha }_1,{\alpha }_2,{\alpha }_3,{\beta }_1,{\beta }_2,{\beta }_3,{\sigma }_{11},{\rho }_{12},{\rho }_{13},{\sigma }_{22},{\rho }_{23},{\sigma }_{33},{\sigma }_1,{\sigma }_2,{\sigma }_3\right\}$ are fixed effect parameters to be estimated.

As the drift and diffusion functions of Equation (1) confirm the conditions of existence and uniqueness of the strong solution, then the SDE in Equation (1) has a unique strong solution for any initial condition x₀ [22,23]. The strategy to solve Equation (1) consists of a specific transform of the type $Y^i\left(t\right)={\left(e^{{\beta }_{1}t}ln\left(X_1^i\left(t\right)\right),e^{{\beta }_{2}t}ln\left(X_2^i\left(t\right)\right),e^{{\beta }_{3}t}{X}_3^i\left(t\right)\right)}^T$. By applying Ito’s formula [24] to the transformation $Y^i\left(t\right)$ the solution takes the following form:

$$\begin{gathered} Y^i\left(t\right)=Y^i\left(t_0\right)+\left(\underset{\_}{e^{\beta t}}-\underset{\_}{e^{\beta t_0}}\right)\underset{\_}{\alpha }+B^{\frac{1}{2}}\cdot {\displaystyle {{\displaystyle \int }}_{t_0}^t}\underset{\_}{e^{\beta u}}dW\left(u\right), \\ P\left(Y_1^i\left(t_0\right)=ln\left(x_{10}\right),Y_2^i\left(t_0\right)=ln\left(x_{20}\right),Y_3^i\left(t_0\right)=e^{\beta t_0}{x}_{30}\right)=1, i=1, \dots , M, \end{gathered}$$

(6)

where $\underset{\_}{\alpha }={\left(\frac{1}{{\beta }_1}\left({\alpha }_1+{\phi }_1^i-\frac{{\sigma }_{11}}{2}\right),\frac{1}{{\beta }_2}\left({\alpha }_2+{\phi }_2^i-\frac{{\sigma }_{22}}{2}\right),{\alpha }_3+{\phi }_3^i\right)}^T$, and $\underset{\_}{e^{\beta t}}=\left(\begin{array}{ccc}{e}^{{\beta }_{1}t}& 0& 0\\ 0& e^{{\beta }_{2}t}& 0\\ 0& 0& e^{{\beta }_{3}t}\end{array}\right)$. Considering the transformation $Y^i\left(t\right)={\left(e^{{\beta }_{1}t}ln\left(X_1^i\left(t\right)\right),e^{{\beta }_{2}t}ln\left(X_2^i\left(t\right)\right),e^{{\beta }_{3}t}{X}_3^i\left(t\right)\right)}^T$ and that the last term in Equation (6) is normally distributed, we deduce that the conditional random vector $\left(X^i\left(t\right)|X^i\left(t_0\right)=x_0\right)$ has a trivariate hybrid lognormal–lognormal–normal distribution $L{N}_2{N}_1\left({\mu }^i\left(t\right); \Sigma \left(t\right)\right)$, with the mean vector ${\mu }^i\left(t\right)$ defined by:

$$\begin{array}{c}\hfill {\mu }^i\left(t\right)={\left({\mu }_1^i\left(t\right),{\mu }_2^i\left(t\right),{\mu }_3^i\left(t\right)\right)}^T\hfill \\ =(\left(e^{-{\beta }_j\left(t-t_0\right)}ln\left(x_{j0}\right)+\frac{1-e^{-{\beta }_j\left(t-t_0\right)}}{{\beta }_j}\left({\alpha }_j+{\phi }_j^i-\frac{{\sigma }_{jj}}{2}\right), j=1,2\right),{\alpha }_3+{\phi }_3^i+\hfill \\ \left(x_{30}-\left({\alpha }_3+{\phi }_3^i\right)\right)e^{-{\beta }_3\left(t-t_0\right)}{)}^T,\hfill \end{array}$$

(7)

the variance-covariance matrix $\Sigma \left(t\right)$:

$$\Sigma \left(t\right)={\left(v_{jk}\left(t\right)\right)}_{j,k=1,\dots ,3}={\left(\frac{{\sigma }_{jk}}{{\beta }_j+{\beta }_k}\left(1-e^{-\left({\beta }_j+{\beta }_k\right)\left(t-t_0\right)}\right)\right)}_{j,k=1,\dots ,3},$$

(8)

and the transition probability density function:

$$f\left(x_1,x_2,x_3,t|{\theta }^1,{\phi }^i\right)=\frac{1}{{\left(2\pi \right)}^{\frac{3}{2}}{\left|\Sigma \left(t\right)\right|}^{\frac{1}{2}}{x}_1{x}_2}exp\left(-\frac{1}{2}\Omega \left(x_1,x_2,x_3,t{|\theta }^1,{\phi }^i\right)\right),$$

(9)

$$\Omega \left(x_1,x_2,x_3,t|{\theta }^1,{\phi }^i\right)={\left(\begin{array}{c}ln\left(x_1\right)-{\mu }_1^i\left(t\right)\\ ln\left(x_2\right)-{\mu }_2^i\left(t\right)\\ x_3-{\mu }_3^i\left(t\right)\hfill \end{array}\right)}^T{\left(\Sigma \left(t\right)\right)}^{-1}\left(\begin{array}{c}ln\left(x_1\right)-{\mu }_1^i\left(t\right)\\ ln\left(x_2\right)-{\mu }_2^i\left(t\right)\\ x_3-{\mu }_3^i\left(t\right)\hfill \end{array}\right),$$

(10)

$${\theta }^1=\left\{{\alpha }_1,{\beta }_1,{\alpha }_2,{\beta }_2,{\alpha }_3,{\beta }_3,{\sigma }_{11},{\rho }_{12},{\rho }_{13},{\sigma }_{22},{\rho }_{23},{\sigma }_{33}\right\}, {\phi }^i=\left({\phi }_1^i,{\phi }_2^i,{\phi }_3^i\right).$$

(11)

The univariate marginal distributions of $\left(X_j^i\left(t\right)|X_j^i\left(t_0\right)=x_{j0}\right)$ j = 1,2 are lognormal $L{N}_1\left({\mu }_j^i\left(t\right); v_{jj}\left(t\right)\right)$, and the marginal distribution of $\left(X_3^i\left(t\right)|X_3^i\left(t_0\right)=x_{30}\right)$ is normal $N_1\left({\mu }_3^i\left(t\right); v_{33}\left(t\right)\right)$. The marginal mean, median, mode, p-quantile (0 < p < 1), and variance trajectories $m_j^i\left(t\right)$, $m{e}_j^i\left(t\right)$, $m{o}_j^i\left(t\right)$, $m{q}_j^i\left(t,p\right)$, and $w_j\left(t\right)$ of the tree density, diameter, and height ($j\in \left\{1,2,3\right\}$) are listed in

Table 1. Marginal mean, median, mode, p-quantile (0 < p < 1), and variance trajectories.

Variable	Trajectory Type	Equation
$X_3^i\left(t\right)$	Mean, median and mode	${\mu }_3^i\left(t\right)$
	Quantile (0 < p < 1)	${\Phi }_p^{-1}\left({\mu }_3^i\left(t\right);v_{33}\left(t\right)\right)$ *
	Variance	$v_{33}\left(t\right)$
$X_j^i\left(t\right),$ $j=1,2$	Mean	$exp\left({\mu }_j^i\left(t\right)+\frac{1}{2}{v}_{jj}\left(t\right)\right)$
	Median	$exp\left({\mu }_j^i\left(t\right)\right)$
	Mode	$exp\left({\mu }_j^i\left(t\right)-v_{jj}\left(t\right)\right)$
	Quantile (0 < p < 1)	$L{\Phi }_p^{-1}\left({\mu }_j^i\left(t\right);v_{jj}\left(t\right)\right)$ *
	Variance	$exp\left(2{\mu }_j^i\left(t\right)+v_{jj}\left(t\right)\right)\cdot \left(exp\left(v_{jj}\left(t\right)\right)-1\right)$

* ${\Phi }_p^{-1}\left(\cdot ;\cdot \right)$ is the inverse of the standard normal distribution function; and $L{\Phi }_p^{-1}\left(\cdot ;\cdot \right)$ is the inverse of the lognormal distribution function.

.

The marginal bivariate distribution of $\left(X_1^i\left(t\right),X_2^i\left(t\right)|X_1^i\left(t_0\right)=x_{10},X_2^i\left(t_0\right)=x_{20}\right)$ is lognormal $L{N}_2\left({\mu }^{i,12}\left(t\right); {\Sigma }^{12}\left(t\right)\right),$ with the mean vector ${\mu }^{i,12}\left(t\right)$ and the covariance matrix ${\Sigma }^{12}\left(t\right)$ defined as

$${\mu }^{i,12}\left(t\right)={\left(e^{-{\beta }_j\left(t-t_0\right)}ln\left(x_{j0}\right)+\frac{1-e^{-{\beta }_j\left(t-t_0\right)}}{{\beta }_j}\left({\alpha }_j+{\phi }_j^i-\frac{{\sigma }_{jj}}{2}\right), j=1,2\right)}^T$$

(12)

$${\Sigma }^{12}\left(t\right)={\left(v_{j,k}\left(t\right)\right)}_{j,k=1,2}.$$

(13)

The marginal bivariate distribution of $\left(X_1^i\left(t\right),X_3^i\left(t\right)|X_1^i\left(t_0\right)=x_{10},X_3^i\left(t_0\right)=x_{30}\right)$, is hybrid lognormal-normal $L{N}_1{N}_1\left({\mu }^{i,13}\left(t\right); {\Sigma }^{13}\left(t\right)\right)$ with mean vector ${\mu }^{i,13}\left(t\right)$ and covariance matrix ${\Sigma }^{13}\left(t\right)$ defined in the following forms:

$$\begin{gathered} {\mu }^{i,13}\left(t\right) \\ ={\left(e^{-{\beta }_1\left(t-t_0\right)}ln\left(x_{10}\right)+\frac{1-e^{-{\beta }_1\left(t-t_0\right)}}{{\beta }_1}\left({\alpha }_1+{\phi }_1^i-\frac{{\sigma }_{11}}{2}\right),{\alpha }_3+{\phi }_3^i \\ +\left(x_{30}-\left({\alpha }_3+{\phi }_3^i\right)\right)e^{-{\beta }_3\left(t-t_0\right)}\right)}^T \end{gathered}$$

(14)

$${\Sigma }^{13}\left(t\right)={\left(v_{j,k}\left(t\right)\right)}_{j,k=1,3}.$$

(15)

The marginal bivariate distribution of $\left(X_2^i\left(t\right),X_3^i\left(t\right)|X_2^i\left(t_0\right)=x_{20},X_3^i\left(t_0\right)=x_{30}\right)$ is hybrid lognormal–normal $L{N}_1{N}_1\left({\mu }^{i,23}\left(t\right); {\Sigma }^{23}\left(t\right)\right),$ with mean vector ${\mu }^{i,23}\left(t\right)$ and covariance matrix ${\Sigma }^{23}\left(t\right)$ defined in the following forms:

$$\begin{gathered} {\mu }^{i,23}\left(t\right) \\ ={\left(e^{-{\beta }_2\left(t-t_0\right)}ln\left(x_{10}\right)+\frac{1-e^{-{\beta }_2\left(t-t_0\right)}}{{\beta }_2}\left({\alpha }_2+{\phi }_2^i-\frac{{\sigma }_{22}}{2}\right),{\alpha }_3+{\phi }_3^i \\ +\left(x_{30}-\left({\alpha }_3+{\phi }_3^i\right)\right)e^{-{\beta }_3\left(t-t_0\right)}\right)}^T \end{gathered}$$

(16)

$${\Sigma }^{23}\left(t\right)={\left(v_{j,k}\left(t\right)\right)}_{j,k=2,3}.$$

(17)

The univariate conditional distribution of $\left(X_j^i\left(t\right)|X_j^i\left(t_0\right)=x_{j0}\right)$ j = 1,2 at a given $\left(X_3^i\left(t\right)=x_3\right)$ is univariate lognormal $L{N}_1\left({\eta }_j^i\left(t,x_3\right); {\lambda }_{j3}\left(t\right)\right)$; the univariate conditional distribution of $\left(X_1^i\left(t\right)|X_1^i\left(t_0\right)=x_{10}\right)$ at a given $\left(X_2^i\left(t\right)=x_2\right)$ is univariate lognormal $L{N}_1\left({\eta }_1^i\left(t,x_2\right); {\lambda }_{12}\left(t\right)\right)$, the univariate conditional distribution of $\left(X_2^i\left(t\right)|X_2^i\left(t_0\right)=x_{20}\right)$ at a given $\left(X_1^i\left(t\right)=x_1\right),$ is univariate lognormal $L{N}_1\left({\eta }_2^i\left(t,x_1\right); {\lambda }_{21}\left(t\right)\right)$, the univariate conditional distribution of $\left(X_3^i\left(t\right)|X_3^i\left(t_0\right)=x_{30}\right)$ at a given $\left(X_1^i\left(t\right)=x_1\right)$ is univariate normal $N_1\left({\eta }_3^i\left(t,x_1\right); {\lambda }_{31}\left(t\right)\right)$, and the univariate conditional distribution of $\left(X_3^i\left(t\right)|X_3^i\left(t_0\right)=x_{30}\right)$ at a given $\left(X_2^i\left(t\right)=x_2\right)$ is univariate normal $N_1\left({\eta }_3^i\left(t,x_2\right); {\lambda }_{32}\left(t\right)\right)$, with the mean and variance defined in

Table 2. Mean and variance of conditional probability density functions with one predictor.

Probability Density	$\mathbf{Mean} {\mathit{\eta }}_{\mathit{j}}^{\mathit{i}}\left(\mathit{t},{\mathit{x}}_{\mathit{k}}\right)$	$\mathbf{Variance} {\mathit{\lambda }}_{\mathit{j}\mathit{k}}\left(\mathit{t}\right)$
$L{N}_1\left({\eta }_1^i\left(t,x_2\right); {\lambda }_{12}\left(t\right)\right)$	${\mu }_1^i\left(t\right)+\frac{v_{12}\left(t\right)}{v_{22}\left(t\right)}\left(ln\left(x_2\right)-{\mu }_2^i\left(t\right)\right)$	$v_{11}\left(t\right)-\frac{{\left(v_{12}\left(t\right)\right)}^2}{v_{22}\left(t\right)}$
$L{N}_1\left({\eta }_1^i\left(t,x_3\right); {\lambda }_{13}\left(t\right)\right)$	${\mu }_1^i\left(t\right)+\frac{v_{13}\left(t\right)}{v_{33}\left(t\right)}\left(x_3-{\mu }_3^i\left(t\right)\right)$	$v_{11}\left(t\right)-\frac{{\left(v_{13}\left(t\right)\right)}^2}{v_{33}\left(t\right)}$
$L{N}_1\left({\eta }_2^i\left(t,x_1\right); {\lambda }_{21}\left(t\right)\right)$	${\mu }_2^i\left(t\right)+\frac{v_{12}\left(t\right)}{v_{11}\left(t\right)}\left(ln\left(x_1\right)-{\mu }_1^i\left(t\right)\right)$	$v_{22}\left(t\right)-\frac{{\left(v_{12}\left(t\right)\right)}^2}{v_{11}\left(t\right)}$
$L{N}_1\left({\eta }_2^i\left(t,x_3\right); {\lambda }_{23}\left(t\right)\right)$	${\mu }_2^i\left(t\right)+\frac{v_{23}\left(t\right)}{v_{33}\left(t\right)}\left(x_3-{\mu }_3^i\left(t\right)\right)$	$v_{22}\left(t\right)-\frac{{\left(v_{23}\left(t\right)\right)}^2}{v_{33}\left(t\right)}$
$N_1\left({\eta }_3^i\left(t,x_1\right); {\lambda }_{31}\left(t\right)\right)$	${\mu }_3^i\left(t\right)+\frac{v_{13}\left(t\right)}{v_{11}\left(t\right)}\left(ln\left(x_1\right)-{\mu }_1^i\left(t\right)\right)$	$v_{33}\left(t\right)-\frac{{\left(v_{13}\left(t\right)\right)}^2}{v_{11}\left(t\right)}$
$N_1\left({\eta }_3^i\left(t,x_2\right); {\lambda }_{32}\left(t\right)\right)$	${\mu }_3^i\left(t\right)+\frac{v_{23}\left(t\right)}{v_{22}\left(t\right)}\left(ln\left(x_2\right)-{\mu }_2^i\left(t\right)\right)$	$v_{33}\left(t\right)-\frac{{\left(v_{23}\left(t\right)\right)}^2}{v_{22}\left(t\right)}$

. All Table 2 univariate conditional distributions are based on one predictor.

The conditional mean, median, mode, p-quantile (0 < p < 1), and variance trajectories $m_j^i\left(t,x_k\right)$, $m{e}_j^i\left(t,x_k\right)$, $m{o}_j^i\left(t,x_k\right)$, $m{q}_j^i\left(t,p,x_k\right),$ and $w_j^i\left(t,x_k\right)$ based on one predictor variable for all scenarios listed in Table 2 are presented in

Table 3. Conditional mean, median, mode, p-quantile (0 < p < 1), and variance trajectories with one predictor.

Probability Density	Trajectory Type	Equation
$N_1\left({\eta }_3^i\left(t,x_j\right); {\lambda }_{3j}\left(t\right)\right)$ $j\in \left\{1,2\right\}$	Mean, median, and mode	${\eta }_3^i\left(t,x_j\right)$
	Quantile (0 < p < 1)	${\Phi }_p^{-1}\left({\eta }_3^i\left(t,x_j\right); {\lambda }_{3j}\left(t\right)\right)$
	Variance	${\lambda }_{3j}\left(t\right)$
$L{N}_1\left({\eta }_j^i\left(t,x_k\right); {\lambda }_{jk}\left(t\right)\right)$ $j,k\in \left\{\left(1,2\right),\left(1,3\right),\left(2,1\right),\left(2,3\right)\right\}$	Mean	$exp\left({\eta }_j^i\left(t,x_k\right)+\frac{1}{2}{\lambda }_{jk}\left(t\right)\right)$
	Median	$exp\left({\eta }_j^i\left(t,x_k\right)\right)$
	Mode	$exp\left({\eta }_j^i\left(t,x_k\right)-{\lambda }_{jk}\left(t\right)\right)$
	Quantile (0 < p < 1)	$L{\Phi }_p^{-1}\left({\eta }_j^i\left(t,x_k\right); {\lambda }_{jk}\left(t\right)\right)$
	Variance	$exp\left(2{\eta }_j^i\left(t,x_k\right)+{\lambda }_{jk}\left(t\right)\right)\cdot \left(exp\left({\lambda }_{jk}\left(t\right)\right)-1\right)$

.

The conditional mean, median, mode, p-quantile (0 < p < 1), and variance trajectories $m_j^i\left(t,x_k,x_m\right)$, $m{e}_j^i\left(t,x_k,x_m\right)$, $m{o}_j^i\left(t,x_k,x_m\right)$, $m{q}_j^i\left(t,p,x_k,x_m\right),$ and $w_j^i\left(t,x_k,x_m\right)$ based on two predictors for all scenarios listed in

Table 4. Mean and variance of conditional probability density functions with two predictors.

Probability Density	$\mathbf{Mean} {\mathit{\Pi }}_{\mathit{j}}^{\mathit{i}}\left(\mathit{t},{\mathit{x}}_{\mathit{k}},{\mathit{x}}_{\mathit{m}}\right)$	$\mathbf{Variance} {\mathit{\Lambda }}_{\mathit{j}\mathit{k}\mathit{m}}\left(\mathit{t}\right)$
$L{N}_1\left({\Pi }_1^i\left(t,x_2,x_3\right); {\Lambda }_{123}\left(t\right)\right)$	${\mu }_1^i\left(t\right)+\left(\begin{array}{cc}{v}_{12}\left(t\right)& v_{13}\left(t\right)\end{array}\right){\left(\begin{array}{cc}{v}_{22}\left(t\right)& v_{23}\left(t\right)\\ v_{23}\left(t\right)& v_{33}\left(t\right)\end{array}\right)}^{-1}\left(\begin{array}{c}ln\left(x_2\right)-{\mu }_2^i\left(t\right)\\ x_3-{\mu }_3^i\left(t\right)\hfill \end{array}\right)$	$v_{11}\left(t\right)-\left(\begin{array}{cc}{v}_{12}\left(t\right)& v_{13}\left(t\right)\end{array}\right){\left(\begin{array}{cc}{v}_{22}\left(t\right)& v_{23}\left(t\right)\\ v_{23}\left(t\right)& v_{33}\left(t\right)\end{array}\right)}^{-1}\left(\begin{array}{c}{v}_{12}\left(t\right)\\ v_{13}\left(t\right)\end{array}\right)$
$L{N}_1\left({\Pi }_2^i\left(t,x_1,x_3\right); {\Lambda }_{213}\left(t\right))\right)$	${\mu }_2^i\left(t\right)+\left(\begin{array}{cc}{v}_{12}\left(t\right)& v_{23}\left(t\right)\end{array}\right){\left(\begin{array}{cc}{v}_{11}\left(t\right)& v_{13}\left(t\right)\\ v_{13}\left(t\right)& v_{33}\left(t\right)\end{array}\right)}^{-1}\left(\begin{array}{c}ln\left(x_1\right)-{\mu }_1^i\left(t\right)\\ x_3-{\mu }_3^i\left(t\right)\hfill \end{array}\right)$	$v_{22}\left(t\right)-\left(\begin{array}{cc}{v}_{12}\left(t\right)& v_{23}\left(t\right)\end{array}\right){\left(\begin{array}{cc}{v}_{11}\left(t\right)& v_{13}\left(t\right)\\ v_{13}\left(t\right)& v_{33}\left(t\right)\end{array}\right)}^{-1}\left(\begin{array}{c}{v}_{12}\left(t\right)\\ v_{23}\left(t\right)\end{array}\right)$
$N_1\left({\Pi }_3^i\left(t,x_1,x_2\right); {\Lambda }_{312}\left(t\right)\right)$	${\mu }_3^i\left(t\right)+\left(\begin{array}{cc}{v}_{13}\left(t\right)& v_{23}\left(t\right)\end{array}\right){\left(\begin{array}{cc}{v}_{11}\left(t\right)& v_{12}\left(t\right)\\ v_{12}\left(t\right)& v_{22}\left(t\right)\end{array}\right)}^{-1}\left(\begin{array}{c}ln\left(x_1\right)-{\mu }_1^i\left(t\right)\\ ln\left(x_2\right)-{\mu }_2^i\left(t\right)\end{array}\right)$	$v_{33}\left(t\right)-\left(\begin{array}{cc}{v}_{13}\left(t\right)& v_{23}\left(t\right)\end{array}\right){\left(\begin{array}{cc}{v}_{11}\left(t\right)& v_{12}\left(t\right)\\ v_{12}\left(t\right)& v_{22}\left(t\right)\end{array}\right)}^{-1}\left(\begin{array}{c}{v}_{13}\left(t\right)\\ v_{23}\left(t\right)\end{array}\right)$

are presented in

Table 5. Conditional mean, median, mode, p-quantile (0 < p < 1), and variance trajectories with two predictors.

Probability Density	Trajectory Type	Equation
$N_1\left({\Pi }_3^i\left(t,x_1,x_2\right); {\Lambda }_{312}\left(t\right)\right)$	Mean, median, and mode	${\Pi }_3^i\left(t,x_1,x_2\right)$
	Quantile (0 < p < 1)	${\Phi }_p^{-1}\left({\Pi }_3^i\left(t,x_1,x_2\right); {\Lambda }_{312}\left(t\right)\right)$
	Variance	${\Lambda }_{312}\left(t\right)$
$L{N}_1\left({\Pi }_j^i\left(t,x_k,x_m\right); {\Lambda }_{jkm}\left(t\right)\right)$ $j,k,m\in \left\{\left(1,2,3\right),\left(2,1,3\right)\right\}$	Mean	$\exp \left({\Pi }_j^i\left(t,x_k,x_m\right)+\frac{1}{2}{\Lambda }_{jkm}\left(t\right)\right)$
	Median	$exp\left({\Pi }_j^i\left(t,x_k,x_m\right)\right)$
	Mode	$exp\left({\Pi }_j^i\left(t,x_k,x_m\right)-{\Lambda }_{jkm}\left(t\right)\right)$
	Quantile (0 < p < 1)	$L{\Phi }_p^{-1}\left({\Pi }_j^i\left(t,x_k,x_m\right); {\Lambda }_{jkm}\left(t\right)\right)$
	Variance	$exp\left(2{\Pi }_j^i\left(t,x_k,x_m\right)+{\Lambda }_{jkm}\left(t\right)\right)\cdot \left(exp\left({\Lambda }_{jkm}\left(t\right)\right)-1\right)$

.

For the Gompertz- and Vasicek-type SDEs, it is possible to derive the stationary univariate marginal and conditional distributions if parameters β₁, β₂, and β₃ are positive [25]. As time t goes to infinity, the trivariate diffusion process defined by Equation (1) assumes a stationary trivariate hybrid lognormal–lognormal–normal distribution $L{N}_2{N}_1\left({\mu }^i; \Sigma \right)$ with the mean vector ${\mu }^i$ defined by

$${\mu }^i={\left({\mu }_1^i,{\mu }_2^i,{\mu }_3^i\right)}^T={\left(\left(\frac{1}{{\beta }_j}\left({\alpha }_j+{\phi }_j^i-\frac{{\sigma }_{jj}}{2}\right), j=1,2\right),{\alpha }_3+{\phi }_3^i\right)}^T,$$

(18)

the variance–covariance matrix $\Sigma$:

$$\Sigma ={\left(v_{jk}\right)}_{j,k=1,\dots ,3}={\left(\frac{{\sigma }_{jk}}{{\beta }_j+{\beta }_k}\right)}_{j,k=1,\dots ,3},$$

(19)

and the stationary probability density function:

$$f\left(x_1,x_2,x_3|{\theta }^1,{\phi }^i\right)=\frac{1}{{\left(2\pi \right)}^{\frac{3}{2}}{\left|\Sigma \right|}^{\frac{1}{2}}{x}_1{x}_2}exp\left(-\frac{1}{2}\Omega \left(x_1,x_2,x_3{|\theta }^1,{\phi }^i\right)\right),$$

(20)

$$\Omega \left(x_1,x_2,x_3|{\theta }^1,{\phi }^i\right)={\left(\begin{array}{c}ln\left(x_1\right)-{\mu }_1^i\\ ln\left(x_2\right)-{\mu }_2^i\\ x_3-{\mu }_3^i\hfill \end{array}\right)}^T{\left(\Sigma \right)}^{-1}\left(\begin{array}{c}ln\left(x_1\right)-{\mu }_1^i)\\ ln\left(x_2\right)-{\mu }_2^i\\ x_3-{\mu }_3^i\hfill \end{array}\right).$$

(21)

The marginal distributions converge to a stationary distribution with the means and variances defined by Equations (18) and (19), and conditional distributions converge to a stationary distribution with the means and variances listed in

Table 6. Mean and variance of conditional probability density functions with one predictor.

Probability Density	$\mathbf{Mean} {\mathit{\eta }}_{\mathit{j}}^{\mathit{i}}\left(\mathit{t},{\mathit{x}}_{\mathit{k}}\right)$	$\mathbf{Variance} {\mathit{\lambda }}_{\mathit{j}\mathit{k}}\left(\mathit{t}\right)$
$L{N}_1\left({\eta }_1^i\left(x_2\right); {\lambda }_{12}\right)$	${\mu }_1^i+\frac{v_{12}}{v_{22}}\left(ln\left(x_2\right)-{\mu }_2^i\right)$	$v_{11}-\frac{{\left(v_{12}\right)}^2}{v_{22}}$
$L{N}_1\left({\eta }_1^i\left(x_3\right); {\lambda }_{13}\right)$	${\mu }_1^i+\frac{v_{13}}{v_{33}}\left(x_3-{\mu }_3^i\right)$	$v_{11}-\frac{{\left(v_{13}\right)}^2}{v_{33}}$
$L{N}_1\left({\eta }_2^i\left(x_1\right); {\lambda }_{21}\right)$	${\mu }_2^i+\frac{v_{12}}{v_{11}}\left(ln\left(x_1\right)-{\mu }_1^i\right)$	$v_{22}-\frac{{\left(v_{12}\right)}^2}{v_{11}}$
$L{N}_1\left({\eta }_2^i\left(x_3\right); {\lambda }_{23}\right)$	${\mu }_2^i+\frac{v_{23}}{v_{33}}\left(x_3-{\mu }_3^i\right)$	$v_{22}\left(t\right)-\frac{{\left(v_{23}\right)}^2}{v_{33}}$
$N_1\left({\eta }_3^i\left(x_1\right); {\lambda }_{31}\right)$	${\mu }_3^i+\frac{v_{13}}{v_{11}}\left(ln\left(x_1\right)-{\mu }_1^i\right)$	$v_{33}\left(t\right)-\frac{{\left(v_{13}\right)}^2}{v_{11}}$
$N_1\left({\eta }_3^i\left(x_2\right); {\lambda }_{32}\right)$	${\mu }_3^i\left(t\right)+\frac{v_{23}}{v_{22}}\left(ln\left(x_2\right)-{\mu }_2^i\right)$	$v_{33}\left(t\right)-\frac{{\left(v_{23}\right)}^2}{v_{22}}$

and

Table 7. Mean and variance of conditional probability density functions with two predictors.

Probability Density	$\mathbf{Mean} {\mathit{\Pi }}_{\mathit{j}}^{\mathit{i}}\left({\mathit{x}}_{\mathit{k}},{\mathit{x}}_{\mathit{m}}\right)$	$\mathbf{Variance} {\mathit{\Lambda }}_{\mathit{j}\mathit{k}\mathit{m}}$
$L{N}_1\left({\Pi }_1^i\left(x_2,x_3\right); {\Lambda }_{123}\right)$	${\mu }_1^i+\left(\begin{array}{cc}{v}_{12}& v_{13}\end{array}\right){\left(\begin{array}{cc}{v}_{22}& v_{23}\\ v_{23}\left(t\right)& v_{33}\end{array}\right)}^{-1}\left(\begin{array}{c}ln\left(x_2\right)-{\mu }_2^i\\ x_3-{\mu }_3^i\hfill \end{array}\right)$	$v_{11}-\left(\begin{array}{cc}{v}_{12}& v_{13}\end{array}\right){\left(\begin{array}{cc}{v}_{22}& v_{23}\\ v_{23}& v_{33}\end{array}\right)}^{-1}\left(\begin{array}{c}{v}_{12}\\ v_{13}\end{array}\right)$
$L{N}_1\left({\Pi }_2^i\left(x_1,x_3\right); {\Lambda }_{213})\right)$	${\mu }_2^i+\left(\begin{array}{cc}{v}_{12}& v_{23}\end{array}\right){\left(\begin{array}{cc}{v}_{11}& v_{13}\\ v_{13}\left(t\right)& v_{33}\end{array}\right)}^{-1}\left(\begin{array}{c}ln\left(x_1\right)-{\mu }_1^i\\ x_3-{\mu }_3^i\hfill \end{array}\right)$	$v_{22}-\left(\begin{array}{cc}{v}_{12}& v_{23}\end{array}\right){\left(\begin{array}{cc}{v}_{11}& v_{13}\\ v_{13}& v_{33}\end{array}\right)}^{-1}\left(\begin{array}{c}{v}_{12}\\ v_{23}\end{array}\right)$
$N_1\left({\Pi }_3^i\left(x_1,x_2\right); {\Lambda }_{312})\right)$	${\mu }_3^i+\left(\begin{array}{cc}{v}_{13}& v_{23}\end{array}\right){\left(\begin{array}{cc}{v}_{11}& v_{12}\\ v_{12}& v_{22}\end{array}\right)}^{-1}\left(\begin{array}{c}ln\left(x_1\right)-{\mu }_1^i\\ ln\left(x_2\right)-{\mu }_2^i\end{array}\right)$	$v_{33}-\left(\begin{array}{cc}{v}_{13}& v_{23}\end{array}\right){\left(\begin{array}{cc}{v}_{11}& v_{12}\\ v_{12}& v_{22}\end{array}\right)}^{-1}\left(\begin{array}{c}{v}_{13}\\ v_{23}\end{array}\right)$

.

The stationary marginal and conditional mean, median, mode, p-quantile (0 < p < 1), and variance trends of the tree diameter and height marginal and conditional processes are listed in

Table 8. Stationary marginal mean, median, mode, p-quantile (0 < p <1), and variance trends.

Variable	Trajectory Type	Equation
Xj j = 1, 2	Mean	$exp\left({\mu }_j^i+\frac{1}{2}{v}_{jj}\right)$
	Median	$exp\left({\mu }_j^i\right)$
	Mode	$exp\left({\mu }_j^i-v_{jj}\right)$
	Quantile (0 < p < 1)	$L{\Phi }_p^{-1}\left({\mu }_j^i;v_{jj}\right)$
	Variance	$exp\left(2{\mu }_j^i+v_{jj}\right)\cdot \left(exp\left(v_{jj}\right)-1\right)$
X3	Mean Median Mode	${\mu }_3^i$
	Quantile (0 < p < 1)	${\Phi }_p^{-1}\left({\mu }_3^i;v_{33}\right)$
	Variance	$v_{33}$

,

Table 9. Stationary conditional mean, median, mode, p-quantile (0 < p < 1), and variance trends with one predictor.

Probability Density	Trajectory Type	Equation
$N_1\left({\eta }_3^i\left(x_j\right); {\lambda }_{3j}\right)$ $j\in \left\{1,2\right\}$	Mean, median, and mode	${\eta }_3^i\left(x_j\right)$
	Quantile (0 < p < 1)	${\Phi }_p^{-1}\left({\eta }_3^i\left(x_j\right); {\lambda }_{3j}\right)$
	Variance	${\lambda }_{3j}$
$L{N}_1\left({\eta }_j^i\left(t,x_k\right); {\lambda }_{jk}\left(t\right)\right)$ $j,k\in \left\{\left(1,2\right),\left(1,3\right),\left(2,1\right),\left(2,3\right)\right\}$	Mean	$exp\left({\eta }_j^i\left(x_k\right)+\frac{1}{2}{\lambda }_{jk}\right)$
	Median	$exp\left({\eta }_j^i\left(x_k\right)\right)$
	Mode	$exp\left({\eta }_j^i\left(x_k\right)-{\lambda }_{jk}\right)$
	Quantile (0 < p < 1)	$L{\Phi }_p^{-1}\left({\eta }_j^i\left(x_k\right); {\lambda }_{jk}\right)$
	Variance	$exp\left(2{\eta }_j^i\left(x_k\right)+{\lambda }_{jk}\right)\cdot \left(exp\left({\lambda }_{jk}\right)-1\right)$

and

Table 10. Stationary conditional mean, median, mode, p-quantile (0 < p < 1), and variance trends with two predictors.

Probability Density	Trajectory Type	Equation
$N_1\left({\Pi }_3^i\left(x_1,x_2\right); {\Lambda }_{312}\right)$	Mean, median, and mode	${\Pi }_3^i\left(x_1,x_2\right)$
	Quantile (0 < p < 1)	${\Phi }_p^{-1}\left({\Pi }_3^i\left(x_1,x_2\right); {\Lambda }_{312}\right)$
	Variance	${\Lambda }_{312}$
$L{N}_1\left({\Pi }_j^i\left(x_k,x_m\right); {\Lambda }_{jkm}\right)$ $j,k,m\in \left\{\left(1,2,3\right),\left(2,1,3\right)\right\}$	Mean	$exp\left({\Pi }_j^i\left(x_k,x_m\right)+\frac{1}{2}{\Lambda }_{jkm}\right)$
	Median	$exp\left({\Pi }_j^i\left(x_k,x_m\right)\right)$
	Mode	$exp\left({\Pi }_j^i\left(x_k,x_m\right)-{\Lambda }_{jkm}\right)$
	Quantile (0 < p < 1)	$L{\Phi }_p^{-1}\left({\Pi }_j^i\left(x_k,x_m\right); {\Lambda }_{jkm}\right)$
	Variance	$exp\left(2{\Pi }_j^i\left(x_k,x_m\right)+{\Lambda }_{jkm}\right)\cdot \left(exp\left({\Lambda }_{jkm}\right)-1\right)$

.

The stationary marginal bivariate distribution of $\left(X_1^i,X_2^i\right)$, is lognormal $L{N}_2\left({\mu }^{i,12}; {\Sigma }^{12}\right)$ with the mean vector ${\mu }^{i,12}$ and the covariance matrix ${\Sigma }^{12}$ given by the following forms:

$${\mu }^{i,12}={\left(\frac{1}{{\beta }_j}\left({\alpha }_j+{\phi }_j^i-\frac{{\sigma }_{jj}}{2}\right), j=1,2\right)}^T$$

(22)

$${\Sigma }^{12}={\left(v_{j,k}\right)}_{j,k=1,2}.$$

(23)

The marginal bivariate distribution of $\left(X_j^i,X_3^i\right)$, j = 1, 2 is hybrid lognormal–normal $L{N}_1{N}_1\left({\mu }^{i,j3}; {\Sigma }^{j3}\right)$ with mean vector ${\mu }^{i,j3}$ and covariance matrix ${\Sigma }^{j3}$ is given by the following forms:

$${\mu }^{i,j3}={\left(\frac{1}{{\beta }_j}\left({\alpha }_j+{\phi }_j^i-\frac{{\sigma }_{jj}}{2}\right),{\alpha }_3+{\phi }_3^i\right)}^T$$

(24)

$${\Sigma }^{j3}=\left(\begin{array}{cc}{v}_{jj}& v_{j3}\\ v_{j3}& v_{33}\end{array}\right).$$

(25)

2.2. Maximum Log-Likelihood Function

The hybrid SDE model defined by Equation (1) can be fitted to the number of tree density (x₁), diameter (x₂), and height (x₃) of samples $\left\{\left(x_{11}^i,x_{21}^i,x_{31}^i\right),\left(x_{12}^i,x_{22}^i,x_{32}^i\right),\dots ,\left(x_{1n_i}^i,x_{2n_i}^i,x_{3n_i}^i\right)\right\}$ at discrete times (ages) $\left\{t_1^i,t_2^i,\dots ,t_{n_i}^i\right\}$ (n_i is the number of observed trees of the ith stand, i = 1, …, M) using the maximum likelihood procedure. The associated maximum likelihood function for the trivariate fixed-effect parameters hybrid SDE model (in this case, the parameters of random effects ${\phi }_1^i$, ${\phi }_2^i$, and ${\phi }_3^i$ are assumed to be equal to the mean value ${\rm E}\left({\phi }_1^i\right)=0$, ${\rm E}\left({\phi }_2^i\right)=0$ and ${\rm E}\left({\phi }_3^i\right)=0$, respectively, takes the following form:

$$L_1\left({\theta }^1\right)={{\displaystyle \prod }}_{i=1}^M{{\displaystyle \prod }}_{j=1}^{n_i}f\left(x_{1j}^i,x_{2j}^i,x_{3j}^i,t_j^i|{\theta }^1,\left(0,0,0\right)\right),$$

(26)

$${\theta }^1=\left\{{\alpha }_1,{\alpha }_2,{\alpha }_3,{\beta }_1,{\beta }_2,{\beta }_3,{\sigma }_{11},{\rho }_{12},{\rho }_{13},{\sigma }_{22},{\rho }_{23},{\sigma }_{33}\right\}$$

and the maximum log-likelihood function as

$$L{L}_1\left({\theta }^1\right)={{\displaystyle \sum }}_{i=1}^M{{\displaystyle \sum }}_{j=1}^{n_i}ln\left(f\left(x_{1j}^i,x_{2j}^i,x_{3j}^i,t_j^i|{\theta }^1,\left(0,0,0\right)\right)\right),$$

(27)

where the probability density function $f\left(x_1,x_2,x_3,t|{\theta }^1,\left(0,0,0\right)\right)$ takes the form defined by Equations (9)–(11).

The maximum likelihood estimator seeks to make the probability density function $f\left(x_1,x_2,x_3,t|{\theta }^1,\left(0,0,0\right)\right)$ the most likely fit the observed dataset $\left\{\left(x_{11}^i,x_{21}^i,x_{31}^i\right),\left(x_{12}^i,x_{22}^i,x_{32}^i\right),\dots ,\left(x_{1n_i}^i,x_{2n_i}^i,x_{3n_i}^i\right)\right\}$ at discrete times (ages) $\left\{t_1^i,t_2^i,\dots ,t_{n_i}^i\right\},$ starting at the initial age t₀ = 0, x₁₀ = 5000, $x_{20}$= 0.1, and $x_{30}$ = 0.1, from which a given set of SDEs (1) evolves.

The associated maximum likelihood function for the trivariate mixed-effect parameters hybrid SDE model takes the following form:

$$L_2\left({\theta }^2,\Psi \right)={{\displaystyle \prod }}_{i=1}^M{{\displaystyle \int }}_{R^3}{{\displaystyle \prod }}_{j=1}^{n_i}f\left(x_{1j}^i,x_{2j}^i,x_{3j}^i,t_j^i|{\theta }^1,{\phi }^i\right){{\displaystyle \prod }}_{k=1}^{3}p\left({\phi }_k^i|{\sigma }_k^2\right)d{\phi }_k^i,$$

(28)

where $p\left({\phi }_k|{\sigma }_k^2\right)$, k = 1, 2, 3 is a normal probability density function with zero mean and constant variance, ${\sigma }_k^2, {\theta }^2=\left\{\delta ,{\alpha }_1,{\beta }_1,{\alpha }_2,{\beta }_2,{\alpha }_3,{\beta }_3,{\sigma }_{11},{\rho }_{12},{\rho }_{13},{\sigma }_{22},{\rho }_{23},{\sigma }_{33},{\sigma }_1,{\sigma }_2,{\sigma }_3\right\}$, and $\Psi =\left({\phi }^1,\dots ,{\phi }^M\right)$. The maximum log-likelihood function is defined as

$$\begin{gathered} L{L}_2\left({\theta }^2,\Psi \right)={{\displaystyle \sum }}_{i=1}^M{{\displaystyle \int }}_{R^3}\left({{\displaystyle \sum }}_{j=1}^{n_i}ln\left(f\left(x_{1j}^i,x_{2j}^i,x_{3j}^i,t_j^i|{\theta }^1,{\phi }^i\right)\right) \\ +{{\displaystyle \sum }}_{k=1}^{3}ln\left(p\left({\phi }_k^i|{\sigma }_k^2\right)\right)\right)d{\phi }_1^{i}d{\phi }_2^{i}d{\phi }_3^i. \end{gathered}$$

(29)

For the mixed-effect parameters hybrid SDE model, the two-step approximated maximum log-likelihood procedure takes the following form:

$$L{L}_2\left({\theta }^2,\stackrel{\wedge }{\Psi }\right)\approx {{\displaystyle \sum }}_{i=1}^M\left(g\left(\stackrel{\wedge }{{\phi }^i}|{\theta }^2\right)+\frac{3}{2}ln\left(2\pi \right)-\frac{1}{2}ln\left(det{\left(\left[-\frac{{\partial }^{2}g\left({\phi }^i|{\theta }^2\right)}{\partial {\left({\phi }^i\right)}^2}\right]\right)}_{{\phi }^i=\stackrel{\wedge }{{\phi }^i}}\right)\right),$$

(30)

where

$$\stackrel{\wedge }{{\phi }^i}=\underset{{\phi }^i}{argmax}g\left({\phi }^i|\stackrel{\wedge }{{\theta }^2}\right), i=1, \dots , M,$$

(31)

$$\begin{gathered} g\left({\phi }^i|\stackrel{\wedge }{{\theta }^2}\right)\equiv {{\displaystyle \sum }}_{j=1}^{n_i}ln\left(f\left(x_{1j}^i,x_{2j}^i,x_{3j}^i,t_j^i|\stackrel{\wedge }{{\theta }^1},{\phi }^i\right)\right)+ln\left(p({\phi }_1^i|{\stackrel{\wedge }{\sigma }}_1\right)+ln\left(p\left({\phi }_2^i|{\stackrel{\wedge }{\sigma }}_2\right)\right)+ \\ ln\left(p({\phi }_3^i|{\stackrel{\wedge }{\sigma }}_3\right). \end{gathered}$$

(32)

The maximization of $L_2\left({\theta }^2,\Psi \right)$ is a two-step optimization problem. The internal optimization step estimates the ${\phi }^i$ for every plot i = 1, …, M with Equation (31). The external optimization step maximizes $L{L}_2\left({\theta }^2,\stackrel{\wedge }{\Psi }\right)$ after plugging $\stackrel{\wedge }{{\phi }^i}$, i = 1, …, M into Equation (30). Parameter estimates using the maximum likelihood procedure are costly processes in terms of time and computational power.

The random effects $\phi =\left({\phi }_1,{\phi }_2,{\phi }_3\right)$ calibrated from the newly observed dataset $\left\{\left(x_{1,1},x_{2,1},x_{3,1}\right),\left(x_{1,2},x_{2,2},x_{3,2}\right),\dots ,\left(x_{1,m},x_{2,m},x_{3,m}\right)\right\}$ at discrete previous times (ages) $\left\{t_1,t_2,\dots ,t_m\right\}$ can be calibrated using the following form:

$$\begin{gathered} \stackrel{\wedge }{\phi }=\underset{\left({\phi }_1,{\phi }_2,{\phi }_3\right)}{argmax}\left({{\displaystyle \sum }}_{j=1}^{m}ln\left(f\left(x_{1,j},x_{2,j},x_{3,j},t_j|\stackrel{\wedge }{{\theta }^1},\phi \right)\right)+ln\left(p({\phi }_1|{\stackrel{\wedge }{\sigma }}_1\right)+ \\ ln\left(p\left({\phi }_2|{\stackrel{\wedge }{\sigma }}_2\right)\right)+ln\left(p({\phi }_3|{\stackrel{\wedge }{\sigma }}_3\right)\right). \end{gathered}$$

(33)

2.3. Data

Mountain pine (Pinus mugo Turra) is dwarf, slow growing plant grown on the coastal dunes of the Curonian Spit (Kuršių Nerija) in Lithuania. The first plantations of Mountain pine were established on the coastal sand dunes nearly 200 years ago [26]. Silvicultural policy in the Curonian Spit, a UNESCO World Heritage Site, started in 2013 when mountain pine trees were cut down in some places to open up the spit’s sandy hills by its gradual replacement with Pine sylvestris. A field study was accomplished in western Lithuania (Kuršių Nerija). Age of sampled stands varies from 53 till 123 years. The radius of the 31 circular plots was 6.9 m. Size of round sample plot was 150 m². In all plots, the diameter at breast height of all trees larger than 1 mm was measured (7005 trees). Diameter was measured to an accuracy of 1 mm. The observed dataset used to develop the model consisted of 702 Mountain pine tree height and diameter pairs of measurements taken from 31 plots, with a wide range of stand ages. A tree was measured once its height exceeded 1.3 m. For every sample tree, the tree diameter over bark at 1.30 m (in centimeters), the tree height (in meters), and the age (in years) were recorded. Height was measured by using clinometer, with precision to the nearest 0.1 meter. Tree age was identified with stand age which provides a general timeframe for when stands first established. The available observed dataset was randomly divided into model estimation (23 plots) and model validation (8 plots) parts. Both estimation and validation datasets are presented in

Table 11. Characterization of the datasets used to fit the SDE (1).

Data	Number of Trees (Plots)	Min	Max	Mean	St. Dev.	Number of Plots	Min	Max	Mean	St. Dev.
Estimation						Validation
N	(23)	2875	25,139	15,244	6345	(8)	6418	19,255	14,717	4371
d (cm)	527 (23)	1.20	13.40	3.96	1.61	175 (8)	1.80	8.00	4.15	1.25
h (m)	527 (23)	1.50	7.90	3.55	1.11	175 (8)	2.07	6.56	3.84	0.91
t (year)	527 (23)	53.00	123.00	92.21	20.82	175 (8)	53.00	103.00	80.60	17.50

.

3. Results

3.1. Parameter Estimators

Parameter estimations of trivariate hybrid SDE (1) use samples where the diameter, height, and age are measured on every tree and tree density per hectare on every stand. The maximum likelihood estimation technique capacities for unbiased minimum variance estimators having approximate normal distributions and approximate sample variances given by the Fisher [27] information matrix will most likely favor the other techniques where the emphasis is placed on predicting particular observations in a given situation. For the fixed- and mixed-effect scenario models, the parameter estimators ${\theta }^1=\left\{{\alpha }_1,{\alpha }_2,{\alpha }_3,{\beta }_1,{\beta }_2,{\beta }_3,{\sigma }_{11},{\rho }_{12},{\rho }_{13},{\sigma }_{22},{\rho }_{23},{\sigma }_{33}\right\}$ and ${\theta }^2=\left\{{\alpha }_1,{\beta }_1,{\alpha }_2,{\beta }_2,{\alpha }_3,{\beta }_3,{\sigma }_{11},{\rho }_{12},{\rho }_{13},{\sigma }_{22},{\rho }_{23},{\sigma }_{33},{\sigma }_1,{\sigma }_2,{\sigma }_3\right\}$, respectively, were calculated by the maximization of the log-likelihood function defined by Equation (27) (fixed effect scenario) and by the two-step maximization technique defined by Equations (30) and (31) (mixed effect scenario). Maximization was performed using the NLPSolve procedure in MAPLE software [28]. The initial condition takes the form if t = 0, then $X^i\left(0\right)=x_0=\left(x_{10},x_{20},x_{30}\right)={\left(5000,0.1,0.1\right)}^T$. The parameter estimators are summarized in

Table 12. Estimates of fixed-effect parameters.

Scenario	Parameters of Drift Term
	α1	β1	α2	β2	α3	β3
Fixed	0.1884	0.0186	0.0833	0.0558	3.5455	0.1548
Mixed	0.5949	0.0598	0.7442	0.4851	3.7622	0.3152
	Parameters of Diffusion Term
	σ11	ρ12	ρ13	σ22	ρ23	σ33	σ1	σ2	σ3
Fixed	0.0047	–0.6443	–0.9999	0.0165	0.8814	0.3845	-	-	-
Mixed	1.9 × 10−5	–0.9999	–0.8711	0.1313	0.7965	0.3284	0.0398	0.1312	1.1268

.

3.2. Bivariate Distributions

The trivariate dynamics of the tree density, diameter, and height denominated by the trivariate probability density function (see Equations (9) and (10) or, for stationary case, Equations (20) and (21)) yield a unified system of forest stand development. The focus of this paper is the methodology of growth and yield modeling using a hybrid trivariate SDE. The derived trivariate hybrid probability density function can be used for calculation of stand volume, marginal bivariate can be used for calculation of stand basal area, and so on. The conditional distributions of the size variable can be used for formalizing a wide spectrum of tree or stand variable relationships. The manner in which tree size variable distributions vary with the others provides information about how these size variables are related, and these distributions can be described in part, as for any univariate distributions, by their expectations and variances. This modeling framework does not introduce a new modern concept, only suggests a new use for the previous trivariate hybrid diffusion process concept.

To visualize the differences between the fixed- and mixed-effects modeling scenarios, we provide standard plots of the bivariate density functions and tolerance regions. We focus on bivariate distributions because they are easier to visualize than higher dimensional distributions. The same ideas could be carried over to three dimensions; however, the visualization would require considerable effort. The estimated bivariate densities use the estimators of the fixed effect parameters $\stackrel{\wedge }{{\theta }^1}$ and $\stackrel{\wedge }{{\theta }^2}$ listed Table 1, and the random effects ${\phi }_1,{\phi }_2, \mathrm{and} {\phi }_3$ calibrated by Equation (33). Figure 1 shows the estimated bivariate probability density functions (tree density and diameter, tree density and height, and diameter and height), as well as fixed- and mixed-effect scenarios for the same randomly selected stand. The fixed-effect scenario densities are considerably flatter than the mixed-effect scenario, which reflects the spatial hierarchy of the observed dataset. The diameter and height bivariate density function appears steeper than the other two distributions.

To illustrate that the observed dataset corresponding to the randomly selected stand from the validation dataset of the Mountain pine trees sufficiently corresponds to the marginal bivariate estimated probability density function, we used a simple graphical technique called a tolerance region. A tolerance region captures a specified proportion or more of size variables, with a given tolerance level β. Conversely, a confidence region provides a method of estimating parameter vectors using a corresponding sample statistic to a given level of confidence γ. The tolerance region is the region pertaining to the entire stand and not just to a specific vector of parameters. It is expected that 100β% of the entire population will lie in the tolerance region. This procedure has two probability values involved: (1) the coverage probability β is the fixed percentage of the data to be covered; (2) the confidence level γ. Hence, with γ% confidence, at least β% of the data fall within the given region. For example, if γ = 0.95 and β = 0.90, we have a 95% confidence region for 90% coverage. Basically, the bivariate normal tolerance region plot can be used for the cases with normally distributed data. For lognormally-distributed data, we can plot a tolerance region using a logarithmic axis. A tolerance region is a region for vector x satisfying the inequality:

$${\left(x-{\mu }^{i,12}\left(t\right)\right)}^T{\left[{\Sigma }^{12}\left(t\right)\right]}^{-1}\left(x-{\mu }^{i,12}\left(t\right)\right)=K$$

(34)

where K is the tolerance factor that is subject to the condition that the region in Equation (34) contains at least β% of the normal distribution with confidence γ% [29]. Figure 2 illustrates the tolerance regions with fixed-effect and mixed-effect scenarios for β = 0.90, 0.95, and confidence level γ = 0.95, which correspond to the randomly selected observed validation dataset. The random effects were calibrated by Equation (33). Figure 2 shows that the tolerance regions of the mixed-effect parameters bivariate probability density functions better centered the observed data points for the randomly selected stand than the tolerance regions of the fixed-effect parameters scenario. For the tolerance region plots, the K values were chosen from Table 1 in [29]: the first setting β = 0.90 and γ = 0.95 produces a K value of 7.74; the second setting β = 0.95 and γ = 0.95 produces a K of 10.02.

Several studies have suggested maintaining a continuous-cover forest stand and gathering particular benefits [30]. A special case of continuous-cover forest stand process could be formulated as a stationary system of SDEs that ensures the equilibrium state of tree size variables. The newly developed stationary bivariate distributions are derived in Equations (22)–(25). For comparison of bivariate distributions in Figure 1 and Figure 2 with their stationary analogues, we present the graphics of stationary distributions in Figure 3 and Figure 4.

4. Discussion

4.1. Models of Tree Density, Diameter, and Height Mean Dynamics

A system of trivariate hybrid SDE models for forest growth and yield estimation was presented in Section 2.1. The developed SDE system is novel in its formation using the different types of diffusion processes and the representatively-observed Mountain pine plot data for modeling tree and stand growth. Since planning and forest management involve forecasting yield over short and long time periods, the diffusion processes relating tree (stand) age with size components produce the most accurate modeling framework. The trivariate hybrid SDE easily connects univariate marginal and conditional probability density functions of tree density, diameter, and height with their multivariate probability density functions.

The marginal probability densities, lognormal $L{N}_1\left({\mu }_j^i\left(t\right); v_{jj}\left(t\right)\right)$, j = 1, 2, and normal $N_1\left({\mu }_3^i\left(t\right); v_{33}\left(t\right)\right)$, i = 1, …, M with mean and variance defined by Equations (7) and (8), and the conditional probability densities lognormal $L{N}_1\left({\eta }_1^i\left(t,x_2\right); {\lambda }_{12}\left(t\right)\right)$, $L{N}_1\left({\eta }_1^i\left(t,x_3\right); {\lambda }_{13}\left(t\right)\right)$, $L{N}_1\left({\eta }_2^i\left(t,x_1\right); {\lambda }_{21}\left(t\right)\right)$, $L{N}_1\left({\eta }_2^i\left(t,x_3\right); {\lambda }_{23}\left(t\right)\right)$, $L{N}_1\left({\Pi }_1^i\left(t,x_2,x_3\right); {\Lambda }_{123}\left(t\right)\right)$, $L{N}_1\left({\Pi }_2^i\left(t,x_1,x_3\right); {\Lambda }_{213}\left(t\right))\right)$, and normal $N_1\left({\eta }_3^i\left(t,x_2\right); {\lambda }_{32}\left(t\right)\right)$, $N_1\left({\eta }_3^i\left(t,x_1\right); {\lambda }_{31}\left(t\right)\right),$, $N_1\left({\Pi }_3^i\left(t,x_1,x_2\right); {\Lambda }_{312}\left(t\right)\right)$ with the mean and variance listed in Table 2 and Table 3, allowed us to predict the mean, median, mode, p-quantile (0 < p < 1), and variance trajectories for all response variables. The mean of the tree density, diameter, and height dynamics in the forestry literature have been formulated using a wide range of mathematical relationships, from linearized fixed-effect parameters regression equations to generalized nonlinear mixed-effect parameters relationships [29,30,31,32]. The linkage between the diameter, height, and the other stand variables such as tree density has many applications in forest inventory and remote sensing for stand basal area and volume calculation [33].

In this paper, after the fitting the SDE (1) to the estimated dataset of Mountain pine trees, the four statistical measures of goodness of fit were calculated and are listed in

Table 13. Statistical measures and p-values of the Student’s t-test for the fixed-effect scenario models.

(Equation): (Predictors)	Estimation Dataset (Prediction)					Validation Dataset (Forecast)
	B (%)	AB (%)	RMSE (%)	R2	T p-Value	B (%)	AB (%)	RMSE (%)	R2	T p-Value
Tree density
(4): (t)	−48.4102 (−15.42)	4332.14 (33.97)	4947.44 (69.67)	0.2564	0.8332	−724.47 (−12.35)	1985.58 (19.43)	3270.42 (46.35)	0.1629	0.0040
(10): (d,h,t)	−80.9793 (−8.12)	3011.79 (21.94)	3751.34 (41.67)	0.5725	0.6204	333.50 (−2.07)	2670.50 (20.51)	3253.55 (30.20)	0.1716	0.1768
(13): (d,t)	−98.8060 (−9.62)	3350.29 (24.37)	4142.50 (47.37)	0.4787	0.5842	123.34 (−4.11)	2735.68 (21.35)	3425.94 (34.64)	0.0814	0.6345
(13): (h,t)	−72.9577 (−8.35)	3014.94 (22.25)	3777.83 (42.36)	0.5664	0.6577	274.36 (−2.56)	2614.51 (20.31)	3248.79 (30.36)	0.1738	0.2655
Diameter
(4): (t)	−0.0088 (−15.54)	1.1923 (34.02)	1.6504 (47.31)	0.0000	0.9025	0.2702 (−2.27)	1.0449 (25.57)	1.3191 (31.87)	0.0000	0.0074
(10): (N,h,t)	−0.0156 (−5.85)	0.7091 (19.46)	0.9528 (26.54)	0.6482	0.7065	0.0662 (−2.98)	0.7928 (19.83)	1.0176 (25.51)	0.3427	0.3906
(13): (N,t)	−0.0271 (−10.60)	0.9516 (26.56)	1.3637 (38.20)	0.2733	0.6501	0.2704 (0.27)	0.9947 (24.07)	1.2880 (30.28)	0.0000	0.0061
(13): (h,t)	−0.0128 (−6.00)	0.7254 (19.81)	0.9705 (26.85)	0.6350	0.7627	0.0300 (−3.89)	0.7931 (20.21)	1.0210 (26.50)	0.3383	0.6977
Height
(4): (t)	0.0047 (−8.54)	0.8118 (23.69)	1.1143 (31.48)	0.0000	0.9229	0.2944 (2.85)	0.6847 (16.66)	0.9155 (21.42)	0.0000	0.00003
(10): (N,d,t)	−0.0014 (−2.76)	0.4862 (14.21)	0.6448 (18.90)	0.6652	0.9605	0.1127 (0.46)	0.5192 (13.13)	0.6964 (16.74)	0.4212	0.0336
(13): (N,t)	−0.0047 (−5.42)	0.6118 (18.10)	0.8411 (25.16)	0.4302	0.8979	0.2745 (3.20)	0.6898 (17.02)	0.9009 (21.42)	0.0313	0.00008
(13): (d,t)	0.0015 (−3.00)	0.5338 (15.56)	0.6993 (20.33)	0.6061	0.9604	0.0693 (−0.50)	0.5541 (14.20)	0.7031 (17.15)	0.4101	0.1937

The mean prediction bias, $B=\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n\left(y_i-\stackrel{\wedge }{y_i}\right)$; the percentage mean prediction bias, $\%B=\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n\frac{y_i-\stackrel{\wedge }{y_i}}{y_i}\ast 100$; the absolute mean prediction bias, $AB=\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n\left|y_i-\stackrel{\wedge }{y_i}\right|$; the percentage mean absolute prediction bias, $\%B=\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n\left|\frac{y_i-\stackrel{\wedge }{y_i}}{y_i}\right|$; the root mean square error, $RMSE=\sqrt{\frac{1}{n-1}{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\stackrel{\wedge }{y_i}\right)}^2}$; the percentage root mean square error, $\%RMSE=\sqrt{\frac{1}{n-1}{{\displaystyle \sum }}_{i=1}^n{\left(\frac{y_i-\stackrel{\wedge }{y_i}}{y_i}\right)}^2}\ast 100$; coefficient of determination, $R^2=1-\frac{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\stackrel{\wedge }{y_i}\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y}\right)}^2},$ where $n={{\displaystyle \sum }}_{i=1}^M{n}_i$ is the total number of observations used to fit the model, M is the number of stands, n_i is the number of measured trees in the ith stand, and $y_i$, $\stackrel{\wedge }{y_i}$, and $\overline{y}$ are the measured, estimated, and average values of the dependent variable (number of trees per hectare, N; diameter, d; height, h; stand basal area, G; stand volume, V_S).

and

Table 14. Statistical measures and p-values of the Student’s t-test for the mixed-effect scenario models.

See Table: (Predictors)	Estimation Dataset					Validation Dataset
	B (%)	AB (%)	RMSE (%)	R2	T p-Value	B (%)	AB (%)	RMSE (%)	R2	T p-Value
Tree density
(1): (t)	1.6709 (−0.01)	6.2768 (0.05)	7.8379 (0.09)	1.0	0.0000	0.2336 (−0.01)	5.4773 (0.04)	6.9829 (0.7)	1.0	0.6587
(5): (d,h,t)	−1.0965 (−0.06)	80.4988 (0.46)	108.7885 (0.60)	0.9996	0.8171	−1.172 (−0.01)	68.4800 (0.43)	88.6621 (0.54)	0.9994	0.8614
(3): (d,t)	−1.0340 (−0.01)	81.9300 (0.46)	111.9989 (0.63)	0.9996	0.8322	−1.260 (−0.01)	71.8156 (0.45)	89.7277 (0.55)	0.9994	0.8528
(3): (h,t)	0.6278 (−0.01)	76.4693 (0.45)	102.8041 (0.60)	0.9997	0.8886	−0.877 (−0.01)	68.7991 (0.43)	91.9456 (0.57)	0.9994	0.8997
Diameter
(1): (t)	–0.0823 (−11.44)	0.8511 (21.83)	1.1336 (32.29)	0.5019	0.0963	−0.1291 (−10.72)	0.8903 (23.79)	1.0774 (31.52)	0.2631	0.1149
(5): (N,h,t)	–0.0070 (−5.05)	0.6436 (17.76)	0.8485 (24.63)	0.7210	0.8501	−0.0202 (−5.26)	0.7311 (18.86)	0.9636 (25.63)	0.4106	0.7821
(3): (N,t)	0.0067 (−8.50)	0.8299 (23.61)	1.1200 (33.33)	0.5138	0.8901	−0.0255 (−7.84)	0.8742 (22.96)	1.0721 (29.81)	0.2704	0.7530
(3): (h,t)	–0.0454 (−5.33)	0.6726 (18.19)	0.8871 (25.32)	0.6950	0.2402	−0.0552 (−5.68)	0.7577 (19.57)	1.0122 (27.08)	0.3497	0.4718
Height
(1): (t)	0.0012 (−2.93)	0.4065 (12.53)	0.5483 (18.03)	0.7578	0.9580	0.0003 (−2.09)	0.3971 (10.96)	0.5182 (15.30)	0.6796	0.9937
(5): (N,d,t)	0.0008 (−1.47)	0.3224 (9.68)	0.4214 (13.07)	0.8570	0.9639	−0.00002 (−1.25)	0.3332 (8.99)	0.4495 (12.07)	0.7588	0.9994
(3): (N,t)	0.0007 (−2.77)	0.4050 (12.47)	0.5451 (17.85)	0.7607	0.9751	−0.0005 (−2.02)	0.3939 (10.87)	0.5165 (15.16)	0.6816	0.9900
(3): (d,t)	0.0010 (−1.19)	0.3306 (9.94)	0.4298 (13.38)	0.8512	0.9584	0.0003 (−1.07)	0.3468 (9.30)	0.4662 (12.27)	0.7406	0.9937

to compare the models concerning the dynamics of mean tree density, mean diameter, and mean height. Additionally, Table 13 and Table 14 provide the p-value of the Student’s t-test [34], which determines whether the mean value of model predictions is statistically different from zero. The used parameter estimates were selected from Table 12, and the random effects for the observed validation dataset were calibrated by Equation (33).

Within comparison problems, the aim of predictive SDE models of the tree density, diameter, and height is to accurately predict an outcome value from a set of predictors (age, tree density, diameter, and height). Common assumptions of these models are nonlinearity; that is, the expected outcome value is modeled by a nonlinear combination of predictors, and the underlying covariance structure is considered to drive changes in the predictor variables (tree density, diameter, and height). The comparison of the models using the four statistical measures, presented in Table 13 and Table 14, highlighted that the mixed-effect models are superior to the fixed-effect models. Consequently, the mixed effect scenario models account unobserved ecological factors such as temperature, solar radiation, precipitation, nutrient condition and much more. The statistical measures computed for the estimation dataset indicated a moderate accuracy level, with lower performance on the validation dataset, as expected. Therefore, as demonstrated in Table 13 and Table 14, the maximum impact on the tree density dynamics demonstrated the tree height. For the impact on diameter dynamics, the height was shown to be the most important predictor variable; for the impact on height dynamics, diameter was the most important predictor variable; for the impact on tree density dynamics, height was the most important predictor variable.

Figure 5 shows the tree density, diameter, and height mean and 5% and 95% quantile trajectories via age for three randomly selected stands from the validation dataset using marginal univariate probability densities. The figure reveals the superiority of the mixed-effects scenario over the fixed-effects scenario. The mixed-effect 5% and 95% quantile curves of the tree density, diameter, and height for all three stands in the validation dataset (Figure 5) closely matched the observed data points in the validation dataset. The fixed-effect scenario modeling technique produced quite a few of the points from the observed datasets that passed the 5% and 95% quantile limits for tree height.

4.2. Models of Stand Basal Area

The stand basal area is a measure of tree density, which is the sum of the basal area of all (living) trees in a stand, expressed in m²/ha, and widely used in forestry to describe the average amount of an area occupied by trees. Stand basal area is a useful measure for understanding forest–wildlife habitat relationships and making timber harvest decisions, and can be used to estimate stand volume per hectare or as a useful measure of the competition in a stand [35]. Basal area production describes the stand development over age and is defined as [36]:

$$G=\frac{\pi \cdot D^2}{40000}\cdot N.$$

(35)

The multitude of existing machine learning models, when applied to stand basal area and volume, can be seen as regression models, and they need the underlying growth processes to be stationary, which means that we assume the mean and the variance of stand size components are constant over time. This might seem obvious to some, but this is seldom mentioned in the scientific literature. All the above regression models can be generalized to handle this non-stationarity by applying SDEs to the dynamics of stand size components such as the stand basal area and volume per hectare. According to the marginal bivariate (tree density and diameter) probability density function defined by Equations (12) and (13), the basal area dynamics takes the following form:

$$\begin{gathered} G\left(t\right)=\frac{\pi }{40000}{{\displaystyle \int }}_{-\infty }^{+\infty }{{\displaystyle \int }}_{-\infty }^{+\infty }{x}_2{}^2\cdot x_1\cdot f\left(x_1,x_2|\stackrel{\wedge }{\Theta }\right)d{x}_{1}d{x}_2, \stackrel{\wedge }{\Theta }=\left(\stackrel{\wedge }{{\theta }^2},\stackrel{\wedge }{\Psi }\right), \\ {\theta }^2=\left\{{\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,{\sigma }_{11},{\rho }_{12},{\sigma }_{22},{\sigma }_1,{\sigma }_2,{\sigma }_3\right\}. \end{gathered}$$

(36)

The aim of developing SDE stand basal area models was to estimate the present and future values of stand basal area at the whole Mountain pine stand level in Lithuania. Figure 6 illustrates the basal area development over time with the observed dataset for the mixed-effect scenario (the random effects were calibrated by Equation (33)). Figure 6 demonstrates that the mixed-effect scenario basal area formula, defined by Equation (36), for three stands provides predicted values close to the observed stand basal area values. The results showed that there were marked differences between the mixed-effect parameters and the fixed-effect parameters simulations of the stand basal area; moreover, the fixed-effect scenario model alone could not explain the full range variation between stands. The fixed-effect scenario stand basal area models should instead be considered as some of the possible dynamics that are not necessarily close to the true trajectory of the stand basal area.

Table 15. Statistical measures and p-values of the Student’s t-test for the stand basal area (m²/ha) (Equation (36)).

Model	B (%)	AB (%)	RMSE (%)	R2	T
Mixed	0.1587 (1.10)	1.1292 (5.16)	1.6082 (8.81)	0.9451	0.5909
Mixed-stationary	0.1587 (1.10)	1.1292 (5.16)	1.6082 (8.81)	0.9451	0.5909

presents the statistical measures for the stand basal area mixed-effect parameters non-stationary and stationary SDE models. Table 15 shows that the mixed-effect scenario SDE model achieved high accuracy in predicting stand basal area. The non-stationary and stationary stand basal area models produced similar goodness results from the numerical simulations. We suggest this conjunction is largely due to the parameter estimates of the SDE (1), which focused on Mountain pine trees growth experiments with the majority of data points from the stationary phase (more than 50 years).

In this study, we used the diameter, height, and tree density as SDE system variables, but the other variables at the stand level (stand basal area and volume per hectare) were analyzed by prioritizing their definition.

4.3. Models of Stand Volume

Typically, stand level models evaluate stand volume as a function of stand level variables such as age, site index, and density. Our developed methodology of stand volume modeling is based on a general trivariate (diameter, height, and density) distribution function that challenges conventional thinking in forest stand modeling. The European national forest inventories traditionally use regression models to estimate tree volume from measurements of tree variables, such as the diameter at breast height and height, which are measured during inventory.

Prior to the stand volume analysis, individual tree volume prediction was examined using a few general nonlinear models using power and q-exponential regression forms [37,38]. We incorporated the q-exponential function, which is defined as

$$V={\beta }_1+{\beta }_2{h}^{{\beta }_3}{\left[-{\beta }_4\left(1-exp(\left(1-{\beta }_5\right)d\right))\right]}_{+}^{\frac{1}{1-{\beta }_5}}$$

(37)

where d is diameter at breast height, h is tree height, β₁–β₅ are the unknown parameters to be estimated, and ${[a]}_{+}=\{\begin{array}{c}a, if a\ge 0,\\ 0,if a<0\end{array}$. The parameters were estimated using observed dataset (217 trees): β₁ = −0.0004; β₂ = 0.0022; β₃ = 0.8422; β₄ = 0.2344; β₅ = −0.3139 [12].

According to the trivariate (tree density, diameter, and height) probability density function defined by Equations (7)–(10), the volume evolution takes the following form:

$$\begin{gathered} V_S\left(t\right)=\frac{\pi }{40000}\cdot {{\displaystyle \int }}_{-\infty }^{+\infty }{{\displaystyle \int }}_{-\infty }^{+\infty }{{\displaystyle \int }}_{-\infty }^{+\infty }V\left(y,z\right)\cdot z\cdot f\left(x,y,z,t|\stackrel{\wedge }{\Theta }\right)\cdot dx\cdot dy\cdot dz, \stackrel{\wedge }{\Theta }\in \left(\stackrel{\wedge }{{\theta }^1},\stackrel{\wedge }{\Psi }\right), \\ {\theta }^1=\left\{{\alpha }_1,{\beta }_1,{\alpha }_2,{\beta }_2,{\alpha }_3,{\beta }_3,{\sigma }_{11},{\rho }_{12},{\rho }_{13},{\sigma }_{22},{\rho }_{23},{\sigma }_{33}\right\}. \end{gathered}$$

(38)

Figure 7 shows the dynamics of the stand volume per hectare as a function of stand age using the mixed-effect scenario.

Table 16. Statistical measures and p-values of the Student’s t-test for the stand volume (m³/ha) (Equation (40)).

Model	B (%)	AB (%)	RMSE (%)	R2	T
Mixed	0.7361 (1.79)	4.8464 (6.40)	7.8747 (9.80)	0.9979	0.5204
Mixed-stationary	0.7361 (1.79)	4.8464 (6.40)	7.8747 (9.80)	0.9979	0.5204

shows the predictive ability for both the newly-developed non-stationary and stationary mixed-effect scenario volume per hectare models defined by Equation (38).

5. Conclusions

In this study, we developed new ideas about how the forest stands of Mountain pine trees evolve from natural or artificial establishment toward old growth forest stands. The old growth forest stand dynamics are unclear, so we proposed a hybrid multivariate diffusion process defined by the multivariate SDE to define changing frames between tree- or stand-state variables over time. First, the marginals we used of the hybrid multivariate SDE in the Gompertz and Vasicek types include the asymmetric and symmetric growth patterns in the biological world. Second, our derived stationary hybrid multivariate probability density function enables the clearer establishment of what constitutes an old growth forest stand (Table 14 and Table 15). The results demonstrated the high accuracy of the newly developed models for stand basal area and volume predictions. Many stand or tree attributes and the different tree species may be examined using stochastic process analogy.