Sums, differences, etc., of p-boxes that are positively or negatively dependent

which consist of the supremal convolution t_P,+ and infimal convolution r_P,+ using the independence copula P as the lower bounding copula and addition as the operation.

We have implemented the Williamson and Downs’ formula for adding positively dependent p-boxes, and have checked that the implementation seems to be computing what the formula says. We would like to double check the implementation, and perhaps produce some arguments about the correctness of the software output based on theoretical considerations or numerical examples. For example, it would be useful pedagogically to have numerical examples with various dependencies (such as those developed on the previous two pages for Need 13) that illustrate why the bounds are what they are.

We expected that it should be fairly straightforward to extend the Williamson and Downs algorithm to handle various related cases. In particular, we want to extend the software to handle subtracting p-boxes that are positively dependent, and to adding p-boxes that are negatively dependent. Negative dependence is when P(X £ x, Y £ y) £ P(X £ x) P(Y £ y). These extensions may be more complicated than we at first expected however. As we see below, the results may be surprising.

We must also extend these operations to multiplication and division. We know that the Williamson and Downs algorithm does not apply for multiplication or division when one or both operands straddle zero. But they should be extended to cases in which one or both operands is strictly negative by juggling the negative sign. For instance, suppose A and B are positively dependent, and A is strictly positive in sign but B is strictly negative. The product A´B should be the same as the negation of the negatively dependent convolution of the A and the -B. This may seem a bit complex. Basically, we’re pulling the negative sign out of the B and reapplying it to the final product. But when we pull the negative sign out of the B we create something that is now no longer positively associated with A but rather negatively associated. See? It makes perfect sense, right? So what are the specific rules that govern the extension to multiplication and division of p-boxes having various signs?

The C++ code below in green is supposed to implement Williamson and Downs’ formula for adding p-boxes that are positively dependent. Keep in mind that this C++ code actually does the calculations in terms of the quantile functions (i.e., the quasi-inverse distribution functions) rather than the distribution functions F and G. And we generally have only bounds on these functions because we start with p-boxes rather than precise distributions for the addends.

// POSITIVE

case RandomNbr::plus:

{

for (i=0; i<n; i++)

{

outlier = toobig;

for (j = i; j<n; j++)

{

here = x.d[j] + y.d[n*(i+1)/(j+1)-1];

if (here<outlier) outlier = here;

}

z.d[i] = outlier;

outlier = -toobig;

for (j = 0; j<=i; j++)

{

int kk = n * (((float) i)/((float) n)-((float) j)/((float) n))/(1-((float) j)/((float) n));

here = x.u[j] + y.u[kk];

if (here>outlier) outlier = here;

}

z.u[i] = outlier;

}

z.mymean = x.mean() + y.mean();

z.myvar = env(x.variance()+y.variance(), x.variance()+y.variance()+2*sqrt(x.variance()*y.variance()));

}

For comparison, the code below in blue does addition for the Fréchet case where no assumption is made about the dependence of the operands. Notice that the only difference (other than the moments) between the algorithms for the positive and Fréchet cases is the index of the y-variable on the lines that compute the value here.

// FRECHET

case RandomNbr::plus:

{

for (i=0; i<n; i++)

{

outlier = toobig;

for (j = i; j<n; j++)

{

here = x.d[j] + y.d[i - j + n - 1];

if (here<outlier) outlier = here;

}

z.d[i] = outlier;

outlier = -toobig;

for (j = 0; j<=i; j++)

{

here = x.u[j] + y.u[i - j];

if (here>outlier) outlier = here;

}

z.u[i] = outlier;

}

z.mymean = x.mean() + y.mean();

z.myvar = env(x.variance()+y.variance()-2* sqrt(x.variance()*y.variance()),

x.variance()+y.variance()+2* sqrt(x.variance()*y.variance()));

}

Below is a third implementation of the Williamson and Downs formula for adding positively dependent p-boxes. It is a script for R that will become part of the ESTpbox.r source library. It was derived from the green C++ code above.

source('ESTpbox.r')

