Phenomenology of stochastic exponential growth

TL;DR

This paper explores a broader class of models for stochastic exponential growth beyond the traditional Geometric Brownian Motion, emphasizing power law noise to better match experimental observations of biological and physical systems.

Contribution

It introduces a generalized phenomenological model with power law multiplicative noise, providing analytical solutions and a method to infer parameters from experimental data.

Findings

Mean-rescaled distributions are approximately stationary at long times.

Power law noise models better fit experimental data than GBM.

Analytical solutions and parameter inference methods are developed.

Abstract

Stochastic exponential growth is observed in a variety of contexts, including molecular autocatalysis, nuclear fission, population growth, inflation of the universe, viral social media posts, and financial markets. Yet literature on modeling the phenomenology of these stochastic dynamics has predominantly focused on one model, Geometric Brownian Motion (GBM), which can be described as the solution of a Langevin equation with linear drift and linear multiplicative noise. Using recent experimental results on stochastic exponential growth of individual bacterial cell sizes, we motivate the need for a more general class of phenomenological models of stochastic exponential growth, which are consistent with the observation that the mean-rescaled distributions are approximately stationary at long times. We show that this behavior is not consistent with GBM, instead it is consistent with power law multiplicative noise with positive fractional powers. Therefore, we consider this general class of phenomenological models for stochastic exponential growth, provide analytical solutions, and identify the important dimensionless combination of model parameters which determines the shape of the mean-rescaled distribution. We also provide a prescription for robustly inferring model parameters from experimentally observed stochastic growth trajectories.

Equations

$$\frac{\text{d}x}{\text{d}t}=\kappa x+\sqrt{D}x\xi (t),$$

$$⟨\xi (t)⟩=0,\text{and}⟨\xi (t)\xi (t^{\prime })⟩=\delta (t-t^{\prime }).$$

$$\frac{\text{d}x}{\text{d}t}=\kappa x+\sqrt{D}{x}^{\gamma }\xi (t).$$

$$\frac{\text{d}x}{\text{d}t}=\kappa x+f(x)\xi (t).$$

$$⟨x(t)⟩=x_0{e}^{\kappa t}.$$

$$\frac{\text{d}y}{\text{d}t}=x_0^{-1}{e}^{-\kappa t}f\left(y{x}_0{e}^{\kappa t}\right)\xi (t).$$

$$\begin{array}{rl}\frac{\text{d}{y}^2}{\text{d}t}=& x_0^{-2}{e}^{-2\kappa t}{f}^2\left(y{x}_0{e}^{\kappa t}\right)\\ & +2y{x}_0^{-1}{e}^{-\kappa t}f\left(y{x}_0{e}^{\kappa t}\right)\xi (t).\end{array}$$

$$\frac{\text{d}}{\text{d}t}\text{var}\left[y(t)\right]=x_0^{-2}{e}^{-2\kappa t}⟨f^2\left(y{x}_0{e}^{\kappa t}\right)⟩,$$

$$\text{var}\left[y(t)\right]=x_0^{-2}{\int }_0^{t}d{t}^{\prime }{e}^{-2\kappa t^{\prime }}⟨f^2(x(t^{\prime }))⟩.$$

$${\int }_0^{\infty }dt{e}^{-2\kappa t}⟨f^2(x(t))⟩<\infty .$$

$$⟨f^2(x)⟩\simeq f^2(x_0{e}^{\kappa t}).$$

$$\underset{x\to \infty }{\lim }\log x\frac{f^2(x)}{x^2}=0.$$

$$\frac{\text{d}x}{\text{d}t}=\kappa x+\sqrt{D}{x}^{\gamma }\xi (t).$$

$$\frac{\text{d}x}{\text{d}t}=\kappa x+\sqrt{D}\xi (t).$$

$${\Sigma }_0^2\equiv \frac{D}{2\kappa x_0^2},$$

$$\text{var}\left[y(t)\right]={\Sigma }_0^2\tau (t),$$

$$\tau (t)\equiv 1-e^{-2\kappa t};$$

$$p(x,t|x_0,0)=\frac{1}{x_0}{e}^{-\kappa t}{q}_{\mathrm{abs}}\left(\frac{1}{x_0}{e}^{-\kappa t}x,\tau (t)\right),$$

$$\begin{array}{rl}{q}_{\mathrm{abs}}(y,\tau )\equiv & \frac{2}{\sqrt{2\pi \tau (t)}{\Sigma }_0}\exp \left[-\frac{y^2+1}{2{\Sigma }_0^2\tau (t)}\right]\\ & \times \sinh \left[\frac{y}{{\Sigma }_0^2\tau (t)}\right].\end{array}$$

$$p_{\mathrm{abs}}(t)$$

$$\begin{array}{rl}⟨y^2(t)⟩=& (1+v(t))\text{mbox}Erf\left[\frac{1}{\sqrt{2v(t)}}\right]\\ & +\sqrt{\frac{2v(t)}{\pi }}\exp \left[-\frac{1}{2v(t)}\right],\end{array}$$

$$\frac{\text{d}x}{\text{d}t}=\kappa x+\sqrt{D}x\xi (t).$$

$$x(t)=x_0{e}^{\sqrt{D}W(t)+(\kappa -D/2)t},$$

$$y(t)=e^{\sqrt{D}W(t)-\frac{1}{2}Dt}.$$

$$q(y,t)=\frac{1}{y\sqrt{2\pi Dt}}\exp \left[-\frac{{\left(\ln (y)+\frac{1}{2}Dt\right)}^2}{2Dt}\right].$$

$$\text{var}\left[y(t)\right]=e^{Dt}-1,$$

$$\mathbb{P}(y(t)>ϵ)={\int }_{ϵ}^{\infty }q(y,t)dy$$

$$={\int }_{\log ϵ+\frac{1}{2}Dt}^{\infty }\frac{dz}{\sqrt{2\pi Dt}}{e}^{-\frac{1}{2Dt}{z}^2}$$

$$=\Phi \left(-\frac{1}{\sqrt{Dt}}(\log ϵ+\frac{1}{2}Dt)\right).$$

$$\frac{\text{d}x}{\text{d}t}=\kappa x+\sqrt{Dx}\xi (t).$$

