Efficient Computation of the Stochastic Behavior of Partial Sum Processes

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

t

-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 $X_{i}$ , $1 \leq i \leq n$ be independent random variables. In general, we consider probabilities of the form

P ((X_{1}, \dots, X_{n}) \in A) =

P ((X_{1}, X_{2}, X_{3}) \in A_{1}, (X_{1} + X_{2}, X_{3}, X_{4}) \in A_{2}, \dots, (X_{1} + \dots + X_{n - 2}, X_{n - 1}, X_{n}) \in A_{n - 2})

(1)

for sets $A_{i} \subseteq R^{3}$ , $1 \leq i \leq n - 2$ . 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 (X_{1} + \dots + X_{n} \in B)

(2)

for a set $B \in R$ 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 $n$ .

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

P (X_{1} \in B_{1}, X_{1} + X_{2} \in B_{2}, \dots, X_{1} + \dots + X_{n} \in B_{n})

(3)

for sets $B_{i} \subseteq R$ , $1 \leq i \leq n$ , 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 $n$ .

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

\int \dots \int (x_{1}, \dots, x_{n}) \in S h_{1} (x_{1}) \dots h_{n} (x_{n}) d x_{1} \dots d x_{n}

(4)

can be evaluated recursively as a series of $2$ -dimensional integral calculations when the set $S \subseteq R^{n}$ is defined by the conditions

(x_{1} + \dots + x_{i}, x_{i + 1}, x_{i + 2}) \in I_{i} \subseteq R^{3}

for $1 \leq i \leq n - 2$ . This is an application of the general discussion of recursive integration given in Hayter (2006) with $d = 2$ . 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

I_{i}

equal to the sets

A_{i}

and the functions

h_{i} (x_{i})

equal to the probability density functions

f_{i} (x_{i})

. In addition, if the random variables

X_{i}

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 $X_{i}$ conditioned on an event $A$ can also be tractable since

P ((X_{1}, \dots, X_{n}) \in C | (X_{1}, \dots, X_{n}) \in A) = \frac{P ((X_{1}, \dots, X_{n}) \in C \cap A)}{P ((X_{1}, \dots, X_{n}) \in A)}

(5)

where the numerator is tractable for certain sets $C$ . In particular, if the set $C \subseteq R^{n}$ can also be defined in terms of the partial sums as

(X_{1} + \dots + X_{i}, X_{i + 1}, X_{i + 2}) \in C_{i} \subseteq R^{3}

for $1 \leq i \leq n - 2$ , then both the numerator and denominator of equation (5) are of the form of equation (4) with the sets $I_{i}$ equal to the sets $A_{i} \cap C_{i}$ or $A_{i}$ and the functions $h_{i} (x_{i})$ equal to the probability density functions $f_{i} (x_{i})$ .

Furthermore, it can be noted that the moments and covariances of the independent random variables $X_{i}$ conditioned on an event $A$ are also tractable since

E [X_{1}^{r_{1}} \dots X_{n}^{r_{n}} ∣ (X_{1}, \dots, X_{n}) \in A] = \frac{D}{P ((X_{1}, \dots, X_{n}) \in A)}

where $D$ is of the form of equation (4) with the sets $I_{i}$ equal to the sets $A_{i}$ and the functions $h_{i} (x_{i})$ equal to $x_{i}^{r_{i}} f_{i} (x_{i})$ .

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 $2$ -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 $X_{i}$ 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 $y_{i} = x_{1} + \dots + x_{i}$ , $1 \leq i \leq n$ , is employed, then equation (4) becomes

\int \dots \int (y_{1}, \dots, y_{n}) \in Ψ h_{1} (y_{1}) h_{2} (y_{2} - y_{1}) \dots h_{n} (y_{n} - y_{n - 1}) d y_{1} \dots d y_{n}

(6)

where the set $Ψ \subseteq R^{n}$ is defined by the conditions

