In this article, we introduce the package Moment Closure, which may be used to generate closure differential equations and closure approximations of the cumulants (moments) of a nonlinear stochastic compartmental model with Markov transitions. Specifically, this package defines the pair of functions MomentClosureSystem and MomentClosurePlots that achieves moment closure through the neglect of high-order cumulants. We demonstrate the application of these functions through the analysis of several test models. In select cases, the resulting cumulant approximations are compared across neglect levels and to exact answers.

Introduction

The complexity of directly estimating the probability distribution of nonlinear stochastic systems with Markov transitions rapidly increases toward intractability with the magnitude of the assumed state space. This increase in complexity is magnified in a network setting when there are multiple nodes (compartments) and numerous intensities. The formulation of a partial differential equation that describes the moment-generating function of such systems is immediate, yet the solution to this is often intractable for even the simplest systems, thereby leaving the direct estimation of the moments unattained. Moment-closure methods, however, may be used to specify a functional relationship between the moments of the system, which thereby allows for the approximation of only a few moments through a closed set of approximating differential equations. These functional relationships are achieved either through the imposition of a known parametric probability distribution on the state space of the system or through the neglect of the high-order cumulants.

The field of moment closure was pioneered by Peter Whittle in the late 1950s [1]. His work consisted of the imposition of a normal distribution on the state space of the system, thereby providing a two-moment (cumulant) closure scheme. While the normal assumption provides computational simplicity in the estimation of the moments, these approximations are often not accurate due to the frequent large deviation of the actual, yet unknown, state distribution from the normal. Since this pioneering work, several researchers have investigated using other parametric distributions and/or high neglect levels with moment-closure methods [2, 3, 4]. Findings suggest that the accuracy of carefully chosen parametric distributions is often superior to that of neglect, yet the difficulty in selecting such a distribution for a given system does not make this approach robust in the general sense. As an alternative, several recent investigations have focused on simply raising the level of cumulant neglect to achieve accurate approximations [5]. Findings suggest that this approach may lead to increased accuracy, yet at the expense of computational effort and possible stability problems. Other related investigations have focused on the accuracy of using moment-closure approximates in moment-based density approximations [6], and on the stability of moment-closure methods [7, 8].

The application of moment closure with neglect has been applied to many diverse areas in science and engineering, including queueing network modeling [9, 10], air traffic modeling [11], ecological models [12], and epidemic modeling [13]. As readily available computational resources increase, applications for moment closure with high-order neglect will become more prevalent across disciplines. The package Moment Closure was created to provide an efficient mechanism to generate the closure differential equations and closure approximations of the cumulants of a stochastic Markov system with up to 25 nodes by neglecting those cumulants that exceed a user-defined level. In the remaining sections of this article, we will provide a general overview of moment closure with cumulant neglect, and demonstrate the use of this package on several test systems. Limitations of the package and areas of future research in moment-closure methods will be discussed.

Overview of Moment Closure

Consider a stochastic Markov -node system with the state space for . The set of instantaneous changes in the state of this network in is given by with the corresponding set of intensities , where is a polynomial intensity function with non-negative powers corresponding to the set in . As shown in Bailey [14], a partial differential equation of the moment-generating function of is given by

(1)

where is a partial differential operator defined for all that replaces each term with the partial derivative of with respect to . Note that equation (1) may be rewritten for the cumulant-generating function as

(2)

An expansion of the partial derivatives and exponentials in equation (2) yields an expression that is polynomial in . An equation for these coefficients yields a set of ordinary differential equations, which is subsequently closed through the definition of a neglect (truncation) level above which all cumulants will be set to zero.

As an example, consider a simple two-compartment birth-death model. The birth intensity is given by , the transition by , and the death by . For and , the partial differential equation of the cumulant-generating function is given by

(3)

An expansion of the partials and exponentials under a neglect level of two yields the following set of ordinary differential equations.

(4)

The solution to this set of equations approximates the first- and second-order cumulants (moments) of the system. Note that a cumulant neglect level of two is equivalent to the imposition of a parametric normal distribution, as was originated by Whittle.

