Statistical analysis and first-passage-time applications of a lognormal diffusion process with multi-sigmoidal logistic mean.

We consider a lognormal diffusion process having a multisigmoidal logistic mean, useful to model the evolution of a population which reaches the maximum level of the growth after many stages. Referring to the problem of statistical inference, two procedures to find the maximum likelihood estimates of the unknown parameters are described. One is based on the resolution of the system of the critical points of the likelihood function, and the other is on the maximization of the likelihood function with the simulated annealing algorithm. A simulation study to validate the described strategies for finding the estimates is also presented, with a real application to epidemiological data. Special attention is also devoted to the first-passage-time problem of the considered diffusion process through a fixed boundary.

Introduction

Growth curves with sigmoidal behavior are widely used in several applied fields including biology (see, for instance, Brauer and Castillo-Chavez (2012)), software reliability (cf. Erto et al. (2020), Inoue and Yamada (2013)) and economics (see, for example, Smirnov and Wang (2020)). During the times different kinds of sigmoidal curves have been introduced such as logistic, Gompertz, Korf, Bertalanffy, etc. Many efforts have been made essentially for two main purposes: (i) unification of classical models (see Chakraborty et al. (2019)), and (ii) generalizations of growth curves (see, for example, Asadi et al. (2020), Di Crescenzo and Spina (2016) and Romero et al. (2016)).

The differential equations which drive the growth of the aforementioned deterministic models are very useful to describe population dynamics. However, in order to make them more realistic, it is necessary to introduce a noise term in the equation. In this way, the differential equations are replaced by stochastic ones. Most of the times, the analysis of the resulting stochastic equation is quite complex, and the transition probability density of the resulting diffusion process cannot be determined (for example, see Campillo et al. (2018), in which the authors propose, for this reason, a new approach to find the maximum likelihood estimates). Models based on diffusion processes are commonly used in various fields of applications, for example plant dynamics (cf. Rupšys et al. (2020), where a hybrid growth is based on Gompertz and Vasicek models), resources consumption (for instance, Nafidi et al. (2019) use the Brennan-Schwartz process to model electricity consumption in Morocco) or particular fish species growth (cf. a stochastic version of the open-ended logistic model considered in Yoshioka et al. (2019)).

In a recent paper, Di Crescenzo et al. (2021) focus on the generalization of the classical logistic growth model introducing more than one inflection point. To this end, firstly, two different birth-death processes, one with linear birth and death rates and the other with quadratic rates, were considered. Then, a diffusive approximation was performed leading to a non-homogeneous lognormal diffusion process with mean of multi-sigmoidal logistic type. Attention was also given to the description of its main features of interest in applied contexts. For instance, the mean of the process is a generalized version of the classical logistic function (see, for instance, Di Crescenzo and Paraggio (2019)) with more than one inflection point. The transition probability density of the process has been obtained explicitly and has been applied to plant dynamics.

Starting from the theoretical results of the previous works, in the present paper we approach the problem of the inference of the stochastic model. This is done by means of the maximum likelihood method, thanks to the availability in closed form of the likelihood function. We also address the treatment of some collateral problems that emerge in the development carried out, such as: (i) obtaining initial solutions to solve the system of likelihood equations, and (ii) bounding the parametric space for addressing the estimation by metaheuristic procedures. All development is supported by simulation examples. Subsequently, in order to provide an example of application to real phenomena, we adopt the proposed model to describe the behavior of the data on the evolution of COVID-19 in different European countries during the two first waves of infection. Indeed, some of the main features of the diffusion process, such as the mean, the mode and the quantiles, may be used for prediction purposes and they are expressed as a function of the parameters of the process.

The problem of parameters estimation has been considered in several papers, for instance in Shimizu and Iwase (1987) and in Tanaka (1987). See also the more recent works of Garcia (2019), in which the author converts the maximization of the likelihood function into an equivalent problem regarding the minimization of a square error, and of Ramos-Ábalos et al. (2020) where maximum likelihood estimates of the parameters of the powers of the homogeneous Gompertz diffusion process are obtained.

Two different strategies to obtain the maximum likelihood estimates of the parameters are introduced. The first is based on the solution of the system of the critical points of the likelihood function, and the other stems from a meta-heuristic optimization method (simulated annealing) to maximize the likelihood function.

This is the outline of the content of the paper. In Sect. 2, the most relevant characteristics of the deterministic and the corresponding stochastic model are recalled. Then, the problem of finding the maximum likelihood estimates of the involved parameters is described in Sect. 3. In several contexts of population dynamics, it may be relevant to know how long the population spends below a certain control threshold. For this reason the first-passage-time (FPT) problem is also addressed. More precisely, in Sect. 4, the R-package fptdApprox (see Román-Román et al. (2020)) is used to determine the approximated FPT density of the lognormal diffusion process through a constant boundary. With the purpose of validating the described procedures for finding the maximum likelihood estimates, a simulation study is presented in Sect. 5. Finally, in Sect. 6 we propose an application of the model to real data concerning the COVID-19 infections in France, Italy, Spain and United Kingdom.

The multisigmoidal logistic model and the corresponding diffusion process

Consider the classical logistic equationwith . If the intrinsic growth rate r is replaced by a polynomial P(t), then the solution of this equation, with the initial condition , is given bywhere . With the hypotesis that when , the carrying capacity of this generalized model is given by , and thus it is independent from the initial condition . In order to obtain a generalized logistic function in which the carrying capacity is dependent on the initial condition, we consider the following equation (cf. in Di Crescenzo et al. (2021))withfor , with , whereand . Under these assumptions, the solution of the ordinary differential equation (1), with initial condition , is the so-called multisigmoidal logistic function given byWe point out that the function may exhibit more than one inflection point, and its carrying capacity iswhere and is defined in Eq. (3). It is easy to note that the function (4) is not monotonous in general, since the monotonicity intervals depend on the coefficients of the polynomial , and the carrying capacity is the maximum value attainable by the function . See Fig. 1 for some plots of the multisigmoidal logistic function.

Fig. 1

The multisigmoidal logistic function for some choices of the parameters: , , , , a and, from bottom to top, ; b and, from bottom to top,

The investigation of the inflection points in the case of multisigmoidal growth curves are of great interest. Unfortunately, since the expression of function (4) is quite complex, these points cannot be obtained explicitly, but it is possible to provide an equation in the unknown t solved by the inflection points, that isIn Fig. 2, the multisigmoidal logistic function and the corresponding inflection points are shown for some choices of the parameters.

Fig. 2

The multisigmoidal logistic function and the corresponding inflection points for , , , , a and ; b and

The corresponding diffusion process