positiveconv.pbox <- function(x,y,op='+') {

# if (op=='-') return(negativeconv.pbox(x,negate.pbox(y),'+'))

# if (op=='/') return(negativeconv.pbox(x,reciprocate.pbox(y),'*'))

x <- makepbox(x)

y <- makepbox(y)

n = pbox.steps+1

zu <- zd <- rep(0,pbox.steps)

for (i in 1:pbox.steps)

{

j <- i:pbox.steps

k = trunc(n*i/j-1)

zd[[i]] <- min(do.call(op,list(x@d[j],y@d[k]))) # here = x.d[j] + y.d[n*(i+1)/(j+1)-1];

j <- 1:i

k = trunc(n*(i/n-j/n)/(1-j/n))+1

zu[[i]] <- max(do.call(op,list(x@u[j],y@u[k]))) # int kk = n * (((float) i)/((float) n)-((float) j)/((float) n))/(1-((float) j)/((float) n)); here = a.u[j] + b.u[kk];

}

# mean

ml <- -Inf

mh <- Inf

if (op %in% '+') { ml <- x@ml + y@ml; mh <- x@mh + y@mh }

if (op %in% '-') { ml <- x@ml - y@mh; mh <- x@mh - y@ml }

# variance

vl <- 0

vh <- Inf

if (op %in% c('+')) {

vl <- min(x@vl+y@vl, x@vl+y@vl+2*sqrt(x@vl*y@vl));

vh <- max(x@vh+y@vh, x@vh+y@vh+2*sqrt(x@vh*y@vh));

}

pbox(u=zu, d=zd, ml = ml, mh = mh, vl=vl, vh=vh, dids=paste(x@dids,y@dids))

}

This R routine is perhaps even more cryptic than the C++ code because of the compact use of arrays that R demands. We can check the routine above with the following direct calculation (not in terms of inverses).

x = 1:200 / 10

p = n = NULL

for (z in x) {

y = z - x

F = punif(x,2,3)

G = punif(y,2,3)

p = c(p,max(F * G)) # inverse of d

n = c(n,min(1-(1-F)*(1-G))) # inverse of u

}

positiveconv.pbox(U(2,3),U(2,3))

points(x,n,col='red')

points(x,p)

This check reveals agreement between the R routine positiveconv.pbox and the direct calculation in terms of the distribution functions F and G. And this result agrees with the p-box corresponding to p = 0 (where T_p = P) that is depicted in Figure 21 of Williamson and Downs (1990). Thus, I believe that the R function positiveconv.pbox is probably correctly implemented. So what about extending it to other operations and to negative dependencies? The examples below try a naked, very simplistic extension to multiplication, subtraction and division:

###############################

# Examples

old.par <- par(mfrow=c(2,2),mar=c(2,2,3,1))

pbox.plotting.every = pbox.plotting = FALSE

# addition

u = U(1,2)

c = positiveconv.pbox(u,u,'+')

c = positiveconv.pbox(u,u)

i = u %|+|% u

p = perfectconv.pbox(u,u)

f = frechetconv.pbox(u,u)

plot(f,col='blue',lwd=3)

lines(c,col='seagreen',lwd=3)

lines(i,col='red')

lines(p,col='tan')

# multiplication

u = U(1,2)

c = positiveconv.pbox(u,u,'*')

i = u %|*|% u

p = perfectconv.pbox(u,u,'*')

f = frechetconv.pbox(u,u,'*')

plot(f,col='blue',lwd=3)

lines(c,col='seagreen',lwd=3)

lines(i,col='red')

lines(p,col='tan')

# subtraction

u = U(1,2)

c = positiveconv.pbox(u,u,'-')

i = u %|-|% u

p = perfectconv.pbox(u,u,'-')

f = frechetconv.pbox(u,u,'-')

plot(f,col='blue',lwd=3)

lines(c,col='seagreen',lwd=3)

lines(i,col='red')

lines(p,col='tan')

# division

u = U(1,2)

c = positiveconv.pbox(u,u,'/')

i = u %|/|% u

p = perfectconv.pbox(u,u,'/')

f = frechetconv.pbox(u,u,'/')

plot(f,col='blue',lwd=3)

lines(c,col='seagreen',lwd=3)

lines(i,col='red')

lines(p,col='tan')

par(old.par)

The results of this script are shown in fours graphs below. The graphs depict the results for addition and multiplication (across the top) and subtraction and division (across the bottom). The Fréchet bounds are shown in each graph in blue. The red distribution is the result of the operation under an independence assumption. The tan distribution is the result of assuming perfect dependence between the operands. The green p-box is the result of the positiveconv.pbox routine which assumes the operands are positively dependent.

It looks like addition and multiplication produce perfectly reasonable results. For subtraction, the result of positive dependence is the same as the result of the Fréchet case. This is not actually untenable. Williamson and Downs (1990, page 99, pages 131ff) say that only the lower copula bound contributes to narrowing the resulting p-box, so it may be correct that the bounds on the sum of negatively dependent p-boxes are exactly the same as the Fréchet bounds on the sum, which make no assumption about their dependence. It is hard to say whether or not the software we have so far handles these cases correctly, or, if not, what corrections need to be made. We need to derive the right algorithms to extend the Williamson and Downs formula.

These examples depicted on the next page used precise distributions, i.e., p-boxes whose edges are coincident. We will have to make sure that the codes we produce also work correctly on nondegenerate p-boxes.