Open Markov Chains: Cumulant Dynamics, Fluctuations and Correlations

In this work we propose a model for open Markov chains that can be interpreted as a system of non-interacting particles evolving according to the rules of a Markov chain. The number of particles in the system is not constant, because we allow the particles to arrive or leave the state space according to prescribed protocols. We describe this system by looking at the population of particles on every state by establishing the rules of time-evolution of the distribution of particles. We show that it is possible to describe the distribution of particles over the state space through the corresponding moment generating function. This description is given through the dynamics ruling the behavior of such a moment generating function and we prove that the system is able to attain the stationarity under some conditions. We also show that it is possible to describe the dynamics of the two first cumulants of the distribution of particles, which in some way is a simpler technique to obtain useful information of the open Markov chain for practical purposes. Finally we also study the behavior of the time-dependent correlation functions of the number of particles present in the system. We give some simple examples of open chains that either, can be fully described through the moment generating function or partially described through the exact solution of the cumulant dynamics.


Introduction
Markov chains are discrete-time models for stochastic evolution, widely used to model systems in physics [1,2], chemistry [2], biology [3,4] among other disciplines [5,6].Roughly speaking, a Markov chain consists of a sequence of random variables {X t ∈ S : t ∈ N 0 } taking values from a (finite or countable) set S, called state space.The jump of the random variable from one state to another in one time step occurs with a prescribed probability, and the probabilities of all the possible jumps are collected in a matrix called Markov matrix or stochastic matrix.A natural way to interpret a Markov chain comes from physics; we can think of the random variable X t as the position at time t of given particle, and this particle moves on the discrete space S. Thus, the stationary probability vector π (provided it exists) is interpreted from a point of view of ensembles: if we have a collection of N non-interacting particles moving according to the rules of the Markov chain, then, the stationary distribution of particles on S is N π if N is large enough.This point of of view clearly shows that a Markov chain is a closed system, since there is no inflow of particles to S nor outflow of particles from S.
In this paper we shall be concerned with the case in which we allow the particles to arrive or leave the state space according to a prescribed protocol.To be precise let us consider a state j ∈ S. On one hand, at every time step we allow a certain number of particles already present in the state j to leave this state to the "outside".On the other hand we also allow a certain number of particles (from the "outside") to arrive at the state j.Both, the number of incoming particles and the number of outcoming particles, are modeled as random variables (or sequences of random variables) with a distribution given a priori.Our main goal in this paper is to describe the population of particles on the state space as well as its fluctuations.Clearly this model for open Markov chains seems to be a kind of non-equilibrium system and we are particularly interested in the behavior of both, the space correlations and time correlations for several possible scenarios for the incoming and outcoming protocols.
At this point it is convenient to mention some works related to our model.For example, in a recent work [7], Floriani et al have considered the case in which some quantity (that they call "mass") moves through the states of the chain according to the Markov transition probabilities.Particularly Floriani et al focus their study to the case where the Markov chain has some absorbing states (in which the mass accumulates) and the mass is supplied by either, a constant source or a periodic source.In contrast, in our model the number of incoming particles at every time step is not constant but a sequence of random variables not necessarily independent and identically distributed (i.i.d.).Moreover, instead of modeling the outflow by absorbing states, we define a protocol of outcoming particles, which allows every particle to leave the chain with a state-dependent probability.Another work which is worth mentioning here is the one of Pollard and co-workers [8][9][10].They consider a class of open Markov process in which the state space is discrete and the time continuous.They assume that the incoming and outcoming fluxes are regulated by means of a set of special states of S (which they call "boundary").One of the main differences of our approach with respect to the one proposed by Pollard is that they suppose that the elements of the boundary has a distribution prescribed a priori (which may be even time-dependent).In contrast, our approach considers every state as a source or sink of particles.In this sense our model can be thought of as a Markov chain in contact with a reservoir of particles.Thus, according to a prescribed protocol the particles goes from the reservoir to the chain and viceversa, the particles go from the chain to the reservoir with a prescribed protocol.This way of modeling the source of particle is, in some way, similar to the grand canonical ensemble in thermodynamics, in which a system is in contact with a reservoir allowing particles to be interchanged.