In Di Crescenzo et al. (2021), a special time-dependent lognormal diffusion process has been considered, with and infinitesimal momentswhere is defined in (2), and . The aforementioned process is determined by the following stochastic differential equation, obtained from Eq. (1) by adding a multiplicative noise term, where means equality in distribution, and where W(t) denotes a Wiener process independent from the (possibly random) initial state , for . We point out that this is not the only way to randomize the growth deterministic equation. Indeed, in the case of random catastrophes, it may be more appropriate to consider as a noise term a Poisson process (see for example Schlomann (2018)). The solution of Eq. (8) iswith The existence and uniqueness of solution of the linear stochastic differential equation (8) is ensured by virtue of the continuity of function (see, for example, Arnold (1974)). Moreover, if either is degenerated at , in the sense that , or follows a lognormal distribution , then the finite dimensional distributions of the process are lognormal. Namely, for any and , the vector follows an n-dimensional lognormal distribution , where the entries of the vector are given byand the components of the matrix are given byFurther, the conditional distribution of the process follows a lognormal distribution, i.e. for From the above mentioned distributions, some characteristics associated to the process can be obtained (cf. Di Crescenzo et al. (2021)). For example, the mean of X(t) conditional on is given byMoreover, if then the mean of X(t) isand the -percentiles for arefor , where is the -percentile of the standard normal random variable. Note that the conditional mean (11) and the mean (12) are multisigmoidal logistic functions of t, in the sense that they solve the multisigmoidal logistic equation (1).

Maximum likelihood estimations

The stochastic model introduced in Sect. 2.1 can be employed in several applications, especially for describing real populations that exhibit a growth pattern with more than one inflection point. Clearly, in order to apply this model to real data, the unknown parameters need to be estimated. In Sect. 2.1 we obtained the distribution of the diffusion process X(t) defined in (9). Now we propose to estimate the parameters by means of the classical maximum likelihood method. The adoption of this strategy is particularly suggested by the availability in closed form of the transition distribution of the process X(t). Hence, we follow the same lines introduced in Román-Román et al. (2018) for general lognormal diffusion processes. We consider a discrete sampling of X(t) based on d independent sample paths, with different observation instants for the i-th sample path, i.e. , , for . For simplicity, assume that the first observation time is identical for any trajectory, i.e. , . Moreover, let the vector contain the variables of the i-th sample path, for , and let . By supposing that follows a one-dimensional lognormal distribution and by considering the transitions of the process X(t), the probability density function of has the following expressionwhere is a vector of dimension , withRecalling (10), the parameters in (14) are given byandwithIn order to obtain a more manageable expression of the density (14), the following change of variables may be considered:Hence, the probability density function of the vector , with , and , isfor , with n defined in (15), where , and , , , with and , for and .

By setting and supposing that and are functionally independent, the log-likelihood function is given bywithandThe maximum likelihood estimations (MLEs) of can be computed easily. Indeed, by differentiating , from (17) we obtainFurther on, in order to find the maximum likelihood estimates of , two different approaches are available: Hereafter, in the Sects. 3.1 and 3.2 we provide a description of the two strategies, whereas in Sect. 5 we present an application to a simulation study that involves the given strategies.

solving the nonlinear system

maximizing the objective function .

The availability of the probability density function of in (14) allows to obtain explicitly the log-likelihood function given in (17). Consequently, following the maximum likelihood estimation procedure, in Sect. 3.1 we obtain the associated system of equations, the final form being reported in Eq. (23) below. However, since such system does not have an explicit solution, its resolution must be obtained by adopting numerical methods.

Solving the nonlinear system

Recalling that , the partial derivatives of are given bywherewithHence, the MLEs are the solutions of the following system of nonlinear equationsBy defining the following quantitieswith , the last equations of the system (19) can be written as followsSubstituting the expression (16) of in the previous equations, one haswhere, for any , one hasHence, until now, the expression of the system solved by the MLEs isThe first equation of system (21) can be further simplified. Indeed, by settingone hasConsequently, the system (21) finally becomesNote that (23) is a system of equations in the unknowns contained in .

Remark 1

For the first equation of the system (23) in the unknown , since , andthe only acceptable solution is

Clearly, since in general system (23) cannot be solved analytically, then a numerical approach is needed. Specifically, we adopt the well-known Newton-Raphson method to solve (23) (for instance, see Dennis and Schnabel (1996)). For such an iterative method, an initial approximation for the solution of the system is needed. It can be obtained by a procedure similar to that used by Román-Román et al. (2019). For the initial solution of the vector , by considering the multisigmoidal logistic function, i.e.it can be supposed, without loss of generality, that (see Remark 2.1 of Di Crescenzo et al. (2021)), so thatThen, considering the sampling defined in Sect. 3, consisting of d independent sample paths of the process X(t), for simplicity we suppose that any sample path of the process has the same number of observations, i.e. for any . However, the following remarks hold even in more general cases. Moreover, let be the values of the mean of the sample paths at the time , for , that iswhere is the value of the i-th sample path at the time .

In general, the carrying capacity is unknown. We suppose that the observations are available over a large time interval, such that the evolution of the population is terminated over such an interval. Hence, the carrying capacity can be approximated with the last value of the sample mean . This approximation can be adopted also in the other cases, since it is used just to construct an initial solution for the parameters of the Newton-Raphson method for the estimate of . Thus, we can consider a polynomial regression for the pairsThe coefficients of the approximating polynomial will be the initial values for the parameters . Thus, the initial solution for is given by .

Finally, in order to construct the initial solution of , let us now recall that for a lognormal distribution , one has , so that the quantity gives an approximation for , where m and are respectively the arithmetic sample mean and the geometric sample mean of a random sample from Y. Hence, one hasSince is estimated by the sample mean m, we haveAs a consequence, in our setting an estimate for is given bywhere and denote respectively the arithmetic and the geometric sample mean of the observations performed at the time . Hence, an initial approximation for can be obtained by performing a simple linear regression of against .

In conclusion, in order to obtain the maximum likelihood estimates of the parameters contained in , the steps of the proposed strategy to solve the system (23) are: The adoption of the above strategy requires to start from good initial solutions for the unknown parameters. Unfortunately, even in this case it is not always possible to guarantee the convergence of this method. For this reason, recently various procedures have been proposed aimed at addressing the maximization of the likelihood function, by viewing this as a direct optimization problem. Indeed, there is a wide range of stochastic metaheuristic methods, which can be classified into two large families: those based on trajectories and those based on swarms. Hereafter, in Sect. 3.2 we employ one of the most widely used, the Simulated Annealing. This method requires necessarily to bound the parametric space, and this matter is the object of Sect. 3.2.2.

finding an initial solution for the parameters and with a polynomial regression of against , for any ;

finding an initial solution for with a simple linear regression of against , with , for any and where can be obtained by means of the second of Eqs. (18);

using the Newton-Raphson method to solve the system (23), with the initial solutions determined at steps (i) and (ii).

Maximizing the log-likelihood function

Let us now illustrate a strategy based on Simulated Annealing (S.A.) and finalized to obtain the MLEs for the parameters of the process (9). We first provide a brief description of this method in Sect. 3.2.1. Then, in Sect. 3.2.2 we describe a suitable criterion to restrict the parametric space, this being essential to apply the S.A. method in the remainder of the paper.

Brief notes on Simulated Annealing