The Moment Closure Package

The Moment Closure package defines the pair of functions MomentClosureSystem and MomentClosurePlots that respectively return closure differential equations and closure approximations of the cumulants of a nonlinear stochastic compartmental model with Markov transitions. The package is compatible with Mathematica 4.0 or higher, and should be placed in a directory that is accessible for loading. It is available from content.wolfram.com/sites/19/2010/10/MomentClosure.m. It may be called from the front end as follows.

The function MomentClosureSystem receives a numerically or symbolically defined system and a level of neglect as input, and returns a set of closure differential equations inclusive of cumulants through the neglect level. The function is entered as where:

• Int is a list of the instantaneous transitions and intensities of the system.
• Neg set to m defines that level above which all cumulants will be neglected.

Note that as a matter of formatting for the input, state variables must be subscripted with the compartment number, Greek letters must be used as parameters in symbolically defined systems, and multiplication is not implied. Further help for the usage of this function may be called from the front end by entering the following command.

The function MomentClosurePlots receives a numerically defined system, initial conditions for the state of the system, the neglect level, the time horizon, a list of those cumulants for which displayed output is desired, and a list of the desired range for the axis for the plotted output of the cumulants. The function is entered as , where:

• Int is a list of the instantaneous transitions and intensities of the system.
• Init is a list of the cumulants that have a nonzero initial value at time zero.
• Neg set to m defines that level above which all cumulants will be neglected.
• TimeSpan set to t defines the upper endpoint of the interval through which evaluation occurs.
• Fns is a list of those cumulants whose solution will be plotted.
• Range is a list of upper and lower endpoints of the axis of one-to-one correspondence with PlotFunctions.

The default of Init is an empty list that denotes all cumulants initially set to zero. The default of Range is an empty list with the range of the axis of all plots set to All. As before, note that state variables must be subscripted with the compartment number and multiplication is not implied. Help may be obtained within the front end by calling the usage of this function as follows.

Internal to the Moment Closure package code, the functions MomentClosureSystem and MomentClosurePlots call the functions LHSPDE, EQUATE, INITIAL, and NUMERICALSOLUTION. In particular:

• LHSPDE defines a moment- and cumulant-generating function of proper dimensionality for the defined system, and generates the left-hand side (LHS) of the partial differential equation of the cumulant-generating function (CGF) given in equation (2).
• EQUATE generates the right-hand side (RHS) of the partial differential equation of the CGF given in equation (2), and equates the LSH and RHS coefficients of the partial differential equation to generate a closed set of ordinary differential equations (ODEs) for the cumulants, an example of which is given in equation (4).
• INITIAL creates a list of initial conditions for the ODEs.
• NUMERICALSOLUTION calls the NDSolve function to numerically solve the ODEs.

Numerical Applications of Moment Closure with Cumulant Neglect

We use the MomentClosureSystem and MomentClosurePlots functions of the Moment Closure package to generate the closure differential equations and cumulant approximations for various stochastic compartmental models with topologies that are frequently encountered. Select cases highlight the accuracy of cumulant approximations and demonstrate the effect of increasing the neglect level. A discussion about stability and other limitations of this package follows.

Two-Compartment Birth-Death Model

Consider the two-compartment birth-death model presented previously. Recall that the birth intensity is given by , the transition by , and the death by . The following command calls the function MomentClosureSystem to generate the moment-closure differential equations for the cumulants of the random state vector under a neglect level of . Note that the state variables of the multiple compartment systems are subscripted, that is, corresponds to the instantaneous number of entities in compartment i, and the symbolic parameters of the intensity functions are entered as Greek letters.

To further this example, let , , and . Assume that the first cumulant of the first node is initially set to 10 and all other cumulants are zero. The following command calls the function MomentClosurePlots to display the closure-based expectation and variance of through time under a neglect level of over the respective ranges and . Note that and are the first and second cumulants of the first compartment, which correspond to the expectation and variance of , respectively.