A model for open Markov chains
The main idea behind our model for an "open" Markov chain is that we allow the particles to arrive and leave the state space S. We have mentioned that the particle can enter the state space according to a protocol which is modeled as a sequence of random variables.Such a sequence of random variables determines the number of particles arriving at certain state every time step.On the other hand, the particle can leave the state space depending on which state they are.The most natural way to model this situation is by defining, for every state, a given probability with which the particle leaves such a state towards the reservoir.This probability must satisfy a compatibility condition, consisting in the fact that a particle in a given state has only two options i) jump to any other state in S or ii) jump to the reservoir.The sum of all these probabilities should be one, in order for the "jump" to be well defined.Notice that due to the compatibility condition, we have that the protocol of outcoming particles is completely determined by means of a non-negative matrix Q with a spectral radius strictly less than one.This is because the "missing probability" in Q (the necessary for Q to be a stochastic matrix) is interpreted as "jump probabilities" towards the outside.Definition 2.1 (Open Markov chain) Let S be a finite set, whose cardinality is denoted by #S = S, and let Q : S × S → [0, 1] ⊂ R be an irreducible and aperiodic matrix with spectral radius strictly less than one.Let {J t : t ∈ N} be a sequence of random vectors taking values in N S 0 .We say that S, Q, {J t : t ∈ N} is an open Markov chain with state space S, jump matrix Q and incoming protocol {J t : t ∈ N}.Now let {N t : t ∈ R} be a sequence of random vectors taking values in N S 0 .Such a sequence is defined as follows.Given the initial random vector N 0 with a given distribution, we define N t recursively as where R t is a random vector taking values in N S 0 , whose components are given by, In the above expression the set ) for all t ∈ N 0 , where q i,j is the (i, j)th component of Q.Then we say that N t is the distribution over the state space at time t with initial condition N 0 .
In the above definition we adopted the notation N 0 to indicate the set {0, 1, 2, 3, . ..} to distinguish it from the set N := {1, 2, 3, . ..}.Notice that Eq. ( 1) establishes the evolution of the number of particles.This equation states that the number of particles at time t + 1 is the number of particles having arrived as the state space (represented by J t ) plus the number of particles having remained in the state space (which is represented by R t ).Observe that the random vector R t has independent components, and every component corresponds to a sum of random variables with binomial distribution.The latter is a consequence of the fact that the number of particles jumping from the state i to the state j is a random variable with binomial distribution with parameters q i,j (the probability to go to the state j from the state i) and N t i (the number of particle in the state i at time t).It is clear that N t+1 depends on N t through the random variable R t which is the responsible of the redistribution of particles among the internal states.
From the above definition we can appreciate that the quantity of foremost interest is N t , the distribution of particles on the state space at time t.Our goal is then to provide a way for determining the probability function of N t and to determine whether or not the process {N t : t ∈ N 0 } is able to reach an stationary distribution.We will see that {N t : t ∈ N 0 } is actually a Markov process and the main goal is to determine its properties.Before establishing a result on this direction let us give some examples of open Markov chains.
Example 2.1 Let us consider the most simple case in which the system has only one state.In this case the matrix Q consists of a single number q (a 1 × 1 matrix), which should be non-negative and strictly less than one, i.e., 0 ≤ q < 1.The most simple case for the incoming protocol consists of a sequence of constant random variables all these taking the same value, which we call J 0 ∈ N 0 .This means that the number of particles arriving at every time-step is J 0 .Every particle arriving at the unique available state has only two options jump to the outside or remain in its state.The probability of remaining in the state is q and the probability of jumping to the outside is 1 − q.This simple example for open Markov chain can be illustrated as a graph with only one vertex and three edges as shown in figure 1.Notice that there are two special edges, one establishing the incoming protocol (labeled by J 0 ) and one defining the outcoming protocol (labeled by 1 − q).
One-vertex open chain.The circle represent the unique state available for the particles.The arrows in the graph stand for the "jump" rules allowed for the particles.Notice that the arrows that do not connect two states represent the incoming and outcoming protocols.
Example 2.2 A less trivial example is provided by giving a matrix larger than 1 × 1.

