# Efficient Computation of the Stochastic Behavior of Partial Sum Processes

In this paper the computational aspects of probability calculations for dynamical partial sum expressions are discussed. Such dynamical partial sum expressions have many important applications, and examples are provided in the fields of reliability, product quality assessment, and stochastic control. While these probability calculations are ostensibly of a high dimension, and consequently intractable in general, it is shown how a recursive integration methodology can be implemented to obtain exact calculations as a series of two-dimensional calculations. The computational aspects of the implementaion of this methodology, with the adoption of Fast Fourier Transforms, are discussed.

## Authors

• 6 publications
• 2 publications
• 2 publications
• 306 publications
09/15/2020

### Asymptotics of the Lebesgue constants for bivariate approximation processes

In this paper asymptotic formulas are given for the Lebesgue constants g...
04/16/2021

### An implementation of an efficient direct Fourier transform of polygonal areas and volumes

Calculations of the Fourier transform of a constant quantity over an are...
03/24/2022

### The SAGEX Review on Scattering Amplitudes, Chapter 4: Multi-loop Feynman Integrals

The analytic integration and simplification of multi-loop Feynman integr...
10/14/2019

### 3rd-order Spectral Representation Method: Part I – Multi-dimensional random fields with fast Fourier transform implementation

This paper introduces a generalised 3rd-order Spectral Representation Me...
01/10/2013

### Toward General Analysis of Recursive Probability Models

There is increasing interest within the research community in the design...
03/30/2017

### Minimum energy path calculations with Gaussian process regression

The calculation of minimum energy paths for transitions such as atomic a...
06/15/2020

### Efficient Ab-Initio Molecular Dynamic Simulations by Offloading Fast Fourier Transformations to FPGAs

A large share of today's HPC workloads is used for Ab-Initio Molecular D...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The tracking of the stochastic behavior of a partial sum process is an important problem with many applications. In general, calculations of the probabilistic properties of such a partial sum process require an ability to compute high-dimensional multivariate probabilities of partial sum variables. Computing these multivariate probabilities is in fact a high-dimensional integration problem which in general cannot be performed efficiently by any numerical method presently available.

Computing the multivariate probability of an event in high dimensions, in its most general form, is an intractable problem. An efficient solution may possibly be devised by exploiting any special structures of the problem under consideration. For example, with a general high dimensional density function, if the event is a convex set, a Markov chain Monte Carlo approach can be devised to efficiently approximate the probability of the event (Smith (1984), Belisle et al. (1993), Lovász (1999) and Kiatsupaibul et.al. (2011)). On the other hand, if the probability distribution is a multivariate standard normal or

-distribution with some special correlation structure, an efficient numerical integration may be constructed to compute a rectangular event, say (Dunnett & Sobel (1955), Soong & Hsu (1997)), or an event based on a complete ordering (Kiatsupaibul et al.(2017)).

The partial sums of independent random variables have a nice structure that can be exploited to devise an efficient numerical algorithm for the calculation of their probabilistic properties. The objective of this paper is to illustrate how such probability calculations for the stochastic behavior of a partial sum process of independent variables can be performed efficiently based upon the adoption of recursive numerical integration techniques.

The specific problem considered in this paper can be described as follows. Let , be independent random variables. In general, we consider probabilities of the form

 P((X1,…,Xn)∈A)=
 P((X1,X2,X3)∈A1,(X1+X2,X3,X4)∈A2,…,(X1+…+Xn−2,Xn−1,Xn)∈An−2) (1)

for sets , . The methodologies discussed in this paper are applicable to the evaulation of this general expression.

However, a special and important case of equation (1) is the sole consideration of the sum of the random variables

 P(X1+…+Xn∈B) (2)

for a set which has many varied applications. When the sum of the random variables does not have an identifiable distribution, the evaluation of this probability is ostensibly challenging, although it is shown in this paper how its evaluation is in fact straightforward for any value of .

More generally, probabilities concerning the stochastic behavior of the partial sum process of the random variables of the form

 P(X1∈B1,X1+X2∈B2,…,X1+…+Xn∈Bn) (3)

for sets , , are also a special case of equation (1). In this paper it is also shown how the evaluation of this expression is in fact also straightforward for any value of .