We continue the example presented previously by raising the neglect level to and and displaying comparative plots for all first- and second-order cumulants, which correspond to the expectations, variances, and covariance of and . Note that the input to the function MomentClosurePlots that was used to generate these plots is omitted due to space considerations, yet consists of repetitively evaluating the function under neglect levels of , , and , and displaying the plots together through the Show command.

Simple Birth-Death Model

Consider a simple birth-death model defined by and . The following command calls the function MomentClosureSystem to generate the moment-closure differential equations for the cumulants of under a neglect level of .

The following command calls the function MomentClosurePlots to display the closure-based first and second cumulants of through time under a neglect level of .

These cumulant approximations are compared to exact measures in the following displayed plots. (The exact measures were obtained through the direct solution of the Kolmogorov equations [15].) Note that the input used to generate these plots is omitted due to space considerations.

Three-Compartment Model with Finite Calling Population

Consider a three-compartment model with a finite calling population defined by and , which is graphically depicted in the following sketch.

Assuming that there are initially 100 units available in the calling population, the following command yields first cumulant approximations for each compartment and for the sum of all three compartments through time under a neglect level of .

Limitations

The package Moment Closure is comprehensive in the sense that it does not place any restrictions on the topology of the system that is entered, other than it must be less than 25 nodes in size. The ensuing set of ordinary differential equations, however, may not have a stable solution for the set of initial conditions under which it is being solved. In particular, stationary saddle points may define regions of instability in the feasible space [7], which is not uncommon for complex multidimensional systems. The Moment Closure package is limited, in that it does not check for the stability of a solution prior to attempting to solve the differential equations numerically. As such, instability is often detected only through error messages indicating a singularity while attempting to solve the differential equations over the range specified. Specifying alternative initial conditions or reparametrizing the system may lead to a stable solution in such cases.

Conclusion

In this article, we introduced the package Moment Closure and have demonstrated its power and simplicity in analyzing Markov systems. The study and application of moment-closure methods with cumulant neglect, however, has been limited historically. It is our conjecture that this was due largely to the lack of computational hardware and software necessary to perform such computations beyond the normal assumption. In our search of the literature, we were unable to find any other documented software contained in a single environment that is able to take general nonlinear Markov systems as input and return moment-closure approximations based on cumulant neglect as output. In this regard, we hope that the package Moment Closure will support the research efforts of others and open avenues for insightful research related to moment closure. Areas of related research that we are presently pursuing include the characterization of moment-closure stability, the development of semi-parametric approaches to moment closure, and applying moment closure to biological systems.

Acknowledgments

The creation of Moment Closure was partially supported by the Air Force Office of Scientific Research under grant #F49620-03-1-0310. The authors would like to acknowledge the work of Karl Adams in the primary coding of Moment Closure and the inspiring work of Qi Zheng, whose package MomCumConvert [16] was central to early versions of this work. We would also like to graciously thank the anonymous reviewers, whose insightful comments greatly improved this article.

References