Let us consider for example the matrix given by
Notice that the above matrix is not stochastic, because some rows does not add to one, but less than one.The latter means that not all the states allow the particles to leave to the outside.As we can see the first row adds to 3/4, which means that every particle on the state 1 has a probability 1/4 to go out to the reservoir.The second row adds to 1/2, this means that a given particle on the state 2 has a probability 1/2 to go out to the outside.Finally, the third row adds to one, meaning that a particle on the state 3 can only jump to the other states 1 or 2 or remains in its current state, but it cannot leave the state space to go to the outside.Now assume that the incoming protocol {J t : t ∈ N} is a set of i.i.d.random vectors, i.e. the protocol is time-independent.Then the number of incoming particles at every time-step can be considered as independent realizations of a single random vector, which we denote by J. To be more precise, we can chose in particular this random vector as J = (J 1 , 0, J 3 ), with J 1 and J 2 two independent random variables with Bernoulli distribution with parameters p 1 and p 3 respectively.If, for instance, the parameters of every Bernoulli distribution were given by p 1 = 0.1 and p 3 = 0.6, then this would mean that the average number of particles arriving at the vertex 1 is lower than the average number of particles arriving at the vertex 3. Figure 2 shows the graphical representation of this open Markov chain.3. The evolution of the particle distribution

Evolution of the moment generating function
In this section we will establish the evolution equation for the process N t , which, as we have anticipated, is a Markov process.This fact can actually be appreciated in equation ( 1), where we have indicated that the number of particles at time t + 1 is uniquely determined by the number of incoming particles and the redistribution of particles that were present at time t.Clearly the last condition states the Markov property for the process {N t : t ∈ N 0 }.In order to obtain the stochastic matrix governing the behavior of N t , let us define p t (n) as the probability vector associated to N t , i.e., p t (n) := P(N t = n).
(4) We will refer to p t (n) as the distribution of particles over the state space or, to simplify, distribution over the state space Let us consider the probability vector at time t + 1.
Notice that, due the fact that N t+1 depends only on N t , which is a consequence of the relation (1), as we have mentioned above, it is clear that Observe that the last expression is a consequence of the Markov property of the process which was assumed in equation (1).Notice that equation ( 5) establishes the evolution equation for p t , which can be written as Now let us call K(k, n) the conditional probability appearing in the above equation, i.e., and observe that this quantity can be rewritten as follows The dependence on k in the above expression is implicit in the random vector R t , since the redistribution of particles depends on the number of particles on every state at time t, which is indeed given by k.Thus, the function K : N S 0 × N S 0 → [0, 1] ⊂ R can be thought as the stochastic matrix corresponding to the process {N t : t ∈ N 0 }, since the evolution equation for the probability vector p t (n) is given in terms of K as follows, In order to solve equation ( 9) for p t (n) it is necessary to make some assumptions on the nature of the random vectors J t and R t for all t.First of all, it is natural to assume that J t and R t are independent.This assumption actually means that the number of particles incoming to the state space does not have any influence on the redistribution of particles already present in the chain.This implies that the joint probability for the random vectors J t and R t can be factorized as the product of its corresponding probability vectors, i.e., P(J t = j; R t = r) = P(J t = j)P(R t = r).
The above equality allows us to express the kernel K(k, n) as, Next, we will solve for p t by using the well-known technique of the moment generating function (m.g.f.).To this end let us introduce some notation.Let G t : R S × R S → R be the m.g.f. of N t , which is defined as At this point it is important to describe our convention for vectors in R S .First of all we should emphasize that we interpret the vectors n, α, N t , etc., as row vectors (i.e.matrices of size 1 × s).Thus, the superscript T means, as usual, matrix transposition implying that the vector α T is a column vector (a matrix of size s × 1).Thus, within this convention, the product nα T should be understood in the sense of the usual matrix product, which in this case results in a single number.Analogously, we also define the moment generating functions for J t and R t as follows, Notice that we omitted the superscript t in the m.g.f. for R t , since two random vectors, say for example R t and R s , are independent and have the same distribution, and consequently, share the same moment generating function.
Next, our objetive is to provide a recurrence relation for G t using the evolution equation ( 9).Thus, let us consider the m.g.f. for N t+1 , and note that where, in the last equality, we used the form of K given in equation (10).We can appreciate that the summation over n together with the summation over the restriction j + r = n results in a double sum over the "indices" j and r without restrictions, i.e., we obtain two sums over independent indices.This observation allows us to write In the above equation we can observe that the summation over j and r results in the m.g.f. for J t and R t respectively.This implies that In Appendix A we show that H(α) can be written as where the function Observe that the equation ( 17) shows explicitly the dependence on k of the m.g.f. for R t .Then, if we substitute the relation ( 17) into (16), we obtain, We can easily see that the summation over k results in the m.g.f. for N t .Thus, Equation ( 20) is a recurrence relation governing the time-dependence of the m.g.f. for N t .This equation can be formally solved to obtain where G 0 stand for the m.g.f for N 0 (the initial distribution on the state space).We should emphasize that the superscript notation H (r) stands for the rth iteration of the function H, i.e., H (r) := H • H • . . .H, r times.Now, let us assume that the process {J t : t ∈ N} is a sequence of i.i.d.random vectors.In this case, the m.g.f. for J t is the same for all t, consequently the formal solution for G t can be expressed as, This result, together with the fact that H (r) (α) → 0 as r → ∞ in an open neighborhood around α = 0 (see Appendix A for a proof), implies that for the case in which the random vectors J t are identically distributed for all t (not necessarily independent), we have that the process {N t : t ∈ N 0 } admits an stationary solution.This is because the m.g.f.G t attains a limit when t → ∞.Such a limit can be written as, whenever the infinite product exist.If this is the case, the m.g.f.G stat (α) is additionally a solution for the evolution equation (20).This means that G stat (α) corresponds to a m.g.f. of a distribution p stat (n) over the state space, which is invariant under the dynamics (9).