The key result of this paper is that the -dimensional integral expression

 ∫⋯∫(x1,…,xn)∈Sh1(x1)⋯hn(xn) dx1⋯dxn (4)

can be evaluated recursively as a series of -dimensional integral calculations when the set is defined by the conditions

 (x1+…+xi,xi+1,xi+2)∈Ii⊆R3

for . This is an application of the general discussion of recursive integration given in Hayter (2006) with . Recursive computational techniques similar to the ones developed in this paper have been applied to the problem of confidence band construction for a distribution function in Kiatsupaibul & Hayter (2015), and to ranked constrained computations in Kiatsupaibul et al. (2017).

Of course, the probability in equation (1

) can be put in this form for continuous random variables with the sets

equal to the sets and the functions

equal to the probability density functions

. In addition, if the random variables have discrete distributions then the results of this paper are still valid with the integrals replaced by sums and the probability density functions replaced by the probability mass functions (see Hayter (2014), for example).

General discussions of stochastic control can be found in Wendell & Rishel (1975) and Øksendal (2014), for example. Moreover, in finance the problem of option pricing is also considered a stochastic control problem. Fusai & Meucci (2008) discuss pricing discretely monitored Asian options, and recursive integration techniques in pricing barrier options have been discussed in Aitsahlia & Lai (1997), Sullivan (2000), Andricopoulos et al. (2003) and Fusai & Recchioni (2007).

The results obtained in this paper can also be used to calculate conditional probabilistic expressions and moments for the stochastic behavior of these partial sum processes. For example, probabilities for the independent random variables conditioned on an event can also be tractable since

 P((X1,…,Xn)∈C|(X1,…,Xn)∈A)=P((X1,…,Xn)∈C∩A)P((X1,…,Xn)∈A) (5)

where the numerator is tractable for certain sets . In particular, if the set can also be defined in terms of the partial sums as

 (X1+…+Xi,Xi+1,Xi+2)∈Ci⊆R3

for , then both the numerator and denominator of equation (5) are of the form of equation (4) with the sets equal to the sets or and the functions equal to the probability density functions .

Furthermore, it can be noted that the moments and covariances of the independent random variables conditioned on an event are also tractable since

 E[Xr11…Xrnn∣(X1,…,Xn)∈A]=DP((X1,…,Xn)∈A)

where is of the form of equation (4) with the sets equal to the sets and the functions equal to .

The layout of this paper is as follows. In section 2 it is shown how the integral in equation (4) can be evaluated recursively as a series of -dimensional integral calculations. Recursive formulas are given for the general case, and also for the special case of equation (2). In section 3 a discussion is provided of the implementation details of the methodology. The adoption of Fast Fourier Transforms is illustrated as a way to improve the computational efficiency of the methodology, and an error analysis of the numerical integrations is provided. Some illustrations of the implementation of the methodology are provided in section 4, with examples in the fields of reliability, product quality assessment, and stochastic control. Finally, a summary is provided in section 5.

## 2 Recursive Integration Methodology

In this section the recursive integration of the integral in equation (4) is discussed. First a change of variables is used to put the expression into a more convenient form, and then the general recursive formulas are provided. The special case of equation (2) is then considered separately. It should be remembered that if the random variables have discrete distributions, then the results of this section are still applicable with the integrals replaced by sums, and the probability density functions replaced by the probability mass functions (see Hayter (2014), for example).

### 2.1 Change of Variables

If the change of variables , , is employed, then equation (4) becomes

 ∫⋯∫(y1,…,yn)∈Ψh1(y1)h2(y2−y1)⋯hn(yn−yn−1) dy1⋯dyn (6)

where the set is defined by the conditions

 (yi,yi+1,yi+2)∈Ji⊆R3

for , and where the set is derived from the set through the relationship

 (x1+…+xi,xi+1,xi+2)∈Ii⇔(yi,yi+1,yi+2)∈Ji.

Notice that in equation (6) the integrand is the product of terms that only involve two adjacent , while the integration region is defined by conditions on only three adjacent . Consequently, equation (6) is of the form given in section 1 of Hayter (2006) with , which implies that it can be evaluated recursively by a series of 2-dimensional integral calculations. Specific formulas for this recursive integration are now provided.

### 2.2 General Recursive Formulas

Let

 Jk(⋅,u,v)={x∈R:(x,u,v)∈Jk}.

