Hyperbolastic Models from a Stochastic Differential Equation Point of View

Abstract

A joint and unified vision of stochastic diffusion models associated with the family of hyperbolastic curves is presented. The motivation behind this approach stems from the fact that all hyperbolastic curves verify a linear differential equation of the Malthusian type. By virtue of this, and by adding a multiplicative noise to said ordinary differential equation, a diffusion process may be associated with each curve whose mean function is said curve. The inference in the resulting processes is presented jointly, as well as the strategies developed to obtain the initial solutions necessary for the numerical resolution of the system of equations resulting from the application of the maximum likelihood method. The common perspective presented is especially useful for the implementation of the necessary procedures for fitting the models to real data. Some examples based on simulated data support the suitability of the development described in the present paper.

1. Introduction

In recent years, widespread interest has arisen concerning the description and analysis of dynamic systems governed by growth curves, owing to the need to explain the evolution of a variety of phenomena in the most precise way possible. Some examples include the study of cell growth (stem cells, tumors, and tissue regeneration), the diffusion of innovations (social networks and the implementation of new technologies), or the exploitation of energy resources (production peaks and length of exploitation periods). In this regard, the latest rage is probably the study of the evolution of epidemics such as COVID-19, in which concepts such as contagion peaks, inflection points in the spread of disease, duration, and onset times of successive waves are of paramount importance for the purposes of planning and policy making. In all cases, the objective is to explain the behavior of the phenomenon and allow for its prediction while keeping in mind that certain external influences may affect its evolution.

The study of growth curves has taken two distinct directions: the unification and the generalization of models (Chakraborty et al. [1]). Unification refers to the compact representation of different models through a single function, which allows for improved mathematical manageability. As for generalization, it refers to the process that starts from a simple equation to which parameters or functions are subsequently incorporated. Its goal is to generate more flexible forms and thus increase applicability to a wider range of research areas. Some outstanding examples of unification are the works of Martínez et al. [2] and Tjørve [3] which propose a variety of approaches for a growth equation that unifies many of the existing sigmoidal growth curves. As for generalization, we can cite the works of Tsoularis and Wallace [4] and that of Koya and Goshu [5]. Although generalized equations offer greater flexibility in form, the resulting models have a harder time of fitting real data sets on account of their complexity. For this reason, the idea of including new parameters has been recently combined with that of introducing flexible functions. In this sense, Tabatabai et al. [6] have constructed the so-called hyperbolic curves and shown their usefulness in the study of the evolution of tumor processes and stem cell growth (Tabatabai et al. [7]). Furthermore, and more precisely, such models have been used to describe the growth of the solid Ehrlich carcinoma under combined treatments (see in [8]). To this end, all three hyperbolastic models are used and compared in different scenarios. Related to this, the type II is employed to model the dynamics of cervical cancer mortality rates [9]. See also in [10], where hyperbolastic curves are used to introduce a model of wound healing. All of these models are the starting point for the development of others, such as the oscillabolastic, which aim to model oscillatory growth (Eby and Tabatabai [11]). Similarly, T-type models (Tabatabai et al. [12]) are capable of representing biphasic sigmoidal growths characterized by stages of growth and decrease, that is, multisigmoidal patterns of behavior. Along these lines, Erto and Lepore [13] have recently defined a new type of sigmoidal curve which, given that certain parametric conditions are met, presents more than one inflection point.

As the variability of the phenomena under study cannot be controlled in deterministic models associated with growth curves, randomness is included in the ordinary differential equations from which they originate. To this end, we employ stochastic differential equations whose solutions are, under certain conditions, diffusion processes. The most widely studied curves in this context are, without a doubt, the logistic and Gompertz curves. For example, Schurz [14] considered a fairly general version of the stochastic differential equation associated with logistic growth, and Rupšys [15] introduced several versions of the differential equation that included delays, whereas Schlomann [16] studied logistic diffusion processes in the presence of catastrophes. This field is unceasingly producing new research, as exemplified by the recent works of Rajasekar et al. [17,18] in which the authors analyzed a stochastic version of the SIR models for the spread of the COVID-19 pandemic. On the other hand, the fact that the Gompertz curve is an excellent model for the description of tumor growth has motivated the introduction of several diffusion processes associated with it (see, for example, Lo [19] and Ferrante [20]). On the basis of these processes, various modifications have subsequently been made with the goal of describing the evolution of tumor growth in the presence of therapeutic treatments (Albano et al. [21,22,23]).

With regard to diffusion processes associated with other growth curves, the existing references are scarce. Perhaps the main difficulty lies in the fact that, although the stochastic differential equation may have a solution (and sometimes not in a totally explicit way), the diffusion process obtained therefrom is not totally manageable. This is a problem even in the logistic case. Skiadas [24] presented an exact solution for the logistic stochastic differential equation with multiplicative noise, although its shape showed that some characteristics cannot possibly be obtained in an explicit fashion, as is the case of transition probability density functions. Therefore, when studying this process one must resort to the discretization of the differential equation. However, convenient reformulations of some growth curves have succeeded in verifying a differential equation of the Malthusian type, which is the way to construct diffusion processes based on modifications of the inhomogeneous lognormal process. Such modeling allows for diffusion processes with transition densities whose functional form is perfectly known. Such is the case of the Bertalanffy (Román-Román et al. [25]), Richards (Román-Román and Torres-Ruiz [26]), and Hubbert curves (Luz Sant’Ana et al. [27]). An alternative way is the modification of ordinary differential equations verified by certain growth curves, which can be done by including in them certain temporal functions. Thus, from the logistic and Weibull equations, Barrera et al. [28,29] introduced the hyperbolastic diffusion processes of type I and III. Likewise, and following this same idea, a multisigmoidal stochastic version of the logistic process has been obtained (Di Crescenzo et al. [30]) as well as a diffusion processes linked to the oscillabolastic curve (Barrera et al. [31]).

The main aim of this paper is to present a common model for diffusion processes associated with the family of hyperbolastic curves. Our proposal not only allows for the results obtained for the type I and type III models to be unified, but also extends them to the hyperbolastic curve of type II. The present paper is structured as follows. Section 2 makes a brief summary of the family of hyperbolastic curves. Additionally, and with the definition given by Tabatabai et al. in [6] as a starting point, it is verified how each of the curves satisfies an ordinary differential equation of the Malthusian type, which extends the deterministic models to stochastic ones that have the non-homogeneous lognormal diffusion process as common point, this being the objective of Section 3. Furthermore, and given that these models are appropriate to describe real phenomena governed by hyperbolastic curves, the inference in them must be developed. Section 4 describes a unified version for the three models presented. Of particular interest is the description of how initial solutions are obtained for the resulting system of maximum likelihood equations. Note that this joint vision enables an efficient computational implementation of the procedures derived from the inferential development herein described. Finally, Section 5 presents some simulation-based examples.

2. An Overview of Hyperbolastic Growth Curves

For the sake of clarity, this section will provide a brief summary of the hyperbolastic curves to date. One of the most interesting aspects of the following description is the fact that all three hyperbolastic models can be seen from a unified perspective.

The hyperbolastic model of type I ($H_1$) appears as an extension of the classic logistic model in which a certain hyperbolic function is introduced. Concretely, the $H_1$ curve is the solution of the initial value problem