[1]	P. Whittle, “On the Use of the Normal Approximation in the Treatment of Stochastic Processes,” Journal of the Royal Statistical Society, Series B (Methodological) 19, 1957 pp. 268-281. www.jstor.org/pss/2983819.
[2]	T. Wehrly, J. Matis, and G. W. Otis, “Approximating Multivariate Distributions in Stochastic Models of Insect Population Dynamics,” in Multivariate Environmental Statistics (G. Patil and C. Rao, eds.), Amsterdam: North-Holland, 1993.
[3]	J. Matis, Q. Zheng, and T. Kiffe, “Describing the Spread of Biological Populations Using Stochastic Compartmental Models with Births,” Mathematical Biosciences, 126(2), 1995 pp. 215-247. library.wolfram.com/infocenter/Articles/3224.
[4]	I. Nåsell, “Extinction and Quasi-Stationarity in the Verhulst Logistic Model,” Journal of Theoretical Biology 211, 2001 pp. 11-27. www.math.kth.se/~ingemar/forsk/verhulst/ver10.pdf.
[5]	S. Wojtkiewicz, B. Spencer, and L. Bergman, “On the Cumulant-Neglect Closure Method in Stochastic Dynamics,” International Journal of Non-Linear Mechanics 31(5), 1996 pp. 657-684. doi:10.1016/0020-7462(96)00029-7.
[6]	E. Renshaw, “Applying the Saddlepoint Approximation to Bivariate Stochastic Processes,” Mathematical Biosciences 168(1), 2000 pp. 57-75. doi:10.1016/S0025-5564(00)00037-7.
[7]	I. Nåsell, “An Extension of the Moment Closure Method,” Theoretical Population Biology 64(2), 2003 pp. 233-239. doi:10.1016/S0040-5809(03)00074-1.
[8]	H. Burchard and E. Deleersnijder, “Stability of Algebraic Non-Equilibrium Second-Order Closure Models,” Ocean Modelling 3(1-2), 2001 pp. 33-50. doi:10.1016/S1463-5003(00)00016-0.
[9]	T. Matis and R. Feldman, “Transient Analysis of State-Dependent Queueing Networks via Cumulant Functions,” Journal of Applied Probability 38, 2001 pp. 841-859. www.jstor.org/pss/3215768.
[10]	T. Matis and R. Feldman, “Correction: Transient Analysis of State-Dependent Queueing Networks via Cumulant Functions,” Journal of Applied Probability 42, 2005 p. 302. www.jstor.org/pss/30040790.
[11]	D. Long, et. al., “Modeling Air Traffic Management Technologies with a Queuing Network Model of the National Airspace System,” National Aeronautics and Space Administration, Contractor Report #1999-208988, 1999. portal.acm.org/citation.cfm?id=888107.
[12]	M. Keeling, “Multiplicative Moments and Measures of Persistence in Ecology,” Journal of Theoretical Biology 205, 2000 pp. 269-281. doi:10.1006/jtbi.2000.2066.
[13]	A. Lloyd, “Estimating Variability in Models for Recurrent Epidemics: Assessing the Use of Moment Closure Techniques,” Theoretical Population Biology 65, 2004 pp. 49-65. www4.ncsu.edu/~allloyd/pdf_files/variability.pdf.
[14]	N. Bailey, The Elements of Stochastic Processes, New York: John Wiley and Sons, 1964.
[15]	S. Ross, Stochastic Processes, 2nd ed., New York: John Wiley and Sons, 1996.
[16]	Q. Zheng. “MomCumConvert: A Package for Conversion between Moments and Cumulants.” Wolfram Library Archive. (Feb 19, 2010) library.wolfram.com/infocenter/MathSource/807.
T. Matis and I. G. Guardiola, “Achieving Moment Closure through Cumulant Neglect,” The Mathematica Journal, 2010. dx.doi.org/doi:10.3888/tmj.12-2.

About the Authors

Timothy I. Matis is an associate professor in the Department of Industrial Engineering and a site director for the Center for Engineering Logistics and Distribution at Texas Tech University. In addition to moment-closure methods, his research interests include queueing theory, random processes, engineering statistics, and the pedagogy of engineering education. He has been the principal investigator on several projects related to the development and application of moment-closure methods.

Ivan G. Guardiola is an assistant professor in the Department of Engineering Management and Systems Engineering at the Missouri University of Science and Technology (formerly known as University of Missouri-Rolla). His research interests lie within stochastic modeling, optimization, and complex engineering system design. He has contributed to various research projects relative to the use of moment-closure methods and stability analysis of nonlinear state-dependent systems.

Timothy I. Matis
Department of Industrial Engineering
Texas Tech University
Box 43061
Lubbock, Texas 79409-3061
timothy.matis@ttu.edu

Ivan G. Guardiola
Department of Engineering Management and Systems Engineering
Missouri University of Science and Technology
204 Engr. Mgmt. Bldg.
600 W. 14th St.
Rolla, MO 65409-0370
guardiolai@mst.edu