The aim of this section is to determine the MLEs by using the S.A. algorithm. The aforementioned method, introduced in Kirkpatrick et al. (1983), is a meta-heuristic optimization algorithm used for problems like finding . It is considered more suitable with respect to other numerical algorithms since it needs less restrictive conditions regarding the regularity of the domain and the analytical properties of the objective function f. The algorithm works such that in every step a random point is chosen in the solution space. If the new solution is better than the previous one, then the latter is replaced. Otherwise, if the new solution is worse than the previous, then the latter may be replaced with a probability rate which depends on the increase of the objective function and on a suitable scale factor T, that is named ‘temperature’ in agreement with the metallurgical process of annealing that inspired this algorithm. We recall that the S.A. is successful because it avoids local minima. In recent years it has been widely used in the context of estimation in diffusion processes (see, for example Luz Sant’Ana et al. (2018) and Román-Román and Torres-Ruiz (2015)).

In this context, the algorithm works in the following way. It begins with an initial choice for the parameters of interest, then is generated from an uniform distribution in a neighborhood of . Then, a new value of is obtained in such a wayConsequently, if , then and therefore is replaced by . Otherwise, if , then may be accepted anyway with probability . The temperature T is defined in such a way that at the beginning the probability of accepting is high, and during the execution of the algorithm the function T decreases. The initial temperature must be sufficiently large so that the algorithm accept the solutions which let the objective function increases with a large probability . In literature, the choices of the initial parameters are usually and , where denotes the average increase of the objective function in an application test where all the solutions which cause an increase are accepted. The cooling process which defines the temperature T is usually chosen of geometric type, i.e. for . Usually the constant is chosen among 0.8 and 0.99 in order to have a slow cooling procedure. In our case, we set . In any iteration of the algorithm, a chain of L new solutions is obtained, for . As required, the algorithm stops when at least one of the following rules is satisfied: (i) the last L obtained values are equal, (ii) the maximum number of iterations (1000, in our case) is attained, (iii) the final temperature is reached.

Bounding the parametric space

S.A. needs a restriction of the solution space , namely the set which contains the parameters . Until now, this space is continuous and unbounded, sinceWe consider so that the simulated sample paths are less variable around the sample mean, and thus the multisigmoidal logistic profile is advisable. For the parameters , we find the confidence intervals by using the data of the polynomial regression performed previously to find the initial solutions. More in detail, it is known that the carrying capacity of the multisigmoidal logistic model with is (see Eq. (5)). The carrying capacity can be approximated with the last value of the sample mean, whereas the initial value with the first value of the sample mean (24), so that one hasFrom Eq. (25), it easily follows and thus an approximation of isConsidering Eqs. (4) and (5), for one hasso thatHence, by replacing with its estimate , we can use the resulting confidence intervals of the parameters of the polynomial regression as intervals of variation for the parameter of the diffusion process. We adopt a confidence level equal to 0.999, to attain a high probability that the true parameter belong to the computed intervals.

In order to approximate the range of variation of , from Eq. (25) we have that the last value of the i-th sample path satisfieswhere with and are the sample data. Hence, for the range of variation of one has , whereIn conclusion, the following bounded intervals are employed:

for we consider the confidence intervals of the coefficients of the polynomial regression of against , for , where ,

for we consider the interval , with a and b defined in (26),

for we consider the interval .

Asymptotic distribution of the MLEs

On the ground of the results given in Sect. 5 of Román-Román et al. (2018), in this section we aim to determine the asymptotic distribution of the MLEs (i) of the parameters of the initial distribution, and (ii) of the parameters of the process. can be expressed aswhere is given bywithand is defined in the third of Eqs. (20). Moreover, is defined asFinally, is given in the second of Eqs. (22). We point out that the matrix (27) will be used in Sect. 5 to determine the asymptotic variances for the estimates of the parameters and an approximation of the confidence intervals. Indeed, by applying the delta method (cf. Oehlert (1992)), any q-parametric function with asymptotically has a q-dimensional normal distribution, i.e. (cf. Román-Román et al. (2018))where is the vector of partial derivatives of with respect to .

The exact distribution of is normal , whereas the exact distribution of is chi-square , cf. Román-Román et al. (2018).

The asymptotic distribution of is a -dimensional normal distribution with mean and covariance matrix , i.e. , where denotes Fisher’s information matrix of . For the diffusion process X(t) with a multisigmoidal logistic mean,

In the following section we address a relevant problem for the applications, namely the FPT problem of the diffusion process X(t) through a continuous boundary. Subsequently, in Sect. 5 we adopt a simulation-based approach as the basis of both computational methods described so far, namely the Newton-Raphson method and the S.A. method. The estimates of the parameters obtained through these methods are then used to perform inference on the FPT density.

First-passage-time problem

The FPT problem of a stochastic process X(t) through a boundary S(t) is a problem of great interest in many fields of application, such as medicine, biology or mathematical finance, since the threshold S(t) may represent a critical value of the modeled population size. Considering a stochastic process , the FPT of the process X(t) through the continuous boundary S(t), given , is defined as the following random variableFinding the expression of the distribution of the variable T is hard in general. However, in literature there are several studies for particular types of processes, for example diffusion processes. It has been shown that if S(t) is a continuous and differentiable function, then the density of T, denoted by , solves a II-kind Volterra equation (cf. Eq. (2.4) of Buonocore et al. (1987)). The aforementioned Volterra equation has an explicit solution only for certain special boundaries (see for example Sects. 2.3 and 4.3 of Giorno and Nobile (2019) in which the FPT density through special boudaries has been obtained for the restricted Gompertz-type diffusion processes). In certain instances, it is appropriate to adopt numerical procedures in order to approximate its solution. To this aim Buonocore et al. (1987) proposed a simple but efficient algorithm, based on the composite trapezoidal formula. More in detail, Theorem 4 of Buonocore et al. (1987) proves the convergence of the approximated FPT density to the theoretical one. However, the application of the proposed numerical procedure requires (i) the choice of a suitable step h of integration which ensures a good approximation of the real solution, (ii) the choice of an initial time instant and (iii) the choice of the final time instant . Román-Román et al. (2008) studied the problems related to the practical application of the numerical procedure. The first problem is linked with a suitable choice of h. Indeed, taking into account the result of Theorem 4 of Buonocore et al. (1987), it is easy to note that the convergence is ensured when . Consequently, the value of h should be small enough, but sufficiently far from 0. Indeed, if h is excessively small, then the computational cost may increase in vain because, with a larger integration step, a similar approximation may be obtained with a smaller number of iterations. On the other hand, if h is excessively large, the approximation may be unsatisfactory. These problems depend on the localization of the FPT T, and may be solved if the range of variation of T is known. For this reason, Román-Román et al. (2008) introduced a function, called ‘FPT location’ (FPTL), finalized to obtain, from a heuristic point of view, the range of variation of T. Specifically, the FPTL function is defined as followswhere is the transition distribution of the process X(t). Referring to the diffusion process with infinitesimal moments given by Eq. (7), and by considering a fixed and constant boundary , given , the FPTL function for the process X(t) is given bywhere is the standard normal distribution andThe information provided by the FPTL function is relevant for an efficient application of the algorithm proposed by Buonocore et al. (1987). Indeed, thanks to the FPTL function, an adaptive step of integration can be obtained. In this way, the execution time of the algorithm is reduced.