Example 3.1 Let us consider the open chain consisting of only one vertex given in example 3.1.
In this case we will assume that the incoming protocol {J t : t ∈ N} consists of a sequence of i.i.d.random variables having Bernoulli distribution with parameter p.Since all J t have identical distribution, then we have only one m.g.f.F characterizing them.This function is given by, On the other hand, due to the fact that Q is a 1 × 1 matrix, we have that H is a real-valued function depending on one variable, α, given by, where we denoted by q the unique element of Q.In this case, the rth iteration of H can be exactly calculated.A straightforward calculation shows that, Notice that H (r) (α) → 0 as r → ∞, as we have anticipated above.Hence, in this case, the stationary m.g.f. is given by The above result shows that the stationary distribution can be interpreted as the convolution of an infinite sequence of Bernoulli distributions with parameters pq r .This particularly means that the random variable N t , when it has reached the stationarity can be written as a sum of Bernoulli random variables X r (with parameter pq r ), The above result allows us to compute, for example, the expected value E[N t ] and the variance, We should observe that the expected value E[N t ] = p/(1 − q) for the number of particles inside the unique vertex shows that an equilibration is attained between the number of incoming particles and the number of outcoming particles.It is clear that the mean number of arriving particles per unit time is p, while, the number of leaving particles per unit time is (1 − q)E[N t ].Once N t has reached the stationarity, we have an equality between these quantities, giving the result stated above.Although we have obtained the mean number of particles by using the argument of equilibration, the same line of reasoning cannot be applied to the variance.We have thus provided a way to compute the fluctuations in the number of particles that are present in the vertex.

Cumulant dynamics
Up to now we have obtained two main results, a recurrence relation of the timedependent m.g.f. of N t and a formal expression for the m.g.f. of N t when the system has attained the stationarity.Now we will focus in two quantities which are of special interest, namely, the two first cumulants for N t , and how do these quantities evolve in time.Let us start by obtaining the first cumulant of N t .Notice that the first cumulant (which coincides with the first moment) can be obtained by taking the first derivative of the m.g.f.G t (α).Let us denote by µ t the first moment of N t , i.e., to which we will refer to as the mean distribution of particles over the state space at time t, or simply, the mean distribution at time t.Now let us notice that, where (µ t ) i is the ith component of µ t and α i is the ith component of α.Observe that the ith component of the mean distribution at time t + 1 can be obtained from equation (20), giving Notice that H(α = 0) = 0 and that any moment generating function evaluated at 0 is one.Hence we have, where ǫ t stands for the expected value of J t , which can be obtained by means of the first derivative of F t , i.e., ∂F t (α) In equation (30) made use of the fact that which is proved in Appendix A. Thus, it is clear that equation (30) can be written as, The above expression states how do the mean distribution µ t evolves in time.This evolution has two components, one involving the internal dynamics (which is given by the term µ t Q giving the internal redistribution of particles) and other one involving the external dynamics (which is given by the time-dependent mean number of incoming particles).Now let us explore the behavior of the second cumulant of N t .The second cumulant corresponds to the variance matrix Var(N t ) which we will denote by Σ t .Notice that this matrix has entries given by Next we will use the dynamics of the m.g.f. of N t to obtain a recurrence for the expected value The above expression together with the evolution equation (20) leads to At this point it necessary to introduce some notations.First let us denote by ∆ the variance matrix of the incoming flux, i.e., The above quantity can be obtained through the second derivative of the m.g.f.F t (α) as follows, Thus, performing the evaluation of equation ( 34) at α = 0 we obtain, where we used the fact that (see Appendix A), Let us simplify equation (37) by noticing that the summations can be written as matrix products, where we defined the matrix Λ t as Equation (39) allows us to obtain the variance matrix Var(N t+1 ), Rearranging terms in the above expression and denoting by Σ t the matrix variance Var(N t ), we obtain Finally, observing that ǫ t + µ t Q = µ t+1 , it is clear that the variance matrix satisfy the evolution equation, Equations ( 31) and ( 43) govern the dynamics of the first and second cumulants and are valid even when the incoming flux has a time-dependent distribution.In the case where the protocol of incoming particles {J t : t ∈ N} is a stationary process (for which the two first cumulants are time-independent) we have that the system can reach the stationarity.Particularly we have that the dynamics equations ( 31) and ( 43) have stationary solutions (proved in Appendix A) given by, where ǫ and ∆ are the mean vector and the variance matrix of the stationary process {J t : t ∈ N} and µ and Σ denote the mean distribution E[N t ] and the variance matrix Var(N t ) when the process {N t : t ∈ N 0 } has reached the stationarity.We also defined Λ as