We assume that can always be represented by a finite union of disjoint closed intervals, so that

 Jk(⋅,u,v)=∪i[ak,i(u,v),bk,i(u,v)]. (7)

 J12k={(x,y)∈R2:∃z∈R∋(x,y,z)∈Jk},

and

 J23k={(y,z)∈R2:∃x∈R∋(x,y,z)∈Jk}.

To compute equation (6), at each first evaluate

 G1(a,u)=∫a−∞h1(x)h2(u−x)dx, (8)

for all , and then compute

 g1(u,v)=∫J1(⋅,u,v)h1(x)h2(u−x)dx=∑i[G1(b1,i(u,v),u)−G1(a1,i(u,v),u)]. (9)

Next, for , at each , and for at each , evaluate

 Gk(a,u)=∫a−∞gk−1(x,u)hk+1(u−x)dx, (10)

for all , and letting for but . Then compute

 gk(u,v)=∫Jk(⋅,u,v)gk−1(x,u)hk+1(u−x)dx=∑i[Gk(bk,i(u,v),u)−Gk(ak,i(u,v),u)]. (11)

Finally, the evaluation of equation (6), and hence of equation (4), is obtained as

 P((X1,…,Xn)∈A)=∬J23n−2gn−2(u,v)hn(v−u)dudv. (12)

Notice that the steps in this evaluation each have the computational intensity of a two-dimensional numerical integration.

### 2.3 Recursive Formulas for the Sum of Independent Random Variables

Now consider the special case where the are independent random variables with probability density functions , respectively. Recursive formulas are now provided for evaluating some probabilistic properties of . First, notice that it follows from equation (6) that

 P(T≤τ)=∫⋯∫yn≤τf1(y1)f2(y2−y1)⋯fn(yn−yn−1)dy1⋯dyn. (13)

This expression can be computed simply by first evaluating

 g1(u)=∫∞−∞f1(x)f2(u−x)dx

at each . Then, sequentially, for , evaluate

 gk(u)=∫∞−∞gk−1(x)fk+1(u−x)dx

at each . Again, notice that the steps in this evaluation each have the computational intensity of a two-dimensional numerical integration. Finally, the required expression is obtained as

 P(T≤τ)=∫τ−∞gn−1(x)dx.

To compute the expectation of conditional on , or , first evaluate

 g11(u)=∫∞−∞w(x)f1(x)f2(u−x)dx.

Then, sequentially, for , evaluate

 g1k(u)=∫∞−∞g1k−1(x)fk+1(u−x)dx

for each . The expectation of conditional on can then be obtained as

 E[w(X1)∣T≤τ]=∫τ−∞g1n−1(x)dxP(T≤τ),

while the expectation of conditional on can be obtained as

 E[w(X1)∣T≥τ]=∫∞τg1n−1(x)dxP(T≥τ).

Notice that expectations for can be obtained from these expressions by reordering the indices of the .

To compute the expectation of conditional on , or , first evaluate

 g21(u)=∫∞−∞w1(x)w2(u−x)f1(x)f2(u−x)dx.

Then, sequentially, for , evaluate

 g2k(u)=∫∞−∞g2k−1(x)fk+1(u−x)dx

for each . The expectation of conditional on can then be obtained as

 E[w1(X1)w2(X2)∣T≤τ]=∫τ−∞g2n−1(x)dxP(T≤τ),

while the expectation of conditional on can be obtained as

 E[w1(X1)w2(X2)∣T≥τ]=∫∞τg2n−1(x)dxP(T≥τ).

Again, expectations for can be obtained from these expressions by reordering the indices of the .

Finally, notice that the expectation of conditional on either or is

 E[T∣T≤τ]=n∑i=1E[Xi∣T≤τ]andE[T∣T≥τ]=n∑i=1E[Xi∣T≥τ],

which becomes

 E[T∣T≤τ]=nE[X1∣T≤τ]andE[T∣T≥τ]=nE[X1∣T≥τ]

when the are identically distributed.

## 3 Implementation details

In this section a discussion is provided of the implementation details of the methodology. The adoption of Fast Fourier Transforms (see Carverhill & Clewlow (1990), for example) is illustrated as a way to improve the computational efficiency of the methodology, and an error analysis of the numerical integrations is provided.