Example 1

Let X(t) be a diffusion process with infinitesimal moments (7), with , , , , , and a.s. See Fig. 3 for the plot of 100 simulated sample paths of the process. Let us study the FPT density through the fixed boundary , by using the information provided by the FPTL function and the R package fptdApprox (for references, see Román-Román et al. (2012), Román-Román et al. (2014) and Román-Román et al. (2020)).

Fig. 3

100 simulated sample paths of the process X(t) with , , , , and . The black line represents the sample mean of the process, while the red line represents the boundary

Figure 4a shows the FPTL function (obtained by means the function FPTL of the package fptdApprox), whereas the approximated FPT density (obtained using the package fptdApprox) is plotted in Fig. 4b. Other useful quantities related to the FPT density are given in Table 1.

Fig. 4

a The FPTL function and b the approximated FPT density of the process X(t) through the constant boundary , for the same assumptions of Fig. 3

Table 1

The mean, the standard deviation, the mode, the first, the fifth and the ninth decile of the FPT of the process X(t) through the boundary

Mean	St. dev.	Mode	1st decile	5th decile	9th decile
40.18765	1.568392	39.92321	39.02346	40.11264	41.58065

Simulation

In Sect. 3, two procedures have been introduced to obtain the MLEs of the parameters involved in the diffusion process (9). The former procedure is based on the numerical resolution of a system of nonlinear equations, whereas the latter is based on the application of S.A. algorithm. In this section, a simulation study is developed to verify the validity of the two aforementioned procedures. We consider the diffusion process X(t) with infinitesimal moments (7), for , and , , , and . These choices of the parameters are performed arbitrarily, to obtain different patterns of the growth curve. For example, the choice , , , refers to the case of a non monotonous multisigmoidal logistic function, whereas the choice , , and to the case of an increasing multisigmoidal logistic curve (see Fig. 5). To estimate the parameters in , we consider the 32 combinations of the values of the parameters listed in Table 2, with in every case. For each case, we simulate 200 sample paths of X(t), by generating 501 simulated points at equidistant times for .

Fig. 5

100 simulated sample paths of the diffusion process for , and a , , and b , , (simulation study)

Table 2

The values of the parameters (simulation study)

Case no.	β_1	β_2	β_3	η	σ	Case no.	β_1	β_2	β_3	η	σ
1	0.1	-0.009	0.0002	e^-1	0.01	17	0.5	-0.009	0.0002	e^-1	0.01
2	^′′	^′′	^′′	^′′	0.05	18	^′′	^′′	^′′	^′′	0.05
3	^′′	^′′	^′′	e^-3	0.01	19	^′′	^′′	^′′	e^-3	0.01
4	^′′	^′′	^′′	^′′	0.05	20	^′′	^′′	^′′	^′′	0.05
5	^′′	^′′	0.0003	e^-1	0.01	21	^′′	^′′	0.0003	e^-1	0.01
6	^′′	^′′	^′′	^′′	0.05	22	^′′	^′′	^′′	^′′	0.05
7	^′′	^′′	^′′	e^-3	0.01	23	^′′	^′′	^′′	e^-3	0.01
8	^′′	^′′	^′′	^′′	0.05	24	^′′	^′′	^′′	^′′	0.05
9	^′′	-0.007	0.0002	e^-1	0.01	25	^′′	-0.007	0.0002	e^-1	0.01
10	^′′	^′′	^′′	^′′	0.05	26	^′′	^′′	^′′	^′′	0.05
11	^′′	^′′	^′′	e^-3	0.01	27	^′′	^′′	^′′	e^-3	0.01
12	^′′	^′′	^′′	^′′	0.05	28	^′′	^′′	^′′	^′′	0.05
13	^′′	^′′	0.0003	e^-1	0.01	29	^′′	^′′	0.0003	e^-1	0.01
14	^′′	^′′	^′′	^′′	0.05	30	^′′	^′′	^′′	^′′	0.05
15	^′′	^′′	^′′	e^-3	0.01	31	^′′	^′′	^′′	e^-3	0.01
16	^′′	^′′	^′′	^′′	0.05	32	^′′	^′′	^′′	^′′	0.05

The remainder of this section is organized as follows: (a) since the degree of the polynomial is unknown a priori, we propose the use of the strategy described in Román-Román et al. (2019), by increasing the degree until the goodness of fit is optimal; (b) considering the degree obtained at the step (a), we use the two procedures described in Sects. 3.1 and 3.2 to find the MLEs of the parameters.

The choice of the best degree of the polynomial is performed under the goodness of fit criteria based on the four following measures:

the absolute relative error (RAE) between the sample mean and the estimated mean, i.e. where denotes the mean of the estimated process considering a polynomial of degree p;

the Akaike information criterion (AIC), which is defined as

the Bayesian information criterion (BIC), which is given by where n represents the number of observations,

the resistor-average distance () between the sample distribution and the p-th estimated distribution , for , which is defined as the following harmonic mean (cf. Johnson and Sinanovic (2001)): where denotes the Kullback-Leibler divergence. Assuming that the sample distribution is lognormal with parameters and that the estimated distribution is lognormal with parameters the Kullback-Leibler divergence between the sample distribution and the p-th estimated distribution for is given by, for any with defined in (10). Clearly, if the theoretical distribution of the process is known, one can alternatively compute the resistor-average distance between the theoretical and the estimated distribution. We consider the expected distance and the median of the distance as reference values for the resistor-average distance.

In cases (ii) and (iii), the stochastic model is characterized by parameters. Moreover, is defined in (17), and and are the MLEs of the parameters and . The best fit is attained for the smallest value of the considered goodness measures. Table 3 shows the estimated parameters for the case no. 1 of Table 2, which is obtained by solving the system (23) for different degrees of the polynomial . Furthermore, the results about the goodness of measures are given in Table 4 and in Fig. 6. It can be noticed that the estimated parameters for and are almost identical, and that is very close to zero in the case . Hence, the results concerning the measures of goodness obtained in these two cases are quite similar. This conclusion is also confirmed by the analysis of the RAE measures (in Table 4), that are often used to measure the fit error of the model in terms of the fit of the mean function. The analysis is performed in terms of the scale of judgment of the model accuracy based on the Mean Absolute Percentage Error (MAPE), cf. Klimberg et al. (2010) and Lewis (1982). Indeed, the judgment suggested by the MAPE shows that and are referred as highly accurate, whereas is evaluated as good forecast, with both and considered as reasonable forecast. Consequently, the choice is taken as the best, since it involves the lowest number of parameters.

Table 3

The estimated parameters obtained by solving system (23) for different degrees of the polynomial (simulation study)

	Value	β_1	β_2	β_3	β_4
p=2	Initial	-0.19227458	0.006134021	–	–
	Estimated	-0.02167632	0.001417851	–	–