Example 3.2 Let us consider a three states open chain with jump matrix given by
with p a parameter restricted to take values in the interval 0 < q < 1/2.Let us also assume that the protocol of incoming particles {J t = (J t 1 , J t 2 , J t 3 ) : t ∈ N} is a stationary process whose joint probability distribution at time t can be written as ). Particularly, we chose f 1,2 and f 3 as Notice that each random variable J t i can only take the values 0 or 1.Moreover, the above choice for the joint probability distribution for the random vector (J t 1 , J t 2 , J t 3 ) implies that the random variables J t 1 and J t 2 are dependent and that J t 3 is independent.It is easy to see that all these random variable have, separately, a Bernoulli distribution with parameter 1/2, i.e., at every time step on every node arrives one particle with probability 1/2 and arrives no particles with probability 1/2.In figure 3 we show a graphical representation of the open chain.A straightforward calculation shows that the first and second moments of J t are given by Now we calculate the mean stationary distribution of particle µ and the stationary variance matrix Σ.The calculation of µ is straightforward.We only need to obtain the inverse of the matrix 1 − Q, which is given by + q 2 q + q 2 q + q 2 1 − q 2 q + q 2 q + q 2 q + q 2 1 − q 2   .

With the above result we can see that the mean stationary distribution µ is given by
The last result implies that the particles on the state space are equidistributed when the system has reached the stationarity.Moreover, we can also appreciate that the number of particles on every state diverges as the parameter p tends to 1/2.This divergence is actually a consequence of the fact that the system does not allow that the particles leave the state space when p = 1/2.Thus, for such a parameter value, the system is still receiving particles but it does not allow that the particles escape to the outside, thus increasing indefinitely the number of particles inside the system.

Now let us compute the stationary variance matrix Σ. Recall that this quantity can be obtained by means of the formula,
First let us obtain the explicit form of Λ.According to equation ( 46) we have that Thus, we have Next we need to compute the nth power of the matrix Q.It is not hard to prove that This result shows that for this particular case because Q is itself a symmetric matrix.These results allows us to see that Thus we have to compute the inverse of 1 − Q 2 and the matrix product Q k ∆Q k .Actually, the above quantities can be exactly computed by using symbolic calculations performed in the software Mathematica.For example, we have that inverse matrix The infinite summation appearing in ( 53) can be also be exactly determined (because all the summands are exponential, thus resulting in geometric series).However such an expression is too long to write down here.Thus, instead of giving a explicit expression for the variance matrix Σ we will write down three correlation functions: the first one between the random variables N 1 and N 2 , a second one between the random variables N 1 and N 3 , and finally a third one between the random variables N 2 and N 3 ‡.These correlation functions are defined as, Then, performing symbolic calculations we obtain that κ 2,3 = κ 1,3 , and that, In order to test the results we obtained for this example we performed numerical simulations.We have simulated the dynamics of the open chain for several parameter values during a time of 5 × 10 5 steps.Then, from the time series obtained we estimated the correlation functions κ 1,2 and κ 2,3 and we compare them against the theoretical prediction given in equations ( 58) and ( 59).In figures 4a and 4b we show the correlation functions obtained by numerical simulations and computed from equations ( 58) and ( 59).We take the parameter value q = 0.25 and plotted κ 1,2 (solid line, filled circles) and κ 2,3 (dashed line, filled squares) as a function of p (figure 4a).The same graph is done but using the parameter value q = 0.45 (figure 4b).As we can see, our the results obtained from numerical simulations are consistent with the formulas that we obtained theoretically.

