Environmental variability and mean-reverting processes

Environmental variability is often incorporated in a mathematical model 
 by modifying the parameters in the model. 
 In the present investigation, two common methods to incorporate 
 the effects of environmental variability in stochastic differential equation 
 models are studied. The first approach 
 hypothesizes that the parameter satisfies a mean-reverting stochastic process. 
 The second approach hypothesizes that the parameter is a linear function of 
 Gaussian white noise. The two approaches are discussed and 
 compared analytically and computationally. Properties of several mean-reverting processes 
 are compared with respect to nonnegativity and their asymptotic stationary 
 behavior. The effects of different environmental 
 variability assumptions on population size and persistence time for 
 simple population models are studied and compared. 
 Furthermore, environmental data are examined for a gold mining stock. It is 
 concluded that mean-reverting processes possess several advantages 
 over linear functions of white noise in modifying parameters for environmental 
 variability.


1.
Introduction. There are numerous fundamental biological parameters, such as per capita birth or death rate, carrying capacity, infection contact rate, and recovery rate, that are constant for a given species in a non-varying environment. However, in a varying environment, these parameters also vary. In the present investigation, two approaches are examined for incorporating environmental variability into the parameters of a mathematical model for a dynamical biological system. To define the problem and to study its nature, consider the simple birth-death process modeled deterministically by with y(0) = y 0 where b(t) and d(t) are per capita birth and death rates respectively. In a laboratory setting with no external influences, b(t) and d(t) are constants and (1) is an accurate model for the average population size y(t). However, demographic randomness in the birth-death process is present even in a well-controlled laboratory setting. Through generalizing the deterministic model to a stochastic differential equation model, equation (1) can be modified to account for the randomness in 2. Modeling environmental variability in SDEs. Two approaches to model environmental variability in the parameters b(t) and d(t) in (1) are considered. The first approach involves applying a mean-reverting SDE process, sometimes referred to as colored noise [11], to the parameter. The second approach is through modeling the parameter with a linear function of Gaussian white noise (WN), specifically, a linear function of the derivative of a Wiener process. A primary difference between the two approaches is that, unlike Gaussian WN, a mean-reverting process is continuously varying with time. For example, for a small time interval ∆t, a meanreverting process X(t) has correlation coefficient ρ(X(t), X(t + ∆t)) = 1 − O(∆t) while ρ(X(t), X(t + ∆t)) = 0 for Gaussian white noise. White noise processes are often used to model random disturbances with a very small correlation period. However, the environments of many biological systems, due to the large number of interacting variables, may be considered to vary in a continuous manner. Thus, mean-reverting processes may be more appropriate for modeling the parameters.
Other differences between a mean-reverting process and a linear function of white noise for modeling environmental variability are examined in the following discussion.
Consider the per capita birth rate b(t). The per capita death rate d(t) can be treated in an analogous manner. A possible model for the birth rate in a randomlyvarying environment is a linear function of Gaussian WN: where W (t) is a standard Wiener process. A conceptual problem immediately occurs in that dW (t)/dt is not defined except in a generalized sense. The per capita birth rate, b(t) in (3), only exists in an abstract sense and, additionally, it is difficult to directly estimate the parameters γ and σ in (3) given data for per capita birth rate. For example, by directly integrating (3), the average per capita birth rate over an That is, the average per capita birth rateb over an interval of length T has a variance which goes to infinity as T → 0. Indeed, successive averages of a linear function of white noise experience large random oscillations. To see this, let be averages of per capita birth rate b(t) over successive time intervals of length ∆t.
and, for example, if σ = 1 and ∆t = 10 −6 , then That is, over 30% of successive average birth rates differ by more than 1400. (If instead ∆t = 10 −8 , then over 30% of successive average birth rates differ by more than 14000.) The question arises of whether it is reasonable for the average value over a time interval of a fundamental parameter, such as per capita birth rate, to become increasingly more variable as the time interval decreases.
A possible way to justify use of a linear function of Gaussian white noise is to study the effects it has on the population dynamics. Let y(t) be the population size at time t. Disregarding demographic variability and assuming no deaths in the population, equation (1) and (3) imply that y(t) satisfies the Itô SDE Using Itô's formula [9,17], it is seen that log(y(t)) satisfies the equation log(y(t)) = log(y(0)) + (γ − σ 2 /2)t + σ W (t) and so, log(y(t)) ∼ N log(y(0)) + (γ − σ 2 /2)t, σ 2 t .
That is, at any time t, the population size satisfies a log-normal distribution with E(log(y(t)) = log(y(0)) + (γ − σ 2 /2)t. Although the average of the logarithm of population size over an interval of length T has variance which goes to infinity as T → 0, a log-normal distribution may appear reasonable for a population undergoing environmental variability. Also, the parameters γ and σ can be estimated given data on the population size with time. Unfortunately, the per capita birth rate is still undefined except in a generalized setting. An important question still remains, specifically, how can (5) be biologically justified and, thus, how can (5) be relied on for providing insight into the dynamics of a population in a randomly varying environment.
A possible way to eliminate the conceptual and practical difficulties associated with linear functions of Gaussian WN is by considering the parameters in a randomly varying environment as mean-reverting processes. In this alternate approach, the parameters b(t) and d(t) satisfy separate SDEs. The SDEs are structured so that the per capita rates cannot drift off to infinity, i.e., they are forced back to asymptotic mean values. These SDEs are called mean-reverting and are commonly used, for example, in financial models. As above, consider in detail the per capita birth rate. Let b e be the asymptotic mean value of the per capita birth rate in the varying environment. When b(t) = b e , then the probability of moving closer to b e is greater than the probability of moving further away from b e . In this way, unrealistic values for the per capita birth and death rates are avoided. Indeed, discrete stochastic models can be derived that lead to mean-reverting processes [4,6]. In the following, several well-known and useful mean-reverting models are examined.
A classic mean-reverting process is the Ornstein-Uhlenbeck (OU) process which, for the per capita birth-rate, has the form: With a similar equation for d(t), a stochastic differential equation system is obtained for (2) of the form for y(t), b(t), d(t) ∈ [0, ∞) × R 1 × R 1 where W i (t), i = 1, 2, 3 are independent standard Wiener processes. (Three Wiener processes are necessary to model the demographic and environment changes.) System (7) represents a stochastic model of a single population experiencing both demographic and environmental variability. The stochastic differential equation (6) for b(t) can be solved exactly to yield Indeed, equation (8) implies that and, asymptotically with time t, it follows that the per capita birth rate b(t) is normally distributed with mean b e and variance σ 2 b /(2γ b ). Thus, it is inherently assumed in this stochastic model that random variations in the environment cause the per capita birth rate to vary about an asymptotic mean value b e for any large time t. An analogous result applies to the per capita death rate. Now consider the mean birth rate over an interval of size T assuming, for convenience, that b(0) = b e . Then, integrating (8) and simplifying, It follows that E(b) = b e and for small T , Thus, as T → 0,b → b e with zero variance. Unlike when the per capita birth rate is a linear function of white noise in (4), the variance goes to 0 rather than ∞ as T → 0. A disadvantage of the Ornstein-Uhlenbeck process is that it is possible for b(t) to have negative values. To avoid negative values, financial mathematicians often use a Cox-Ingersoll-Ross (CIR) [1,2] mean-reverting SDE model which, for birth rate b(t), has the form with b(0) = b 0 > 0. Applying Itô's formula, the mean and mean square in the birth rate for this model satisfy the ordinary differential equations . An important feature of the CIR model is that the birth rate, b(t), is nonnegative with probability one for any t ≥ 0. The probability density of the stochastic differential equation (10) satisfies [12,13]: and I q (z) is the modified Bessel function of order q. That is, where Γ is the gamma function. Using the identity ≥ 0 with probability one for any time t ≥ 0 for the CIR SDE model.
The OU and CIR models in equations (6) and (10), respectively, are important special cases of a general family of mean-reverting processes of the form where 0 ≤ ν ≤ 1. It is known, for example, that SDE (11) has nonnegative solutions when 1/2 ≤ ν ≤ 1 [1,2]. A third mean-reverting process, in addition to the OU and the CIR models, that is commonly applied but is not a member of SDE family (11) is a log-normal process. For the per capita birth rate, log(b(t)) satisfies the mean-reverting SDE This model is referred to as the Black-Karasinski (BK) model in mathematical finance [7]. Letting y(t) = log(b(t)) and y e = log(b e ), then y(t) satisfies the OU SDE dy(t) = γ b (y e − y(t))) dt + σ b dW 2 (t).
(13) Assuming that y(0) = y e and as equations (6) and (13) are identical in form, it follows from the argument for equation (6) that Since b(t) = exp(y(t)), then b(t) is log-normal for all t. For large t, the probability density of b(t) approaches a stationary log-normal density with mean and variance equal, respectively, to b e exp( 3. Properties of several mean-reverting processes. Some additional properties of the three mean-reverting processes, discussed in the previous section, are presented in this section. Denoting the three random variables for these processes as r i (t), i = 1, 2, 3, they satisfy the Itô SDEs: dr 2 (t) = γ 2 (r e,2 − r 2 (t)) dt + σ 2 (r 2 (t)) 1/2 dW 2 (t) (15) d log((r 3 (t)) = γ 3 (log(r e,3 ) − log(r 3 (t))) dt + σ 3 dW 3 (t) (16) where r i (0) = r i,0 , W i (t) are standard Wiener processes, and γ i and σ i are constants for i = 1, 2, 3. SDEs (14), (15), and (16) are the OU, CIR, and log-normal BK models, respectively. Several useful properties of these three mean-reverting processes are listed in Table 1. They possess several relevant properties for a biological parameter including continuity, reversion to a mean value, and approach to an asymptotic stationary distribution. In addition, the CIR and BK models are nonnegative. The stationary probability densities for these three SDEs can be derived from the formulas presented in Section 2 and are given, respectively, by For example, for the CIR model (15) noting from Section 2 that p(t, b) can be written as , v →v =ĉb, and which is the desired result. The stationary probability densities are illustrated in Figure 1 each with a mean of 10 and a variance of 12.5. Notice that the BK and CIR densities are both nonnegative.
The sample paths of these mean-reverting processes are influenced by the stationary means and variances as well as the correlation coefficient ρ(X(t), X(t + ∆t)) which is a measure of the dependence of the sample path on the previous position. Assuming at time t that X(t) satisfies the asymptotic probability density and estimating X(t + ∆t) for small ∆t using the Euler-Maruyama [14] approximation, the correlation coefficient satisfies 1 − ρ(X(t), X(t + ∆t)) ∝ γ i ∆t. This implies that the dependence of the sample path on the previous position decreases as γ i increases. This property is illustrated in Figure 2 for the CIR process (15) where the stationary mean and variance are kept constant at 100 and 10, respectively, but γ 2 increases from 5 to 20 to 80 in the three plots. (In the calculations, r e,2 is fixed and σ 2 changes with γ 2 to keep the stationary variance constant.) This figure indicates the flexibility of mean-reverting processes for modeling the variability in parameters as a result of a varying environment. 4. Comparison between environmental variability hypotheses. In this section, it is investigated how the population dynamics are affected by the assumption made about the parameters in a varying environment. Specifically, the dynamics of three simple population models are compared where, in these models, the parameters are assumed to be either linear functions of white noise or mean-reverting  Table 1. Several Properties of the Mean-Reverting SDEs (14), (15), (16) processes. The first model is a simple death model where the variance in population size is compared for the two assumptions on environmental variability. In the second model, carrying capacity varies due to environmental variability and population size is compared. In the third model, growth rate varies due to environmental variability and persistence time is compared. These examples illustrate that the population dynamics can differ significantly for WN or mean-reverting hypotheses. 4.1. Death process. The simplest example, which illustrates a difference in the population dynamics between the two hypotheses, is the elementary death process where per capita death rate d(t) is randomly varying due to the environment. Assuming that d(t) is a linear function of white noise, specifically, then population size y(t) satisfies the SDE dy(t) = −d e y(t) dt − σy(t) dW (t), y(0) = y 0 .
Next, a CIR mean-reverting process for d(t) is assumed of the form For model (17), with d(t) satisfying the mean-reverting SDE (21), population size y(t) satisfies y(t) = y 0 exp − t 0 d(s) ds . For this model, the mean E(y(t)) and variance Var(y(t)) are bounded for all time t as d(t) ≥ 0 with probability one.

EDWARD ALLEN
To illustrate a difference in the population dynamics between the two hypotheses with respect to the form chosen for d(t), notice that if the mean time-averaged death rate d e is less than σ 2 /2, then the variance Var(y(t)) → ∞ as t → ∞ for the linear function of white noise while the variance is bounded for the mean reverting process.

4.2.
Population carrying capacity. In this subsection, the simple population model dy(t) = a(N (t) − y(t)) dt (22) is studied where N (t) is considered the carrying capacity of the population. It is assumed in this subsection that N (t) is randomly varying due to the environment and parameter a is constant. To see how environmental variability affects the dynamics, demographic variability is not treated. First, it is assumed that N (t) satisfies a linear function of Gaussian WN: where α and β are constants and W (t) is a standard Wiener process. Substituting equation (23) directly into equation (22) results in the Itô SDE for population size y(t) of the form dy(t) = a(α − y(t)) dt + aβ dW (t).
Second, it is assumed for the same population problem (22) that N (t) satisfies an OU mean-reverting process. (The OU mean-reverting process is considered here but computations indicate that a CIR or BK mean-reverting process gives similar results.) In this case, the population model with environmental variability is the where N (0) = α, y(0) = y 0 , and a, α, r 1 , and r 2 are constants. System (28) is readily solved for N (t) and y(t) to obtain and N (t) = α + exp(−r 1 t) t 0 r 2 exp(r 1 τ ) dW (τ ).
The two different hypotheses about how N (t) randomly varies, either as an OU or a WN process, makes surprisingly little difference in the dynamics of the population size. Both models for environmental variability in N (t) give a stochastic differential equation system for y(t) where y(t) is continuous with mean and variance approaching constant values as t increases. That is, normal stationary distributions are asymptotically approached with time t under either hypothesis on the variability in N (t).
There are some very important differences, though, between the two hypotheses. The first difference is in the interpretation of the parameter a in model (22). In particular, the variances of the two hypotheses are much different with respect to parameter a. As a increases, the variance of y(t) for any time t increases without bound in the WN model while in the OU model, the variance approaches a constant value as a increases. In other words, how parameter a is viewed in (22) is much different for the two hypotheses.
Secondly, suppose now that both a and N (t) are varying in the environment. Then, it is impossible to model both a and N using linear functions of white noise. If a and N are similarly treated with linear functions of white noise, the equation for y(t) would have undefined squared terms of white noise. However, the meanreverting models have no such restriction and a can be readily modeled in the same manner as N (t), specifically, using a separate SDE for parameter a.

4.3.
Population growth rate and persistence time. In this subsection, a simple persistence time problem is examined where the population size y(t) satisfies with y(0) = y 0 and where r(t) is the per capita population growth rate. Per capita population growth rate is a fundamental biological parameter which depends on the species and the environment. In a randomly varying environment, r(t) varies with time t and experiences random changes from the environment. As in the previous subsection, to study how environmental variability affects the dynamics, demographic variability is not considered. The persistence time, defined as the time until the population size decreases below unity, is studied. Two forms for the variability in r(t) are hypothesized. First, a linear function of WN is assumed where r(t) has the form and where W (t) is a standard Wiener process with constants α < 0 and β > 0 and so, y(t) satisfies the Itô SDE dy(t) = αy(t) dt + βy(t) dW (t).
Second, an OU mean-reverting Itô SDE is assumed for r(t). Specifically, r(t) and y(t) satisfy    dr(t) = γ(r e − r(t)) dt + β dW (t), with r(0) = r e where W (t) is a standard Wiener process and γ > 0, r e < 0, and β > 0 are constants. The two models for r(t) differ in that the mean-reverting model (34) is continuous in time t with probability one. However, the two models for r(t) have several similar features. For example, lettingr = 1 T T 0 r(s) ds be the time-averaged growth rate, then E(r) = α and Var(r) = β 2 /T for (32) while E(r) = r e and Var(r) = for (34). Indeed, the mean time-averaged growth rate E(r) is constant and Var(r) decreases to zero as T increases for either the OU or WN hypothesis.
The mean persistence time is examined for these two models. Let T (z) be the persistence time and u 1 (z) = E(T (z)) be the mean persistence time where the population's initial size is z. For the WN model, where population size y(t) satisfies equation (33), the mean persistence time satisfies the boundary-value problem [3,5]: with u 1 (1)=0 and u 1 (∞) = 0. Solving the differential equation gives which is finite for α < 0. The second problem (34) is more difficult to analyze for persistence time. However, it can be simplified by letting z(t) = t 0 r(s) ds. Then, and z(t) satisfies the Itô SDE where α(t) = β 2 γ 2 (1 − exp(−γt)) 2 . As exit occurs when population size y(t) reaches unity, or y 0 exp( t 0 r(s) ds) = 1, then exit occurs at the identical time when z(t) = − log(y 0 ). Thus, an equivalent problem for the persistence time for the SDE system (34) when y(0) = y 0 and exit occurs at y(t) = 1 is to study the persistence time for the scalar SDE (37) when z(0) = 0 and exit occurs at z(t) = − log(y 0 ).
Persistence-time problem (37) is a diffusion process with drift for one absorbing boundary. For example, with zero drift, i.e. r e = 0, the mean persistence time is infinite. Also, for example, with α(t) equal to a constant, the mean persistence time is exactly equal to − log(y 0 )/r e . For problem (37), it is useful to consider the reliability function R(x, t) which is equal to the probability that the exit time is greater than t assuming that z(0) = x. The mean persistence time is equal to ∞ 0 R(x, t) dt and the reliability function R(x, t) satisfies the equation where R(x, 0) = 1 for x > − log(y 0 ) and R(− log(y 0 ), t) = 0. By examining R(x, t), a lower bound on the mean persistence time can be shown to be proportional to − log(y 0 )/r e . (Another way to estimate a lower bound is by considering the mean persistence time for a problem with two absorbing boundaries, for example, at − log(y 0 ) and −1/r e .) In particular, as r e → 0 and a pure diffusion process is approached, the mean persistence time for (37) goes to infinity. Correspondingly, the mean persistence time for (34) goes to infinity as r e → 0. Summarizing, it is observed that the growth rate models OU and WN give different interpretations with respect to how the parameters influence the mean persistence time. Let E(r) equal the average growth rate over a large time interval T . Then, as noted earlier, E(r) = α for WN and E(r) = r e for OU. For WN, the mean persistence time goes to a constant value, 2 log(y 0 )/β 2 , as E(r) → 0, whereas for OU, the mean persistence time goes to ∞ as E(r) → 0. The interpretation of mean persistence is different for the two different hypotheses, WN or OU.
5. An example of environmental variability. To test and compare the SDE models in the present investigation, data directly related to a biological problem would be useful. Unfortunately, it is very difficult to find appropiate data of sufficient quality and quantity to compare SDE models for biological populations undergoing evironmental variability. However, there is a plethora of accurate financial data and this data is useful in studying how populations (companies) interact in the financial environment [16,22,24]. Indeed, the financial market is a highlyvarying environment where companies undergo population-type processes. Several researchers have demonstated the applicability and usefulness of modeling companies as populations undergoing competition, growth, and death [16,22,24]. In the present section, the OU and WN hypotheses are tested on data obtained from the financial market.
Gold mining company stock prices are considered as competitive populations in an environment where gold price acts as an important environmental variable. Consider, for example, Barrick Gold Corporation (NYSE:ABX). Barrick Gold Corporation is one of the world's largest gold mining companies with mines and development projects on four continents. The stock price of Barrick Gold Corporation is correlated strongly with gold price. Consider, for example, annual ABX stock price along with gold price over the nine-year period 2002-2010 summarized in Table 2. A simple model for Barrick Gold Corporation stock price as a function of gold price is B ≈ 1.4G 1/2 where B is the gold corporation stock price. Thus, gold price acts as an environmental variable for ABX stock price and a possible stochastic model for ABX stock price B(t) is the Itô SDE where G(t) is gold price and α, β, ν and γ are constants.

Date
Gold Consider how gold price behaves dynamically as an environmental variable. Specifically, does gold price behave as a mean-reverting stochastic process or as a linear function of white noise? Two models for gold price are studied for the London afternoon gold prices on 497 trading days in 1994 and 1995, specifically, an OU mean-reverting process and a linear function of WN. These gold prices are plotted in Figure 3. The OU mean-reverting stochastic process for gold price is The linear function of white-noise is By integrating each equation over a time interval of one day, the equations are fit to the daily price data set. Specifically, letting t i+1 − t i = h = 1, integrating equation (39) over [t i , t i+1 ], and approximating the integral on the right-hand side leads to the Euler-Maruyama approximation [14,15]: where η i are normally distributed numbers with zero mean and unit variance. The gold price, G i , for i = 1, . . . , 497 are the data values. The transition probability of G i+1 given G i is normally distributed, i.e., p(G i+1 |G i ) ∼ N(G i + ha(G e − G i ), b 2 h).
In this case, p(G i+1 |G i ) is normally distributed with mean a and variance b 2 /h and the maximum likelihood estimates for a and b 2 /h in equation (40)   A problem with the linear function of white noise (42) is the sensitivity and consistency of the time interval interval h to the estimation of parameters a and b. For example, if h = 4, over an average of four days rather than one day, then b increases from b = 4.88 to b = 9.31. A second problem with the linear function of white noise is the lack of fit with the data. One sample path of the OU and WN models, (41) and (42) respectively, were computationally simulated using the fitted parameter values. The sample paths are presented in Figure 3 along with the actual gold price data. Figure 3 clearly shows that the continuous mean-reverting process is a more realistic model for gold price than a linear function of white noise. 6. Summary and conclusions. The properties of three mean-reverting stochastic processes, defined by Itô SDEs, are summarized and compared to a linear function of Gaussian WN. The mean-reverting processes possess several important features that better characterize environmental variability in biological systems than does a linear function of white noise. Most importantly, mean-reverting processes are conceptually and biologically realistic. Other advantages of mean-reverting processes over a linear function of Gaussian WN are continuity, nonnegativity, practicality, possession of asymptotic distributions, and ease of fitting the parameters to environmental data.