$$\begin{array}{rl}p(x,t|x_0;0)=& \frac{\kappa }{D\sinh (\kappa t/2)}\exp \left[-\frac{2\kappa (x+x_0{e}^{\kappa t})}{D(e^{\kappa t}-1)}\right]\\ & \times \sqrt{\frac{x_0}{x}}{I}_1\left(\frac{2\kappa }{D\sinh (\kappa t/2)}\sqrt{x_{0}x}\right),\end{array}$$

Full text

I Introduction

Stochastic exponential growth occurs when a quantity of interest multiplies geometrically, with a random multiplier and/or with stochastic waiting times between growth spurts. Such processes have been observed in various contexts, across a range of length and timescales. Familiar examples include microbial population growth Ginzburg et al. (1982); Keiding (1975); Lewontin and Cohen (1969); Iyer-Biswas et al. (2016), nuclear fission Everett and Ulam (1948); Birkhoff and Varga (1958), copy number dynamics of key biochemicals Iyer-Biswas et al. (2014a); Tsuru et al. (2009); Friedman et al. (2006); Cai et al. (2006), increase of individual cell sizes during interdivision periods Iyer-Biswas et al. (2014b, a); Pfeuty and Kaneko (2007); Iyer-Biswas and Zilman (2016), tissue growth in developmental and cancer biology Eden (1961); Furusawa and Kaneko (2000); Kohandel et al. (2007); Halpin-Healy and Zhang (1995), spreading of epidemics Girvan et al. (2002), growth of internet and propagation of viral social media posts Huberman and Adamic (1999); Kumar et al. (2000), financial markets Samuelson (1965); Black and Scholes (1973), stochastic divergence of neighboring phase space trajectories Stoddard and Ford (1973), and cosmological expansion of the universe Diamandis et al. (1999).

Despite the variety of contexts in which the phenomenon occurs, literature on modeling stochastic exponential growth has predominantly focused on one model, Geometric Brownian Motion (GBM). Such processes appear naturally in the context of proportionate growth processes; for a general introduction to these processes and their applications see Mitzenmacher (2004). A stochastic process, $x(t)$, is said to undergo GBM if its time evolution obeys a Langevin equation with linear drift and linear multiplicative noise van Kampen (2007):

[TABLE]

in which $\xi(t)$ is Gaussian white noise, traditionally interpreted in the Itô sense van Kampen (2007), with

[TABLE]

It is well-known that this process yields a log-normal distribution for $x$, and the ratio of its standard deviation to mean diverges with time van Kampen (2007).

Recent studies have revealed a scaling law governing fluctuations during stochastic exponential growth of individual microbial cell sizes in between divisions: the mean-rescaled distribution is stationary and has finite variance at long times Iyer-Biswas et al. (2014a, b). This observation is inconsistent with the behavior predicted from GBM, in which the coefficient of variation (ratio of standard deviation to mean) diverges with time, and so the mean-rescaled distribution is not stationary van Kampen (2007).

Motivated by these observations, we seek a general phenomenological description of stochastic exponential growth, consistent with the time independence of the mean-rescaled distribution at long times. To this end we consider a stochastic Markovian growth process which is a generalization of GBM, in which the multiplicative noise term is given by a power-law with exponent $\gamma$:

[TABLE]

This model is well-studied in finance literature in the context of applications in option pricing Cox (1996); Lo et al. (2000). We find constraints on $\gamma$ for which these dynamics are consistent with stationarity of the asymptotic mean-rescaled distribution. We show that the probability density of the mean-rescaled variable, $y(t)\equiv x(t)/\left\langle x(t)\right\rangle$, asymptotically approaches a time-independent distribution with finite variance provided $0\leq\gamma<1$. Thus we focus on the phenomenology of stochastic exponential growth for the class of models given by Eq. (3), with $0\leq\gamma<1$. For this class of models we provide analytical solutions for the time dependent dynamics, and use them to motivate a robust prescription for inferring model parameters ($\gamma$, $\kappa$ , and $D$) from observed stochastic growth trajectories.

In Section II we comment on generic properties of a more general class of models than Eq. (3), namely, ones with a general form of the multiplicative noise term. Remarkably, a stable distribution for the mean-rescaled variable, $y(t)$, is found to exist in the long time limit for this general class of models. We derive a sufficient condition for finiteness of the variance of this asymptotic distribution, and use this constraint to motivate why we focus on the phenomenology of the class of models in Eq. (3). In Section III, we discuss analytical solutions and general properties of Eq. (3). We also examine in detail the limiting cases corresponding to $\gamma=0$ and $\gamma=1$, and a special case, $\gamma=1/2$, which arises naturally from a microscopic model for stochastic exponential growth Iyer-Biswas et al. (2014a). Using the analytical solutions for Eq. (3), we find that a specific combination of the second and third central moments of the distributions is particularly sensitive to the exponent $\gamma$, and not to other model parameters. In Section IV we utilize this idea to provide a robust prescription for inferring the exponent $\gamma$ from time series data, and test the prescription using synthetic data from numerical simulations of Eq. (3). Once $\gamma$ has been determined, estimation of the parameters $\kappa$ and $D$ is a straightforward statistical inference problem. Here we use the Maximization of Likelihood Estimation (MLE) approach of parametric statistical regression, and present its application to the models considered.

II Condition for Finite Variance of the Mean-Rescaled Distribution

To begin with, we consider properties of the mean-rescaled distribution for a general class of stochastic exponential growth models with multiplicative noise:

[TABLE]

We will assume the initial condition, $x(t=0)=x_{0}$, and an absorbing boundary at the origin, $x=0$, to ensure that $x(t)$ remains non-negative. The phenomenological model introduced in Eq. (3) corresponds to the subset of models with multiplicative noise, $f(x)$, of the form $\sqrt{D}x^{\gamma}$.

In the absence of noise , $x$ grows exponentially with time: $x(t)=x_{0}\,e^{\kappa t}$. When fluctuations are present, the ensemble averaged value undergoes exponential growth with growth rate $\kappa$:

[TABLE]

This follows directly from Eq. (4) upon taking the expectation value and integrating over time. We remind the reader that Eq. (4) is to be interpreted in the Itô sense.