With variables a direct implementation of the methodology requires a calculation with a computational intensity that is equivalent to a sequence of two-dimensional numerical integrations. This is already efficient considering that the original problem is ostensibly an -dimensional numerical integration. However, in the case when the limits of integration in section 2.2 are invariant over pairs ,the computation can be accelerated even further using a Fast Fourier Transform convolution.

### 3.1 Fast Fourier Transform Convolution

It can be observed that the recursive integration formulas given in section 2 involve the convolution of two functions. Consequently, in some cases the speed of the computation can be increased with a Fast Fourier Transform technique (FFT). As the well-known Convolution Theorem states (see, for example, Smith (2007)), a convolution with respect to the variable in the original domain is equivalent to multiplication with respect to the variable in the transformed domain.

More formally, letting denote the Discrete Fourier Transform (DFT) and its inverse, convolutions between two functions and can be computed as

 f∗g=F−1(F(f)⋅F(g)).

The functions are decomposed into the transformed domain using the DFT, multiplied in the transformed domain, and then transformed back into the original domain using the inverse DFT.

Notice that the DFT and its inverse can be calculated by the FFT algorithm. Using a grid size of , the overall computational intensity of conducting the convolution in this way using the FFT is (see, for example, Smith (2007)), which is lower than the computational intensity obtained with the direct computation of the convolution in the time domain. The comparative accuracies and efficiencies of the two methods are now demonstrated.

### 3.2 Accuracy and Efficiency

In order to illustrate and compare the accuracies and efficiencies of the implementations of the recursive integration formula introduced in section 2.2

, the formula is applied to the calculation of the cumulative distribution function, the conditional cumulative distribution function, and the conditional expectation of the sum of 10 independent identically distributed exponential random variables with parameter

= 1. In this case the sum of these random variables has a known gamma distribution, so that the exact values of the calculated quantities are known.

The formulas in section 2.2 are implemented with a truncation of the support at 30. Table 1 shows the computed values and the errors of the required quantities, together with their computational times, obtained from implementations with the direct convolution and the FFT convolution. Different grid sizes are used, and both methods are implemented with SciPy’s Python library (see Jones et al. (2001)) with an Intel Core i5 CPU.

It can be seen from Table 1 that the two implementations have similar errors, but the implementation with the FFT convolution is significantly faster. Consequently, it is useful to apply the FFT technique to the evaluation of the recursive formulas given in section 2.2 when the limits of integration are invariant over pairs .

## 4 Examples and Illustrations

In this section the methodology presented in this paper is illustrated through applications to problems in the fields of reliability, product quality assessment, and stochastic control that require probability calculations for partial sums of independent random variables. The first example concerns a reliability problem where failed components are successively replaced with new components, while the second example concerns a product quality assessment problem where batches are evaluated based on a measurement of the sum of their individual items. Finally, the third problem concerns discrete time stochastic control.

### 4.1 Reliability Example

Suppose that a machine contains “identical” components which are deployed successively. Thus, the first component is deployed until it fails, whereupon the second component is deployed, and so on. The machine operates until the th component has failed. Furthermore, suppose that an observer can tell whether or not the machine is operating, but not how many components have failed if the machine is still operating.

If the component lifetimes are taken to be independent with specified distributions, then the methodology presented in this paper can be used to investigate the probabilistic properties of the lifetime of the machine. Some illustrative calculations are provided when the component lifetimes are taken to be independent identically distributed Weibull distributions. Without the methodologies presented here, calculations on the sum of Weibull distributions are generally intractable and would usually be assessed with simulations.

The following are examples of the kinds of probability calculations that can be performed using the recursive integration methodology presented in this paper. If the component lifetimes are with distributions , so that the machine lifetime is , then an obvious quantity of interest is the machine survival function

 P(T≥t).

If the machine is observed at time , then if the machine is still operating the conditional survival function is

 P(T≥t∣T≥τ)=P(T≥t)P(T≥τ).

If the machine has failed at time then the conditional survival function is

 P(T≥t∣T≤τ)=P(t≤T≤τ)P(T≤τ).

The expected failure time of the machine is simply times the individual expected component failure time, but if the machine is observed to be still operating at time , then the conditional expected failure time is

 ∫∞t=τtf(t∣t≥τ)dt=GP(T≥τ)