Distribution of particles leaving the state space
Up to now we have given an expression for the m.g.f. for the number of particles found at every state of the system.Since the system is open, at every time step there is a number of particles arriving to the system, which is determined by the random vector J t .The total number of particles per step arriving to the system is then When the system reaches the stationarity the mean number of particles within the system attains a constant value, meaning that there is a equilibration between the  number of incoming and outcoming particles.To be precise, if we denote by O t the total number of particles leaving the state space, then, under stationarity we should have that a fact that can be inferred by the "conservation" of the mean number of particles at stationarity.Our goal here is to go beyond the above expression, we would like to characterize how the random variable O t evolves in time and how much it is influenced by the incoming number of particles and the "jumping" rules of the Markov chain.To this end, let us define some quantities which will allow to describe the random variable O t explicitly.
Definition 3.1 Let S, Q, {J t : t ∈ N} be an open Markov chain with state space S, jump matrix Q (with components (Q) i,j = q i,j ) and incoming protocol {J t : t ∈ N}.Let e i be defined as the escape probability from the state i, i.e. the probability with which a particle in the state i leaves the system to the outside, Next, let U t = (U t 1 , U t 2 , . . ., U t S ) be a random vector whose components have binomial distribution as follows, U t i ∼ Binom(e i , N t i ) where N t is time-dependent distribution over the state space.Then we say the U t i is the number of particles leaving the state i to the outside at time t.The total number of particles O t leaving the system is then the random variable given by Finally, we define the vector e, with components e i = (e) i , which will be referred to as the escape probability vector of the chain.Additionally let E be a diagonal matrix with components (E) i,j = E i,j defined as which will be referred to as the e escape probability matrix of the chain.
Our main goal is now to characterize the random vector U t giving the number of particles leaving the chain.We should notice that the distribution of U t depends on N t which is also a random vector.This implies that, given the value of N t , we can specify the conditional distribution for U t , i.e., which is given by Once we know T (n; m), we can write an expression for the probability distribution for the random vector U t , If we denote by r t (m) the probability distribution of the random vector U t , i.e., then it is clear that equation (65) can be rewritten as The next steps consists in obtaining the moment generating function of U t .This is because, as we saw above, the expression for the m.g.f. of N t can be written explicitly.Thus, let R t : R S → R be the m.g.f. of Thus we can see that, using expression (66), R t can be written as follows Since T (n; m) is a product of binomial distributions, it is clear that the sum over m results in the product of moment generating functions of binomial random variables, Moreover, if we define the function C = (C 1 , C 2 , . . ., C S ) : R S → R S as follows, it is clear that we can write The above expression, together with equation ( 68), gives us The last relation states that the moment generating function of U t can be obtained by means of the moment generating function of N t , which is mediated by the transformation C.
Observe that R t (α) allows us to obtain the first and the second moment of number of leaving particles.Taking the first derivative (gradient) of R t (α) and evaluating it in α = 0 we obtain, The above result allows us to calculate the mean number of leaving particles Now, when the system reaches stationarity we have to replace µ t by µ, Thus, from the definition of e i we obtain and recalling that µ satisfy the equation which is the relation we have anticipated by invoking the particle number conservation principle.
The variance matrix Var(U t ) can also be obtained by means of the second derivative of R t (α).A calculation achieved in Appendix A shows that where D is a diagonal matrix with components where δ i,j is the well-known Krönecker delta.The variance matrix U t allows us in turn to obtain a closed formula for the variance of the total number of leaving particles, where 1 stands for the S × S identity matrix.We observe that the variance matrix of the number of outcoming particles is not the same as the variance matrix of the incoming particles nor the variance matrix of the particles in the system.Thus we have that the fluctuations and correlations are modulated by the internal dynamics of the system.This is important because measuring the correlations and fluctuations of the outcoming flux gives information on the internal dynamics of the system.
Example 3.3 Let us consider again the one-vertex model given in example 3.1.We have seen that the m.g.f. of N t is given by (1 − pq r + pq r e α ) . (80) Notice that the probability escape "vector" consists of a single number, given by e 0 = 1 − q (recall that the jump matrix is also a single number).Thus we have that the transformation C : R → R can be written as C(α) = log (1 − e 0 + e 0 e α ) = log (q + (1 − q)e α ) .
Now we can obtain the m.g.f. of U , the number of particles leaving the system, α) .
Notice that e C(α) can be written as e C(α) = q + (1 − q)e α , thus we have, The last expression establishes that the distribution of U t (at stationarity) can be seen as an infinite convolution of Bernoulli distributions, with parameters pq r (1 − q) for r ∈ N 0 .Thus, the random variable U t can be written as an infinite sum of i.i.d.random variables Y r (with the above-mentioned Bernoulli distribution), whenever U t has reached stationarity.Particularly the mean number of leaving particles as well as its variance can be exactly determined, Var We should emphasize that the above expressions can be obtained through the formulas for E[O t ] and Var(O t ), given in equations ( 74) and ( 78).Since S = 1 (because we have only one state) we have that where we used the expressions for E[N t ] and Var(N t ) given in equations ( 25) and (26).