What can be said about the spread of stochastic trajectories, $x(t)$, around the expectation value given in Eq. (5), in the presence of the noise term in Eq. (4)? We are especially interested in the mean-rescaled variable $y(t)\equiv x(t)/\left\langle x(t)\right\rangle$, and its behavior in the long time limit, $t\to\infty$. Under general conditions111Assuming that the noise term $f(x)$ satisfies appropriate regularity conditions which are required for the existence and uniqueness of the solution of the stochastic differential equation Eq. (4), it can be shown that $y(t)$ is a continuous time, non-negative martingale. It is known that any non-negative supermartingale (martingale) is $L^{1}$ integrable. By the Doob’s martingale convergence theorem [37], it follows that the limit $\lim_{t\to\infty}y(t)=y_{\infty}$ exists, in an almost sure sense. Furthermore, $y_{\infty}$ is also $L^{1}$ integrable. This implies that $y(t)$ converges to a well-defined random variable in long time limit., $y(t)$ approaches a limiting distribution in the long time limit. However, this result does not place constraints on the finiteness of the variance of the asymptotic distribution of $y(t)$. Motivated by experimental observations in Iyer-Biswas et al. (2014b, a), we are interested in possible functional forms of the noise term $f(x)$, for which the distribution of $y(t)$ has a finite variance as $t\to\infty$.

To this end, we first find an expression for the variance of $y(t)$. Changing variables in Eq. (4),

[TABLE]

Using Itô’s lemma McKean (1969),

[TABLE]

The variance, $\text{var}\left[y(t)\right]$, becomes equal to $\left\langle y(t)^{2}\right\rangle-1$, since $y(t)$ is a mean-rescaled variable and its average is $\left\langle y(t)\right\rangle=1$. Thus, using Eq. (7), we have

[TABLE]

which implies

[TABLE]

Therefore, the mean rescaled distribution has a finite variance only if

[TABLE]

What are the possible functional forms of $f(x)$ which satisfy the condition Eq. (10)? Evidently, any bounded function $f(x)<M$ will satisfy this condition. Thus, any model of multiplicative noise for which the amplitude of the noise saturates for large $x$ will produce a finite variance distribution for $y$ in the large time limit.

The small-noise limit. In this limit, $f(x)\ll\sqrt{\kappa}\,x$, we can approximate the expectation with the deterministic value,

[TABLE]

Therefore the condition on $f(x)$, for finiteness of variance of the asymptotic mean-rescaled distribution, can be simply stated as $f^{2}(x)/x^{2}$ must fall off faster than $1/\log x$ in the limit $x\to\infty$, such that the integral in Eq. (10) converges as $t\to\infty$. This analysis suggests a more explicit condition on $f(x)$:

[TABLE]

Thus if the leading behavior of $f(x)$ in the large $x$ limit is a power law, i.e., $f(x\rightarrow\infty)\sim x^{\gamma}$, then the mean-rescaled distribution can become stationary only for $0\leq\gamma<1$. Hence our focus on the phenomenological models of stochastic exponential growth given by Eq. (3).

In the following Section, before proceeding to give a general solution for models of the form in Eq. (3), we will contextualize the results by carefully considering the limiting cases, $\gamma=0$ (constant additive noise) and $\gamma=1$ (linear multiplicative noise, corresponding to GBM). In addition, we will also consider the $\gamma=1/2$ (square-root multiplicative noise) case in detail, since it arises naturally in discrete stochastic models (for e.g., the stochastic Hinshelwood cycle Iyer-Biswas et al. (2014a), birth-death and Galton-Watson branching processes Ethier and Kurtz (1986)).

In relation to Eq. (12), the expected behavior of the asymptotic mean-rescaled distributions for these cases are summarized below.

• •

Constant additive noise, $f(x)=\sqrt{D}$: since the noise has bounded intensity for this case, the asymptotic distribution of $y(t)$ has finite variance. This distribution can be found exactly and is Gaussian (see Section III).

• •

Linear multiplicative noise, $f(x)=\sqrt{D}x$: this corresponds to geometric Brownian motion (GBM). Here the condition Eq. (12) is not satisfied, and thus, as explicitly derived in the following Section, the variance of $y(t)$ increases without limit.

• •

Power-law multiplicative noise, $f(x)=\sqrt{D}x^{\gamma}$, with $0<\gamma<1$: since the condition in Eq. (12) is satisfied for this choice of $f(x)$, the variance of the asymptotic mean-rescaled distribution is finite. We will show this explicitly in Section III.

III Power-law multiplicative noise: solutions and properties

As motivated in the preceding section, we now focus on the phenomenology and analytical solutions of stochastic exponential growth process with power-law multiplicative noise:

[TABLE]

To contextualize the general results, we first review the properties of the familiar special cases of this equation, i.e., $\gamma\in\{0,1,1/2\}$ in Section IIIA — C, before proceeding to the general solution in Section IIID. The constant additive noise case, $\gamma=0$, corresponds to a simple stochastic differential equation which is known to have a Gaussian solution Uhlenbeck and Ornstein (1930); van Kampen (2007). Linear multiplicative noise, $\gamma=1$, corresponds to geometric Brownian motion (GBM). This process was proposed as a model for risky asset prices by Samuelson Samuelson (1965), and was used by Black and Scholes in their pricing theory for equity derivatives  Black and Scholes (1973). Square-root multiplicative noise, $\gamma=1/2$, corresponds to the square-root process Feller (1951); this is obtained as the continuous limit of discrete stochastic models such as the the stochastic Hinshelwood cycle Iyer-Biswas et al. (2014a). More generally, square-root processes arise in the context of birth-death and Galton-Watson branching processes Ethier and Kurtz (1986); Iyer-Biswas et al. (2009); Iyer-Biswas and Jayaprakash (2014); Iyer-Biswas (2009).

III.1 Constant additive noise, $\gamma=0$

The simplest model of stochastic exponential growth has constant amplitude noise, $f(x)=\sqrt{D}$, and the corresponding stochastic differential equation has a simple Gaussian noise term:

[TABLE]

This immediately suggests that the dimensionless ratio,

[TABLE]

which can be identified as the variance of the mean-rescaled distribution, should be a key determinant of its shape (the subscript [math] for $\Sigma^{2}$ indicates that this corresponds to $\gamma=0$). This dimensionless variable is the inverse Péclet number of the process, and as we will see, we will need to generalize the Péclet number (or its inverse) for different values of $\gamma$.

In the absence of a boundary condition at the origin, this process is a mean-repelling Ornstein-Uhlenbeck process Uhlenbeck and Ornstein (1930). We remind the reader that the Ornstein-Uhlenbeck process has the form ${dx}/{dt}=\kappa(\mu-x)+\sqrt{D}\xi(t)$, where $\kappa>0$ corresponds to the mean-reverting case, and $\kappa<0$ is the mean-repelling case. For both cases this is known to yield a Gaussian distribution. Therefore, the mean-rescaled distribution, $q(y,t)$, is also Gaussian, with variance given by