where is the conditional distribution of the failure time and

 G=∫…x1+…+xn≥τ∫(x1+…+xn)f1(x1)…fn(xn) dx1…dxn.

This can be evaluated as the sum of

separate integrals which are identical if the component lifetimes are identically distributed. The variance of the conditional failure time can be obtained by having

in place of in the integrand, so that can be found from terms with and in the integrand.

Finally, if the machine is observed to be still operating at time , then the distribution of the number of failed components at time can be obtained, for , as

 P(no more than i−1 components have failed by time τ)=
 P(X1+…+Xi≥τ∣T≥τ)=P(X1+…+Xi≥τ)P(T≥τ).

Table 2 shows the computed results (with computation times using the Fast Fourier Transform technique) of these probabilities when

are independent, identically distributed Weibull random variables with shape parameter equal to 2 and scale parameter equal to 1. These random variables have an expectation of 0.886 and a standard deviation of 0.463.

### 4.2 Product Quality Example

Consider a product quality assessment problem where a measureable property of an item is satisfactory if it is no smaller than a specified level . Let , , represent the values of these properties for a batch of items, and suppose that they can be modelled as being independent with an identical probability density function .

Suppose that instead of the costly approach of testing each item in the batch, it is possible and simple to obtain information about the sum . This is the case, say, if the weight of the item is of interest or the radiation emitted from the item. It is useful to be able to make probability statements about the number of satisfactory items in the batch based upon the information obtained about . In practice, the exact value of may be observed, or a lower or an upper bound may be obtained.

