Table Of Content

Jump Locations of Jump-Diffusion Processes with State-Dependent Rates Christopher E. Miles∗ and James P. Keener Department of Mathematics, University of Utah 155 South 1400 East, Salt Lake City, Utah 84112, USA (Dated: January 10, 2017) We propose a general framework for studying jump-diffusion systems driven by both Gaussian noiseandajumpprocesswithstate-dependentintensity. Ofparticularnaturalinterestarethejump locations: thesystemevaluatedatthejumptimes. However,thestate-dependenceofthejumprate provides direct coupling between the diffusion and jump components, making disentangling the two to study individually difficult. We provide an iterative map formulation of the sequence of 7 distributions of jump locations. Computation of these distributions allows for the extraction of the 1 interjump time statistics. These quantities reveal a relationship between the long-time distribution 0 of jump location and the stationary density of the full process. We provide a few examples to 2 demonstrate the analytical and numerical tools stemming from the results proposed in the paper, n including an application that shows a non-monotonic dependence on the strength of diffusion. a J 8 I. INTRODUCTION modifying the state, which jump-diffusion allows. Diffu- sion with switching, which bares resemblance to jump- h] Stochastic models driven by both Gaussian noise (dif- diffusion, has also been studied [18] but is distinctly dif- c fusion)andadiscretejumpcomponent,orsocalledjump- ferent in the role of discrete noise and not addressed in e diffusion models have been used in a wide variety of ap- this work. m plications. Perhapsmostcommonly,thesejump-diffusion It is worth mentioning the issue of blowup in these - models are used in mathematical finance [2–5] to de- models. That is, for some choices of the diffusion and t a scribe frequent small transactions with occasional larger jump components, the system will experience an infinite st ones. Jump-diffusiondescriptionshavealsoserveduseful number of jumps or blowup (or both) in finite time, or . in neuroscience, typically describing the noisy buildup is said to be unstable. Although thoroughly studied[19– t a of voltage in a neuron leading to stochastic firing [6– 21] and shown subtle even for relatively simple models, m 9]. Other applications include biophsyical descriptions stability is deliberately not focused upon in this work. - of movement of chromosomes [10] and the interaction However, in the appendix of this work, we present some d between soil moisture and rainfall events [11]. Also rele- thoughtsonhowtheproposedframeworkcouldprovidea n vant to the results described in this paper is a subset of unique perspective in studying these phenomena, as this o c jump-diffusion models that neglect Gaussian noise and is primarily a partial differential equations (PDEs) for- [ are driven by the dynamics of the jump process alone. mulation rather than the existing path-wise approaches. This type of models appears in an equally eclectic vari- If a jump process is state-dependent, then the distri- 1 v ety of applications [12–15]. bution (or statistics) of the jump locations can be stud- 0 Of particular interest in this work are jump-diffusion ied. That is, we refer to the locations as the state of the 2 processeswithstate-dependentjumprates,whichappear system evaluated at the jump times. For instance, this 9 in some previous work [13, 14, 16]. This is a natural coulddescribethestatisticsofthevoltageatwhichaneu- 1 generalization, as any realistic system is unlikely to have ron fires or the price at which a jump occurs in finance. 0 a constant jump rate, but rather one that depends on This seemingly natural quantity of interest is seemingly . 1 thecurrentstateofthesystem. Forinstance,aneuronis unstudied and the main focus of this work. Although 0 more likely to fire as the voltage increases or, in finance, efficient Monte Carlo methods exist for jump-diffusion 7 an asset to crash as it rises in price. and Lévy processes [22, 23], the jump component may 1 Previous work has thoroughly investigated jump- be rare and therefore difficult to find an accurate distri- : v diffusion models for constant jump rates [1, 11]. The butionof. AfullPDEdescriptionoftheprocesscanalso Xi interjump time distributions for discrete noise alone has becomputed, butitisseeminglydifficult(orimpossible) alsobeenstudied[13,16]primarilyusingrenewaltheory. to disentangle the jump component alone. r a However, in general, state-dependent jump-diffusion is In this work, we present a general (accommodating a not arenewalprocess,assequentialdistributionsarenot wide variety of models) framework for studying jump- independent. The theory of Cox processes [17] is also diffusion systems and particularly focus on studying the well-established and refers to a diffusion process driv- sequenceofjumplocations. Wepresentaniterativemap ing a state-dependent Poisson process, but this does not formulation explicitly describing the distribution of the complete the feedback loop of the jump process further ith jump location, which allows for statistics to be ex- tracted about the interjump times even with diffusion, a limitation of previous works [13]. By assuming the processreachesstationarity,wefindidentifytermbehav- ∗ [email protected] ior of the distribution of jump locations and statistics 2 of the interjump times. The stationary jump locations B. Diffusion Component turn out to be closely related to the full stationary process, being equal if and only if the system is trivially Consider the classical one-dimensional SDE driven by state dependent. A few examples are provided, illustrat- Brownian motion described by ing the usefulness of the theory toward analytical and (cid:112) numerical results. One example provides an interesting dY =A(Y )+ 2D(Y )dW , (2) t t T t non-monotonic effect of diffusion on the interjump time statistics in a particular model. where W is a Brownian motion. The Fokker-Planck t (forwardChapman-Kolmogorov)equationdescribingthe evolution of the probability density function p(y,t) corresponding to the path-wise description (2) is II. FORMULATION ∂ p(y,t)=Lp:=−∂ {A(y)p}+D(y)∂ p. t y yy WeusetheFokker-Planck operator Ltocharacterizethe A. Jump-Diffusion diffusion process for the remainder of the paper. Let X denote a state-dependent jump-diffusion pro- t cess and p˜(x,t) be the probability density (PDF) of this C. Jump Component process. The evolution of the density is described by the forward Champman-Kolmorgov equation [24, 25] LetZ beaninhomogeneous(state-dependent)Poisson t process with intensity (rate) λ(X ). Also necessary is a t description of what actually occurs at the jumps, some- ∂ p˜(x,t)=Lp˜−λ(x)p˜+Jλ(x)p˜, (1) t times referred to as the reset map [20]. The behavior of the reset map is characterized by the jump operator J, where the operator L characterizes the diffusion compo- which is a probability density flux, ensuring that (1) in- nent, λ(x) is the state-dependent intensity at which the deeddescribestheevolutionofaprobabilitydensity. We jump process occurs, and the operator J describes the now describe a few possible choices of this operator. jump. Afeatureweseektoemphasizeisthestatedepen- dence of the jump-rate λ(x), resulting in direct coupling between the diffusion and jump components of the pro- 1. Jump Operator Examples cess. 1. Constantjumpsize. Atthejumptimes,theprocess increments by a fixed amount ∆, so 0.8 10 X →X+∆. 0.6 8 ] 1− Then, the corresponding jump operator J is then st] 0.4 me a shift by that fixed quantity X [dit 0.2 46 X) [tit Jp:=p(x−∆). ( λ 0.0 2 2. Reset. Fundamentally different than the previous example, rather than jumping by a fixed displace- 0.2 0 ment,theprocessresetstoadeterministicquantity, 0.0 0.2 0.4 0.6 0.8 1.0 η, so t [time] X →η. FIG. 1. An example realization of a jump-diffusion process X with Lp = ∂ {axp}+Dp , λ(x) = αexp(cid:8)−x2/β(cid:9) and t x xx The corresponding J operator is Jp=p(x−∆). Alsoshownisthestatedependentjumprate λ(x). Inthisparticularrealization,3jumpsoccur. Parameter (cid:90) ∞ values used a=.1,D=1,α=10,β =10,∆=1. Jp:=δ(x−η) p(x,t)dx, −∞ whereδ(x)istheDiracdeltafunction,againnoting Although the particular change choices of L and J that the scaling ensures flux preservation. significantly change the behavior and characteristics of the process, the description (1) is flexible enough to ac- 3. Random jump size. A generalization of the first commodate a wide variety of models, which we will now example, the particle can also jump a random dis- elaborate upon. placement, where the size of the jump is described 3 by the probability density µ(∆,x), which can also which means the first jump time density p (t) is then τ be state-dependent. The J operator is ds (cid:90) ∞ (cid:90) ∞ pτ(t)=−dt =− ∂tpdx. Jp:= p(x−∆,t)µ(∆,x)d∆. −∞ −∞ Then, using (3) and the fact that p→0 as x→±∞, we obtain 4. Other examples. Although not directly considered in this paper, reset maps of the form (cid:90) ∞ (cid:90) ∞ p (t)= λpdx= ∂ qdx. τ t X →γX =⇒ Jp:=p(γx) ∞ ∞ Duetothelackofspatialfluxinq(x,t),thedensityofexit have been studied in other work [20] and can be locations is q(x,∞), which can obtained by integrating formulatedinourframeworkinthedescribedman- over all possible jump times ner. A generalization of the reset case could also be made to allow for resets to a random location, (cid:90) ∞ p (x)= ∂ qdt=q(x,∞)−q(x,0)=q(x,∞) but this does not seem to provide additionally in- j t 0 teresting structure. by noting that q(x,0)≡0. It is worth noting that in the reset case, the system is indeed a renewal process, as the distributions between Although this provides useful information about the jumps are entirely independent, meaning that the tools distributions (in time and space) at which a particular of [11, 13, 16] can be utilized. This reset case has also (thefirst)jumpoccurs, weareinterestedinstudyingthe been studied more explicitly in [26, 27]. distributions of all jumps. III. RESULTS B. Jump Location Sequential Mapping A. Jump Distributions We instead study the interjump dynamics by noting thatbetweenjumps,thefullprocess(1)canbedescribed by the absorbing process (3). Define t to be the ith A more convenient process to study is what we refer i jump time and let X := X be the ith jump location. to as the absorbing process, described by i ti Consider the iterative definition of p (x) for i=2,3,... i (cid:40) ∂ p(x,t)=Lp−λp  ∂tq(x,t)=λp. (3) ∂tpî(x,t)=Lpî−λpî t pˆ(x,0)=Jp (4) i i−1 p (x)=(cid:82)∞λpˆ dt. Notethatwedistinguishbetweenthisandthefullprocess i+1 0 i p˜(x,t),whichwerefertohenceforthasthereinjectedpro- Thereisaslightclumsinesswiththefirstiterationofthis cess. Whiletheabsorbingprocessdoesnotcaptureallof procedure. Let p be some known starting distribution the behavior of the full reinjected process, the absorbing 0 of the process, then the first iterate becomes process fixes the distribution (in both space and time) of jump events and consequently serves more fruitful in  ∂ pˆ (x,t)=Lpˆ −λpˆ disentangling the jump component from the full process.  t 1 1 1 pˆ (x,0)=p 1 0 Theorem 1. The densities of the first jump location p (x)=(cid:82)∞λpˆ dt, 1 0 1 p (x) and first jump time p (t) of (3) are described by j τ noting the difference in the starting conditions. (cid:90) ∞ p (x)=q(x,∞)= λ(x)p(x,t)dt j Theorem 2. The distribution of the ith jump location 0 X is p (x) and is described iteratively by (4). ti i and ThisconstructionfollowsfromdirectlyfromTheorem (cid:90) ∞ 1. Effectively,pˆ servesasintermediatequantitytracking i p (t)= λ(x)p(x,t)dx. τ thedistributionofallpossiblejumplocations(andtimes) −∞ for that particular iterate. Then, to start the next iter- Proof. Define the survival probability s(t) be the proba- ate, the distribution of jump locations must be modified bility that the first jump has not occurred by time τ jump procedure characterized by J. Although tracking the jump locations rather than the jump times seems (cid:90) ∞ counterintuitive at first, no natural analogous formula- s(t):=P[τ >t]= p(x,t)dx, tion is apparent. −∞ 4 Define a more convenient quantity to study Elaborating a bit on this relationship: while (9) and (8)mayappeartoproducethesameresult,thesolutions u :=p (x)/λ(x). (5) i i toeachrequiredifferentscalings. Sincepˆ isaprobability s If λ(x) = 0 for some x, then necessarily p (x) = 0 since density, it must be that i λ(x) = 0 implies no jump can occur at this location, (cid:90) ∞ hence this quotient causes no difficulties. pˆ dx=1. (10) s Theorem 3. The description (4) is equivalent to −∞ Tui+1 :=[λ−L]ui+1 =Jλui, However, u(cid:63)λ=p(cid:63), where p(cid:63) is a probability density, so this means that or more explicitly (cid:90) ∞ ui+1 =T−1Jλui. (6) u(cid:63)(x)λ(x)dx=1. (11) −∞ Proof. Integratingbothsidesof(4)withrespecttotand noting that pˆ(x,∞)=0 we are left with i C. Moments (cid:90) ∞ (cid:90) ∞ pˆ(x,0)= Lpˆ dt−λ(x) pˆ. i i i 0 0 Define τ to be the ith interjump time. That is, i The linear operator L is a differential operator in x and τ := t −t . Also, let τ be the interjump time for i i i−1 (cid:63) consequently commutes with the time integral. Using the stationary process. (cid:82)∞ the initial condition and p /λ = pˆ dt, we obtain the i 0 i Theorem 5. The mean interjump time τ can be recov- desired result. i ered from the distribution of the ith jump location and The map (6) provides an explicit description for the is sequence (cid:90) ∞ (cid:90) ∞ p (x) u1(x)→u2(x)→···→u(cid:63)(x)→u(cid:63)(x)→··· , (7) (cid:104)τi(cid:105)= uidx= λi(x) dx. (12) −∞ −∞ where, we are assuming the process reaches stationarity. Also, the mean stationary interjump time can be com- Convergence of the sequence (7) corresponds to a fixed puted from the stationary distribution of the jump loca- point of the map (6). Although we seek not to focus on tions p this aspect, a natural question to ask is under what con- (cid:63) ditions the process reaches stationarity, or equivalently, (cid:90) ∞ (cid:90) ∞ p (x) when does u converge to some u . In Appendix B, we (cid:104)τ (cid:105)= u dx= (cid:63) dx. (13) i (cid:63) (cid:63) (cid:63) λ(x) provide a brief commentary about how the relationship −∞ −∞ (6) can be thought of as an iterated linear non-negative Proof. We integrate both sides of (4) with respect to x, integral operator which can be studied as such. Results again noting that pˆ →0 as x→±∞, resulting in from theoretical ecology literature are then cited which i provideaheuristicanalysisofconditionsforconvergence. (cid:90) ∞ (cid:90) ∞ Proceeding assuming the iterations converge to sta- ∂t pîdx=− λpîdx. tionarity,wecanassumethat(6)hasafixedpoint,which −∞ −∞ must be of the following form. However, from Theorem 1, we see that the right-hand Theorem 4. The stationary distribution of the jump lo- sideisexactlythedistributionoftheinterjumptimepτi, cations p (x) is described by so we have (cid:63) (cid:90) ∞ 0=Lu −λu +Jλu , (8) (cid:63) (cid:63) (cid:63) ∂ pˆ dx=−p (t) (14) t i τi where u :=p /λ. −∞ (cid:63) (cid:63) Corollary 1. The stationary distribution of the jump Taking the mean on both sides with respect to τi, locations p is the same as the stationary distribution of (cid:63) (cid:90) ∞(cid:90) ∞ the full process if and only if λ is constant. (cid:104)τ (cid:105)=− t∂ pˆ dxdt, i t i Proof. This is an immediate consequence of taking the 0 −∞ full process to be in stationarity, so dpˆ/dt = 0 in (1), (cid:82)∞ after integrating by parts and noting that pˆ dt = which satisfies 0 i p /λ, we get the desired result. i 0=Lpˆ −λpˆ +Jλpˆ . (9) s s s While the first moment (mean) of the interjump time This is exactly the same relationship as (8). Recalling can be computed directly with a quadrature of the pre- that λu = p and that both p and pˆ are both prob- sumed known p , the higher order moments are less (cid:63) (cid:63) (cid:63) s i ability densities, the only way that p = pˆ is if u is a straightforward. Appendix A describes how, in theory, (cid:63) s (cid:63) rescaling of p or, equivalently, λ is constant. knowledge of the p can be used to extract higher order (cid:63) i 5 moments of τ . Solving for higher order moments using where I,K are Bessel functions [28] such that I → 0 i this approach requires solving a hierarchy of differential as x → −∞ and K → 0 as x → ∞ (and are linearly equations, meaning in practicality, this may not be so independent). The jump condition (18) then becomes feasible. However, a more practical numerical approach may be to solve for pi and then solve the absorbing (3) D{c2u(cid:48)R(0)−c1u(cid:48)L(0)}=−1. (19) process run to extract explicitly p using the relation- τi However, we also need to impose continuity of u , so we ships from Theorem 1. (cid:63) also have the requirement c u (0)=c u (0). (20) IV. EXAMPLES 1 L 2 R Theconditions(19),(20)provideustwoequationsfortwo In this section, we briefly provide some examples to il- unknowns, which yield lustrate the usefulness of the proposed results. The first, a relatively simple case demonstrates the ability to de- 2 √ 2 √ c = 2K (2/ D), c = 2I (2/ D). rive analytical quantities for models which were difficult 1 D a/D 2 D a/D tostudypreviously. Thesecondexampleprovidesabrief Thus far, our solution is then defined only up to a con- analysis identifying that certain models may not be fea- stant, butweknowthatu mustsatisfythescaling(11), sible to study with Fourier analysis directly but become (cid:63) which our choice of c ,c serendipitously already satisfy. so when formulated in our framework. The third and 1 2 final example, which is studied numerically from the re- 0.5 0.7 lationships proposed in this work, provides interesting p (x) pˆ (x) s interactions between the state-dependence of the jump 0.4 Xti 0.6 Xt process and the strength of diffusion. 0.5 y y nc0.3 nc0.4 e e u u A. Constant drift with reset to origin freq0.2 freq0.3 0.2 0.1 0.1 Theresetcasehasbeenstudiedinotherworks[26,27], butwillserveheretoillustratetheabilitytoderiveexact 0.0 0.0 4 2 0 2 4 6 4 2 0 2 4 quantities utilizing the proposed framework. Consider a x [dist] x [dist] jump-diffusion model with constant drift and diffusion withjumpcomponenttoberesettozerowithrateλ(x)= FIG. 2. A comparison of Monte Carlo simulations of the ex- exp{x}. Note that the reset structure means that the ampledescribedby(15)andanalyticallyderivedcorrespond- ing quantities. Left: The stationary jump location density processisarenewalprocess,andconsequentlyp =p for i (cid:63) p (x)iscomparedwithahistogramofrealizationsofthestate all i. That is, the process is in stationarity after a single (cid:63) oftheprocessX evaluatedatthejumptimest . Right: The jump. Then, in stationarity, the fundamental relation in t i derived stationary distribution of the full process pˆ is com- this work (8) becomes s pared with the histogram of a realization of the process X . t (cid:90) ∞ Parameter values used are a=1,D=1. 0=−∂ {au }+D∂ xu −exu +δ(x) exu dx. x (cid:63) x (cid:63) (cid:63) (cid:63) −∞ By computing u explicitly, we can use (5) to immedi- (15) (cid:63) ately obtain the distribution of jump locations p . From From (11), the integral scaling the δ function is simply (cid:63) (13),wecouldalso,inprinciple,computestatisticsabout equal to 1, so the interjump times. Finally, from Corollary 1, we also δ(x)=−∂ {au }+D∂ xu −exu . (16) have the stationary density of the full process, pˆ , which x (cid:63) x (cid:63) (cid:63) s is a rescaling of q such that (10) is satisfied. We seek to This can be computed in a similar manner to a Green’s (cid:63) emphasizethatbeingabletocomputethedistributionof function, noting that for x(cid:54)=0, we have stationary jump locations or the stationary density itself 0=−∂ {au }+D∂ xu −exu , (17) (either analytically as done in this example or numeri- x (cid:63) x (cid:63) (cid:63) cally) allows for the extraction of the other, illustrating and integrating (16) from (−ε,ε) and using the fact that the usefulness and practicality of our results. u must be continuous, we get the matching condition (cid:63) In Figure 2, a comparison is shown between Monte D(cid:8)u(cid:48)(0+)−u(cid:48)(0−)(cid:9)=−1. (18) CarlosimulationsoftheprocessXtcorrespondingto(15) and the analytically derived distributions in this section. By solving (17), we get that our solution is of the form Specifically,anEuler-Maruyama[29]schemecorrespond-  (cid:16) √ (cid:17) ingtothepath-wisedescriptionwassimulated,recording c1uL(x):=c1ea2xbIa/b 2 bbex x<0 all of the values on the realization of the path as well as u(cid:63) =c2uR(x):=c2ea2xbK−a/b(cid:16)2√bbex(cid:17) x>0, tluhsetrvaatleusesthoaftXbotthatththeeanjuamlyptictaimllyesdetir.iveTdheexpfirgeusrseionils- 6 for the stationary jump locations p (x) and the station- partnershipwithFourieranalysis. Forthisexample,take (cid:63) ary distribution of the full process pˆ (x) agree with the aconstantdriftanddiffusionwithrandomjumpsizedis- s Monte Carlo simulations. tribution µ(∆). Then, in stationarity, (8) becomes B. Constant drift with random jump size Thisexampleispresentedonlybrieflytoillustratethat this framework can provide some insight to problems in (cid:90) ∞ 0=−∂ {au }+D∂ u −λ(x)u + µ(∆)λ(x−∆)u (x−∆)d∆. (21) x (cid:63) xx (cid:63) (cid:63) (cid:63) −∞ One possible quantity of interest is the first moment of However, recall from (12) that U(0) = (cid:104)τ (cid:105), the mean (cid:63) jump locations, so interjump time, so (cid:90) ∞ (cid:90) ∞ (cid:104)µ(cid:105) (cid:104)j(cid:63)(cid:105)= xp(cid:63)dx= xu(cid:63)λdx. (cid:104)τ(cid:63)(cid:105)=− a . (24) −∞ −∞ This result is quite interesting for two reasons. For one, Althoughλ(x)isleftarbitrary,theFouriertransformcan we have an explicit expression for the mean interjump surprisingly be used to gain some insight. Define the time revealed by the Fourier transform, independent of Fourier transform of u(cid:63) to be the choice of λ(x). More interestingly is the sign of the (cid:90) ∞ quantity. Since (24) represents an interjump time, it U(k):=F[u ]= e−ikxu (x)dx. must be non-negative, meaning this quantity only exists (cid:63) (cid:63) −∞ if a,(cid:104)µ(cid:105) differ in sign. This is intuitive as the drift and jump process must oppose each other to have a chance Also define the transformed quantities of reaching stationarity, which this demonstrates. Λ(k):=F[λu ], M(k):=F[µ] (cid:63) C. Ornstein-Uhlenbeck with Gaussian jump rate from which, we note that Λ(cid:48)(0)=−i(cid:104)j(cid:105), M(cid:48)(0)=−i(cid:104)µ(cid:105). This example serves to identify curious subtle behaviors of the coupling between diffusion and state- Taking the transform of (21), dependent jumps illuminated by our results. Consider, for this example, an Ornstein-Uhlenbeck diffusion with 0=−aikU −Dk2U −Λ(k)+M(k)Λ(k), (22) drift −ax and diffusion coefficient D. The jump rate is taken to be Gaussian which, evaluated at k =0 yields (cid:114) α2β λ(x)= exp{−βx2/2} (25) [1−M(0)]Λ(0)=0. (23) 2π However,weknowµ(∆)andu λ=p arebothprobabil- and the jump is a fixed distance ∆. This could perhaps (cid:63) (cid:63) ity densities and consequently M(0) = 1 and Λ(0) = 1. be thought of as a molecular motor [30] attached via a Thus, (23) is trivially satisfied. This is not seemingly Hookean linker, causing a relaxation to an un-stretched useful on its own, but we can couple the higher order x = 0 state. As the motor becomes less stretched, it moments by taking a k derivative of (22) to yield preferentially takes a step of size ∆ at a faster rate. Of possible interest would be to study how the cargo diffu- 0=−aikU(cid:48)−aiU −2DkU −2k2U(cid:48)−Λ(cid:48)+M(cid:48)Λ+MΛ(cid:48), sion coefficient D effects the effective rate at which the motor steps. The resulting stationary relationship (8) which, at k =0 and using M(0)=1,Λ(0)=1 yields becomes 0=−aiU −Λ(cid:48)(0)+M(cid:48)(0)+Λ(cid:48)(0), 0=∂x{axu(cid:63)}+D∂xxu(cid:63)−λ(x)u(cid:63)+λ(x−∆)u(cid:63)(x−∆). (26) which says that necessarily Theredoesnotappeartobeafruitfulapproachtostudy- ing(26)analytically. However,wesolveitnumericallyus- M(cid:48)(0)=−i(cid:104)µ(cid:105)=−aiU(0) ing a Crank-Nicolson upwind scheme [31]. From this, u (cid:63) 7 andconsequentlyp isobtained. Thefirstmomentofthe V. CONCLUSION (cid:63) interjumptime(cid:104)τ (cid:105)isimmediateafteraquadrature(13), (cid:63) but higher moments require solving differential equation Inthiswork,wehavepresentedamathematicalframe- solving, asdescribedby Appendix A.Insteadoftaking work for studying jump-diffusion systems with state- thisapproach,oncep isobtained,Jp isusedasastart- (cid:63) (cid:63) dependentjumprates. Theformulationisflexibleenough ingconditionfortheabsorbingprocess3,whichissolved to accommodate a variety of behaviors, providing rele- numerically using the same PDE solver. The advantage vance to a wide range of applications. The particular ofthisisthatthefulldistributionofτ isobtainedusing (cid:63) objectsoffocusinthisworkareisthedistributionofthe the relationship described by Theorem 1. jump locations: the values of X at the jump times t . t i We first reformulate the full process into another that isequivalentbetweenjumps. Usingthisformulation, the x10-1 explicit distribution of the jump locations can be ex- tracted and related to the previous jump location. This relation allows us to derive an iterated map description of the sequence of jump locations. With explicit knowledge of the distribution of jump locations, moments of the interjump times can be computed in more generality than previous work [13]. Assuming the process reaches stationarity, the distribution of the full process and the stationaryjumplocationsareshowntobecloselyrelated. A brief discussion of possible conditions on convergence x10-2 tostationarityinthelensofthisframeworkareprovided, utilizing tools from theoretical ecology motivated by the observation that the iterated map is effectively a non- negative integral linear operator. Three examples of particular choices of the diffusion and jump components are briefly discussed. The first showsthatmodelswithrelativelysimplefunctionalforms can be solved exactly using results from our work and likely would not be possible or difficult to derive other- wise. The second example shows that although studying full jump-diffusion models directly using Fourier space FIG.3. Thefirstandsecondmomentsofthestationaryjump may be infeasible, translating the problem to our pro- locations j and interjump times τ of the process described posed framework and then studying via the Fourier (cid:63) (cid:63) by (25),(26) as the strength of diffusion is varied, computed transform may serve useful. The last example is com- by solving (26) numerically. A non-monotonic dependence puted numerically from results derived in the paper and on diffusion is seen in the mean interjump time. Parameter shows that diffusion can have a non-monotonic influence values used are α=1,β =100. on when coupled with state-dependent jump processes, showing that our results are useful in finding subtle behaviors even when only numerical approaches are feasi- The results of this procedure as a function of the dif- ble. fusionstrengthD areshowninFigure3. Plottedarethe first and second moments of the stationary jump locations j and interjump times τ . A noteworthy behav- (cid:63) (cid:63) ior is observed: for certain parameter values, the mean VI. ACKNOWLEDGEMENTS interjump time is a non-monotonic function of the diffusion strength, D. Although this is perhaps counterin- This research was partially supported by NSF grant tuitive, discrete events with non-monotonic dependence DMS 1122297 and DMS-RTG 1148230. on the strength of diffusion have been observed in other biophysical phenomena [32]. The interjump times seem to change significantly for different values of the drift, the jump locations seem to change relatively little. In- Appendix A: Higher order interjump time moments tuitively, the variance of both quantities increases with diffusion strength and lower drift. Intheory,higherordermomentsoftheinterjumptimes This example provides an illustration that even when canbecomputedwithknowledgeofthesequenceofjump analyticalsolutionscannotbefound,theframeworkpro- distributions pi as was discussed with the first moment vides relationships that can be computed numerically to in (12). illuminatenewbehaviorsnotidentifiedinpreviouswork. To illustrate this, consider (14) and take the second 8 moment with respect to τ on both sides, resulting in Now,iterationsofourmapcorrespondtoiteratingthe i integral operator A with kernel G˜. This formulation is (cid:90) ∞ (cid:90) ∞ exactly that of the so-called integral projection models (cid:104)τ2(cid:105)=− t2∂ pˆ dtdx. i t i (IPM) in theoretical ecology. [33] Stability results from −∞ 0 the IPM literature can be used to understand the con- (cid:82)∞ vergence of our iterative procedure. Again, noting that u = p /λ = pˆ dt, after integrat- i i 0 i The main result comes from [34] and is effectively a ing by parts, we get statement of the Krein-Rutman theorem [35, 36], the in- (cid:90) ∞ finite dimensional analog of the Perron-Frobenius theo- (cid:104)τ2(cid:105)=− v dx, rem. [25, 37] i i −∞ Although L1 seems like the natural function space to study the spectral properties of these operators, since which looks similar to (12), however, vi satisfies they must preserve probability, L2 turns out to be far more accessible due to issues establishing compactness [λ−L]vi =ui. (A1) inL1. InAppendixCof[33],theauthorsprovideamore thorough commentary on these complications. In L2, we Thus, computation of the second moment requires solv- have the property ing the differential relationship (A1). Continuing to Theorem B.1 ([36]). If higher order moments using same approach yields a fur- therhierarchycoupledinadifferentialmanner. Interest- (cid:90)(cid:90) (cid:12) (cid:12)2 ingly, the differential operator (left hand side) of (A1) is (cid:12)G˜(x,ξ)(cid:12) dxdξ <∞, (cid:12) (cid:12) exactlyT,thesameas(6),meaningtheGreen’sfunction R2 (effectively T−1, discussed more in Appendix B) could then the integral operator with kernel G˜ maps from L2 to perhaps be reused to compute these quantities. L2. We then cite the main theorem which establishes the existence of a dominant eigenvalue for A. Appendix B: Spectral properties of iterated map Theorem B.2 (Easterling1998,[34]). SupposethatG˜ ∈ Note that (6) involves the inverse of a differential lin- L2 andisnon-negative, thenifthereexistsanα>0,β > ear operator T, which is exactly determined by its cor- 0,u such that 0 responding Green’s function. Let G(x,ξ) be the corresponding Green’s function to T, meaning that α(x)u (ξ)≤G˜(x,ξ)≤β(x)u (ξ) 0 0 for all x,ξ, then the integral operator A has a dominant T G(x,ξ)=δ(x−ξ). x eigenvalue with associated eigenfunction. Then, (6) can be rewritten as This establishes conditions for the existence of a dominant eigenvalue and eigenvector, which establish the (cid:90) ∞ long-term behavior. u = G(x,ξ)Jλu (ξ)dξ. i+1 i −∞ Theorem B.3 (Easterling1998,[34]). Assuming that A satisfies the previous condition, then the stationary dis- After a change of variables, define G˜ by tribution u is given by the eigenfunction φ associated (cid:63) 1 with the dominant eigenvalue λ , (cid:90) ∞ (cid:90) ∞ 1 u = G˜(x,ζ)u (ζ)dζ := G(x,ξ)Jλu (ξ)dξ. i+1 i i u −∞ −∞ lim i =κφ , i→∞λi 1 Forexample,forfixedjumpsizes: Jp=p(x−∆),then 1 where κ is a scaling parameter. G˜(x,ζ)=G(x,ζ+∆)λ(ζ+∆). The previous theorem suggests that the sequence u i converges if and only if |λ | = 1. Thus, this summary Abbreviate this linear non-negative integral operator 1 providesaroughheuristicandnovelangleindetermining (cid:90) ∞ whether the sequence ui (and the full process) approach Aq := G˜(x,ζ)q(ζ)dζ. stationarity, although this is not the focus of this work. −∞ [1] F. Hanson, Stochastic Processes and Control for Jump- (SIAM, Philadelphia, 2007). Diffusions, Advances in Design and Control, Vol. 13 9 [2] G.YanandF.B.Hanson,Am.ControlConf.2006,2989 arXiv:1303.6999v1. (2006). [22] Y. Xia and M. B. Giles, in Springer Proc. Math. Stat., [3] S. G. Kou, Manage. Sci. 48, 1086 (2002). Vol. 23 (2012) pp. 695–708. [4] S.G.Kou,inHandbooksOper.Res.Manag.Sci.,Vol.15 [23] P. Glasserman and N. Merener, Proc. R. Soc. A Math. (2007) pp. 73–116. Phys. Eng. Sci. 460, 111 (2004). [5] J. Figueroa-Lo´pez, Handb. Comput. Financ., edited by [24] C. Gardiner, Stochastic Methods: A Handbook for the J.-C. Duan, W. K. Ha¨rdle, and J. E. Gentle (Springer Natural and Social Sciences, 4th ed., Springer Series in Berlin Heidelberg, Berlin, Heidelberg, 2012) pp. 1–28. Synergetics, Vol. 13 (Springer Berlin Heidelberg, 2009) [6] W.GerstnerandW.M.Kistler,Spiking Neuron Models: p. 447. Single Neurons, Populations, Plasticity (2002) p. 494. [25] P. C. Bressloff, Stochastic Processes in Cell Biology, In- [7] L. Sacerdote and M. T. Giraudo, in Stoch. Biomath. terdisciplinary Applied Mathematics, Vol. 41 (Springer Model., Lecture Notes in Mathematics, Vol. 2058 International Publishing, 2014) p. 679. (Springer Berlin Heidelberg, 2013) pp. 99–148, [26] M.R.EvansandS.N.Majumdar,Phys.Rev.Lett.106, arXiv:1101.5539v1. 1 (2011), arXiv:1102.2704. [8] B. Lindner, B. Doiron, and A. Longtin, Phys. Rev. E [27] J.M.Meylahn,S.Sabhapandit, andH.Touchette,Phys. 72, 061919 (2005). Rev. E 92, 1 (2015), arXiv:1510.02431. [9] V.DiMaio,P.Lańsky´, andR.Rodriguez,Gen.Physiol. [28] M. Abramowitz and I. A. Stegun, Handbook of Mathe- Biophys. 23, 21 (2004). maticalFunctions,DoverBooksonMathematics,Vol.55 [10] B. Shtylla and J. P. Keener, J. Theor. Biol. 263, 455 (Dover Publications, 1966). (2010). [29] P. P. E. Kloeden and E. Platen, Stochastics An Int. J. [11] E. Daly and A. Porporato, Phys. Rev. E 73, 026108 Probab. Stoch. Process., Stochastic Modelling and Ap- (2006). plied Probability, Vol. 23 (Springer, 1992) pp. 1–7. [12] Y. Wu and W. Q. Zhu, Phys. Rev. E 77, 1 (2008). [30] J. Howard, Mechanics of Motor Proteins and the Cy- [13] E. Daly and A. Porporato, Phys. Rev. E 75, 1 (2007). toskeleton (Sinauer Associates, 2001) p. 367. [14] P. D’Odorico, F. Laio, and L. Ridolfi, Am. Nat. 167, [31] R. J. LeVeque, Finite Difference Methods for Ordinary E79 (2006). and Partial Differential Equations (SIAM, Philadelphia, [15] M. S. Bartlett, E. Daly, J. J. McDonnell, A. J. Parolari, 2007). andA.Porporato,Proc.R.Soc.AMath.Phys.Eng.Sci. [32] C.E.MilesandJ.P.Keener, (2016),arXiv:1605.06488. 471, 20150389 (2015). [33] S. P. Ellner and M. Rees, Am. Nat. 167, 410 (2006). [16] E. Daly and A. Porporato, Phys. Rev. E 74, 041112 [34] M.R.Easterling,The integral projection model: Theory, (2006). analysis and application, Ph.D. thesis, North Carolina [17] D. R. Cox, J. R. Stat. Soc. Ser. B 17, 129 (1955). State University (1998). [18] S. D. Lawley, J. C. Mattingly, and M. C. Reed, SIAM [35] M.G.KreinandM.A.Rutman,Am.Math.Soc.Transl. J. Math. Anal. 47, 3035 (2015), arXiv:1407.2264. 26, 128 (1950). [19] F. Xi, Stoch. Process. their Appl. 119, 2198 (2009). [36] M. A. Krasnosel’skii, Positive Solutions of Operator [20] L. DeVille, S. Dhople, A. D. Dom´ınguez-Garc´ıa, and Equations (P. Noordhoff, 1964). J. Zhang, SIAM J. Appl. Dyn. Syst. 15, 526 (2016), [37] J. P. Keener, Principles of Applied Mathematics: Tran- arXiv:1411.6600. formation and Approximation (WestviewPress,2000)p. [21] B. Cloez and M. Hairer, Bernoulli 21, 505 (2015), 603.