[TABLE]

with

[TABLE]

$\tau(t)$ identifies the characteristic timescale of approaching stationary-state. Thus, the mean-rescaled distribution, $q(y,t)$, asymptotically approaches a time-independent Gaussian distribution with variance $\Sigma_{0}^{2}$, for $t\gg 1/(2\kappa)$. By definition, the mean of the distribution is unity at all times.

In proximity of the origin this stochastic process can become negative. However, in real-world applications of the stochastic exponential growth process (e.g., the growth dynamics of bacterial cell sizes), the dynamical variable (e.g., the cell’s size) may be constrained to always remain positive. This can be addressed in the model by imposing an absorbing boundary condition at the origin, $x=0$, to avoid negative values for the dynamical variable. We can compute the solution for $p(x,t|x_{0},0)$, the probability density of $x$ at time $t$ given the initial condition $x_{0}$ at time $t=0$, by using the method of images (see Appendix A for details of derivation). We obtain

[TABLE]

where

[TABLE]

Since the modified time, $\tau(t)$, approaches a finite limit ($=1$) in the long time limit, it follows that $y(t)$ approaches a stationary distribution asymptotically; this distribution is given by using $\tau(t)\to 1$ in Eq. (18).

The probability of absorption at origin is

[TABLE]

The mean of $y(t)$ is equal to $1$, as expected, and the second moment is given by

[TABLE]

where $v(t)\equiv\Sigma_{0}^{2}\,\tau(t)$. For $v(t)\ll 1$, one has $\text{var}\left[y(t)\right]\simeq v(t)$. This variance approaches a limiting value in the infinite time limit: $\text{var}\left[y(t)\right]\simeq{D}/({2\kappa\,x_{0}^{2}})\equiv\Sigma_{0}^{2}$.

III.2 Linear multiplicative noise (geometric Brownian Motion), $\gamma=1$

The case of $\gamma=1$ corresponds to another simple model. Here the noise amplitude is proportional to the magnitude of the process, $x(t)$ and so

[TABLE]

This process has been well studied and is known as geometric Brownian motion (GBM) van Kampen (2007). The solution for $x(t)$ is

[TABLE]

where $W(t)$ denotes standard Brownian motion.

Using Itô’s lemma, we find that the mean-rescaled process also follows geometric Brownian motion:

[TABLE]

The probability density of the mean-rescaled variable, $q(y,t)$, is given by

[TABLE]

The variance of $y(t)$ is given by

[TABLE]

while the mean is always unity.

In the small time limit, $t\ll 1/D$, the variance increases linearly as $\text{var}\left[y(t)\right]\sim t$, while in long time limit, $t\gg 1/D$, it increases exponentially as $\text{var}\left[y(t)\right]\sim e^{Dt}$. Thus, for $\gamma=1$, the variance of $y(t)$ does not approach a finite asymptotic value in the infinite time limit. The distribution of $y(t)$ becomes increasingly concentrated near zero, as can be seen by computing the probability that $y(t)$ is larger than an arbitrary value, $\epsilon$.

[TABLE]

Here $\Phi(z)=\int_{z}^{\infty}\frac{dt}{\sqrt{2\pi}}e^{-\frac{1}{2}t^{2}}$ is the normal cumulative distribution function. This probability goes to zero as $t\to\infty$, for any $\epsilon>0$. At the same time, the variance of $y(t)$ increases without limit. The seemingly paradoxical conclusion that the distribution of $y(t)$ becomes more and more concentrated near the origin while its variance grows without limit can be understood using the argument provided by Lewontin and Cohen Lewontin and Cohen (1969), who formulated it in the general context of the proportionate growth models Mitzenmacher (2004).

III.3 Square-root multiplicative noise, $\gamma=\frac{1}{2}$

Next we consider the special case, $\gamma=\frac{1}{2}$, which arises naturally in the context of phenomenological models corresponding to discrete stochastic exponential growth process Feller (2008). Rewriting Eq. (3) with a square-root multiplicative noise term, one has

[TABLE]

Processes with multiplicative square-root noise often appear in models of population genetics and in the Langevin description of birth-death processes Feller (2008). The properties of the stochastic differential equations with multiplicative square-root noise have been studied in detail by Feller in Ref. Feller (1951).

In the absence of the noise term this process describes exponential growth with a constant growth rate, $\kappa>0$. The noise term introduces fluctuations which have variance proportional to $x$. The transition density $p(x,t|x_{0},0)$ can be read off from Ref. Feller (1951) (also see Equation (1.6) in Ref. Mazzon (2012)), and is given by

[TABLE]

where $I_{1}(x)$ is the modified Bessel function of the first kind.

One can check by direct calculation that the transition density satisfies the following relations.

[TABLE]

The relation Eq. (29) encodes normalization of the probability while accounting for the non-zero probability of absorption at the boundary. As expected, in Eq. (30) we find that the mean grows exponentially with time, with growth rate, $\kappa$.

Using Itô’s lemma, we find that the second moment, $\left\langle x^{2}(t)\right\rangle$, satisfies the following time evolution equation:

[TABLE]

Solving it with initial condition, $\left\langle x^{2}(0)\right\rangle\equiv x_{0}^{2}$, one has

[TABLE]

Thus the variance of $x(t)$ is

[TABLE]

Once again, we can identify a dimensionless variance, a combination of parameters which determines the shape of the mean rescaled distribution, from the preceding equation:

[TABLE]

We can find the density of the mean rescaled variable $y(t)$ directly from the distribution of $x(t)$ given in Eq. (28). We present next an alternative approach, based on the method of the time change. This has the advantage that it is applicable also to the case of general exponent $\gamma$, which will be discussed next.

An application of Itô’s lemma yields the following stochastic differential equation for $y(t)$:

[TABLE]

with the initial condition $y_{0}=1$. The time-dependent factor can be absorbed into a deterministic rescaling of time. To this end, we define the modified dimensionless time, $\tau(t)$, by

[TABLE]

which implies $d\tau=\text{e}^{-\kappa t}\kappa\,dt$. In terms of $\tau$, Eq. (35) can be written as

[TABLE]

The probability density of $y(t)$ can be read off from Eq. (28):