$$\begin{array}{c}{\displaystyle \frac{d}{dt}}x\left(t\right)=\left(\rho +{\displaystyle \frac{\tau }{\sqrt{1+t^2}}}\right){\displaystyle \frac{M-x\left(t\right)}{M}}x\left(t\right)\hfill \\ \\ x\left(t_0\right)=x_0,x_0>0,\hfill \end{array}$$

(1)

resulting in

$$x\left(t\right)=\frac{M}{1+\left(\frac{M}{x_0}-1\right){ϵ}_{{\nu }_1}{\left(t_0\right)}^{-1}{ϵ}_{{\nu }_1}\left(t\right)},{\nu }_1={\left(\rho ,\tau \right)}^T,$$

(2)

where ${ϵ}_{{\nu }_1}\left(t\right):=exp(-\rho t-\tau arcsinh\left(t\right))$. Here, M represents the carrying capacity of the system under consideration whereas real parameters $\tau$ and $\rho$ jointly determine growth rate. We note that parameter $\rho$ gives the intrinsic growth rate, whereas $\tau$ controls the distance from the $H_1$ curve to the logistic one. Regarding function ${ϵ}_{{\nu }_1}$, it will play an important role in the unification of the strategies required to work with this and similar models, as it encodes the action of the temporal function that distinguishes models I and III, and which also characterizes the type-II model. Another interesting aspect is the value of the asymptote. As it is verified that ${\lim }_{t\to \infty }{ϵ}_{{\nu }_1}\left(t\right)=0$, then ${\lim }_{t\to \infty }x\left(t\right)=M$ if $\rho >0$, and ${\lim }_{t\to \infty }x\left(t\right)=0$ if $\rho <0$. In any case, the value of M is independent from the initial value, which can significantly restrict some applications. In practice, there are situations in which the growth phenomenon under study shows an $H_1$-type sigmoidal growth and several sample paths are available, each with the same growth pattern but with different initial values and upper bounds. This can be avoided by re-parameterizing the curve. Indeed, after replacing ${ϵ}_{{\nu }_1}\left(t_0\right)/\left(M/x_0-1\right)$ with a new parameter $\eta$, and M with $x_0(\eta +{ϵ}_{{\nu }_1}\left(t_0\right))/\eta$, expression (2) becomes

$$x_{{\theta }_1}\left(t\right)=x_0\frac{\eta +{ϵ}_{{\nu }_1}\left(t_0\right)}{\eta +{ϵ}_{{\nu }_1}\left(t\right)},$$

for which it is verified that ${\lim }_{t\to \infty }{x}_{{\theta }_1}\left(t\right)=x_0(\eta +{ϵ}_{{\nu }_1}\left(t_0\right))/\eta$. Here, ${\theta }_1$ represents parametric vector ${\theta }_1={\left(\eta ,{\nu }_1^T\right)}^T$. This last reformulation of the model allows us to describe trajectories with different values of asymptotic stabilization in terms of the initial state $x_0$.

Regarding the hyperbolastic curve of type II ($H_2$), its introduction presents an increase in complexity with respect to the type-I model. This curve is defined as

$$y\left(t\right)=\frac{M}{1+\left(\frac{M}{y_0}-1\right){\displaystyle \frac{arcsinh{ϵ}_{{\nu }_2}\left(t\right)}{arcsinh{ϵ}_{{\nu }_2}\left(t_0\right)}}},{\nu }_2={\left(\beta ,\gamma \right)}^T,$$

(3)

being the solution of the ordinary differential equation

$$\frac{d}{dt}y\left(t\right)=\alpha \beta \gamma t^{\gamma -1}{y}^2\left(t\right)tanh\left(\frac{M-y\left(t\right)}{\alpha y\left(t\right)}\right)M^{-1},$$

with initial condition $y\left(t_0\right)=y_0$. Here, ${ϵ}_{{\nu }_2}\left(t\right):=exp(-\beta t^{\gamma })$, parameters M, $\beta$, and $\gamma$ can take any real value and

$$\alpha =\left(\frac{M-y_0}{y_0}\right)/arcsinh{ϵ}_{{\nu }_2}\left(t_0\right).$$

In this case, $\beta$ represents the intrinsic growth rate whereas $\gamma$ is an allometric constant. Note also that (3) tends to M when $\beta >0$, and to extinction ($y\left(t\right)\to 0$) when $\beta <0$. This leads $\beta$ to assume a key property in the dynamic evolution of the model. Again, this goes to show how the limit value does not depend on the initial condition. By applying the same reasoning as in the previous case, we arrive at the following reformulation of curve

$$y_{{\theta }_2}\left(t\right)=y_0\frac{\eta +arcsinh{ϵ}_{{\nu }_2}\left(t_0\right)}{\eta +arcsinh{ϵ}_{{\nu }_2}\left(t\right)},$$

verifying that ${lim}_{t\to \infty }{y}_{{\theta }_2}\left(t\right)=y_0(\eta +{ϵ}_{{\nu }_2}\left(t_0\right))/\eta$, being now ${\theta }_2={\left(\eta ,{\nu }_2^T\right)}^T$.

Finally, the hyperbolastic model of type III ($H_3$) follows the structure established for the $H_1$ curve, but this time using the classic Weibull ordinary differential equation as a starting point. In this way, the $H_3$ ordinary differential equation arises, resulting in

$$\frac{d}{dt}z\left(t\right)=\left(M-z\left(t\right)\right)\left(\beta \gamma t^{\gamma -1}+\frac{\tau }{\sqrt{1+{\tau }^2{t}^2}}\right).$$

(4)

Considering $z\left(t_0\right)=z_0>0$, the solution of (4) is

$$z\left(t\right)=M-\left(M-z_0\right)\frac{{ϵ}_{{\nu }_3}\left(t\right)}{{ϵ}_{{\nu }_3}\left(t_0\right)},{\nu }_3={\left(\beta ,\gamma ,\tau \right)}^T,$$

being ${ϵ}_{{\nu }_3}\left(t\right):=exp\left(-\beta t^{\gamma }-arcsinh\left(\tau t\right)\right)$. The real parameters $\beta$ and $\gamma$ are an intrinsic growth rate and an allometric value, respectively, whereas $\tau$ represents the distance from the $H_3$ curve to the classical Weibull model. By applying to the limit value of the curve a reasoning analogous to the ones described above, we arrive at the following reformulation of the curve

$$z_{{\theta }_3}\left(t\right)=z_0\frac{\eta -{ϵ}_{{\nu }_3}\left(t\right)}{\eta -{ϵ}_{{\nu }_3}\left(t_0\right)},{\mathit{\theta }}_3={\left(\eta ,{\nu }_3^T\right)}^T,$$

for which it is verified that ${lim}_{t\to \infty }{z}_{{\theta }_3}\left(t\right)=z_0\eta /(\eta -{ϵ}_{{\nu }_3}\left(t_0\right))$.

As it has been shown, each of the hyperbolastic curves is the solution to a certain ordinary differential equation. However, it is not difficult to verify that each curve also verifies a Malthusian-type differential equation. Specifically, one has

$$\frac{d}{dt}{j}_{\ell }\left(t\right)=h_{{\theta }_{\ell }}\left(t\right)j_{\ell }\left(t\right),\ell =1,2,3,$$

(5)

where $j_{\ell }={\delta }_{1\ell }{x}_{{\theta }_1}+{\delta }_{2\ell }{y}_{{\theta }_2}+{\delta }_{3\ell }{z}_{{\theta }_3}$, being ${\delta }_{··}$ the Kronecker delta and