p=3	Initial	-0.02332665	-0.005246427	0.0001724780	–
	Estimated	0.10200076	-0.009115730	0.0002019617	–
p=4	Initial	-0.16488586	0.011922169	-0.0004348227	6.54664e-03
	Estimated	0.10313862	-0.009310721	0.0002108884	-1.240487e-07
p=5	Initial	-0.23209031	0.024605250	-0.0012041705	2.521331e-05
	Estimated	-0.23209031	0.024605250	-0.0012041705	2.521331e-05
p=6	Initial	-0.36570406	0.060641889	-0.0044843917	1.579124e-04
	Estimated	-0.36570406	0.060641889	-0.0044843917	1.579124e-04
p=2	Initial	–	–	0.3688831	0.0092356136
	Estimated	–	–	0.2342132	0.0013732663
p=3	Initial	–	–	0.3688831	0.0092356136
	Estimated	–	–	0.3649943	0.000103298
p=4	Initial	–	–	0.3688831	0.0092356136
	Estimated	–	–	0.3649325	0.0001032229
p=5	Initial	-1.572548e-07	–	0.3688831	0.0092356136
	estimated	-1.572548e-07	–	0.3688831	0.0000859656
p=6	Initial	-2.617632e-06	1.706174e-08	0.3688831	0.0092356136
	Estimated	-2.617632e-06	1.706174e-08	0.3688831	0.0000859656

Table 4

The goodness measures for different degree

Degree	RAE	BIC	AIC	Median of D_RA	Mean of D_RA
2	0.163198755	-6986.46	-7006.983	1.076721237	1.874847040
3	0.001667815	-10261.07	-10286.721	0.001035228	0.001668560
4	0.001390662	-10254.13	-10284.916	0.001149273	0.001916131
5	0.413978453	34449.77	34449.77	54.690728902	84.305882377
6	0.416847669	47306.94	47306.94	69.316641926	94.836470468

For the resistor-average distance , the estimated and the theoretical distributions are considered (simulation study)

Fig. 6

The resistor average distance between a–b the sample and the estimated distribution, and c–d the theoretical and estimated distribution for the case 1 of Table 2, for different degrees of the polynomial (simulation study)

The same result can be obtained for the other parameters choices, but it is omitted for brevity. Here, we limit to mention that the AIC and its Bayesian version, the BIC, provide a global measure of the adjustment to the model in terms of the likelihood that the model itself gives to the observed sample, so that these measures also allow for model selection criteria. The AIC and the BIC are seen often as complementary measures to the use of the Resistor Average Distance between the sample and the theoretical distributions of the model. However, there is no criterion that indicates that one measure is better than another, and thus in general the use of several alternative measures is recommended, as usual in practical applications. In our analysis, the coincidence of the conclusions suggested by these measures supports the final decision. Hence, from now on, a polynomial of degree will be considered.

Table 5 shows the estimated values of the parameters obtained by solving the nonlinear system (23) by means of the Newton-Raphson method. These values provide good parameters estimates, especially when is small. The last column of the Table 5 contains the RAE. In this case, it is defined aswhere are the values of the sample mean and are the values of the estimated mean at the time considering a polynomial of degree . For a comparison between or and the RAE, see Fig. 7a–b: it can be noticed that the value of the RAE shows an increasing trend with respect to the parameter , whereas it shows a constant trend with respect to . In Fig. 8a–b the theoretical, sample and estimated sample means for the parameters choices number 1 and 2 of the Table 2 are shown. Clearly, the best estimation is obtained when is small.

Table 5

The estimated values of the parameters obtained by solving the system (simulation study)

Case no.	β^_1	β^_2	β^_3	η^	σ^	RAE
1	0.1008708	-0.009083141	0.0002016053	0.3748345	0.009948509	0.009948509
2	0.09647959	-0.009128214	0.0002035713	0.4174093	0.04861361	0.00968444
3	0.1006515	-0.009032199	0.0002004328	0.04941139	0.01006651	0.001394514
4	0.09813563	-0.008927315	0.0001989717	0.05311648	0.05045509	0.004833166
5	0.09917466	-0.008943618	0.0002994268	0.3684543	0.009948515	0.0007958334
6	0.1014043	-0.008998974	0.0002981807	0.3892023	0.05016611	0.004735055
7	0.1000035	-0.008977959	0.000299122	0.0499255	0.01019311	0.001671998
8	0.0998637	-0.009123996	0.0003010265	0.05329735	0.04990423	0.006741047
9	0.0996273	-0.006921284	0.000198009	0.3691476	0.01007512	0.001185446
10	0.09453137	-0.006522038	0.0001887209	0.3961056	0.05014104	0.008249555
11	0.0994893	-0.006961563	0.000199408	0.05006629	0.01005502	0.001227369
12	0.09991228	-0.007300655	0.0002073627	0.05440441	0.04977635	0.008537779
13	0.09871006	-0.006934799	0.0002995149	0.3732344	0.009929474	0.001391576
14	0.1089433	-0.00769021	0.0003170773	0.3827415	0.05020807	0.007503426
15	0.09915895	-0.006970525	0.0003002155	0.05005078	0.009887142	0.001247189
16	0.099621	-0.007308948	0.0003089962	0.05292541	0.05039523	0.01051539
17	0.5027283	-0.009975521	0.0002718814	0.3675535	0.01005231	0.001725124
18	0.4994427	-0.0113279	0.0008357429	0.3763182	0.04958661	0.01836523
19	0.5009121	-0.009399375	0.0002318097	0.04989883	0.01014064	0.001231152
20	0.4889807	-0.007392059	0.0001686769	0.05133072	0.05002574	0.02180217
21	0.4966796	-0.005064475	9.070535e-05	0.3645112	0.008760448	0.00882512
22	0.5169986	-0.008982219	0.0007673324	0.3842942	0.05027921	0.01926412
23	0.5007486	-0.00918498	0.0003122469	0.05015442	0.009886917	0.001718233
24	0.5004829	-0.009757746	0.0004111398	0.05113107	0.0501877	0.03498688
25	0.4992278	-0.007530001	0.0002736762	0.367121	0.01017701	0.003303379
26	0.522767	-0.01616745	0.001661617	0.3978223	0.0504873	0.02065117
27	0.4947951	-0.004930811	-3.751225e-05	0.04872577	0.0100036	0.01431122
28	0.4973394	-0.004424543	-5.398909e-05	0.0509033	0.04989332	0.008392656
29	0.5007914	-0.005683407	0.0001988298	0.3686014	0.01002725	0.001054943
30	0.6589188	-0.01472611	7.767629e-05	0.3884501	0.05398134	0.02949292
31	0.5005393	-0.007046342	0.0003130846	0.04990925	0.009948212	0.002032919
32	0.4944903	-0.006618593	0.0004977852	0.0531368	0.05019473	0.01046167

Fig. 7

The RAE for a and , b and and c , and with respect of the number of replications. In all the cases (simulation study)

Fig. 8