(y_{i}, y_{i + 1}, y_{i + 2}) \in J_{i} \subseteq R^{3}

for $1 \leq i \leq n - 2$ , and where the set $J_{i}$ is derived from the set $I_{i}$ through the relationship

(x_{1} + \dots + x_{i}, x_{i + 1}, x_{i + 2}) \in I_{i} \Leftrightarrow (y_{i}, y_{i + 1}, y_{i + 2}) \in J_{i} .

Notice that in equation (6) the integrand is the product of terms that only involve two adjacent $y_{i}$ , while the integration region is defined by conditions on only three adjacent $y_{i}$ . Consequently, equation (6) is of the form given in section 1 of Hayter (2006) with $d = 2$ , 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

J_{k} (\cdot, u, v) = {x \in R : (x, u, v) \in J_{k}} .

We assume that $J_{k} (\cdot, u, v)$ can always be represented by a finite union of disjoint closed intervals, so that

J_{k} (\cdot, u, v) = \cup_{i} [a_{k, i} (u, v), b_{k, i} (u, v)] .

(7)

In addition, let

J_{k}^{12} = {(x, y) \in R^{2} : \exists z \in R ∋ (x, y, z) \in J_{k}},

and

J_{k}^{23} = {(y, z) \in R^{2} : \exists x \in R ∋ (x, y, z) \in J_{k}} .

To compute equation (6), at each $(u, v) \in J_{2}^{12}$ first evaluate

G_{1} (a, u) = \int_{- \infty}^{a} h_{1} (x) h_{2} (u - x) d x,

(8)

for all $a \in \cup_{v} \cup_{i} [a_{1, i} (u, v), b_{1, i} (u, v)]$ , and then compute

g_{1} (u, v) = \int_{J_{1} (\cdot, u, v)} h_{1} (x) h_{2} (u - x) d x = \sum i [G_{1} (b_{1, i} (u, v), u) - G_{1} (a_{1, i} (u, v), u)] .

(9)

Next, for $k = 2, \dots, n - 3$ , at each $(u, v) \in J_{k + 1}^{12}$ , and for $k = n - 2$ at each $(u, v) \in J_{k}^{23}$ , evaluate

G_{k} (a, u) = \int_{- \infty}^{a} g_{k - 1} (x, u) h_{k + 1} (u - x) d x,

(10)

for all $a \in \cup_{v} \cup_{i} [a_{k, i} (u, v), b_{k, i} (u, v)]$ , and letting $g_{k} (u, v) = 0$ for $(u, v) \in J_{k}^{23}$ but $(u, v) \notin J_{k + 1}^{12}, k = 2, \dots, n - 3$ . Then compute

g_{k} (u, v) = \int_{J_{k} (\cdot, u, v)} g_{k - 1} (x, u) h_{k + 1} (u - x) d x = \sum i [G_{k} (b_{k, i} (u, v), u) - G_{k} (a_{k, i} (u, v), u)] .

(11)

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

P ((X_{1}, \dots, X_{n}) \in A) = \iint J_{n - 2}^{23} g_{n - 2} (u, v) h_{n} (v - u) d u d v .

(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 $X_{1}, \dots, X_{n}$ are independent random variables with probability density functions $f_{1}, \dots, f_{n}$ , respectively. Recursive formulas are now provided for evaluating some probabilistic properties of $T = X_{1} + \dots + X_{n}$ . First, notice that it follows from equation (6) that

P (T \leq τ) = \int \dots \int y_{n} \leq τ f_{1} (y_{1}) f_{2} (y_{2} - y_{1}) \dots f_{n} (y_{n} - y_{n - 1}) d y_{1} \dots d y_{n} .

(13)

This expression can be computed simply by first evaluating

g_{1} (u) = \int_{- \infty}^{\infty} f_{1} (x) f_{2} (u - x) d x

at each $u \in R$ . Then, sequentially, for $k = 2, \dots, n - 1$ , evaluate

g_{k} (u) = \int_{- \infty}^{\infty} g_{k - 1} (x) f_{k + 1} (u - x) d x

at each $u \in R$ . 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 \leq τ) = \int_{- \infty}^{τ} g_{n - 1} (x) d x .