$$h_{{\theta }_{\ell }}\left(t\right)=\left\{\begin{array}{cc}{\displaystyle \frac{{ϵ}_{{\nu }_1}\left(t\right)}{\eta +{ϵ}_{{\nu }_1}\left(t\right)}}\left(\rho +{\displaystyle \frac{\tau }{\sqrt{1+t^2}}}\right)\hfill & \mathrm{if}\ell =1,\hfill \\ \\ {\displaystyle \frac{\rho {ϵ}_{{\nu }_2}\left(t\right)}{\eta +arcsinh{ϵ}_{{\nu }_2}\left(t\right)}}{\left(1+{ϵ}_{{\nu }_2}{\left(t\right)}^2\right)}^{-\frac{1}{2}}\hfill & \mathrm{if}\ell =2,\hfill \\ \\ {\displaystyle \frac{{ϵ}_{{\nu }_3}\left(t\right)}{\eta -{ϵ}_{{\nu }_3}\left(t\right)}}\left(\beta \gamma t^{\gamma -1}+{\displaystyle \frac{\tau }{\sqrt{1+{\tau }^2{t}^2}}}\right)\hfill & \mathrm{if}\ell =3.\hfill \end{array}\right.$$

Note that the $h_{{\theta }_{\ell }}$ functions do not depend on the initial values taken by the curves. This unified formulation allows us to extend the hyperbolastic deterministic procedure to stochastic versions under the prism of the inhomogeneous lognormal diffusion process, which is the objective of the next section.

Finally, note that an important characteristic of this type of sigmoidal models is the presence of inflection points, as well as the expressions that they must verify. Note that the number of inflections generally increases the flexibility of the curve and its adaptation to a variety of observations. From (5), the inflection time instants are the solutions of the equation

$$\frac{d}{dt}{h}_{{\theta }_{\ell }}\left(t\right)+h_{{\theta }_{\ell }}^2\left(t\right)=0,\ell =1,2,3.$$

(6)

In general, the resulting equations are quite complex and in practical cases strategies for their estimation will be developed from the observed data. On the other hand, such an expression will allow us to integrate a procedure for calculating initial solutions for the different parameters, as we will see in later sections.

3. Hyperbolastic Diffusion Processes

The deterministic hyperbolastic models presented in the previous section enable the description of growth phenomena in contexts in which detailed knowledge of their leading variables exists. However, such models do not account for the effects of external factors influencing the development of phenomena. For this reason, we must address the development of stochastic models which, while based on the previous ones, attempt to overcome their shortcomings.

Such models are commonly built by adding a multiplicative noise to the original ordinary differential equations. The multiplicative noise is considered more suitable than the additive one in the context of growth phenomena. Thus, random influences are made dependent on the state of the system at each instant.

However, this methodology can lead to stochastic differential equations whose solution is not easy to handle, mainly due to lack of explicit characteristics such as transition densities, which may affect practical applications. For example, consider the case of the $H_1$ hyperbolastic model. From (1) we consider the stochastic differential equation

$$dX\left(t\right)=\left(\rho +\frac{\tau }{\sqrt{1+t^2}}\right)\frac{M-X\left(t\right)}{M}X\left(t\right)dt+\sigma X\left(t\right)dW\left(t\right),$$

where $W\left(t\right)$ represents the standard Wiener process, independent of initial condition $X\left(t_0\right)=X_0$ for $t\ge t_0$, being $\sigma >0$ the volatility parameter.

This equation can be solved by using the following integrating factor (see [32])

$$F\left(t\right)=exp\left(\frac{{\sigma }^2}{2}(t-t_0)-\sigma \left(W\left(t\right)-W\left(t_0\right)\right)\right),$$

from which the solution can be expressed as

$$X\left(t\right)=\frac{X_0{ϵ}_{{\nu }_1}\left(t_0\right)/\left({ϵ}_{{\nu }_1}\left(t\right)F\left(t\right)\right)}{{\displaystyle 1+{\displaystyle \frac{X_0}{M}}{\int }_{t_0}^t\left(\rho +\frac{\tau }{\sqrt{1+s^2}}\right)\frac{{ϵ}_{{\nu }_1}\left(t_0\right)}{{ϵ}_{{\nu }_1}\left(s\right)}F{\left(s\right)}^{-1}ds}},$$

where ${ϵ}_{{\nu }_1}\left(t\right)=exp(-\rho t-\tau arcsinh\left(t\right))$.

Note that, although the stochastic differential equation can be solved, the form of the solution does not provide an explicit distribution of $X\left(t\right)$, which renders the analytic treatment of the process difficult. In fact, to this purpose the literature discretizes the stochastic equation and obtains numerical approximations. This is the rationale behind the next subsection, in which we will consider stochastic, more manageable versions of the models.

Hyperbolastic Diffusion Processes from the Non-Homogeneous Lognormal Diffusion Process

Based on the previous comments, this section proposes the construction of hyperbolastic stochastic models which are analytically manageable and in which the mean of the process can be associated with the original curve. To this end, we will employ the lognormal process with exogenous factors, which has the aforementioned properties and allows us to work extensively with the resulting diffusions. As is known, the inhomogeneous lognormal process is built on the basis of a differential equation of the Malthusian type with a time-dependent growth rate, to which a multiplicative noise is added as described above (see Roman et al. [33]). Therefore, and given that the hyperbolic curves all verify a Malthusian-type differential equation, it is therefore possible to associate each of them with an inhomogeneous lognormal diffusion process whose mean is precisely the curve in question. In addition, the properties of this type of process make it perfectly treatable, since its finite-dimensional distributions are explicitly available.

Briefly, the non-homogeneous lognormal diffusion process (or with exogenous factors) is the solution to the linear stochastic differential equation

$$dX\left(t\right)=h_{{\theta }_{\ell }}\left(t\right)X\left(t\right)dt+\sigma X\left(t\right)dW\left(t\right),\ell =1,2,3,$$

with initial condition $X\left(t_0\right)=X_0$. $h_{{\theta }_{\ell }}$ is a continuous function depending on a vectorial parameter ${\mathit{\theta }}_{\ell }\in {\Theta }_{\ell }\subseteq {\mathbb{R}}^k$ whereas $W\left(t\right)$ is the standard Wiener process, independent of $X_0$ for $t\ge t_0$. The solution of this equation is

$$X\left(t\right)=X_{0}exp\left(H_{{\xi }_{\ell }}(t,t_0)+\sigma \left(W\left(t\right)-W\left(t_0\right)\right)\right),{\xi }_{\ell }={\left({\mathit{\theta }}_{\ell }^T,{\sigma }^2\right)}^T,$$

(7)

where

$$H_{{\xi }_{\ell }}(t_0,t)={\int }_{t_0}^t{h}_{{\theta }_{\ell }}\left(u\right)du-\frac{{\sigma }^2}{2}(t-t_0).$$

Indeed, if $Z\left(t\right)=logX\left(t\right)$, then, by Itô’s formula, it follows that

$$dZ\left(t\right)=\left(h_{{\theta }_{\ell }}\left(t\right)-\frac{{\sigma }^2}{2}\right)dt+\sigma dW\left(t\right),$$

which leads to a Gaussian process Z, with X then being lognormal.

Given its Markovian character, the distribution of the process follows from that of $X_0$. Indeed, if $X_0$ is a lognormal distribution ${\mathsf{\Lambda }}_1\left[{\mu }_0,{\sigma }_0^2\right]$, or $X_0$ is degenerate at a value $x_0$ ($X_0=x_0$ a.s.), it can be shown that the finite-dimensional distributions of the process follow a lognormal distribution. Specifically, vector ${\left(X\left(t_1\right),\dots ,X\left(t_n\right)\right)}^T$, for $t_1<\dots <t_n$ with $n\in \mathbb{N}$, follows a n-dimensional lognormal distribution ${\mathsf{\Lambda }}_n\left[\mathit{ϵ},\mathbf{\Sigma }\right]$ where the components of vector $\mathit{\epsilon }$ and matrix $\mathbf{\Sigma }$ are

$${\epsilon }_i={\mu }_0+H_{{\xi }_{\ell }}(t_0,t_i),i=1,\dots ,n$$

and

$${\sigma }_{ij}={\sigma }_0^2+{\sigma }^2(min(t_i,t_j)-t_0),i,j=1,\dots ,n,$$

respectively. In particular, from the two-dimensional distributions $\left(X\left(t\right)\right|X\left(s\right))=x_s$ it follows that

$$X\left(t\right)|X\left(s\right)=x_s\sim {\mathsf{\Lambda }}_1\left[log{x}_s+H_{{\xi }_{\ell }}(s,t);{\sigma }^2(t-s)\right],s<t$$

so the transition probability density function takes the form

$$f(x,t|x_s,s)=\frac{1}{x\sqrt{2\pi {\sigma }^2(t-s)}}exp\left(-\frac{log(x/x_s)-H_{{\xi }_{\ell }}(s,t)}{2{\sigma }^2(t-s)}\right).$$

(8)

Regarding the main characteristics of the process, we highlight the moments of order n, whose expression is

$$E\left[X{\left(t\right)}^n\right]=E\left[X_0^n\right]exp\left(n{\int }_{t_0}^t{h}_{{\theta }_{\ell }}\left(u\right)du\right).$$

From this expression, it is evident that for $n=1$ a convenient choice of function $h_{{\theta }_{\ell }}$ allows many types of curves to be modeled by this process.

Note that by using this methodology in the construction of the processes, certain issues are addressed which could not be if the classical method of introducing a multiplicative noise in the original deterministic differential equation had been followed. Ultimately, we must not forget that our goal is to model phenomena which exhibit behaviors that are governed by these types of curves.

Figure 1 shows, for each hyperbolastic process, a comparative graph of the means obtained from simulated trajectories using each stochastic differential equation (the original and the Malthusian). This illustrates how the behavior is the same in both cases.

4. Inference

As shown above, one of the main strength of the models presented here is their flexibility, which makes them suitable for describing and fitting growth patterns in real cases.

There are, at least, two natural ways to identify the parameters in our stochastic hyperbolastic models. The first is the direct use of likelihood methods, and the second is particle filtering (PF). Herein, we will employ the former due to its widespread use. However, we also want to point out that our unification and linearization work applies equally well to PF (see, e.g., in [34,35]). In this method, the lognormality assumption on $X_0$ is not required. The parameters to be identified would become the signal (with trivial dynamics) and the logarithm of the linearized hyperbolastic SDEs would become the classical observations. Finally, one could use PF on the non-linearized SDEs, which would require a little manipulation of the observation function coming from the integrating factor solution (shown in Section 3 for the type-I case).

Despite the above, in this section we address the estimation of the parameters of the model by means of the maximum likelihood (ML) method. However, it must be taken into account that the functional and parametric complexity of the models makes the estimation procedure remarkably difficult. In order to tackle this issue, the estimation by maximum likelihood will require the use of numerical procedures, many of which will demand the computation of initial solutions for the parameters.

4.1. Maximum Likelihood Equations

The starting point for the estimation procedure is the observation of d sample-paths at time instants $t_{ij}$, $(i=1,\dots ,d,j=1,\dots ,n_i)$. Note that neither the sample sizes nor the times of observation have to be the same. Moreover, we will suppose that $t_{i1}=t_0$, for all $i=1,\dots ,d$. Let ${\mathbf{X}}_i^T$ be the vector containing the random variables of the i-th sample-path, that is ${\mathbf{X}}_i={(X\left(t_0\right),\dots ,X\left(t_{i,n_i}\right))}^T$, $i=1,\dots ,d$, and denote $\mathbf{X}={\left({\mathbf{X}}_1^T|\cdots |{\mathbf{X}}_d^T\right)}^T$. Furthermore, we will consider that $X_0$ follows a lognormal distribution, concretely $X_0\sim {\mathsf{\Lambda }}_1[{\mu }_1,{\sigma }_1^2]$ (as it has been said before, this would not be a requirement if the particle filtering approach were to be used).

From (8), and considering a fixed value $\mathbf{x}$ of $\mathbf{X}$, the log-likelihood function for each model, defined in a unified way for $\ell =1,2,3$, is

$$\begin{array}{cc}\hfill log{\mathbb{L}}_{\mathbf{x}}^{\ell }\left(\mathit{\varsigma },{\mathit{\xi }}_{\ell }\right)=& -\frac{(n+d)log2\pi }{2}-\frac{dlog{\sigma }_1^2}{2}-\frac{nlog{\sigma }^2}{2}\hfill \\ \hfill & -\sum _{i=1}^{d}log{x}_{i1}-\frac{1}{2{\sigma }_1^2}\sum _{i=1}^d{\left(log{x}_{i1}-{\mu }_1\right)}^2-\frac{Z_1+{\mathsf{\Phi }}_{{\xi }_{\ell }}-2{\mathsf{\Gamma }}_{{\xi }_{\ell }}}{2{\sigma }^2}\hfill \end{array}$$

where $n={\sum }_{i=1}^d(n_i-1)$, $\mathit{\varsigma }={({\mu }_1,{\sigma }_1^2)}^T$ is the parametric vector of the initial distribution, and

$$\begin{array}{cc}\hfill Z_1& =\sum _{i=1}^d\sum _{j=2}^{n_i}\frac{1}{{\Delta }_{ij}}{\left(log\frac{x_{ij}}{x_{ij-1}}\right)}^2,\hfill \\ \hfill {\mathsf{\Phi }}_{{\xi }_{\ell }}& =\sum _{i=1}^d\sum _{j=2}^{n_i}\frac{1}{{\Delta }_{ij}}{H}_{{\xi }_{\ell }}^2(t_{ij-1},t_{ij}),\hfill \\ \hfill {\mathsf{\Gamma }}_{{\xi }_{\ell }}& =\sum _{i=1}^d\sum _{j=2}^{n_i}\frac{1}{{\Delta }_{ij}}log\frac{x_{ij}}{x_{ij-1}}{H}_{{\xi }_{\ell }}(t_{ij-1},t_{ij}),\hfill \end{array}$$

being ${\Delta }_{ij}=t_{ij}-t_{i,j-1}$. The expression of function $H_{{\xi }_{\ell }}$ will depend on the hyperbolastic case under study. Nevertheless, it is possible to find a unified formulation by considering $H_{{\xi }_{\ell }}(s,t)={\mathcal{H}}_{{\theta }_{\ell }}(s,t)-\frac{{\sigma }^2}{2}(t-s)$, where

$${\mathcal{H}}_{{\theta }_{\ell }}(s,t):={\zeta }_{\ell }log\frac{\eta +v_{\ell }\left({ϵ}_{{\nu }_{\ell }}\left(s\right)\right)}{\eta +v_{\ell }\left({ϵ}_{{\nu }_{\ell }}\left(t\right)\right)},\ell =1,2,3,$$

with

$${\zeta }_{\ell }=\left\{\begin{array}{cc}1\hfill & \mathrm{if}\ell =1,2,\hfill \\ -1\hfill & \mathrm{if}\ell =3,\hfill \end{array}\right.$$

and

$$v_{\ell }\left(u\right)=\left\{\begin{array}{cc}u\hfill & \mathrm{if}\ell =1,\hfill \\ arcsinh\left(u\right)\hfill & \mathrm{if}\ell =2,\hfill \\ -u\hfill & \mathrm{if}\ell =3.\hfill \end{array}\right.$$

Under the assumption that $\mathit{\varsigma }$ and ${\mathit{\xi }}_{\ell }$ are functionally independent, the estimate of $\mathit{\varsigma }$ is

$${\displaystyle {\widehat{\mu }}_1=\frac{1}{d}\sum _{i=1}^{d}log{x}_{i1}\mathrm{and}{\widehat{\sigma }}_1^2=\frac{1}{d}\sum _{i=1}^d{(log{x}_{i1}-{\widehat{\mu }}_1)}^2,}$$

(9)

while that of ${\mathit{\xi }}_{\ell }$ is obtained (see [33] for details) from the system of equations

$$\begin{array}{cc}\hfill {\Xi }_{{\theta }_{\ell }}+\frac{{\sigma }^2}{2}{\mathsf{\Omega }}_{{\theta }_{\ell }}& =0,\hfill \\ \hfill {\sigma }^2\left(n+\frac{{\sigma }^2}{4}{Z}_2\right)+2X_2^{{\theta }_{\ell }}-X_1^{{\theta }_{\ell }}-Z_1& =0,\hfill \end{array}$$

(10)

where the grouped variables are defined as follows:

$$\begin{array}{cc}\hfill {\Xi }_{{\theta }_{\ell }}& =\sum _{i=1}^d\sum _{j=2}^{n_i}{\Delta }_{ij}^{-1}\left(log\frac{x_{ij}}{x_{ij-1}}-{\mathcal{H}}_{{\theta }_{\ell }}(t_{ij-1},t_{ij})\right)\frac{\partial }{\partial {\mathit{\theta }}_{\ell }}{\mathcal{H}}_{{\theta }_{\ell }}(t_{ij-1},t_{ij}),\hfill \\ \hfill {\mathsf{\Omega }}_{{\theta }_{\ell }}& =\sum _{i=1}^d\sum _{j=2}^{n_i}\frac{\partial }{\partial {\mathit{\theta }}_{\ell }}{\mathcal{H}}_{{\theta }_{\ell }}(t_{ij-1},t_{ij}),\hfill \\ \hfill X_1^{{\theta }_{\ell }}& =\sum _{i=1}^d\sum _{j=2}^{n_i}{\Delta }_{ij}^{-1}{\mathcal{H}}_{{\theta }_{\ell }}^2(t_{ij-1},t_{ij}),\hfill \\ \hfill X_2^{{\theta }_{\ell }}& =\sum _{i=1}^d\sum _{j=2}^{n_i}{\Delta }_{ij}^{-1}\left(log\frac{x_{ij}}{x_{ij-1}}{\mathcal{H}}_{{\theta }_{\ell }}(t_{ij-1},t_{ij})\right),\hfill \\ \hfill Z_2& =\sum _{i=1}^d\sum _{j=2}^{n_i}{\Delta }_{ij}.\hfill \end{array}$$

The partial derivatives of ${\mathcal{H}}_{{\theta }_{\ell }}$ are calculated with respect to the parameters of every model, and therefore they depend on ℓ. Indeed, these parameters involve a common one, say $\eta$, as well as the particular ones of each model: ${\nu }_1={(\rho ,\tau )}^T$ for the $H_1$ model, ${\nu }_2={(\beta ,\gamma )}^T$ for the $H_2$ model, and ${\nu }_3={(\beta ,\gamma ,\tau )}^T$ in the case of the $H_3$ model.

All these expressions can be unified as follows:

$$\begin{array}{cc}\hfill \frac{\partial }{\partial \eta }{\mathcal{H}}_{{\theta }_{\ell }}(s,t)& =\frac{{\zeta }_{\ell }}{\eta +v_{\ell }\left({ϵ}_{{\nu }_{\ell }}\left(s\right)\right)}-\frac{{\zeta }_{\ell }}{\eta +v_{\ell }\left({ϵ}_{{\nu }_{\ell }}\left(t\right)\right)},\hfill \\ \hfill \frac{\partial }{\partial {\nu }_{\ell }}{\mathcal{H}}_{{\theta }_{\ell }}(s,t)& =\frac{{\zeta }_{\ell }{v}_{\ell }^{\prime }\left({ϵ}_{{\nu }_{\ell }}\left(t\right)\right)}{\eta +v_{\ell }\left({ϵ}_{{\nu }_{\ell }}\left(t\right)\right)}\frac{\partial }{\partial {\nu }_{\ell }}{ϵ}_{{\nu }_{\ell }}\left(t\right)-\frac{{\zeta }_{\ell }{v}_{\ell }^{\prime }\left({ϵ}_{{\nu }_{\ell }}\left(s\right)\right)}{\eta +v_{\ell }\left({ϵ}_{{\nu }_{\ell }}\left(s\right)\right)}\frac{\partial }{\partial {\nu }_{\ell }}{ϵ}_{{\nu }_{\ell }}\left(s\right),\hfill \end{array}$$

where $v_{\ell }^{\prime }\left(u\right)={\displaystyle \frac{d{v}_{\ell }\left(u\right)}{du}}$ and $\ell =1,2,3$.

Despite the few parameters involved, this system remains too complex to be solved with analytical methods. Therefore, numerical methods, such as the Newton–Raphson procedure, are the proper way to solve the system. Nevertheless, initial solutions for the parameters are required in order to apply such methods. In this case, such solutions can be obtained from the information provided by the data and taking into account the properties of the curve.

4.2. Initial Solutions

As we mentioned before, solving the system of Equation (10) by means of numerical methods requires finding initial solutions first. This task may prove daunting according to the degree of complexity of the hyperbolastic models. In this section, strategies for the search of initial solutions are proposed for each hyperbolastic stochastic model.

The vector of parameters to be estimated is ${\mathit{\xi }}_{\ell }={({\mathit{\theta }}_{\ell }^T,{\sigma }^2)}^T$, where ${\mathit{\theta }}_{\ell }$ comprises $\eta$ and the rest of the parameters of the model under study. All of these parameters are related within the expressions derived from each model, and therefore the initial value of one parameter may depend on the initial value of another. A strategy is then required which will provide a sequence of steps that guarantee the final computation of all values.

As shown in Section 2, type-I and -II models have two parameters each, besides parameter $\eta$, which is common to both, and diffusion parameter ${\sigma }^2$. On the other hand, type-III models have three parameters: $\beta ,\gamma$, and $\tau$. This model will require a particular step in the procedure. We propose the following search schema.

• Initial solution for ${\sigma }^2$: it is known that, for a lognormal distribution ${\mathsf{\Lambda }}_1[\epsilon ,\chi ]$, an estimation of $\chi$ is provided by the quotient between the arithmetic and geometric means. By taking advantage of this result, the initial value of ${\sigma }^2$, in the case of a degenerate initial variable $X_0$, is obtained by linear regression of ${\sigma }_i^2=2log(m_i/m_i^g)$ against $t_i$, for $i=1,\dots ,n$. Here, $m_i$ and $m_i^g$ are, respectively, the arithmetic mean and the geometric mean of the sample paths at time instant $t_i$. If $X_0$ is not degenerate, i.e., ${\sigma }_0^2>0$, then ${\sigma }_i^2$ provides an estimate of ${\sigma }_0^2+{\sigma }^2(t_i-t_0)$ for each $t_i$, thus a previous estimation of ${\sigma }_0^2$ is required. This can be accomplished from the values of the sample paths at the first observation time $t_0$ by using (9).
• Initial solution for the rest of the parameters: the general procedure involves considering the two parameters at once (in the type-III model, one of the parameters must be fixed) and using expressions derived from the curve. Indeed, we can see how, for each of the reformulations made for the curves, parameter $\eta$ can be expressed as a function of the rest of the parameters and the limit value of the curve. Then $\eta$ can be viewed as a function ${\eta }_{\ell }\left({\mathit{\nu }}_{\ell }\right)$ of the parameters of every curve. Note that this can be accomplished because ${lim}_{t\to \infty }{ϵ}_{{\nu }_{\ell }}\left(t\right)=0$ for all $\ell =1,2,3$. The expressions for every model are then

  $${\eta }_{\ell }\left({\mathit{\nu }}_{\ell }\right)=\left\{\begin{array}{cc}{ϵ}_{{\nu }_1}\left(t_0\right){\left(k_1/x_0-1\right)}^{-1}\hfill & \mathrm{if}\ell =1,\hfill \\ arcsinh{ϵ}_{{\nu }_2}\left(t_0\right){\left(k_2/x_0-1\right)}^{-1}\hfill & \mathrm{if}\ell =2,\hfill \\ {ϵ}_{{\nu }_3}\left(t_0\right){\left(1-x_0/k_3\right)}^{-1}\hfill & \mathrm{if}\ell =3,\hfill \end{array}\right.$$

  where $k_1,k_2$, and $k_3$ are the limit values of every curve, respectively. Given that these values are, as a general rule, not known, an approximation that uses the last values of the mean observed data can be considered.
  The goal then is to use the information provided by the observed inflection time of the curve. Let $t_{*}$ be an inflection time instant, which can be obtained by fitting a spline function to the observed mean values and finding either the maximum of its first derivative or a zero value of its second derivative. Note that $t_{*}$ must verify Equation (6). This fact will allow us to obtain the initial values of the rest of the parameters. To this end, we define function

  $$I\left({\mathit{\theta }}_{\ell }\right)=\left|\frac{d}{dt}{h}_{{\theta }_{\ell }}\left(t_{*}\right)+h_{{\theta }_{\ell }}^2\left(t_{*}\right)\right|.$$
  Here, notation $h_{{\theta }_{\ell }}\left(t\right)$ alludes to the fact that we have substituted functions ${\eta }_{\ell }\left({\nu }_{\ell }\right)$ in the expressions of functions $h_{{\theta }_{\ell }}\left(t\right)$, so that function I depends only on the parameters of interest. Finally, the selected values will be those that minimize function I in a bounded region. In particular, candidates can be obtained from such bounded region satisfying the condition

  $$I\left({\mathit{\theta }}_{\ell }\right)\le \epsilon ,$$

  where $\epsilon >0$ is a prescribed error threshold. We recommend varying its order of magnitude between 0.0001 and 0.1, depending on the sample data. At this step, and in order to achieve a good implementation, it is worth differentiating the following cases:

  (a)

  Type-I and -II models: initial values for pairs of parameters $(\rho ,\tau )$ (model I) and $(\beta ,\gamma )$ (model II) can be obtained by applying the aforementioned methodology where the bounded regions are two-dimensional in every case. The limits of such regions can be proposed following the nature of the parameters. $\rho$ and $\beta$ are growth parameters in their respective models, so their magnitude would be related with the observed growth. On the other side, $\tau$ is the “distance” to the logistic curve and $\gamma$ an allometric constant, and therefore their bounds can also be proposed according to the phenomena observed.

  (b)

  Type-III model: there are three parameters in vector ${\nu }_3$. In order to avoid a three-dimensional region, we propose fixing parameter $\gamma$ provided that $\tau$ is the “distance” to the Weibull curve. That is, by fitting a Weibull model to the observed data, an estimate of $\gamma$ could be considered the initial value of the parameter. Next, the methodology described above could be applied for the pair $(\beta ,\tau )$. Simulation studies may be carried out in order to verify the stability of the solutions according to fluctuations in the fixed value of $\gamma$. On the other hand, considerations about the bounded region are then the same as in the other two models.

• Initial solution for $\eta$: once solutions for the other parameters have been obtained, an initial value of $\eta$ follows from applying ${\eta }_{\ell }\left({\nu }_{\ell }\right)$ to such initial values. As it has been said, limiting values $k_{\ell }$ can be obtained by taking the last value of the observed mean data. Nevertheless, if $t_0=0$, then ${ϵ}_{{\nu }_{\ell }}\left(t_0\right)=1$ for all $\ell =1,2,3$, and the initial value of $\eta$ is independent of the rest of the parameters.

5. Simulation

In this section, the proposed methodology is applied to simulated data. For the three hyperbolastic models, a simulation study is carried out concerning the variability of the diffusion parameter, named $\sigma$, and the number n of observations. For every model, 20 sample paths were simulated with the same final time instant at $t=5$ and time step $0.1$. Thus, 501 points for each path were simulated in every iteration. In all three cases, an initial degenerate distribution at $x_0=1$ were considered, in the sense $P(X\left(t_0\right)=x_0)=1$. Furthermore, three different diffusion coefficient values were used in every simulation, being $\sigma =0.01,0.05$ and $0.075$. Then, three subsets of equidistant points were extracted according to values $n=11,31$ and 51. There were, therefore, 9 sets of data for combinations of n and $\sigma$. In every case, 100 replications were carried out. In every one of these, the sample paths were simulated following (7) and the proposed methodology was applied. In particular, initial solutions were considered according to Section 4.2. Then, the corresponding system of ML equations was solved following Section 4.1. In order to obtain good performance and improved results, after every replication the mean value between the current solution and the mean of the solutions from previous iterations were computed. The particular parameters for every hyperbolastic models were:

• for the $H_1$ model: $t_0=0,\eta =0.05,\rho =1,\tau =0.25$,
• for the $H_2$ model: $t_0={10}^{-3},\eta =1.5,\beta =0.4,\gamma =2.5$,
• for the $H_3$ model: $t_0={10}^{-3},\eta =2,\beta =0.1,\gamma =2.3,\tau =0.05$.

Finally, the relative absolute error (RAE) from the observed mean to a new simulated mean with estimated parameters was computed. The expression for RAE is done by

$$\mathrm{RAE}=\frac{1}{n}\sum _{i=1}^n\frac{\left|m_i-\widehat{E}\left(X\left(t_i\right)\right)\right|}{m_i},$$

where $m_i$ is the mean observed value at time $t_i$, $\widehat{E}\left(X\left(t_i\right)\right)$ is the mean value of the estimated process at time $t_i$, and n is the sample size.

A total number of 2700 computations were performed. All of them were performed in R. The solutions of ML equations were obtained with the packages “nleqslv” [36] and “BB” [37].

Results are shown in the following tables for initial and final solutions. In the latter case, an additional column showing the RAE has been included.

After a more in-depth analysis, it is shown that the RAE is low in all the three hyperbolastic models. Note that this error measures the difference between the observed mean and the mean of a new diffusion process simulated with the estimated parameters. Furthermore, with $\sigma$ being the parameter concerning the variability of the paths, it is expected that the error growths as $\sigma$ does. Nevertheless, the combination with different sample sizes (n) can lead to particular behaviors for each hyperbolastic model.

In the case of the $H_1$ model,

Table 1. Final solutions for simulated data of the $H_1$ process and RAE to estimated mean. Original parameters are $\eta =0.05,\rho =1$, and $\tau =0.25$.

Combinations		$\widehat{\mathit{\eta }}$	$\widehat{\mathit{\rho }}$	$\widehat{\mathit{\tau }}$	${\widehat{\mathit{\sigma }}}^2$	RAE
$\sigma =0.01$	$n=11$	$0.0500$	$0.9994$	$0.2508$	$0.0082$	$0.0017$
	$n=31$	$0.0501$	$1.0099$	$0.2337$	$0.0121$	$0.0019$
	$n=51$	$0.0500$	$0.9991$	$0.2511$	$0.0095$	$0.0019$
$\sigma =0.05$	$n=11$	$0.0500$	$0.9978$	$0.2528$	$0.0408$	$0.0087$
	$n=31$	$0.0501$	$1.0070$	$0.2375$	$0.0472$	$0.0090$
	$n=51$	$0.0499$	$0.9961$	$0.2550$	$0.0477$	$0.0092$
$\sigma =0.075$	$n=11$	$0.0502$	$0.9977$	$0.2529$	$0.0612$	$0.0131$
	$n=31$	$0.0501$	$1.0056$	$0.2392$	$0.0702$	$0.0136$
	$n=51$	$0.0499$	$0.9946$	$0.2570$	$0.0715$	$0.0138$

shows that the RAE increases with $\sigma$ and, for every value of this parameter, it remains approximately the same for all sizes n. This leads to the conclusion that the $H_1$ model behaves well regardless of sample size. Let us remember that such model can be viewed as an extension of the classic logistic growth model. Furthermore, such invariant behavior for different values of n could be expected seeing the contents of

Table 2. Initial solutions for simulated data of the $H_1$ process.

Combinations		${\mathit{\eta }}_0$	${\mathit{\rho }}_0$	${\mathit{\tau }}_0$	${\mathit{\sigma }}_0^2$
$\sigma =0.01$	$n=11$	$0.0540$	$0.7327$	$0.2901$	$0.0097$
	$n=31$	$0.0540$	$0.8327$	$0.3301$	$0.0097$
	$n=51$	$0.0540$	$0.8578$	$0.3342$	$0.0097$
$\sigma =0.05$	$n=11$	$0.0540$	$0.7331$	$0.2901$	$0.0485$
	$n=31$	$0.0540$	$0.8232$	$0.3251$	$0.0485$
	$n=51$	$0.0540$	$0.8482$	$0.3304$	$0.0485$
$\sigma =0.075$	$n=11$	$0.0540$	$0.7302$	$0.2891$	$0.0727$
	$n=31$	$0.0540$	$0.8165$	$0.3211$	$0.0728$
	$n=51$	$0.0540$	$0.8455$	$0.3266$	$0.0728$

, which shows that initial values for parameters $\eta$ and $\sigma$ are exactly the same for every n. In the case of $\eta$, they are also the same for all $\sigma$. Figure 2 shows a comparison between the observed sample mean and the mean of the simulated new process after extracting the corresponding number of equidistant points. It can be seen that the fitness of the model is good, even when variability is high. For that case, the last row of figures shows a slight growing deviation from the observed mean. Nevertheless, the overall behavior of the model is appropriate. On the other hand, the evolution of RAE alongside $\sigma$ and independently of n for the $H_1$ model can be observed clearly in the corresponding table. As an addition to this, and in order to illustrate the procedure, graphs related with the search of initial values for the $H_1$ model in the last combination, corresponding to $n=51$ and $\sigma =0.075$, are shown in Figure 3. Following the methodology proposed in Section 4.2, a search for initial solutions ${\rho }_0$ and ${\tau }_0$ had been performed in square $[0.25,1.5]\times [0.1,0.75]$ divided into 100 points at every dimension, amounting to ${100}^2$ points in the whole region. At every pair of points, the inflection function derived from (6) was tested for an error $ϵ=0.1$. The histograms for all suitable values of ${\rho }_0$ and ${\tau }_0$ are shown in Figure 3, respectively, as well as red lines at the selected values (the median). Furthermore, a plot showing levels of the inflection condition function is shown in the last image. The region at which the error is below of $ϵ$ is filled in blue. Finally, red lines from final ${\rho }_0$ and ${\tau }_0$ are crossed at each selected pair.

The final results for the model $H_2$ are shown in

Table 3. Final solutions for simulated data of the $H_2$ process and RAE to estimated mean. Original parameters are $\eta =1.5,\beta =0.4$ and $\gamma =2.5$.

Combinations		$\widehat{\mathit{\eta }}$	$\widehat{\mathit{\beta }}$	$\widehat{\mathit{\gamma }}$	${\widehat{\mathit{\sigma }}}^2$	RAE
$\sigma =0.01$	$n=11$	$1.4993$	$0.4002$	$2.5021$	$0.0099$	$0.0012$
	$n=31$	$1.5000$	$0.3966$	$2.5105$	$0.0101$	$0.0014$
	$n=51$	$1.5000$	$0.3998$	$2.5040$	$0.0100$	$0.0014$
$\sigma =0.05$	$n=11$	$1.4991$	$0.4006$	$2.5191$	$0.0495$	$0.0059$
	$n=31$	$1.5027$	$0.3951$	$2.5363$	$0.0500$	$0.0065$
	$n=51$	$1.5020$	$0.3988$	$2.5256$	$0.0501$	$0.0066$
$\sigma =0.075$	$n=11$	$1.5010$	$0.4006$	$2.5370$	$0.0742$	$0.0088$
	$n=31$	$1.5059$	$0.3941$	$2.5572$	$0.0749$	$0.0096$
	$n=51$	$1.5047$	$0.3981$	$2.5438$	$0.0751$	$0.0097$

. They shown the expected behavior, with an increment of the error proportional to that of related with $\sigma$. Nevertheless, this model exhibits the lower error among the three hyperbolastics. Initial solutions ${\beta }_0$ and ${\gamma }_0$ (see

Table 4. Initial solutions for simulated data of the $H_2$ process.

Combinations		${\mathit{\eta }}_0$	${\mathit{\beta }}_0$	${\mathit{\gamma }}_0$	${\mathit{\sigma }}_0^2$
$\sigma =0.01$	$n=11$	$1.4991$	$0.1780$	$2.3030$	$0.0096$
	$n=31$	$1.4991$	$0.3012$	$2.3071$	$0.0096$
	$n=51$	$1.4991$	$0.3720$	$2.3192$	$0.0096$
$\sigma =0.05$	$n=11$	$1.5020$	$0.1760$	$2.3038$	$0.0477$
	$n=31$	$1.5020$	$0.3230$	$2.3113$	$0.0477$
	$n=51$	$1.5020$	$0.3588$	$2.3186$	$0.0477$
$\sigma =0.075$	$n=11$	$1.5088$	$0.1759$	$2.3040$	$0.0716$
	$n=31$	$1.5088$	$0.3244$	$2.3113$	$0.0715$
	$n=51$	$1.5088$	$0.3591$	$2.3179$	$0.0716$

) were found in a region $[0.001,1]\times [2,4]$, where bounds for ${\gamma }_0$ are chosen according to the fact that it appears as the exponent of time, that is, $t^{\gamma -2}$, in the expression of the $H_2$ model. The comparisons between observed and estimated mean values for every combination are shown in Figure 4, where the influence of $\sigma$ can be observed in the variability of observed values. Such influence is particularly noticeable in this example because the curve reaches its limit behavior earlier than in other models. This implies that more observations are therefore subject to growing variability. Even in that case, the estimated mean points are a good fit for the observed ones, for all $\sigma$ and all n.

Finally, estimated values and initial solutions for the $H_3$ model are shown in

Table 5. Final solutions for simulated data of the $H_3$ process and RAE to estimated mean. Original parameters are $\eta =2,\beta =0.1,\gamma =2.3$, and $\tau =0.05$.

Combinations		$\widehat{\mathit{\eta }}$	$\widehat{\mathit{\beta }}$	$\widehat{\mathit{\gamma }}$	$\widehat{\mathit{\tau }}$	${\widehat{\mathit{\sigma }}}^2$	RAE
$\sigma =0.01$	$n=11$	$2.0137$	$0.0831$	$2.1378$	$0.0495$	$0.0098$	$0.0394$
	$n=31$	$2.0137$	$0.0926$	$2.1387$	$0.0485$	$0.0099$	$0.0293$
	$n=51$	$2.0137$	$0.0972$	$2.1347$	$0.0485$	$0.0099$	$0.0246$
$\sigma =0.05$	$n=11$	$2.0154$	$0.0862$	$2.1373$	$0.0502$	$0.0492$	$0.0359$
	$n=31$	$2.0154$	$0.0939$	$2.1331$	$0.0488$	$0.0493$	$0.0294$
	$n=51$	$2.0154$	$0.0987$	$2.1244$	$0.0487$	$0.0493$	$0.0256$
$\sigma =0.075$	$n=11$	$2.0167$	$0.0915$	$2.1369$	$0.0507$	$0.0740$	$0.0317$
	$n=31$	$2.0179$	$0.0955$	$2.1258$	$0.0480$	$0.0741$	$0.0309$
	$n=51$	$2.0193$	$0.0999$	$2.1157$	$0.0478$	$0.0741$	$0.0282$

and

Table 6. Initial solutions for simulated data of the $H_3$ process.

Combinations		${\mathit{\eta }}_0$	${\mathit{\beta }}_0$	${\mathit{\gamma }}_0$	${\mathit{\tau }}_0$	${\mathit{\sigma }}_0^2$
$\sigma =0.01$	$n=11$	$2.0136$	$0.0831$	$2.1381$	$0.0495$	$0.0098$
	$n=31$	$2.0136$	$0.0925$	$2.1391$	$0.0485$	$0.0098$
	$n=51$	$2.0136$	$0.0970$	$2.1351$	$0.0485$	$0.0098$
$\sigma =0.05$	$n=11$	$2.0151$	$0.0866$	$2.1394$	$0.0503$	$0.0489$
	$n=31$	$2.0151$	$0.0933$	$2.1357$	$0.0489$	$0.0489$
	$n=51$	$2.0151$	$0.0973$	$2.1285$	$0.0487$	$0.0489$
$\sigma =0.075$	$n=11$	$2.0178$	$0.0911$	$2.1393$	$0.0508$	$0.0733$
	$n=31$	$2.0178$	$0.0942$	$2.1312$	$0.0491$	$0.0733$
	$n=51$	$2.0178$	$0.0977$	$2.1224$	$0.0489$	$0.0733$

, respectively. According to the former, it can be observed that errors (RAE) follow the same pattern for every value of $\sigma$ and their combination with n. Indeed, the estimated values lead to the same errors for all $\sigma$. The reduction of the error in this model appears when n increases, as expected, but being more noticeable that in the other models. In any case, RAE remains practically the same for all combinations. On the other side, initial values of $\eta$ and $\sigma$ are the same when n is varying. In the case of the $H_3$ model, the search region for initial values ${\beta }_0$ and ${\tau }_0$ was the rectangle $[0.005,0.25]\times [0,0.1]$. The initial value ${\gamma }_0$ was obtained, in every combination, after fitting a Weibull curve to the observed mean by nonlinear least squares (a study of the stability of such approach concerning ${\gamma }_0$ was carried out in [29]). This behavior is illustrated by Figure 5. In particular, the most noticeable differences between observed and estimated mean points appear for low values of $\sigma$. When variability increases the model appears to behave better, showing that it would be suitable to describe complex growth dynamics with a high degree of variability.

A comparison between the results of the three models shows that the lowest RAE is achieved with the hyperbolastic model of type II. The type-I model also yields low RAE values, but mainly for low variability. The type-III model is the one showing the highest RAE for $\sigma =0.01$. Despite such error being higher than in the other two models for high values of $\sigma$, it remains stable and does not vary too much. The first conclusion to be drawn from this comparison is that the type-II model displays the best behavior, but the type-III model seems to be more flexible as it can be more successfully adapted to a variety of scenarios.

6. Conclusions

Associated with each of the hyperbolastic functions ($H_1$, $H_2$, $H_3$), a stochastic model has been presented based on a diffusion process whose mean function is the curve under consideration. A previous analysis of the hyperbolastic curves determines that all of them can be considered as the solution of a linear ordinary differential equation, more specifically of the Malthusian type. This provides a joint vision of all the curves and allows us to extend the deterministic models to stochastic ones, which is done through stochastic differential equations obtained by adding a multiplicative noise to the previously mentioned Malthusian equation. Thanks to this unification of models, global inference procedures can be presented for all of them, which include both obtaining the system of maximum likelihood equations and describing the common procedure by which the initial solutions required to solve said system are obtained. This unified perspective is particularly useful in simplifying the computational implementation of numerical procedures to fit data from real phenomena.

A simulation study has been carried out concerning initial and final solutions for different combinations of variability and sample size. Results show that every model behaves as expected but with subtle differences according to the combination of $\sigma$ and n. As a matter of fact, the type-I model is affected by increasing variability, with this being higher as the sample size increases. The type-II model shows the lower errors and, finally, the error in the case of the type-III model remains the same for every combination of variability and sample size, and its reduction with n becomes more noticeable than in other models.

Future research lines might concern the inclusion of the processes herein described in richer models able to describe more complex behaviors. Indeed, hyperbolastic growth models could be useful in addressing random growth in several other models (e.g., epidemiological models such as SIR or SEIR, or any other describing biological phenomena such as tumor growth or tumor suppression). To this end, it must be noted that the power of hyperbolastic models, as is the case of many others, relies mainly on their number of parameters as well as the flexibility of their intrinsic functions. This implies a growing, data-dependent complexity in the estimation of parameters and the performance of other statistical procedures. This limitation must be acknowledged and considered in order to build reliable models.