[TABLE]

where

[TABLE]

The variance of $y(t)$ is given by $\text{var}\left[y(t)\right]=v(t)$. As $t\to\infty$ limit, $\tau$ approaches $1$, thus the mean-rescaled distribution tends to a stationary distribution with variance equal to $\Sigma_{\frac{1}{2}}^{2}$.

In the small noise limit, i.e., $D\ll\kappa\,x_{0}$, the asymptotic probability density, $q^{*}(y)$, is sharply peaked around $1$, and is well approximated as

[TABLE]

This relation follows from the general result, Eq. (38), upon using the asymptotic expansion of the Bessel function, $I_{\alpha}(x)$: in the limit $x\to\infty$, it approaches $\frac{e^{x}}{\sqrt{2\pi x}}(1+O(x^{-1}))$.

III.4 General exponent $\gamma$

We have seen that for $\gamma=0$ and $1/2$, there are characteristic timescales equal to ${1/(2\kappa)}$ and ${1}/{\kappa}$, respectively, after which the mean-rescaled distributions become stationary. After this time, the shape of the distributions only depends on the dimensionless variables $\Sigma_{0}$ and $\Sigma_{1/2}$. In this section, we first find the corresponding timescale for a general $\gamma<1$. Then, we find a dimensionless variable, $\Sigma_{\gamma}$, that determines the shape of the mean rescaled distribution at steady state. This dimensionless variable is a generalization of the inverse Péclet number, which is often used in transport literature Patankar (1980).

Let us start with the case of a power-law multiplicative noise of the form $f(x)=\sqrt{D}x^{\gamma}$, with a general exponent, $\gamma>0$.

[TABLE]

This process appears in financial literature in the context of models known as Constant Elasticity of Variance (CEV) Cox (1996). We will show that its properties can be understood by relating it to a generalization of the square-root process ($\gamma=1/2$) that we considered in the preceding section. The small-noise approximation result, Eq. (12), suggests that $\gamma$ must be restricted to the range $\gamma<1$, in order to have a finite variance for $y(t)$ in the large time limit. We will prove this constraint rigorously.

We would like to derive the probability density, $p(x,t|x_{0},0)$, for the process defined by Eq. (3). It will prove more convenient to first derive the probability density of the mean-rescaled variable, $q(y,t)$, and subsequently obtain $p(x,t|x_{0},0)$ using

[TABLE]

Note that the properties of the mean rescaled distribution should depend on the two dimensionless parameters $\gamma$ and $D/(\kappa\,x_{0}^{2(1-\gamma)})$. It would however simplify the algebra if we define the dimensionless variance as the following combination of these parameters:

[TABLE]

To obtain the probability density, $q(y,t)$, of the mean-rescaled variable, $y$, we start by finding an SDE for the variable $y$ using Itô’s lemma:

[TABLE]

with the initial condition $y_{0}=1$. Following the same approach as before, we absorb the time dependent factor into a deterministic time change. We define the modified dimensionless time, $\tau\equiv\tau(t)$, such that $d\tau=2\kappa(1-\gamma)e^{-2\kappa(1-\gamma)t}dt$. Thus,

[TABLE]

$\tau(t)$ approaches $1$ as $t\to\infty$, provided that $\gamma<1$. The time scale for the mean rescaled distribution to approach its stationary solution is set by the exponent $1/(2\kappa(1-\gamma))$. This timescale itself approaches infinity as $\gamma\to 1$, hence a stationary mean-rescaled distribution is never attained in the GBM case.

With this change of time variable Eq. (43) reads:

[TABLE]

Equation (45) can be reduced to a square root process with constant drift by the following change of variable:

[TABLE]

Using Itô’s lemma, one has the process that $z$ follows:

[TABLE]

with the initial condition $z_{0}=1$. Stochastic differential equations with square-root multiplicative noise have been studied in detail by Feller in Feller (1951); this paper considered the general case of SDEs of the form

[TABLE]

This is a generalization of the square-root process with an additional constant drift term. The process in Eq. (47) is obtained by the following substitutions:

[TABLE]

Feller Feller (1951) showed that the stochastic differential equation Eq. (48) has different qualitative behaviors depending on the relative strength and sign of the coefficients. The following three cases are important to distinguish between: (i) $0<c<a$, (ii) $c<0$, and (iii) $c>a>0$.

For the problem we are interested in, only two of the cases in the Feller classification of solutions are relevant:

(i) $0<c<a$: this corresponds to $\gamma\in(0,\frac{1}{2})$. The fundamental solution of the forward Kolmogorov equation corresponding to Eq. (47) is not unique. There are two independent fundamental solutions, and the problem is well-posed only if we add an additional boundary condition at $z=0$, for example absorbing or reflecting boundary. With an absorbing boundary condition, the transition density is given by Eq. (50) below.

(ii) $c<0$: this corresponds to $\gamma\in[\frac{1}{2},1)$. There is only one fundamental solution, given by

[TABLE]

This can be extracted from equation (1.6) of Mazzon (2012) by taking the $b\to 0$ limit. Substituting for $a$ and $c$ from Eq. (49), we obtain

[TABLE]

The probability of absorption at origin is found to be

[TABLE]

where $\Gamma(a,z)$ is the upper incomplete Gamma function, defined as $\Gamma(a,z)=\int_{z}^{\infty}t^{a-1}e^{-t}dt$. The moments of the $y(t)$ can be computed in closed form using the probability distribution function $q(y,t)$ given in Eq. (51). As expected, the mean, $\left\langle y(t)\right\rangle$, evaluates to unity. The next two moments are

[TABLE]

Here ${}_{1}F_{1}(a,b,z)$ is the confluent hypergeometric function, $\tau(t)$ is the function given in Eq. (44) and $a=2{\Sigma}_{\gamma}^{2}(1-\gamma)^{2}$. It is straightforward to show that these equations reduce to the results derived in the previous section for the $\gamma=\frac{1}{2}$ case. The moment relations Eq. (III.4) and Eq. (III.4) also hold for $\gamma\in(0,1)$, provided that an absorbing boundary condition is imposed at $x=0$.