To compute the expectation of $w (X_{1})$ conditional on $T \leq τ$ , or $T \geq τ$ , first evaluate

g_{1}^{1} (u) = \int_{- \infty}^{\infty} w (x) f_{1} (x) f_{2} (u - x) d x .

Then, sequentially, for $k = 2, \dots, n - 1$ , evaluate

g_{k}^{1} (u) = \int_{- \infty}^{\infty} g_{k - 1}^{1} (x) f_{k + 1} (u - x) d x

for each $u \in R$ . The expectation of $w (X_{1})$ conditional on $T \leq τ$ can then be obtained as

E [w (X_{1}) ∣ T \leq τ] = \frac{\int_{- \infty}^{τ} g_{n - 1}^{1} (x) d x}{P (T \leq τ)},

while the expectation of $X_{1}$ conditional on $T \geq τ$ can be obtained as

E [w (X_{1}) ∣ T \geq τ] = \frac{\int_{τ}^{\infty} g_{n - 1}^{1} (x) d x}{P (T \geq τ)} .

Notice that expectations for $w (X_{i})$ can be obtained from these expressions by reordering the indices of the $X_{i}$ .

To compute the expectation of $w_{1} (X_{1}) w_{2} (X_{2})$ conditional on $T \leq τ$ , or $T \geq τ$ , first evaluate

g_{1}^{2} (u) = \int_{- \infty}^{\infty} w_{1} (x) w_{2} (u - x) f_{1} (x) f_{2} (u - x) d x .

Then, sequentially, for $k = 2, \dots, n - 1$ , evaluate

g_{k}^{2} (u) = \int_{- \infty}^{\infty} g_{k - 1}^{2} (x) f_{k + 1} (u - x) d x

for each $u \in R$ . The expectation of $w_{1} (X_{1}) w_{2} (X_{2})$ conditional on $T \leq τ$ can then be obtained as

E [w_{1} (X_{1}) w_{2} (X_{2}) ∣ T \leq τ] = \frac{\int_{- \infty}^{τ} g_{n - 1}^{2} (x) d x}{P (T \leq τ)},

while the expectation of $w_{1} (X_{1}) w_{2} (X_{2})$ conditional on $T \geq τ$ can be obtained as

E [w_{1} (X_{1}) w_{2} (X_{2}) ∣ T \geq τ] = \frac{\int_{τ}^{\infty} g_{n - 1}^{2} (x) d x}{P (T \geq τ)} .

Again, expectations for $w_{i} (X_{i}) w_{i} (X_{j})$ can be obtained from these expressions by reordering the indices of the $X_{i}$ .

Finally, notice that the expectation of $T$ conditional on either $T \leq τ$ or $T \geq τ$ is

E [T ∣ T \leq τ] = n \sum i = 1 E [X_{i} ∣ T \leq τ] and E [T ∣ T \geq τ] = n \sum i = 1 E [X_{i} ∣ T \geq τ],

which becomes

E [T ∣ T \leq τ] = n E [X_{1} ∣ T \leq τ] and E [T ∣ T \geq τ] = n E [X_{1} ∣ T \geq τ]

when the $X_{i}$ 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 $n$ variables $X_{i}$ a direct implementation of the methodology requires a calculation with a computational intensity that is equivalent to a sequence of $n$ two-dimensional numerical integrations. This is already efficient considering that the original problem is ostensibly an $n$ -dimensional numerical integration. However, in the case when the limits of integration $\cup_{i} {a_{k, i} (u, v), b_{k, i} (u, v)}$ in section 2.2 are invariant over pairs $(u, v)$ ,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 $F$ denote the Discrete Fourier Transform (DFT) and $F^{- 1}$ its inverse, convolutions between two functions $f$ and $g$ can be computed as