Time correlations for the open Markov chain
Up to now we have seen that it is possible to find an explicit expression for the timedependent distribution on the state space.This evolution is fully characterized by the two first cumulants, the mean distribution over the state space µ t and the variance matrix Σ t .It is important to emphasize that we made no assumptions on the time correlations of the sequence of random vectors {J t : t ∈ N}.This is because to obtain the distribution over the state space, N t , for a given time t, it was enough to know the number of incoming particles at time t.On the other hand, if we would like to compute the two-times correlation function for certain observable we necessarily have to known the number of incoming particles a two different times.This information unavoidably will be related to the two-times covariance matrix of the process {J t : t ∈ N}, i.e., the covariance matrix between the random vector J t and J t+s for s, t ∈ N. Our goal in this section is to obtain an expression for the covariance C i,j (t, t + s) between the ith coordinate of N t and the jth coordinate of N t+s , i.e., where µ i is the ith coordinate of the mean stationary distribution.
In order to compute the expected value E[N t i N t+s j ] it is necessary to have an expression for the stationary joint distribution P t,t+s (n, m).This quantity is defined as, Our goal here is to establish a method to obtain the joint distribution P t (n, m).
Actually, we will first determine an expression for a more general quantity.Let P (n 0 , n 1 , n s ) denote the joint probability function of the random vectors N t , N t+1 , . . ., N t+s §, First of all, let us introduce some notation that will be useful to perform further calculations.Let f (j 0 , j 1 , . . ., j s−1 ) be the joint probability function of the random vectors J t , J t+1 , . . ., J t+s−1 , i.e., f (j 0 , j 1 , . . ., j s−1 ) := P J t = j 0 ; J t+1 = j 1 ; . . ., J t+s−1 = j s−1 . (87) Let us also denote by h(r; k) the probability function of the random vector R t , which, as we saw in section 3, depends on the value taken by the random vector N t (a value which we denote by k in the probability function h).
With the above-introduced notation it is possible to write the joint probability function, given in equation ( 86), in terms of the probability functions f and h, Notice that the random vector R t+j depends on n j for 1 ≤ j ≤ s − 1, values which are given a priori.Thus, the random vectors R t , R t+1 , . . ., R t+s−1 are all independent (because the values taken by the random vectors N t+j for 0 ≤ j ≤ s are all fixed), which allows us to write In terms of the probability functions defined above we have, The above expression states that the joint distribution P (n 0 , n 1 , . . ., n s ) can be written in terms of the stationary distribution p stat and the probability functions of the random vector R t and the joint distribution of the random vectors § Notice that, to be strict, the probability function P depends on t, t + 1, . . ., t + s, and we should denote this dependence explicitly by using subscripts, i.e., P = P t,t+1,...,t+1 .However, we will not use such a notation by the sake of simplicity in further calculations.The same convention will be adopted for other "multiple-times" joint probability functions or its corresponding moment generating functions.
The above expression can be used to obtain the moment generating function of (N t , N t+s ) and then the corresponding two-times covariance matrix C i,j (t, t + s) = (Cov(N t , N t+s )) i,j .Those calculations are performed in Appendix A, here we only write down the result, Notice that equation ( 90) is valid even if the system has not necessarily reached stationarity.If we assume that the system has attained the stationarity (which means that Σ t no longer depend on time), we obtain Due to stationarity, it is clear that the covariance matrix depends only on s, the difference between the times t and t + s.
Example 4.1 Let us consider the system introduced in example 3.2.We should notice that we were able to obtain an exact expression for the stationary variance matrix Σ.Thus, computing the covariance matrix involves only the product of two matrices.The resulting expression for the covariance matrix is too long to write down here.Thus, instead of giving explicitly the expression for the covariance matrix, we will show the theoretically computed time-dependent correlation functions defined as Ci,j (s) := (Cov(N t , N t+s )) i,j In figure 5 we show the behavior of the correlation functions C1,1 (s) and C1,2 (s) for the three-states open Markov chain studied in example 3.2.For the parameter values we display the theoretically computed correlation function using equation ( 91) and ( 92).The figure also shows the correlation functions obtained by means of numerical simulations of the system.The total time-steps performed to obtain the data from simulations was 5 × 10 5 .We appreciate that the theoretical results agree with the simulations showing the consistency of our results.