Figure 1 compares the qualitative behavior (on a log-linear scale) of trajectories of the stochastic exponential growth processes described by Eq. (3), for different values of the two dimensionless variables which determine the dynamics, namely, $\gamma$ and $\Sigma_{\gamma}$. The variance of the process is set by $\Sigma_{\gamma}^{2}$, while the exponent $\gamma$ determines the relative stochasticity at short and long times. Small $\gamma$ values correspond to more stochastic/deterministic behavior at short/long times, while for $\gamma$ values around $1$, the log-scale trajectories stay stochastic over a longer time.

The limit of small noise

In the large Péclet number limit, or small noise regime, i.e., when $\Sigma_{\gamma}\ll 1$, the distribution of $y(t)$ is sharply peaked around 1. The probability density, $q(y,t)$, given by Eq. (51), can be approximated using the approximation of the Bessel function for large argument:

[TABLE]

in which

[TABLE]

The exact result for the density Eq. (51) depends on the model parameters $(\kappa,D,x_{0},t)$ through the same combination $v(t)$ which approximate the variance of the distribution:

[TABLE]

This is obtained from Eq. (III.4) using the asymptotic expansion of the confluent hypergeometric function of large positive argument

[TABLE]

For the small noise regime, the $t\to\infty$ limit of the variance is simply our dimensionless variance, $\Sigma_{\gamma}^{2}$:

[TABLE]

which is finite for $\gamma<1$.

Itô vs. Stratonovich interpretations

For many practical applications, the asymptotic variance of the mean-rescaled distribution, which can also be interpreted as the square of coefficient of variation of the original variable, is expected to be small  Iyer-Biswas et al. (2014b). Here, we show that in this region, for $\gamma<1$, there is no difference between the behavior of Eq. (3), and a similar equation where the noise is interpreted in the Stratonovich sense, i.e. ,

[TABLE]

We start by writing the Itô equivalent of Eq. (59):

[TABLE]

Let us first examine the limiting cases: For $\gamma=0$, the Itô and Stratonovich equations are exactly the same. For $\gamma=1$, the Stratonovich equation is the same as the Itô equation with $\kappa\to\kappa+D/2$. For any $\gamma$ in between, the extra term added to the deterministic part of the equation decays as $x$ grows, and therefore, does not affect the dynamics of $x$ at long times. However, the short time dynamics can affect the long time distribution. This is where the small noise limit plays a role.

From the behavior of the two limiting cases, we know that the extra term is more problematic as we get closer to $\gamma=1$. But for $\gamma<1$, the additional term starts with a value strictly smaller than $\Sigma_{\gamma}^{2}$, and its expected value decays exponentially over time. Therefore, for $\Sigma_{\gamma}^{2}\ll 1$, the dynamics of Eq. (59) and Eq. (3) are the same. In summary, one does not need to worry about the choice of Itô vs. Stratonovich schemes when phenomenologically modeling a stochastic exponential growth process.

Let us further quantify this statement by defining $\varepsilon=\left\langle x_{\text{It\^{o}}}-x_{\text{Str}}\right\rangle/\left\langle x_{\text{It\^{o}}}\right\rangle$ to be the relative error of choosing the incorrect prescription. Then we have

[TABLE]

In particular, the relative error is a bounded function of time and is proportional to $\Sigma_{\gamma}^{2}$.

When the system is not in the small noise limit, Itô and Stratonovich schemes are no longer equivalent. In this regime, we can still solve the Stratonovich problem by mapping it to an Itô problem. The general $\gamma$ model in Stratonovich interpretation can be reduced by a change of variable to the $\gamma=0$ case in Itô formalism, discussed in Section III.A. If we change of variables to $y=x^{1-\gamma}/(1-\gamma)$, then we have

[TABLE]

Now the noise is additive and there is no distinction between Itô and Stratonovich interpretations, and so we can use the previously found solution to this model.

Possible generalizations of the models

We comment briefly on the assumptions of our model, and possible generalizations. The assumption of Markovian dynamics is a standard assumption in models of stochastic kinetics van Kampen (2007), and simplifies the analytical treatment. We also assume that the noise is modulated by the amplitude of the process being modeled $x(t)$. A more sophisticated treatment could take into account the stochasticity of the growth rate, $\kappa$, which could represent a stochastic environmental variable. A similar approach is taken in population dynamics modeling, where the analog of $\kappa$ is an environmental variable, which might represent factors such as food supply or temperature  Cohen (1976); Tuljapurkar and Orzack (1980); Cohen (2014). One model of this type in Ref. Pirjol (2015) considered a proportionate growth process where the growth multipliers are modulated by GBM. Interestingly, in this model the growth rate of the mean behaves discontinuously as model parameters cross a critical curve, similar to what is observed in a phase transition.

IV Statistical Inference of the Model Parameters

Our goal in this section is to use the analytical results provided in Section III to find a robust prescription for the determining the exponent $\gamma$ and the other model parameters, in our case $\kappa$ and $D$ from experimentally observed stochastic growth trajectories. We first, in Section IV.1, show that there is a ratio of the third moment and the second moment of the mean rescaled distribution that is only sensitive to the exponent $\gamma$ and can be used to determine the $\gamma$. Then in Section IV.2, we discuss the use of maximum likelihood estimation to find the other parameters of the stochastic dynamics.

IV.1 Determination of $\gamma$

The results for the second and third moments of $y(t)$ in Eq. (III.4) and Eq. (III.4) depend only on the two combinations of variables

[TABLE]

Experimental measurements of the two moments $\text{var}\left[y(t)\right]$ and the central third moment $\left\langle(y(t)-1)^{3}\right\rangle$ can be thus translated into a constraint for the two parameters in Eq. (63). This offers a possibility of determining the noise exponent parameter $\gamma$ from the shape of the $y(t)$ distribution.

As shown in Eq. (56), in the small-noise limit the variance of this distribution gives the parameter

[TABLE]

If the time $t\gg 1/(2\kappa(1-\gamma))$ (much larger than the characteristic growth time), then the variance determines ${\Sigma}_{\gamma}^{2}\tau_{\infty}={\Sigma}_{\gamma}^{2}$.

Using Eq. (51) we computed the shape of the distribution $q(y)$ for several values of $\gamma$, for several values of the variance $\text{var}\left[y(t)\right]=0.01,0.1,0.2$. The results are shown in Fig. 2, where for each case we plot $q(y)$ for $\gamma=0.5$ (black), 0.75 (blue) and 0.9 (red). These plots show that the shape of the distribution is sensitive to $\gamma$.

In order to determine $\gamma$ we propose the measurement of the ratio

[TABLE]