If the exact value of is observed then

 P(exactly i items are satisfactory∣T)=
 (ni)P(X1≥c,…,Xi≥c,Xi+1

where

 H2=∫…x1+…+xn=T∫f(x1)…f(xn) dx1…dxn

and

 H1=∫…x1+…+xn=Tx1≥c,…,xi≥cxi+1

As an illustration, some calculations are shown when , and

is taken to be a Laplace (double exponential) distribution with parameter

, so that

 f(x)=12e−|x|,x∈R.

Table 3 shows the computed values of at different and . The computational time of each entry using the Fast Fourier Transform technique was about 0.3 seconds.

If the bounds or are observed rather than the exact value of , then the expressions for and can be modified so that the integration regions depend on the conditions or . In either case and can be again be evaluated using the recursive integration methodologies presented in this paper.

### 4.3 Discrete Time Stochastic Control Example.

This section illustrates the application of the methodology developed in this paper to a discrete time stochastic control problem. Let , , be the performance measurement of a process at discrete times , where the are non-negative and assumed to be independent and identically distributed when the process is operating correctly. The objective is to dynamically track the partial means of the over time, and to detect any increase in the mean of the by a certain decision rule.

For , denote the partial means up to by

 ¯Xn=∑ni=1Xin.

Suppose that for each , the process is stopped when both and are greater than , for a certain control limit . If the process is not stopped prior to , then the process is deemed to have been operating correctly throughout the time horizon . For a specified distribution of the , it is required to calculate the value of that provides a probability of of not incorrectly stopping the process within the horizon .

The control limit can be obtained by searching for the value of that is the solution to the equation

 P(Xk+1≤¯Xk+c∗ or Xk+2≤¯Xk+c∗, % for k=1,…,N−2)=1−α. (14)

The event in equation (14) is the event that the process is not terminated within the time horizon. This event is in the form of equation (1), which can be computed by the formula in equation (6).

In order to compute equation (6), the in equation(7) have and as the transformed variables

 Yk+1=k+1∑i=1Xi and Yk+2=k+2∑i=1Xi.

Furthermore, given and , the process is in control at time if

 Yk≥kk+1(u−c∗),

or

 Yk≥k(v−u−c∗).

In addition, since the are non-negative random variables it follows that

 Yk≤u,for k=1,…,N−1.

Therefore,

 Jk(⋅,u,v)=[ak(u,v),bk(u,v)],

where

 ak(u,v)=min{kk+1(u−c∗), k(v−u−c∗)}

and

 bk(u,v)=u.

Notice that in this case the Fast Fourier Transform technique cannot be used because the limits of the integrals and vary over and .

To obtain the required control limit , the probability in equation (14) has to be computed at several values of in order to search for the solution. Consequently, for a large time horizon it is essential that an efficient computation methodology, as developed in this paper, is available in order to obtain in practice.

Table 4 shows the control limit for different values of and , together with computational times using the recursive integration methodology developed in this paper, for the case where the are independent, identically distributed exponential random variables with scale parameter equal to 1.

## 5 Summary

The tracking of the stochastic behavior of a partial sum process is an important problem. There are many applications of partial sum processes, and in this paper examples have been provided in the fields of reliability, product quality assessment, and stochastic control.

It has been shown how calculations of the probabilistic properties of such a partial sum process, which ostensibly require an ability to compute high-dimensional multivariate probabilities, and so are consequently intractable in general, can in fact be solved as a sequence of two dimensional computations, with each computation being the convolution of two functions.

Finally, it has been shown how the Fast Fourier Transform technique can be utilized for the evaluation of these convolutions in some cases. The results of this paper allow the efficient computation of the probabilistic properties of many important partial sum processes.

References

Aitsahlia, F. and Lai, T. L., 1997. “Valuation of discrete barrier and hindsight options,” The Journal of Financial Engineering, 6 (2), 169-177.

Andricopoulos, A.D., Widdicks, M., Duck, P.W. and Newton, D.P., 2003. “Universal option valuation using quadrature methods,” Journal of Financial Economics, 67, 447-471.

Belisle, C. J. P., Romeijn, H. E. and Smith, R.L., 1993. “Hit-and-run algorithm for generating multivariate distribution,” Mathematics of Operations Research, 18, 255-266.

Carverhill, A.P. and Clewlow, L.J., 1990. “Flexible convolution”, Risk, 3, 25- 29.

Dunnett, C. W. and Sobel, M., 1955. “Approximations to the probability integral and certain percentage points of a multivariate analogue of Student’s -distribution,” Biometrika, 42, 258-260.

Fleming, W. H. and Rishel, R. W., 1975. “ Deterministic and Stochastic Optimal Control,” Springer, New York.

Fusai, G. and Meucci, A., 2008. “Pricing discretely monitored Asian options under Levy processes,” Journal of Banking and Finance, 32, 2076-2088.

Fusai, G. and Recchioni, M.C., 2007. “Analysis of quadrature methods for pricing discrete barrier options,” Journal of Economic Dynamics and Control, 31 (3), 826-860.

Hayter, A. J., 2006. “Recursive integration methodologies with statistical applications,” Journal of Statistical Planning and Inference, 136, 2284-2296.

Hayter, A. J., 2014. “Recursive formulas for multinomial probabilities with applications,” Computational Statistics, 29 (5), 1207-1219.

Jones E., Oliphant E., Peterson P., et al., 2001. “SciPy: Open Source Scientific Tools for Python,” http://www.scipy.org/.

Kiatsupaibul, S., Hayter, A. J. and Wei, L., 2017. “Rank constrained distribution and moment computations,” Computational Statistics and Data Analysis, 105, 229-242.

Kiatsupaibul, S. and Hayter, A. J., 2015. “Recursive confidence band construction for an unknown distribution function,” Biometrical Journal, 57 (1), 39-51.

Kiatsupaibul, S., Smith, R. L. and Zabinsky, Z. B., 2011. “An analysis of a variation of hit-and-run for uniform sampling from general region,” ACM Transactions on Modeling and Computer Simulation, 21 (3), Article 16.

Lovász, L., 1999. “Hit-and-run mixes fast,” Mathematical Programming, 86, 443-461.

Øksendal, B., 2014. “Stochastic Differential Equations”, 6th Edition, Springer, Heidelberg.

Smith, J. O., 2007. “Mathematics of the Discrete Fourier Transform (DFT), with Audio Applications”, Second Edition, W3K Publishing.

Smith, R. L., 1984. “Efficient Monte Carlo procedures for generating points uniformly distributed over bounded regions,”

Operations Research, 32, 1296-1308.

Soong, W. C., and Hsu, J. C., 1997. “Using complex integration to compute multivariate normal probabilities,” Journal of Computational and Graphical Statistics, 6 (4), 397-415.

Sullivan, M. A., 2000. “Pricing discretely monitored barrier options”, Journal of Computational Finance, 3 (4), 35-52.