f * g = F^{- 1} (F (f) \cdot 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 $N$ , the overall computational intensity of conducting the convolution in this way using the FFT is $O (N log N)$ (see, for example, Smith (2007)), which is lower than the computational intensity $O (N^{2})$ 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.

		Direct Convolution			FFT Convolution
Grid size		Value	Error	Time (sec)	Value	Error	Time (sec)
0.01	$P (T \geq 12)$	0.24275	0.00036	0.030	0.24275	0.00036	0.017
	$P (T \geq 12 ∣ T \geq 10)$	0.52945	0.00013	0.028	0.52945	0.00013	0.017
	$E [X_{1} ∣ T \geq 10]$	1.27289	0.00032	0.055	1.27289	0.00032	0.033
0.001	$P (T \geq 12)$	0.24239	0.00000	2.376	0.24239	0.00000	0.098
	$P (T \geq 12 ∣ T \geq 10)$	0.52932	0.00000	2.342	0.52932	0.00000	0.092
	$E [X_{1} ∣ T \geq 10]$	1.27320	0.00001	4.719	1.27320	0.00001	0.185
0.0001	$P (T \geq 12)$	0.24240	0.00000	388.769	0.24240	0.00000	1.235
	$P (T \geq 12 ∣ T \geq 10)$	0.52932	0.00000	408.467	0.52932	0.00000	1.191
	$E [X_{1} ∣ T \geq 10]$	1.27320	0.00001	838.920	1.27320	0.00001	2.285

Table 1: Comparisons of the implementation methods for the methodology for a sum of 10 independent identically distributed exponential random variables with parameter

λ = 1

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 $\cup_{i} {a_{k, i} (u, v), b_{k, i} (u, v)}$ are invariant over pairs $(u, v)$ .

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 $n$ “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 $n$ 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 $X_{i}$ with distributions $f_{i} (x_{i})$ , so that the machine lifetime is $T = X_{1} + \dots + X_{n}$ , then an obvious quantity of interest is the machine survival function

P (T \geq t) .

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

P (T \geq t ∣ T \geq τ) = \frac{P (T \geq t)}{P (T \geq τ)} .

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

P (T \geq t ∣ T \leq τ) = \frac{P (t \leq T \leq τ)}{P (T \leq τ)} .

The expected failure time of the machine is simply $n$ 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

\int_{t = τ}^{\infty} t f (t ∣ t \geq τ) d t = \frac{G}{P (T \geq τ)}

where $f (t | t \geq τ)$ is the conditional distribution of the failure time and

G = \int \begin{matrix} \dots x_{1} + \dots + x_{n} \geq τ \end{matrix} \int (x_{1} + \dots + x_{n}) f_{1} (x_{1}) \dots f_{n} (x_{n}) d x_{1} \dots d x_{n} .

This can be evaluated as the sum of $n$

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

t^{2}

in place of

t

in the integrand, so that

G

can be found from terms with

x_{i}^{2}

and

x_{i} x_{j}

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 $1 \leq i \leq n - 1$ , as

P (no more than i - 1 components have failed by time τ) =

P (X_{1} + \dots + X_{i} \geq τ ∣ T \geq τ) = \frac{P (X_{1} + \dots + X_{i} \geq τ)}{P (T \geq τ)} .

Table 2 shows the computed results (with computation times using the Fast Fourier Transform technique) of these probabilities when $X_{1}, \dots, X_{10}$

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.

	Computed value	Computational time (sec)
$P (T \geq 8)$	$0.7139490$	$1.876$
$P (T \geq 10)$	$0.2154629$	$2.040$
$P (T \geq 12)$	$0.0206421$	$1.926$
$P (T \geq 12 ∣ T \geq 10)$	$0.0958036$	$1.701$
$P (8 \leq T < 10 ∣ T \leq 10)$	$0.6353888$	$1.617$
$P (X_{1} + \dots + X_{7} \geq 10 ∣ T \geq 10)$	$0.0104016$	$1.776$
$E[T]\textsuperscript\textdagger$	$8.8627912$	$1.797$
$E [T ∣ T \geq 10]$	$12.3020396$	$3.664$
^† $E [T]$ the exact value is equal to $10 Γ (1.5)$

Table 2: Computed results and computation times for reliability example.

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 $c$ . Let $X_{i}$ , $1 \leq i \leq n$ , represent the values of these properties for a batch of $n$ items, and suppose that they can be modelled as being independent with an identical probability density function $f (x)$ .

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 $T = X_{1} + \dots + X_{n}$ . 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 $T$ . In practice, the exact value of $T$ may be observed, or a lower or an upper bound may be obtained.

If the exact value of $T$ is observed then

P (exactly i items are satisfactory ∣ T) =

(\frac{n}{i}) P (X_{1} \geq c, \dots, X_{i} \geq c, X_{i + 1} < c, \dots, X_{n} < c | T) = (\frac{n}{i}) \frac{H_{1}}{H_{2}}

where

H_{2} = \int \begin{matrix} \dots x_{1} + \dots + x_{n} = T \end{matrix} \int f (x_{1}) \dots f (x_{n}) d x_{1} \dots d x_{n}

and

H_{1} = \int \begin{matrix} \dots x_{1} + \dots + x_{n} = T x_{1} \geq c, \dots, x_{i} \geq c x_{i + 1} < c, \dots, x_{n} < c \end{matrix} \int f (x_{1}) \dots f (x_{n}) d x_{1} \dots d x_{n} .

As an illustration, some calculations are shown when $n = 10$ , $c = 1$ and $f (x)$

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

λ = 1

, so that

f (x) = \frac{1}{2} e^{- | x |}, x \in R .

Table 3 shows the computed values of $P (exactly i items are satisfactory ∣ T)$ at different $i = 0, 1, \dots, 10$ and $T = 0, 5, 10, 15, 20$ . The computational time of each entry using the Fast Fourier Transform technique was about 0.3 seconds.

If the bounds $T \leq t$ or $T \geq t$ are observed rather than the exact value of $T$ , then the expressions for $H_{1}$ and $H_{2}$ can be modified so that the integration regions depend on the conditions $x_{1} + \dots + x_{n} \leq t$ or $x_{1} + \dots + x_{n} \geq t$ . In either case $H_{1}$ and $H_{2}$ can be again be evaluated using the recursive integration methodologies presented in this paper.

	$T$
i	0	5	10	15	20
0	0.0774	0.0004	0.0000	0.0000	0.0000
1	0.3629	0.0518	0.0024	0.0002	0.0000
2	0.3896	0.2960	0.0477	0.0076	0.0016
3	0.1461	0.4176	0.2315	0.0688	0.0213
4	0.0225	0.1971	0.3905	0.2374	0.1135
5	0.0015	0.0347	0.2560	0.3568	0.2771
6	0.0000	0.0023	0.0656	0.2443	0.3310
7	0.0000	0.0001	0.0061	0.0751	0.1948
8	0.0000	0.0000	0.0002	0.0094	0.0542
9	0.0000	0.0000	0.0000	0.0004	0.0063
10	0.0000	0.0000	0.0000	0.0000	0.0002

Table 3: The probability of exactly

i

items having a satisfactory weight (weight greater than

c = 1

) given an observed total weight

T

n = 10

items. The items are assumed to have independent and identically distributed weights with a Laplace distribution with parameter

λ = 1

. The computational time of each entry using the Fast Fourier Transform technique was about 0.7 seconds.

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 $X_{i}$ , $i = 1, \dots, N$ , be the performance measurement of a process at discrete times $i$ , where the $X_{i}$ 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 $X_{i}$ over time, and to detect any increase in the mean of the $X_{i}$ by a certain decision rule.

For $n = 1, \dots, N$ , denote the partial means up to $n$ by

{¯ X}_{n} = \frac{\sum_{i = 1}^{n} X_{i}}{n} .

Suppose that for each $n = 3, \dots, N$ , the process is stopped when both $X_{n}$ and $X_{n - 1}$ are greater than ${¯ X}_{n - 2} + c (α, N)$ , for a certain control limit $c (α, N)$ . If the process is not stopped prior to $N$ , then the process is deemed to have been operating correctly throughout the time horizon $N$ . For a specified distribution of the $X_{i}$ , it is required to calculate the value of $c (α, N)$ that provides a probability of $1 - α$ of not incorrectly stopping the process within the horizon $N$ .

The control limit $c (α, N)$ can be obtained by searching for the value of $c^{*}$ that is the solution to the equation

P (X_{k + 1} \leq {¯ X}_{k} + c^{*} or X_{k + 2} \leq {¯ X}_{k} + c^{*}, % for k = 1, \dots, 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 $J_{k} (\cdot, u, v)$ in equation(7) have $u$ and $v$ as the transformed variables

Y_{k + 1} = k + 1 \sum i = 1 X_{i} and Y_{k + 2} = k + 2 \sum i = 1 X_{i} .

Furthermore, given $Y_{k + 1} = u$ and $Y_{k + 2} = v$ , the process is in control at time $k$ if

Y_{k} \geq \frac{k}{k + 1} (u - c^{*}),

Y_{k} \geq k (v - u - c^{*}) .

In addition, since the $X_{i}$ are non-negative random variables it follows that

Y_{k} \leq u, for k = 1, \dots, N - 1.

Therefore,

J_{k} (\cdot, u, v) = [a_{k} (u, v), b_{k} (u, v)],

where

a_{k} (u, v) = min {\frac{k}{k + 1} (u - c^{*}), k (v - u - c^{*})}

and

b_{k} (u, v) = u .

Notice that in this case the Fast Fourier Transform technique cannot be used because the limits of the integrals $a_{k} (u, v)$ and $b_{k} (u, v)$ vary over $u$ and $v$ .

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

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

	$N$
$α$	8	Time (sec)	10	Time (sec)	12	Time (sec)
0.10	1.96	2300	2.28	2402	2.55	3980
0.05	2.65	1964	3.08	2339	3.47	2804

Table 4: The control limit

c (α, N)

α = 0.05, 0.10

and

N = 8, 10, 12

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 $t$ -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.

Efficient Computation of the Stochastic Behavior of Partial Sum Processes

Authors

Asymptotics of the Lebesgue constants for bivariate approximation processes

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

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

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

Toward General Analysis of Recursive Probability Models

Minimum energy path calculations with Gaussian process regression

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

1 Introduction

2 Recursive Integration Methodology

2.1 Change of Variables

2.2 General Recursive Formulas

2.3 Recursive Formulas for the Sum of Independent Random Variables

3 Implementation details

3.1 Fast Fourier Transform Convolution

3.2 Accuracy and Efficiency

4 Examples and Illustrations

4.1 Reliability Example

4.2 Product Quality Example

4.3 Discrete Time Stochastic Control Example.

5 Summary

Efficient Computation of the Stochastic Behavior of Partial Sum Processes

Authors

Related Research

Asymptotics of the Lebesgue constants for bivariate approximation processes

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

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

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

Toward General Analysis of Recursive Probability Models

Minimum energy path calculations with Gaussian process regression

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

This week in AI

1 Introduction

2 Recursive Integration Methodology

2.1 Change of Variables

2.2 General Recursive Formulas

2.3 Recursive Formulas for the Sum of Independent Random Variables

3 Implementation details

3.1 Fast Fourier Transform Convolution

3.2 Accuracy and Efficiency

4 Examples and Illustrations

4.1 Reliability Example

4.2 Product Quality Example

4.3 Discrete Time Stochastic Control Example.

5 Summary