This is defined such that it is equal to 1 for $\gamma=1/2$, and its deviation from 1 is a measure of the difference of $\gamma$ from $1/2$.

In Figure 3 we show the theoretical prediction for $r_{3}(\gamma)$ for three values of the variance $\text{var}\left[y(t)\right]=0.01,0.1,0.2$. This ratio can be computed explicitly as function of $\gamma$ and $v(t)$ as

[TABLE]

where the second and third moments are given explicitly in Eq. (III.4) and Eq. (III.4). The results are shown in Figure 3, which shows good discriminating value in $\gamma$, and little sensitivity to the variance of $y(t)$. Although the theoretical results in Eq. (III.4) and Eq. (III.4) are derived only for $\gamma\in[\frac{1}{2},1)$, we used them also for $\gamma<\frac{1}{2}$, in order to show that the ratio $r_{3}$ varies monotonically with $\gamma$ through $\gamma=\frac{1}{2}$. As mentioned above, the theoretical results of Eq. (III.4) and Eq. (III.4) hold for all $\gamma\in(0,1)$, provided that for $\gamma<1/2$, an absorbing boundary condition is imposed at $x=0.$

The insensitivity of the ratio $r_{3}$ to $v(t)$ at low noise limit can be understood by expanding Eq. (III.4) and Eq. (III.4) using Eq. (57) in $v(t)$ and calculating $r_{3}$:

[TABLE]

showing why a linear relationship is observed in the small $v(t)$ limit. Note that the exact results from Eq. (66) should be used to infer $\gamma$ from distributions with non-zero variance, instead of the approximate linear relationship.

IV.2 Maximum-likelihood estimation

The method for determining $\gamma$ described in the previous section uses only the distributional properties of the $y(t)$ variable at a fixed time $t$. It is also possible to constrain and determine the parameters of the stochastic process, namely, $\kappa,D,$ and $\gamma$, by using time-series data for the $y(t)$ process.

The simplest estimation method for the parameters of a stochastic diffusion is the Maximum-Likelihood Estimation (MLE) method Bartlett (1955); Billingsley (1961). This has been also applied to the estimation of the parameters for a stochastic kinetic model Erdi and Lente (2014).

Here we provide a brief summary of the method, applied to the model at hand. Assume that we observe the stochastic process $x(t)$ on a sequence of times $\{t_{i}\}_{i=1}^{n}$. We would like to use the time series of observations $\{x_{i}\}_{i=1}^{n}$ to determine or constrain the parameters of the continuous-time diffusion which is assumed to have generated the observed data

[TABLE]

The coefficients of the process depend on certain parameters $\theta$.

Define the log-likelihood function

[TABLE]

Here $p_{x}(x_{i},t_{i}|x_{i-1},t_{i-1};\theta)$ is the conditional density of $x_{i}$ given $x_{i-1}$, obtained from the diffusion Eq. (68). The log-likelihood function depends on the parameters $\theta$.

The MLE method determines the parameters $\theta$ by minimizing the log-likelihood function $L_{n}(\theta)$. This method is easy to use in cases where the conditional density $p_{x}$ appearing in Eq. (69) is known in closed form. This is the case for many popular one-dimensional diffusions, like geometric Brownian motion, the square root model Feller (1951), and the Ornstein-Uhlenbeck process Uhlenbeck and Ornstein (1930).

Explicit results in Eq. (28) and Eq. (41) for the transition density in the processes Eq. (27) and Eq. (3), respectively, given in the previous section, can be used to determine the parameters $\theta=\{\kappa,D,\gamma\}$ of these processes from time series of observations of the stochastic variable $\{x_{i}\}_{i=1}^{n}$. We propose combining the MLE method with the approach based on a measurement of $r_{3}(\gamma)$, defined in Eq. (65). Once $\gamma$ has been determined according to the prescription in Sec. IV.1, the MLE will be used to infer the remaining parameters, $\{\kappa,D\}$. As an illustrative example, Appendix B gives a detailed application of this method to the estimation of the parameters of the exponential growth process with square-root noise.

V Summary

To characterize the phenomenology of stochastic exponential growth, we have considered a general class of Markovian stochastic models with power-law multiplicative noise described by Eq. (3). In this model trajectories start at $x=x_{0}$ at time $t=0$ and undergo exponential growth with rate constant $\kappa$ subject to power-law multiplicative noise with strength $D$ and exponent $\gamma$. Since the expected value of this process is a simple exponential growth, $x_{0}e^{\kappa t}$, to understand the dynamics of the distribution, we only need to understand the dynamics of the distribution of the mean rescale variable $y=x/\left\langle x\right\rangle$. This paper provides an exhaustive study of the dynamics of this mean-rescaled distribution.

Starting with Eq. (3), we have shown that the dynamics of the mean-rescaled distributions depend only on the (dimensionless) exponent $\gamma$, and a time-dependent dimensionless factor (see Eq. (39)):

[TABLE]

The latter combines the effect of all the other variables, namely $\kappa$, $D$, $x_{0}$, and, $t$. In the small noise regime, $v(t)$ is approximately the time-dependent variance of the distribution. The timescale given by $1/(2\kappa(1-\gamma))$ sets the time for the distribution to approach steady state. This timescale approaches infinity as $\gamma\to 1$, thus stationarity of the mean-rescaled distribution is never attained in the geometric Brownian motion (GBM) case (the case of $\gamma=1$); this is a unique case for which the mean-rescaled distribution becomes independent of the initial condition $x_{0}$. On the other hand, for $\gamma<1$, the mean-rescaled distributions approach stationary distributions with a finite variances at long times.

For $\gamma>0$, both the deterministic and the stochastic terms in Eq. (3) go to zero as $x\to 0$ which seems to suggest the presence of an absorbing boundary at $x=0$. However, due to singularity of the multiplicative noise at the origin, for $\gamma<1/2$, the trajectories can leak through the origin, and therefore, one needs to artificially impose an additional boundary condition at the origin. The choice of boundary condition (absorbing or reflecting) affects the dynamics of the distribution in $x>0$ region. This is not the case for $\gamma\geq 1/2$, for which $x=0$ is a natural absorbing boundary.