The theoretical, sample and estimated means of the process X(t) for the parameters of the cases number 1 and 2 of Table 2 (from left to right). The results are obtained via Newton-Raphson method in (a) and (b) and via S.A. in (c) and (d). (Simulation study)

Further on, the estimated values obtained via S.A. are given in Table 6, whose last column contains the value of the RAE defined in Eq. (28). Since S.A. is a heuristic algorithm, the MLEs have been computed as the average of the results obtained by 10 uses of the procedure. Figure 8c–d provide the theoretical, the sample and the estimated (via S.A.) means for the cases no. 1 and no. 2 of the Table 2. In addition, in Fig. 7c, the trend of the RAE is plotted as a function of the number of replications: clearly, the goodness of the results improves as the number of replications increases.

Table 6

The estimated values of the parameters obtained via S.A. (simulation study)

Case no.	β^_1	β^_2	β^_3	η^	σ^^2	RAE
1	0.101561	-0.009082552	0.0002018522	0.3625206	9.989558e-05	0.01145807
2	0.09954966	-0.009144748	0.0002033654	0.4041588	0.002549043	0.04471631
3	0.1131058	-0.00993448	0.0002154949	0.05040806	0.0001547924	0.02141325
4	0.1091467	-0.009814184	0.0002144671	0.05347671	0.002515041	0.0371569
5	0.1180621	-0.01109313	0.0003614234	0.3760972	0.0001248723	0.01352716
6	0.1278859	-0.01283086	0.0004142981	0.4122873	0.00251747	0.01633801
7	0.1060342	-0.009544155	0.00031209	0.0501014	0.0001160831	0.007163973
8	0.1233499	-0.01139476	0.000355587	0.05196034	0.002510716	0.02342116
9	0.05999024	-0.003213622	0.0001156288	0.3619285	0.00019005	0.02598964
10	0.1281187	-0.01053716	0.000303384	0.4643666	0.002569706	0.04634721
11	0.07833233	-0.005157051	0.0001608136	0.05258961	0.0001810726	0.04867231
12	0.118225	-0.008582805	0.0002342187	0.05132632	0.002638776	0.035553561
13	0.1177451	-0.00947922	0.0003863382	0.3979672	0.0002070029	0.02261541
14	0.1155341	-0.009715339	0.0004058845	0.4262232	0.002535962	0.01010711
15	0.09977635	-0.007003481	0.0003003812	0.04986716	0.0001068239	0.001787891
16	0.1240548	-0.009570475	0.0003701909	0.05082801	0.002543112	0.02769778
17	0.4996174	-0.007431545	1.808038e-05	0.3697584	9.898729e-05	0.003974179
18	0.5116561	-0.01727471	0.0002780168	0.3700344	0.002477853	0.01019564
19	0.4908564	-0.007392949	5.594539e-05	0.04968137	0.0001481544	0.008720396
20	0.5013114	-0.00809688	0.0001788676	0.0530825	0.00273835	0.009876328
21	0.6249235	-0.02166049	0.000398154	0.4112333	0.000224436	0.06819835
22	0.5686069	-0.02374947	0.000553332	0.3847588	0.002566421	0.01345391
23	0.4865293	-0.006591622	5.45598e-05	0.05038445	0.0001598932	0.01251157
24	0.4833729	-0.00633273	0.000290948	0.05788054	0.002589166	0.03791299
25	0.4994063	-0.008634032	9.26605e-05	0.3633383	0.0001204934	0.003431125
26	0.4987191	-0.004809727	6.356344e-05	0.3925613	0.002584155	0.02507006
27	0.5200495	-0.01434754	0.001020145	0.05181291	0.0001540365	0.02603051
28	0.4965627	-0.004536852	0.0001049973	0.05176822	0.002566167	0.01328192
29	0.4859707	-0.005580455	9.861833e-05	0.3638371	0.0001267092	0.006797028
31	0.4917976	-0.004818751	0.0001038265	0.04923269	0.0001633201	0.01411135
32	0.5002401	-0.008964678	0.0006969377	0.05423012	0.002528847	0.02389477

Moreover, Table 7 contains the estimated values of the parameters (obtained by solving the system (23)), as well as their real values and the asymptotic estimation error. Finally, Table 7 provides various confidence intervals obtained by applying the delta method and using the distribution given in Sect. 3.3 for the case no. 1 of Table 2.

Table 7

The estimated values of the parameters (obtained by solving the system (23)), their real values, their asymptotic estimation error and their , and confidence intervals (simulation study)

Parametric function	η	β_1	β_2
Estimated value	0.3748345	0.1008708	-0.0090831
Real value	0.3678794	0.1000000	-0.0090000
Standard error	1.218018e-03	5.355428e-05	2.683426e-07
95% confidence interval	(0.3064315, 0.4432375)	(0.0865276, 0.1152140)	(-0.0100984,-0.0099352)
90% confidence interval	(0.3174289, 0.4322401)	(0.0888336, 0.1129080)	(-0.0099352,-0.0082311)
75% confidence interval	(0.3346872, 0.4149818)	(0.0925245, 0.1092891)	(-0.0096790,-0.0084872)

Approximation of FPT density

In this section, the FPT problem is analyzed. With reference to a diffusion process X(t) with a multisigmoidal logistic mean and , and , we construct 50 simulated sample paths (see Fig. 9a), each one being formed by 361 data simulating for , . As in Sect. 5, we first chose the optimal polynomial degree (which corresponds to the best fit), and then we found the MLEs of the parameters by solving the system (23). Further on, the R package fptdApprox is used to approximate the FPT density of the process through a constant threshold . Table 8 provides the estimated parameters, whereas Fig. 9b shows the theoretical, the sample and the estimated means, for .

Fig. 9

a 50 simulated sample paths of the diffusion process X(t) for , and . b Theoretical, sample and estimated means of the process X(t) for the FPT density approximation (simulation study—FPT problem)

Table 8

The estimated values of the parameters considering for the FPT density approximation (simulation study - FPT problem)

Degree	η	β_1	β_2	β_3
p=2	2.5287214	-0.09696011	0.002712063	–
p=3	0.3425411	0.09735434	-0.008730983	0.0001936518
p=4	0.9005909	0.08465459	-0.003953699	-0.0001821386
p=5	0.9005909	0.12376611	-0.014259773	0.0006907176
p=6	0.9005909	0.04089375	0.016958673	-0.0032783145

Table 9 provides the four goodness measures for the considered degrees p of the polynomial . Figure 10 shows the resistor-average distances between the theoretical and the estimated distributions, and also between the sample and the estimated distributions. From the given results it follows that the best degree is .

Table 9

The RAE, the AIC, the BIC the median and the mean of the resistor-average distance of the parameters for

Measure of goodness	p=2	p=3	p=4	p=5	p=6
RAE	0.313009817	0.001985486	0.046745173	0.044842133	0.046064145
AIC	-9288.882	-14314.158	-10687.536	-11569.711	-12001.874
BIC	-9267.012	-14286.821	-10654.731	-11531.439	-11958.135
Median of D_RA	3.767780356	0.003568808	0.40969122	0.347181473	0.404933011
Mean of D_RA	4.276620795	0.009886785	0.617753008	0.442544796	0.582836790

