Trajectory composition of Poisson time changes and Markov counting systems

T rajectory composition of Poisson time changes and Marko v counting systems Carles Bret ´ o 1 Departamento de Estad´ ıstica and Instituto Flor es de Lemus, Universidad Carlos III de Madrid, C / Madrid 126, Getafe, 28903, Madrid, Spain Abstract Changing time of simple continuous-time Markov counting processes by independent unit-rate Poisson processes results in Markov counting processes for which we pro vide closed-form transition rates via composition of trajectories and with which we construct nov el, simpler infinitesimally over -dispersed processes. K eywor ds: Compartmental models, Subordination, Random time change, Compound process, Infinitesimal dispersion 1. Introduction The statistical analysis of dynamical systems plays an important role in scientific research. When these systems in volve counts, such analysis may be carried out using continuous-time Markov processes, which are often approxi- mations to real systems and hence fail to capture some features of real data. Luckily , some of these features may be better captured after replacing time in those processes by a Markov random time. Such a time randomization approach to improving statistical modeling was recently proposed and studied in detail in Bret ´ o and Ionides (2011), where the resulting Mark ov time-changed processes are defined via transition rates. Unfortunately , transition rates of such time- changed processes are in general una vailable and can be di ffi cult to obtain in closed form. This lack of closed-form transition rates limits the appeal of this time randomization approach and may e ven discourage applied researchers from using it at all. T o help make this approach more appealing, this paper considers changing by a Poisson process the time of a lar ge family that includes man y continuous-time Markov counting processes used in applications (for example from epidemiology , biochemistry or sociology). For the resulting time-changed process, transition rates are provided in the required closed form. These closed-form transition rates constitute the main result of this paper , which we obtain by composing trajectories of the counting process with those of the random time (instead of by integrating out the random time). Our choice of a Poisson time change seems to be unusual in the applied literature and produces time-changed models simpler than those previously considered, as we illustrate by constructing se veral novel ov er- dispersed counting processes, which can be used as building blocks to construct multiv ariate ov er-dispersed Markov counting systems. Dynamical systems that in volv e counts ha ve been studied in many disciplines by considering Mark ov counting systems without simultaneous ev ents, although compound systems (which allow simultaneity) hav e also receiv ed some attention. Fields where counting systems have been modeled as continuous-time Markov chains include epi- demiology and ecology (K ermack and McK endrick, 1927; Shrestha et al., 2011), pharmacokinetics (Haseltine and Rawlings, 2002; Matis and W ehrly, 1979; Sri vasta va, 2002) and engineering and operations research (Doig, 1957; Jackson, 2002). In these fields, most such processes are in fact Marko v counting systems (Bret ´ o and Ionides, 2011), mainly networks of queues (Br ´ emaud, 1999) or compartmental models (Jacquez, 1996; Matis and Ki ff e, 2000) that Email addr ess: carles.breto@uc3m.es (Carles Bret ´ o) 1 T el. + 34916245855; Fax: + 34916249848 Pr eprint submitted to Elsevier November 1, 2021 rule out the possibility of simultaneous transitions or events. When simultaneous e vents are possible, these counting systems are called compound (Bret ´ o and Ionides, 2011). Compound counting systems can capture better the v ariability in real data thanks to being infinitesimally o ver- dispersed and ha ve been constructed relying on random time changes and defined via closed-form transition rates, which is what this paper is concerned with. Compound Markov counting systems hav e been considered as a means to increase compatibility of theoretical models with real data, for example in the context of DN A sequence alignment and genomic data (Thorne et al., 1992) and of en vironmental stochasticity and epidemiological data (Bret ´ o et al., 2009). The y are also infinitesimally o ver-dispersed (Bret ´ o, 2012a), which is a model feature fav ored by infectious- disease data (Bret ´ o et al., 2009; He et al., 2009; Ionides et al., 2006; Shrestha et al., 2011). Such o ver-dispersion can be modeled with compound processes, which can be constructed under mild conditions (Bret ´ o, 2012b) via the well-known operation of random change of time. Such time randomization approach w as considered in detail in Bret ´ o and Ionides (2011), after being first illustrated in Bret ´ o et al. (2009) where a compound compartmental model was constructed and defined by transition rates e xpressed in closed form. In vestigating such closed-form rates for other models in general is our main concern. The problem that this paper addresses is the di ffi culty deriving closed-form transition rates of time-changed pro- cesses, which lies in the non-linearity of expected values of transition probabilities and which limits the appeal of time randomization in applications. T ransition rates of Markov counting processes can be understood as appropri- ate limits of transition probabilities (Br ´ emaud, 1999). These probabilities are most likely non-linear in time. After time is randomized, transition rates are instead determined by the expected transition probabilities (with respect to the randomized time), but such e xpected v alues need not be readily a vailable in closed form due to the non-linearity . Consider a unit-rate Poisson process whose time index t is changed by random time R ( t ). Its expected probability of k transitions o ver time interval [0 , l ] is the left hand side of (1) belo w , which is an analytic expression. A corresponding closed-form expression can be obtained, for example, assuming that { R ( t ) } is a gamma process with E  R ( t )  = t and V  R ( t )  = t /τ . Such expression is (Bret ´ o and Ionides, 2011; K ozubowski and Podg ´ orski, 2009), if Γ is the gamma function, E R h R ( l ) k e − R ( l ) i / k ! = Γ  l /τ + k  /  k ! Γ  l /τ  1 + τ  l /τ  1 + τ − 1  k  . (1) Howe ver , this closed-form e xpression is based on deri vations specific to the Poisson gamma process of this example and it need not extend straightforwardly to other random times or counting processes (like the non-linear death pro- cesses considered in Section 3). Seemingly technical di ffi culties like this one can pre vent applied researchers from randomizing time. The discouraging limitations imposed on time randomization by una vailable closed-form expressions stem from the necessity to define time-changed models as implicit hierarchies and the resulting complications on model interpre- tation, which we seek to alleviate in this paper . Implicit definitions of a model are those given in terms of numerical procedures to generate realizations or sample paths (Bret ´ o et al., 2009), e.g., for our Poisson gamma example abov e, a realization at time t w ould come from the following hierarchy of random draws: use a value drawn from a gamma random v ariable with mean t (and variance t /τ ) as the mean of a Poisson random v ariable from which to dra w the de- sired process realization. Implicit definitions like this one are all that is needed to do inference using “plug-and-play” methods (Bret ´ o et al., 2009), without the need to work out any closed-form expressions like (1). Ho wev er, implicit models can be harder to interpret. Consider what an applied researcher might ask when deciding what to make of and ho w to interpret results obtained from an implicit model: Is the time-changed model well-beha ved? What as- pects of the original model vary after changing time? Should interpretation of the original parameters change? How does the choice of time change a ff ect the answers to these questions? Interpretation issues like these might be tack- led considering implicit definitions only b ut answers may reflect numerical artifacts and may not be as apparent as with closed-form e xpressions, as illustrated in next section. Helping mitigate such interpretation issues to mak e time randomization more attractiv e is the ultimate goal of this paper . The main contribution of this paper is to provide closed-form transition rates for a large family of Markov counting processes time-changed by Poisson times via composition of trajectories and to illustrate how these closed-form rates facilitate the use of time randomization to improve Markov counting systems used in applications. The sought closed-form e xpressions are pro vided in Section 2 under mild requirements satisfied by many well-behav ed processes considered in the applied literature. These expressions provide the desired details about the time-changed process 2 to help address interpretation issues, promoting the use of time randomization. They are obtained by focusing on process trajectories to get around non-linear expectations like (1), for which the unusual Poisson time change turns out to be con venient. Not only does our Poisson time choice facilitate the deri vation of the expressions, it also a voids increasing the number of parameters, which results in a simpler interpretation of the parameters of the time-changed models. This is illustrated in Section 3, where time randomization is considered for sev eral processes of interest in the biological and social sciences and applied in detail to the widespread Poisson, linear birth and linear death processes to construct nov el, infinitesimally o ver-dispersed processes and multi variate systems without the interpretation issues of time-changed models defined only implicitly . 2. Closed-form transition rates of simple Mark ov counting processes time-changed by Poisson processes Before giving the closed-form time-changed rates in Theorem 1 below , we introduce our notation as we progress through some fundamental aspects of counting processes, random change of time and Poisson processes, which pro- vides a context for the theorem. Continuous-time Markov counting processes are fully characterized by transition rates, which are intimately re- lated to the process compoundness and intensity , and which will also be used to define our time-changed processes. W e denote a time-homogeneous continuous-time Mark ov counting process by { X ( t ) } and define it via its transition rates or transition semigroup local characteristics (Br ´ emaud, 1999), which we write as q X ( x , k ) ≡ lim h ↓ 0 P  X ( t + h ) = x + k | X ( t ) = x  h (2) where h , t ∈ R ≥ 0 , x ∈ N ≥ 0 and k ∈ N ≥ 1 . Transition rates determine whether { X ( t ) } is infinitesimally over -dispersed, which occurs if and only if simultaneous events are possible (Bret ´ o and Ionides, 2011)—i.e., if and only if { X ( t ) } is compound so that there exists at least one x and k ? > 1 with q X  x , k ?  > 0. If { X ( t ) } is not compound, then it is simple (Daley and V ere-Jones, 2003) and q X ( x , 1 ) is the only non-zero transition rate. T ransition rates are also related to the process rate function or intensity (Daley and V ere-Jones, 2003), which we write as λ X ( x ) ≡ lim h ↓ 0 1 − P  X ( t + h ) = x | X ( t ) = x  h . If the rate function summarizes the overall acti vity implied by the transition rates and is finite—i.e., if it satisfies λ X ( x ) = P k ≥ 1 q X ( x , k ) < ∞ for all x , we say that { X ( t ) } is conserv ativ e and stable (Br ´ emaud, 1999). Conserv ation, sta- bility and simpleness are satisfied by many well-beha ved processes used in applications (including those in Section 3 below) and are the only conditions that Theorem 1 imposes on the counting process to be time-changed. (Actually , time homogeneity is also imposed but it can be relax ed at the cost of complicating notation and deriv ations.) Random change of time only needs to be sketched at this point and details are postponed to Appendix A, b ut we do elaborate here on our reasons for choosing a Poisson time change over other time changes. In a nutshell, random time change inv olves a base process, say { X ( t ) } , whose time index t is replaced by a random time process, say { N ( t ) } , giving the time-changed process, say { S ( t ) } ≡ { X ( N ( t )) } . While the base process considered in this paper is rather general, i.e., simple conservati ve stable processes, our time change is v ery specific: a unit-rate Poisson process. The choice of Poisson time change is unusual in the literature but key if we wish to a void directly inte grating out the random time from the hierarchy X  N ( t )  by follo wing our alternati ve approach based on trajectory composition to prov e Theorem 1 belo w . In the time randomization literature, L ´ evy processes (Sato, 1999) are routinely considered as time changes. One such process that has recurrently been chosen to obtain closed-form expressions is the gamma process. Examples include time-changed di ff usions (Madan et al., 1998) or, in the context of counting processes, time-changed Poisson processes (Hougaard et al., 1997; K ozubo wski and Podg ´ orski, 2009; Lee and Whitmore, 1993) and time-changed pure-birth and pure-death processes (Bret ´ o and Ionides, 2011). In spite of processes other than gamma ha ving been considered (K umar et al., 2011; Marion and Renshaw, 2000; V arughese and Fatti, 2008), Poisson processes seem to ha ve been dispensed with. This literature emphasizes the random-variable hierarchy perspecti ve, which drags us to the problematic non-linear expected transition probabilities analogous to (1). This hierarchical 3 perspectiv e giv es the following analytical expression for the time-changed transition rates using (2) directly lim h ↓ 0 E N ( t + h ) " P  B  N ( t + h )  = s + k    B  N ( t )  = s  # h . (3) This expression is not in closed form. T o obtain one, the limit could be handled using T aylor expansions but the expectation is likely to be non-linear (as argued in the introduction). Such non-linear expectations can be circumvented by composing the trajectories of the processes (more on this in Appendix A). Trajectories are easier to compose, as we argue in Appendix A, if the time change is a Poisson process lik e in Theorem 1. Theorem 1 ( T ransition rates of simple Markov counting processes time-changed by Poisson processes ) . Let { X ( t ) } be a continuous-time Markov counting pr ocess that is time-homogeneous, simple, conservative, stable, and defined by r ate function λ X ( x ) = q X ( x , 1 ) < ∞ . Let { N ( t ) } be a P oisson pr ocess defined by λ N ( n ) = 1 . Let also both pr ocesses be independent. Then, the transition r ates and rate function of the time-c hanged pr ocess { S ( t ) } ≡ n X  N ( t )  o ar e q S ( s , k ) = P  X (1) = s + k    X (0) = s  (4) λ S ( s ) = 1 − e − λ X ( s ) . (5) Pr oof. Theorem 1 is prov ed in Appendix A to preserve the flo w of the paper . 3. Applications: making time change more appealing to model o ver -dispersed systems The closed-form expressions in Theorem 1 make time randomization more attracti ve to applied researchers by providing them with details about their univ ariate time-changed process { S ( t ) } that refer to the ov erall process behavior and to the role of specific parameters, as we illustrate below with di ff erent time changes and base processes (ev en when the distribution of X (1) in (4) is di ffi cult to obtain). This distribution only needs to be considered for univ ariate processes, which can then be used as building blocks to construct infinitesimally ov er-dispersed multiv ariate systems of counts. The overall behavior of a process is a subjective concept, which could include features like whether or not it is well-behav ed, compound or whether its transition rates and rate function have any distinctive characteristics, all of which can be complemented by additional information on the role of specific parameters. The Poisson time-changed process { S ( t ) } of Theorem 1 is well-beha ved in the sense that its transition rates q S ( s , k ) ∈ [0 , 1] < ∞ and that it remains conserv ativ e and stable, since λ S ( s ) = 1 − P  X (1) = s | X (0) = s  = P k ≥ 1 q S ( s , k ) . This assures modelers that their time-changed model is not pathological and that it fits in the frame work of most existing Markov process theory . { S ( t ) } is also compound, which implies it is infinitesimally over -dispersed. Such over -dispersion informs applied scientists of the additional variability demanded by the data, which may also be interpreted as the base model { X ( t ) } being misspecified (Bret ´ o et al., 2009; Bret ´ o and Ionides, 2011). Compoundness also has implications for describing and summarizing { S ( t ) } . It suggests complementing the information giv en by the more con ventional infinitesimal mean and v ariance with that giv en by the rate function and the mean and variance of the transition rates (or jump sizes), which ha ve the follo wing distincti ve features in the case of { S ( t ) } . One on hand, moments of the transition rates in (4) are giv en directly by the distribution of random v ariable X (1). On the other hand, the rate function in (5) is necessarily restricted to the interval [0 , 1]. These two features clue in modelers that e vent size moments will inherit any constraint present in the distribution of X (1) and that the e vent rate function of { S ( t ) } is bounded from abov e. Moreov er, the closed form for the rate function giv en in (5) permits quantitatively assessing the impact on λ S ( s ) of parameters appearing in the original λ X ( x ). The role of specific parameters can be analyzed in more detail by using additional information on the distribution of X (1) gi ven X (0) to measure the impact of a parameter on transition rates, on the rate function and on moments, which in addition facilitates connecting and comparing the time-changed process to other processes, as we illustrate 4 in se veral examples below . The required distrib ution of X (1) is a ubiquitous result for Poisson base processes { X ( t ) } with rate α (independent of Poisson time change { N ( t ) } ), for which X (1) ∼ Poisson( α ). This distribution is all that is needed to apply Theorem 1 (from which we borrow the notation) to pro ve the follo wing corollary . Corollary 2 ( P oisson Poisson process ) . Consider Theor em 1 and let λ X ( x ) = α , i.e., let { X ( t ) } be a P oisson pr ocess with rate α ∈ R > 0 . The transition rates of { S ( t ) } are the pr obability mass function of a P oisson random variable parameterized by r ate α , i.e., q S ( s , k ) = α k e − α / k ! and its rate function is λ S ( s ) = 1 − e − α . Note 3 ( On Poisson Poisson pr ocesses ) . The name “P oisson P oisson pr ocess” has been used in the literatur e to r efer to a di ff er ent pr ocess, i.e., to r efer to compound P oisson pr ocesses with P oisson event size distrib ution (see for example Hanson, 2007). Corollary 2 gi ves a rate function and transition rates for { S ( t ) } parameterized solely by α . The impact of α on transition rates, the rate function and moments for { S ( t ) } is di ff erent than for { X ( t ) } . Before the time change, an increase in α linearly increases the rate λ X at which ev ents occur in { X ( t ) } and ev ents hav e in variably size one. After the time change, its impact on the ev ent rate λ S becomes non-linear and, more importantly , it changes the mean and variance of ev ent sizes of { S ( t ) } . These e vent size moments are simply those of X (1), i.e., both α and, although knowing them is informativ e, they di ff er from the infinitesimal moments of { S ( t ) } . These infinitesimal moments and their ratio (the dispersion index D dS ) can be deriv ed by noticing that { S ( t ) } is a compound Poisson process, so that µ dS ( s ) = lim h ↓ 0 E  S ( t + h ) | S ( t ) = s  h = lim h ↓ 0 α  1 − e − α  h h = α  1 − e − α  (6) σ 2 dS ( s ) = lim h ↓ 0 V  S ( t + h ) | S ( t ) = s  h = lim h ↓ 0 α (1 + α )  1 − e − α  h h = α (1 + α )  1 − e − α  (7) D dS ( s ) = σ 2 dS ( s ) µ dS ( s ) = 1 + α (8) The infinitesimal quantities (6)–(8) allo w further analysis of the di ff erence in the role of α in { S ( t ) } and in { X ( t ) } . In addition, they are also important for applied modelers interested in bounding simulation computational load, which is a k ey element when applying the plug-and-play inference methods pointed to in the introduction. Computational load can be bound by approximating { S ( t ) } at certain times using deterministic ordinary di ff erential equations (ODEs), as explored for example in Haseltine and Rawlings (2002). For e xample, the unit-rate Poisson time { N ( t ) } of Corollary 2 could be approximated by function n ( t ) = t , giving the approximating ODE dn ( t ) / d t = 1. Other approximating ODEs regarding Corollary 2 could be d x ( t ) / d t = α and d s ( t ) / d t = α  1 − e − α  , d x  n ( t )  / d t = α . The approximating ODE d s ( t ) / dt = α is correct ne vertheless if the time change used is instead the g amma process used to deriv e (1) (Bret ´ o and Ionides, 2011). The gamma time change used to derive (1) has disadvantages when compared to our Poisson time change, re- gardless of the base process. A fundamental disadvantage is that our approach to finding the closed-form e xpressions in Theorem 1 does not work, basically because the gamma trajectories are more complex (more details in Note 6 in Appendix A). A minor disadv antage is that it introduces an additional parameter , τ . This additional parameter complicates the preceding non-trivial comparison of the role of α in { S ( t ) } with its role in { X ( t ) } . Of course, τ makes { S ( t ) } more flexible, e.g., D dS = 1 + ατ (Bret ´ o and Ionides, 2011), and could always be fixed at some arbitrary v alue, say τ = 1, which would gi ve the same D dS for both the Poisson and the gamma time change. In this sense, the Poisson choice makes { S ( t ) } simpler and still infinitesimally ov er-dispersed. The Poisson time change could be replaced by a compound Poisson time change. This extension would allow for additional parameters that could play the role of τ , providing additional flexibility . It would also still permit proving Theorem 1 with some minor modifications (which are be yond the scope of this paper). The abov e analysis comparing the role of parameters before and after changing time and with di ff erent time changes is also possible for base processes other than Poisson. Base processes for which the distrib ution of X (1) is well-known also include linear birth and linear death processes (see for example Bharucha-Reid, 1960), for which we next provide closed-form rates and moments before considering applying Theorem 1 to processes for which this distribution is more di ffi cult to obtain. Corollary 4 ( Negative-binomial Poisson process ) . Let λ X ( x ) = β x I { x > 0 } , i.e ., let { X ( t ) } be a linear birth pr ocess with individual birth rate β ∈ R > 0 . The transition rates of { S ( t ) } ar e the pr obability mass function of a negative 5 binomial random variable parameterized by success pr obability 1 − e − β and number of failures until stopping s, i.e., q S ( s , k ) =  s + k − 1 k   e − β  s  1 − e − β  k for s ∈ N > 0 and its rate function is λ S ( s ) = 1 − e − β s I { s > 0 } . Corollary 5 ( Binomial Poisson pr ocess ) . Let λ X ( x ) = δ ( d 0 − x ) I { x < d 0 } , i.e ., let { X ( t ) } be the counting pr ocess associated with a linear death pr ocess with individual death rate δ ∈ R > 0 and initial population size d 0 ∈ N > 0 . The transition r ates of { S ( t ) } ar e the pr obability mass function of a binomial r andom variable parameterized by success pr obability 1 − e − δ and number of trials ( d 0 − s ) , i.e ., q S ( s , k ) =  d 0 − s k   1 − e − δ  k  e − δ  ( d 0 − s ) − k with k = 1 , . . . , d 0 − s for each s = 0 , . . . , d 0 − 1 and its rate function is λ S ( s ) = 1 − e − δ ( d 0 − s ) I { s < d 0 } . In both Corollary 4 and 5, the ev ent size moments are the moments of the corresponding ne gative binomial and binomial distrib utions. Infinitesimal moments are pro vided below for the linear death process of Corollary 5 only . (Providing them for Corollary 4 would require a justification to pass the h -limit inside the expectation, a technicality that at this point we prefer leaving for future research.) If { S ( t ) } is the process from Corollary 5, then µ dS ( s ) = d 0 − s X k = 1 k lim h ↓ 0 h − 1 P  S ( t + h ) = s + k    S ( t ) = s  = d 0 − s X k = 1 k P  X (1) = s + k    X (0) = s  = ( d 0 − s )  1 − e − δ  σ 2 dS ( s ) = d 0 − s X k = 1 k 2 P  X (1) = s + k    X (0) = s  + lim h ↓ 0 µ 2 dS ( s ) h 2 + o ( h ) h = µ dS ( s ) " 1 + ( d 0 − s − 1)  1 − e − δ  # D dS ( s ) = 1 + ( d 0 − s − 1)  1 − e − δ  These infinitesimal moments are again di ff erent from those of a binomial gamma process, but in both cases { S ( t ) } is infinitesimally equi-dispersed for s = d 0 − 1 (Bret ´ o and Ionides, 2011). Regarding approximating ODEs corresponding to Corollary 5, the y are similar to those of Corollary 2: d x ( t ) / dt = δ ( d 0 − x ) I { x < d 0 } and d s ( t ) / d t = ( d 0 − s )  1 − e − δ  I { s < d 0 } . All of the preceding analysis of the role of parameters has been greatly facilitated by nice and tractable distributions of X (1). When the distrib ution of X (1) is di ffi cult to obtain, the closed-form results provided by Theorem 1 inherit this lack of tractability , which, in the most extreme case, can force the use of approximations in both univ ariate and multiv ariate settings. An example of a distrib ution that is di ffi cult to handle is that of non-linear death processes with λ X ( x ) = x ( d 0 − x ) I { x < d 0 } , which is of current interest in the context of bacterial disinfection (Chou et al., 2005) and for which Theorem 3 in Billard et al. (1979) gi ves intricate closed-form expressions for P  X (1) = s + k | X (0) = s  . In this case, results can be obtained without the need of any approximation. Approximation of the distribution of X (1) is an option if this distribution is unav ailable in closed form. Such approximation can be done numerically but numerical approximations can be a voided at least for death processes with general death rate, including the follo wing non-linear death processes: logistic and exponential death processes (which are of interest to model in vasion of larv ae, Faddy and Fenlon, 1999); Bass death processes (to model innovation di ff usion in social groups, e.g., Karmershu and Sharma, 2007); and po wer death process (to model chemical formation of activ ated carbons, Fan and Argoti, 2011). For such processes, it is possible to truncate the infinite sum for P  X (1) = s + k | X (0) = s  provided in Theorem 2 in Billard et al. (1980). These numerical or truncation approximations should in any case be e xpected to gi ve more reliable results than empirically attempting to approximate the limit and expectation in (3) to obtain approximate time-changed transition rates. The transition rates considered so far refer to uni variate counting processes only b ut can also be used to construct multiv ariate Markov counting systems. Over -dispersed Marko v counting systems hav e been constructed combining the transition rates of Poisson g amma and binomial gamma processes as b uilding blocks of a multi variate counting system in the conte xt of epidemiological susceptible-infectious-recov ered (SIR) compartmental models (Bret ´ o, 2012a; Bret ´ o and Ionides, 2011). These blocks based on gamma time changes can be replaced by the Poisson-based blocks presented in this paper , resulting in alternativ e SIR compartmental models. Such alternati ve SIR models can then be considered as a more parsimonious alternativ e to the g amma-based models to test for o ver-dispersion (i.e., for model misspecification) and to improve the fit of the base model to data. 6 Acknowledgements This work was supported by Spanish Government Project ECO2012-32401 and Spanish Program J uan de la Cierva (JCI-2010-06898). References Barndor ff -Nielsen, O. E., Shiryaev , A., 2010. Change of time and change of measure. W orld Scientific Publishing. Bharucha-Reid, A. T ., 1960. Elements of the Theory of Marko v Processes and their Applications. McGraw-Hill. Billard, L., Lacayo, H., Langberg, N. A., 1979. A ne w look at the simple epidemic process. Journal of Applied Probability 16, 1, 198–202. Billard, L., Lacayo, H., Langberg, N. A., 1980. Generalizations of the simple epidemic process. Journal of Applied Probability 17, 4, 1072–1078. Br ´ emaud, P ., 1999. Markov Chains: Gibbs Fields, Monte Carlo Simulation, and Queues. Springer, Ne w Y ork. Bret ´ o, C., 2012a. On infinitesimal dispersion of multiv ariate Markov counting systems. Statistics and Probability Letters 82, 720–725. Bret ´ o, C., 2012b. Time changes that result in multiple points in continuous-time Markov counting processes. Statistics and Probability Letters 82, 2229–2234. Bret ´ o, C., He, D., Ionides, E., King, A., 2009. T ime series analysis via mechanistic models. Annals of Applied Statistics 3, 319–348. Bret ´ o, C., Ionides, E., 2011. Compound Markov counting processes and their applications to modeling infinitesimally o ver-dispersed systems. Stochastic Processes and their Applications 121, 2571–2591. Chou, S. T ., fan, L. T ., Argoti, A., V idal-Michel, R., More, A., 2005. Stochastic modeling of thermal disinfection of bacteria according to the logistic law . American Institute of Chemical Engineers Journal 51, 9, 2615–2618. Daley , D., V ere-Jones, D., 2003. An Introduction to the Theory of Point Processes. V olume I: Elementary Theory and Methods. Springer . Doig, A., 1957. A bibliography on the theory of queues. Biometrika 44, 490–514. Faddy , M. J., Fenlon, J. S., 1999. Stochastic modelling of the inv asion process of nematodes in fly larvae. Journal of Applied Statistics 48, 31–37. Fan, L. T ., Argoti, A., 2011. Stochastic modeling for the formation of activated carbons: Nonlinear approach. Industrial & Engineering Chemistry Research 50, 8836–8841. Hanson, F . B., 2007. Applied Stochastic Processes and Control for Jump-Di ff usions: Modeling, Analysis, and Computation (Advances in Design and Control). Society for Industrial and Applied Mathematics. Haseltine, E. L., Rawlings, J. B., 2002. Approximate simulation of coupled fast and slow reactions for stochastic chemical kinetics. Journal of Chemical Physics 117, 6959–6969. He, D., Ionides, E., King, A., 2009. Plug-and-play inference for disease dynamics: Measles in large and small populations as a case study . Journal of the Royal Society Interface 7, 271–283. Hougaard, P ., Lee, M.-L. T ., Whitmore, G. A., 1997. Analysis of overdispersed count data by mixtures of poisson v ariables and poisson processes. Biometrics 53, 4, 1225–1238. Ionides, E. L., Bret ´ o, C., King, A. A., 2006. Inference for nonlinear dynamical systems. Proceedings of the National Academy of Sciences of the USA 103, 18438–18443. Jackson, J., 2002. How netw orks of queues came about. Operations Research 50 (1), 112–113. Jacquez, J. A., 1996. Compartmental Analysis in Biology and Medicine. 3rd edition, BioMedware, Ann Arbor , MI. Karmershu, Sharma, P ., 2007. Truncating the hierarchy of moment equations based on point distrib ution—applications to innovation di ff usion. Mathematical and Computer Modelling 45, 233–240. Kermack, W . O., McKendrick, A. G., 1927. A contribution to the mathematical theory of epidemics. Proceedings of the Royal Society of London, Ser . A 115, 700–721. K ozubowski, T . J., Podg ´ orski, K., 2009. Distributional properties of the negati ve binomial L ´ evy process. Probability and Mathematical Statistics 29, 43–71. Kumar , A., Nane, E., V ellaisamy , P ., 2011. T ime-changed Poisson processes. Statistics and Probability Letters 81, 1899–1910. Lee, M.-L. T ., Whitmore, G. A., 1993. Stochastic processes directed by randomized time. Journal of Applied Probability 30, 2, 302–314. Madan, D. B., Carr , P ., Chang, E. C., 1998. The variance gamma process and option pricing. European Finance Re view 2, 79–105. Marion, G., Renshaw , E., 2000. Stochastic modelling of en vironmental v ariation for biological populations. Theoretical Population Biology 57, 197–217. Matis, J. H., Ki ff e, T . R., 2000. Stochastic Population Models. A Compartmental Perpecti ve. Springer . Matis, J. H., W ehrly , T . E., 1979. Models of compartmental systems. Biometrics 35 (1), 199–220. Nicolai, R., Frenk, J., Dekker, R., 2009. Modelling and optimizing imperfect maintenance of coatings on steel structures. Structural Safety 31 (3), 234–244. Sato, K., 1999. L ´ evy Processes and Infinitely Divisible Distrib utions. Cambridge University Press. Shrestha, S., King, A., Rohani, P ., 2011. Statistical inference for multi-pathogen systems. PLoS Computational Biology 7(8): e1002135. Sriv astav a, R., 2002. Stochastic vs. deterministic modeling of intracellular viral kinetics. Journal of Theoretical Biology 218, 309–321. Thorne, J. L., Kishino, H., Felsenstein, J., 1992. Inching to ward reality: An improv ed likelihood model of sequence evolution. Journal of Molecular Evolution 34, 3–16. V arughese, M., Fatti, L., 2008. Incorporating en vironmental stochasticity within a biological population model. Theoretical Population Biology 74, 115–129. 7 Appendix A. Before proving Theorem 1, we gi ve a formal definition of time change. Time change can be approached rigorously by defining an underlying probability space and by paying attention to the measure theory in volved (as in for example Barndor ff -Nielsen and Shiryae v, 2010; Sato, 1999). Such rigor is not needed for our proof but beginning with a probability space ( Ω , F , P ) clarifies it. On this probability space, the collection (indexed by t ∈ R ≥ 0 ) of random variables X t ( ω ) : Ω → N ≥ 0 defines the base counting process { X ( t ) } of Theorem 1. If instead of indexing by t , we index by ω (which we stress by switching to lowercase), x ω ( t ) is an N ≥ 0 -valued deterministic function called the process trajectory or sample path. Composing trajectory x ω ( t ) with that of Poisson process { N ( t ) } is equiv alent to defining the time-changed process { S ( t ) } as the collection of random variables such that, for all t , S t ( ω ) = X N t ( ω ) ( ω ) (A.1) (as in Sato, 1999) with trajectories s ω ( t ) = x ω  n ω ( t )  . Proof of Theor em 1 . The result is prov ed by inferring, from the properties of the composed trajectory , the con- ditional distribution of inter-e vent times in { S ( t ) } , which determines the rate function in (5), and the conditional distribution of e vent sizes in { S ( t ) } , which determines the transition rates in (4). A key property of the trajectories of the processes in volv ed is that they are flat with step or jump discontinuities, which carry information about the random inter-e vent times. Let these discontinuity points of trajectories x ω ( t ), n ω ( t ) and s ω ( t ) be denoted by { δ x 1 , δ x 2 , . . . } , { δ n 1 , δ n 2 , . . . } and { δ s 1 , δ s 2 , . . . } respecti vely . These three collections are related, as can be seen from trajectory s ω ( t ) implied in (A.2) by definition (A.1): s ω ( t ) ≡ x ω ( n ω ( t )) =                  x ω (0) for 0 ≤ t < δ n 1 x ω (1) for δ n 1 ≤ t < δ n 2 x ω (2) for δ n 2 ≤ t < δ n 3 . . . . . . (A.2) From (A.2),  δ s 1 , δ s 2 , . . .  will be the subsequence  δ n g 1 , δ n g 2 , . . .  of  δ n 1 , δ n 2 , . . .  that skips δ n i if and only if x ω ( i ) = x ω ( i − 1), e.g., g 1 = 1 (or equi valently δ s 1 = δ n 1 ) if and only if x ω (1) , x ω (0), otherwise δ n 1 will be skipped and δ s 1 will be some later discontinuity point of n ω ( t ). Points  δ s 1 , δ s 2 , . . .  are in fact a set of realizations of the random ev ent times in { S ( t ) } , providing a means to characterizing the distrib ution of the i th inter-e vent time of { S ( t ) } . This random inter-e vent time is denoted by T S i and its realizations by t S i . (In addition, we artificially set the first e vent time to 0, so that T S 1 coincides with the first non-zero ev ent time.) The conditional distribution of inter -event times T S i is easiest described via an example of trajectories, which later facilitates describing the conditional distribution of ev ent sizes. Consider for e xample the trajectories in Figure A.1 below: the realized first inter-e vent time t S 1 = t N 1 + t N 2 + t N 3 . In general, T S 1 = T N 1 + . . . + T N G 1 , where G 1 follows a geometric distribution with succes probability π ( X (0)) = 1 − e λ X ( X (0)) . Since X (0) = S (0) by definition of N (0), conditionally on S (0) = s , T S 1 ∼ exponential( π ( s )). Returning to Figure A.1, t S 2 = t N 4 + t N 5 . In general, T S 2 = T N G 1 + . . . + T N G 1 + G 2 , where G 2 follows again a geometric distribution now with success probability π ( X ( G 1 )). Since X ( G 1 ) = S ( T S 1 ), conditionally on S ( T S 1 ) = s , T S 2 ∼ exponential( π ( s )). It follo ws that in general, conditionally on S  P i − 1 i = 1 T S i  = s , T S i ∼ e xponential( π ( s )), which shows (5). Regarding jump or e vent sizes, define the first jump size J 1 ≡ S ( T S 1 ) − S (0). In Figure A.1, J 1 = X (3) − X (0) = X (3) − X (2). In general, J 1 = X ( G 1 ) − X ( G 1 − 1). Since X ( G 1 − 1) = X (0) = S (0), conditionally on S (0) = s , the probability of a first jump of size k > 0 is P  J 1 = k    S (0) = s  = P  X ( G 1 ) = s + k    X ( G 1 − 1) = s  π ( s ) = P  X (1) = s + k    X (0) = s  π ( s ) 8 where the last equality follo ws by time homogeneity of { X ( t ) } . This result generalizes to the i th jump J i by an argument analogous to the preceding one (used to generalize the result for T S 1 to T S i ), which shows (4). Note 6 ( On extending the pr oof to other time changes ) . The distrib utions of T S i and J S i ar e not as straightforwar d to derive using the strate gy used in this pr oof if { N ( t ) } is not a P oisson pr ocess. F or e xample, if it is a gamma pr ocess, T S i is still e xponential but the distribution of J S i depends on that of the over-shoot of the hitting time of { N ( t ) } to the corr esponding δ X i . However , this over-shoot distribution seems to have r emained elusive and we ar e only awar e of appr oximations to its moments (see for example Nicolai et al., 2009). If instead of a gamma pr ocess, { N ( t ) } is some counting pr ocess (of finite activity) other than P oisson, the strate gy followed in this pr oof r emains e ff ective, although the pr oof becomes mor e involved. 0 x ω ( t ) t 1 2 ×××× δ x 1 ×××× δ x 2 3 4 ×××× δ x 3 5 0 n ω ( t ) t ×××× δ n 1 s ω ( δ n 1 )= x (0) ��� ×××× δ n 2 s ω ( δ n 2 )= x (0) ��� ×××× δ n 3 s ω ( δ s 1 )= x (3) s ω ( δ s 1 )= x (3) s ω ( δ s 1 )= x (3) ��������� ×××× δ n 4 s ω ( δ n 4 )= x (3) ��� ×××× δ n 5 s ω ( δ s 2 )= x (5) s ω ( δ s 2 )= x (5) s ω ( δ s 2 )= x (5) ��������� 0 s ω ( t ) t ×××× δ s 1 δ s 1 δ s 1 ×××× δ s 2 δ s 2 δ s 2 Figure A.1: Example of composition of trajectories. The three horizontal lines represent time and show (marked with × ) realizations of several ev ent times (i.e., trajectory discontinuities) for each process. Since the first e vent in { X ( t ) } does not occur until time interval (2 , 3], ev ent time realizations δ n 1 and δ n 2 hav e no e ff ect on { S ( t ) } and δ n 3 triggers the first event time realization of { S ( t ) } (i.e., δ s 1 = δ n 3 ), turning the two events of size one in { X ( t ) } occurred at δ x 1 , δ x 2 ∈ (2 , 3] into a simultaneous event of size two in { S ( t ) } . After δ S 1 , there is δ S 2 = δ N 5 , at which the single event (again of size one) occurred at δ X 3 ∈ (4 , 5] is shifted to δ S 2 and remains of size one. 9