For $\gamma<1$, the dimensionless combination of parameters, $\Sigma_{\gamma}^{2}$, defined by $\Sigma_{\gamma}^{2}\equiv\lim_{t\to\infty}v(t)=Dx_{0}^{2(\gamma-1)}/(2\kappa(1-\gamma))$, sets the asymptotic variance of the mean rescaled distribution and is the analogue of the inverse Péclet number. When $\Sigma_{\gamma}\ll 1$, the dynamics are almost deterministic, and the distributions for all $\gamma$’s are well-approximated with a narrow Gaussian around the mean. However, when $\Sigma_{\gamma}\sim 1$ or larger, the exponent $\gamma$ determines the skewness of the distribution. We have identified a particular metric of skewness, namely, the ratio of the third central moment to the square of variance, which is sensitive to $\gamma$ and very insensitive to $\Sigma_{\gamma}$. Thus, this measure can be used to infer the value of $\gamma$ from experimentally measured distributions. Additionally, we have provided analytical expressions for transition probabilities $p(x,t|x_{0},0)$ (the probability density of $x$ at time $t$, given $x(0)=x_{0}$), and a prescription that incorporates these analytical expressions to infer the remaining model parameters from observed growth trajectories by using MLE methods.

VI Acknowledgments

We thank Rudro Biswas and Sasha Grosberg for insightful discussions. We acknowledge financial support from Purdue University Startup Funds and Purdue Research Foundation.

VII Author Contributions

SI-B conceived of and designed research; DP, FJ and SI-B performed calculations and simulations, and wrote the paper.

Appendix A Constant Noise with Absorbing Boundary: Method of Images

We start with Eq. (13). Recalling that $x(t)=x_{0}e^{\kappa t}y(t)$, the application of Itô’s lemma gives that $y(t)$ follows the diffusion:

[TABLE]

The time dependent factor can be absorbed into the time change

[TABLE]

simplifying Eq. (70) to

[TABLE]

In the absence of the boundary condition at $x=0$, the solution of the SDE in Eq. (72) is a time-changed Brownian motion, shifted by a constant amount. Its distribution is a Gaussian which is centered at $1$ and has variance

[TABLE]

The absorbing condition requires that the probability density of $y$ at $x=0$, $q_{\rm abs}(0,\tau)=0$ for all $\tau\geq 0$. The solution of the Eq. (70) can be found by the method of images and can be constructed by superposition of two Gaussian distributions, centered around $y=1$ and $y=-1$ respectively:

[TABLE]

where

[TABLE]

It is straightforward to show that the solution, Eq. (74), satisfies the absorbing boundary condition at $y=0$ for all $t\geq 0$. Figure 4 shows a plot of the density $q_{\rm abs}(x,\tau)$, and of its construction by the method of images as the difference of two Gaussian distributions. For an alternative derivation of this result, see Feller (1951).

Appendix B Details of application of MLE

We illustrate here an example of the application of the MLE method for determining the parameters $\kappa$ and $D$ of the process

[TABLE]

from a time series of observations $\{x_{i}\}_{i=1}^{n}$ of the process on the sequence of times $\{t_{i}\}_{i=1}^{n}$.

The log-likelihood function is defined in Eq. (69) and requires the conditional transition density function $p_{x}(x_{i},t_{i}|x_{i-1},t_{i-1};\theta)$. For the square root noise process Eq. (76) this is given in Eq. (28), and for the general $\gamma$ model it is given in (41).

The MLE method determines the model parameters by minimization of the log-likelihood function, $L_{n}(\theta)$, over the values of the model parameters $\theta=\{\kappa,D\}$. A judicious choice of the initial starting point for these parameters can speed up the minimization.

We start by considering the simpler case of the square root process $\gamma=\frac{1}{2}$. For this case we propose to choose the initial values $(\kappa_{0},D_{0})$ such that the first two conditional moments are matched exactly. They are given by Eq. (30) and Eq. (32)

[TABLE]

This gives

[TABLE]

For general $\gamma\neq\frac{1}{2}$, the MLE method can be reduced to the estimation of the parameters of a square-root model by working with the time series for $z_{i}=x_{i}^{2(1-\gamma)}$. We assume that $\gamma$ has been determined from a measurement of the ratio $r_{3}(\gamma)$, as explained in Sec. IV.1. The MLE method will be applied to determine the growth rate $\kappa$ and noise parameter $D$, for known $\gamma$.

For this case it is convenient to write the log-likelihood function in terms of the $z_{i}$ variables, for which the transition density is known from Eq. (50) and has a simpler form.

A first estimate for $\kappa$ can be obtained for general $\gamma$ again from matching the first conditional moment (77), which gives the estimator $\kappa_{0}$ in (79). The analog of the relation (78) is given by

[TABLE]

with $\langle y_{\tau}^{2}\rangle$ given by (III.4) with the replacement $t\to\tau$. This gives a complicated non-linear equation for $D$, which has to be solved to obtain $D_{0}$.

A simpler approach starts with the relation

[TABLE]

with $\Delta=1-e^{-2\kappa(1-\gamma)\tau}$. This follows from the diffusion (47) for $z_{t}=y_{t}^{2(1-\gamma)}$. From this equation we obtain the conditional moment

[TABLE]

If $\kappa\tau\ll 1$, we have to a good approximation $\Delta\simeq 2\kappa(1-\gamma)\tau$. We obtain the following simple estimator for $D$

[TABLE]

This estimator is expected to be useful when $\gamma$ is not too close to $\frac{1}{2}$. For this case the numerator and denominator are small, and the estimator may be unstable numerically.

In order to illustrate the application of this method, we generated sample data by simulating the SDE of the square-root process $\gamma=\frac{1}{2}$ using an Euler discretization in time. For this simulation we generated $n=100$ values on a time grid with time step $\tau=0.01$. The model parameters assumed are $\kappa=0.5,\sqrt{D}=0.01$ and starting point $x_{0}=0.01$. Figure 5 shows the sample path for this data.

For the example data set considered here the initial estimators Eq. (79), Eq. (80) are $\kappa_{0}=0.47749$ and $\sqrt{D_{0}}=0.01918$.

Minimization of the log-likelihood function over $\kappa$ and $D$ with the starting point $(\kappa_{0},D_{0})$ gives $\sqrt{\hat{D}}=0.01005$ and $\hat{\kappa}=0.47734$. These values are close to the actual parameter values $\sqrt{D}=0.01$ and $\kappa=0.5$. The error of the determination decreases with $n$, the number of the data points. See Bartlett (1955); Billingsley (1961) for error estimates, and the $n\to\infty$ asymptotics of the MLE errors.