For the resistor-average distance, the estimated and the theoretical distributions are considered (simulation study—FPT problem)

Fig. 10

Resistor-average distance between a the theoretical and the estimated distributions and b between the sample and the estimated distributions for the FPT density approximation (simulation study - FPT problem)

Using the estimated model obtained so far, we now focus on the approximation of the FPT density through the boundary . Figure 11 shows the approximated FPT density and the FPTL function realized with the package fptdApprox. Finally, in Table 10 other useful quantities related to the FPT density are provided. It is worth noting that the results obtained in this section are in agreement with those given in Example 4.1.

Fig. 11

The approximated FPT density and the FPTL function of the process X(t) through the boundary (simulation study—FPT problem)

Table 10

The mean, the standard deviation, the mode, the first, the fifth and the ninth decile of the FPT of the process X(t) through the boundary (simulation study—FPT problem)

Mean	St. dev.	Mode	1st decile	5th decile	9th decile
39.88883	1.034443	39.70804	38.83282	39.83612	41.0529

Application to real data

Multisigmoidal functions are suitable to model several special growth phenomena in which the carrying capacity is reached after various stages. In any of these stages a linear growth trend is followed by an explosion of exponential type which finally flattens to a specific value. A growth of this kind is typical of some fruit species, such as peaches or coffee berries (see, for instance the application given in Sect. 3 of Di Crescenzo et al. (2021)). But also some population diseases follow an expansion with a multisigmoidal trend. In this section we apply the considered stochastic model to data concerning the COVID-19 infections in four different European countries, taken from Worldometers (2020). This is just an example finalized to show an application of the multisigmoidal logistic model, without taking into account specific more sophisticated models that describe epidemiological phenomena with greater precision. First of all, we note that the trend of infections in France, Italy, Spain and United Kingdom is similar (see Fig. 12a).

Fig. 12

a Number of infections in France, Italy, Spain and United Kingdom, the black line represents the sample mean. b Sample and the estimated means obtained by solving the system (23) (real application)

This suggests to view these data as different trajectories of the diffusion process X(t) defined on , having a multisigmoidal logistic mean (cf. Sect. 2.1). Hence, in order to find the MLEs of the parameters, we apply the procedure described in Sect. 3.1. For each country, the initial time corresponds to the 30-th day after the one in which the number of infections exceeded 100 (March 30th for France, March 24th for Italy, March 21st for Spain, April 5th for UK), and the final time is chosen as . For any path, the data are scaled as divided by their maximum value, so to be interpreted as a percentage of the last and therefore the maximum value of the growth curve. The estimated means obtained for different degrees are plotted in Fig. 12b. Table 11 provides the initial and the estimated values of the parameters, whereas Table 12 shows the four measures of goodness, for different degrees of the polynomial. Regarding the RAE, every time the degree increases, the approximation improves, whereas the AIC, the BIC and resistor-average distance show that the best choice is .

Table 11

The initial and the estimated values of the parameters considering different degrees

	Value	β_1	β_2	β_3	β_4
p=2	Initial	0.00644346	5.375995e-05	–	–
	Estimated	0.05356653	0.0011042774	–	–
p=3	Initial	0.04480694	-4.587843e-04	1.540731e-06	–
	Estimated	0.04774851	-0.0004685118	1.506227e-06	–
p=4	Initial	0.04090832	-3.650280e-4	8.831071e-07	1.405770e-09
	Estimated	0.04090832	-0.0003650280	8.831071e-07	1.405770e-09
p=5	Initial	0.06509648	-1.269895e-03	1.176365e-05	-5.092724e-08
	Estimated	0.06509648	-0.0012698948	1.176365e-05	-5.092724e-08
p=6	Initial	0.05891632	-9.396086e-04	5.806280e-06	-3.170420e-09
	Estimated	0.05891632	-0.0009396086	5.806280e-06	-3.170420e-09

The results have been obtained by solving the system (real application)

Table 12

RAE, AIC, BIC, the median and the mean of the resistor-average distance considering different degrees

p	2	3	4	5	6
RAE	0.39085660	0.10404805	0.06853298	0.02087067	0.02030226
AIC	-7225.192	-8175.819	-7582.180	-7677.528	-7664.849
BIC	-7205.561	-8151.281	-7552.734	-7643.173	-7625.587
Median of D_RA	0.2250753	0.1301707	0.1548300	0.1615877	0.1642668
Mean of D_RA	1.0586904	0.2741844	0.3218632	0.3103206	0.3106608

For the resistor-average distance, the estimated and the sample distributions are considered (real application)

Regarding the last measure of goodness, see also Fig. 13a, in which the resistor-average distances between the sample and the estimated distributions are provided. Hence, in view of the results obtained for the measures of goodness, the degree is considered.

Fig. 13

a The resistor-average distances between the sample and the estimated distributions considering different degrees of the polynomial. b The -percentiles of the estimated diffusion process X(t) obtained for a degree and for (real application)

Table 13 shows the estimated values of the parameters, the estimation of their standard error and the , and percentiles.

Table 13

The estimated values, the standard error and the , and confidence intervals for the parameters, considering (real application)

Parametric function	η	β_1	β_2
Estimated value	0.03605835	0.04774851	-0.0004685118
Standard error	6.086478e-03	2.342236e-03	2.825993e-05
95% confidence interval	(0.0241291, 0.0479876)	(0.0431578, 0.0523392)	(-0.0005239,-0.0004131)
90% confidence interval	(0.0260470, 0.0460697)	(0.0438959, 0.0516011)	(-0.0005150,-0.0004220)
75% confidence interval	(0.0290568, 0.0430599)	(0.0450541, 0.0504429)	(-0.0005010,-0.0004360)

Moreover, the -percentiles (13) of the estimated diffusion process with are provided in Fig. 13b with . Let us now consider a restricted time range from to , in order to predict the trend of the growth curve in a short-term prediction analysis. Indeed, forecasting the number of infections during a disease in progress is interesting also in the case of short terms, especially for the goodness of estimation (better in this case than in the long term analysis) and for the timeliness of the results. The considered procedure is the same of the one used above, so (i) the best degree p for the polynomial is chosen by considering various measures of goodness, and (ii) the estimated values of the parameters are used to construct a diffusion process X(t) defined on . The estimated values of the parameters are given in Table 14, the values of the four measures of goodness are given in Table 15, and finally in Fig. 14a we provide the resistor-average distances between the restricted sample and the estimated distributions. See also Fig. 14b for the plots of the estimated means for different degrees of the polynomial.

Table 14

The estimated values of the parameters considering (real application with )

Degree	η	β_1	β_2	β_3
p=2	0.22079965	0.04165889	-0.0002129137	–
p=3	0.03650929	0.04792200	-0.0004714958	1.518106e-06
p=4	0.04101207	0.03985192	-0.0003277795	5.314482e-07
p=5	0.04101207	0.06778144	-0.0013896339	1.350773e-05
p=6	0.04101207	0.05704750	-0.0008066323	2.820787e-06