Conclusions
We have introduced a simple model for open Markov chains by interpreting the state space of a usual Markov chain as physical "sites" where non-interacting particles can be placed and moving throughout it according to "jumping rules" given by a kind of stochastic matrix.The conditions for the chain to be open are given as a protocol of incoming particles, defined by a discrete-time stochastic process, and by a protocol of outcoming particles, which is implicitly defined by the condition that the "stochastic matrix" (called here jump matrix) has a spectral radius strictly less than one.These conditions establish the rules by means of which the particles arrive and leave the state space to the outside.We have shown that this model can be treated by means of the moment generating function technique, allowing us to obtain, in a closed form, the moment generating function of the distribution of particles over the state space.We  have also shown that the system can be partially described by the dynamics of the two first cumulants of the the distribution of particles over the state space.Actually, we have given closed formulas for the two first cumulants when the system is able to reach the stationarity.We have also studied how the correlations in the incoming protocol of particles are processed by the open chain.We have obtained closed formulas allowing compute the two-times covariance matrix for the random vector defined as the number of particles on the states.Our main result is that the stationary two-times covariance matrix does not depend on the correlations of the particles arriving at the state space.This means that the stationary correlation functions essentially behaves as a closed Markov chain, i.e., that the correlations vanishes exponentially in time.The nonstationary correlations might probably content some information on the correlations of the incoming particles, but it would be necessary a more exhaustive study in this direction for a better understanding of such a process.= e i (1 − e i ). (A.25) The above results, together with the fact that C(0) = 0, imply that the variance matrix can be written as Var(U t ) i,j = E[N t i N t j ]e i e j + (µ t ) j e j (1 − e j )δ i,j − (µ t E) i (µ t E) j .= e i (Σ t ) i,j e j + (µ t ) i (µ t ) i e i e j + (µ t ) j e j (1 − e j )δ i,j − (µ t E) i (µ t E) j , where we used the fact that E[N t i N t j ]e i e j = (Σ t ) i,j e j + (µ t ) i (µ t ) i .Now recalling that that (E) i,j = e i δ i,j it is clear that Var(U t ) i,j = (EΣ t E) i,j + (µ t ) j e j (1 − e j )δ i,j , (A.26) which is the expression anticipated in equation (77).Next, our goal consists in obtaining the two-times covariance matrix Cov(N t , N t+s ).To this end we need to obtain an expression for the m.g.f. of the joint distribution P t,t+s (n, m).We have seen that the latter can be expressed in terms of the multiple-times joint distribution as follows

Figure 2 .
Figure 2.An example of a three-states open Markov chain.

Figure 3 .
Figure 3.An example of a three-states open Markov chain.

Figure 4 .
Figure 4.The correlation functions κ 1,2 and κ 2,3 .In panel (a) we plot the correlation functions for the parameter value q = 0.25 fixed and varying the parameter p.The solid line corresponds to κ 1,2 computed from equation (58) and the filled circles corresponds to κ 1,2 but numerically obtained from the simulations of the stochastic dynamics during 5 × 5 time steps.In panel (b) we do the same as in panel (a) but using the parameter value q = 0.45.

Figure 5 .
Figure 5. Time-dependent correlation function for the three-states open chain.Panel (a) shows the correlation function C1,1 (s) as a function of s for the parameter values p = 0.40 and q = 0.45.The theoretically computed correlation is represented by the solid line and the numerically obtained from simulations are represented by the filled cicles.Panel (b) shows the correlation function C1,2 (s) for the same parameter values as in panel (a).The theoretically computed correlation function corresponds to the solid line while the numerically estimated from simulations corresponds to filled squares.