Table 15

The RAE, the AIC, the BIC, the median and the mean of the resistor-average distance of the parameters for

Measure of goodness	p=2	p=3	p=4	p=5	p=6
RAE	0.41031587	0.10421156	0.07598928	0.02665995	0.02500666
AIC	-2722.134	-8035.566	-7476.246	-7591.347	-7571.929
BIC	-2702.568	-8011.108	-7446.896	-7557.105	-7532.796
Median of D_RA	0.2275183	0.1294749	0.1557215	0.1593042	0.1626228
Mean of D_RA	6.2370324	0.2212618	0.2666561	0.2532703	0.2536295

For the resistor-average distance, the estimanted and the sample distributions are considered (real application with )

Fig. 14

a Resistor-average distances between the sample and the estimated distributions in the restricted time range. b Sample and the estimated means obtained by solving the system (23) in the restricted time range (real application with )

Also in this case, the RAE does not provide a good measure of goodness, since every time the degree increases, the approximation improves. From the remaining results, the best choice is , which corresponds to the lowest value of the BIC and of the AIC, and to the lowest resistor-average distance. Hence, considering and the corresponding estimated values of the parameters, Fig. 15a provides the sample and the estimated means in the complete time range, i.e. in . The relative errors between the values of the sample and the estimated means are given in Table 16; note that in all cases they are less than .

Fig. 15

a The sample and the forecasted means in the complete time range considering a degree . b Resistor-average distances between the sample and the estimated distributions in the restricted time range . c Approximated FPT density and d FPTL function in the restricted time range through the boundary (real application—FPT problem)

Table 16

The sample, the forecasted mean and relative error for the last 4 time instants of the complete range , considering a degree (real application)

Time	247	248	249	250
Sample mean	0.97179233	0.98096101	0.98949697	1.00000000
Forecasted mean	1.00554923	1.01157710	1.01720499	1.02244759
Relative error	0.03473674	0.0312103	0.02800212	0.02244759

This section is devoted to the FPT problem. Considering an initial portion of the available data and setting a constant boundary, an estimate of the FPT density is constructed using the numerical procedures recalled in Sect. 4. The resulting FPT density is then compared to the approximated FPT density obtained by using the whole data set.

More in detail, we consider only the first 220 data of COVID-19 infections in the restricted time range and we investigate the best model to fit them. The choice of the optimal degree p of the polynomial is based on the measures of goodness (i)–(iv) described in Sect. 5. The estimated parameters (given in Table 17) are obtained by solving the system (23).

Table 17

The estimated values of the parameters considering in the restricted time range (real application—FPT problem)

Degree	η	β_1	β_2	β_3
p=2	0.09416137	0.02514812	-4.247417e-05	–
p=3	0.04858678	0.05004022	-5.125970e-04	1.693253e-06
p=4	0.06278231	0.04592707	-4.462889e-04	1.183743e-06
p=5	0.06278231	0.07198803	-1.559548e-03	1.646946e-05
p=6	0.06278231	0.05899141	-7.664128e-04	1.337821e-07

By comparing the results given in Table 18 and in Fig. 15b, we choose as the optimal degree.

Table 18

The RAE, the AIC, the BIC, the median and the mean of the resistor-average distance of the parameters for

Measure of goodness	p=2	p=3	p=4	p=5	p=6
RAE	0.28460545	0.08631377	0.06598901	0.02791701	0.03225994
AIC	-6330.849	-7150.963	-6876.159	-6999.085	-6971.142
BIC	-6311.748	-7127.086	-6847.507	-6965.657	-6932.939
Median of D_RA	0.2448219	0.1144287	0.1238069	0.1284333	0.1355944
Mean of D_RA	0.2635491	0.1260033	0.1405463	0.1325032	0.1328700

For the resistor-average distance, the estimated and the sample distributions are considered (real application—FPT problem)

Then, we fix a constant threshold which corresponds to the of the last and maximum data in the complete time range and we use the R package fptdApprox to obtain an estimation of the FPT density. The choice of the constant boundary is not random. Indeed, it is worth observing that the descendent inflection points correspond to the peaks of the function representing the daily increments of the infections. More in detail, by means of Eq. (6), the function representing the sample mean of the infections shows two descendent inflection points, one at the time and the other at the time . The population sizes corresponding to the inflection time instants are and . In the time interval , the mean function has a logistic trend, hence the FPT problem through the boundary is beyond the scope of the present work. Instead, since the mean in the time interval has a multisigmoidal logistic profile, we focus our attention to the FPT problem through the threshold .

The approximated FPT density and the FPTL function of the estimated process X(t) through the boundary are plotted in Fig. 15c–d. In order to validate the predicted results concerning the FPT, we consider also the same problem in the complete time range . The forecasted results for the restricted time range and the approximated results for the complete time range are given in Table 19. We note that the most meaningful index is the mode, since it corresponds to the peak of the FPT density and the two modes (namely, the mode in the restricted and in the complete time ranges) are quite close to each other (the relative error between the two values is about ).

Table 19

The mean, the mode, the 1st and the 5th deciles and the standard deviation of the FPT in the complete time range and in the restricted time range (real application—FPT problem)

Time range	Mean	Mode	1st decile	5th decile	St. deviation
Complete I_C	221.7308	219.8435	212.9876	221.1406	12.21116
Restricted I_R	197.6721	222.2305	214.8875	228.3496	80.0505

Conclusions

During the recent years, many sigmoidal stochastic models have been introduced to study phenomena of interest in various different scientific areas. In order to model more complex population dynamics in which the maximum level of the growth is reached after many stages, we referred to the multisigmoidal logistic stochastic growth model. More in detail, the present work has been devoted to the analysis of the corresponding statistical inference and of the FPT problem. Two procedures useful to find the MLEs of the parameters have been described, one based on the resolution of the system of the critical points of the likelihood function, and the other one based on the maximization of the likelihood function by means of the S.A. algorithm. Then, the described strategies have been validated with a simulation study. The last section of the paper has been devoted to a real application concerning COVID-19 infections in four different European countries (France, Italy, Spain and United Kingdom). The data have been fitted using a suitable multisigmoidal logistic stochastic model. Finally, a study regarding the FPT through a fixed boundary has been also performed.

Future developments can be oriented to find the MLEs of the parameters with other meta-heuristic optimization procedures (such as Variable Neighborhood Search or other swarm-based algorithms) in order to obtain nice estimates in a short computational time. We aim also to introduce a more sophisticated model suitable to describe better epidemiological dynamics with multiple waves, starting from the multisigmoidal logistic equation. Moreover, aiming at a thorough analysis of the convergence speed for parameter estimation in stochastic differential equations, these approaches will be compared with applications of the recent method called ‘covariance matrix adaptation evolution strategy’. Indeed, the latter is used often in the presence of several parameters (cf., for instance, Ghosh et al. (2012) and Willjuice and Baskar (2010)).