A General Method for Construction of Bivariate Stochastic Processes Given two Marginal Processes

I. INTRODUCTION

Suppose two random variables Y, Z are given, with some "real life" interpretation.

Mainly we consider a pair Y, Z to have bio-medical, reliability or econometric interpretation. For example, Y, Z may be (stochastically dependent) levels of some chemicals (hormones, for example) in a human or an animal body. In particular, one may consider Y to describe random level of "bad" cholesterol and Z an accompanying sugar level (both measured at the same time). Apparently, the two (random) quantities are stochastically dependent on each other. Another tandem of random quantities Y, Z can be daily salt consumption Y and blood pressure Z. There are a countless number of such pairs of random quantities that have biomedical meanings. They mostly are stochastically dependent. This means that a measured value of Z has an influence on the probabilities (if, for a fixed value z, the random event occurs and , such that , is its elementary realization). In particular, may be considered a 'critical value' of Y (for example, the critical cholesterol level ) so that occurrence of the random event means occurrence of a disease.

Again, the probability of occurring of the so defined disease may depend on the measured value of the random quantity when, for a given measurement , we observe that , for some fixed value . Here notice that, in this framework, one may disregard an obtained by a measurement particular value and only notice if or otherwise. The more essential value, in this case, is .

In the general case of two random variables Y, Z we will be interested in finding analytical formulas for all the conditional probabilities for any arbitrary value y of Y and not only for a critical one . To achieve this goal we will outline methodology for construction of the joint probability distribution of the random variables Y, Z under various situations.

A similar but different development of methodology of construction of bivariate probability distributions based on Cox and Aalen ideas (see, [2] and [1]) can be found in [3].

The constructions so far rely on building bivariate probability distributions of random vectors. However, the main goal of this work is to extend the obtained methods to methods where the random variables Y, Z are replaced by t-dependent stochastic processes . Thus, given two marginal (originally independent) univariate stochastic processes, we provide tools for constructing a class of bivariate processes such that the "initial" processes and their probability distributions, for each t, remain as the marginals within the constructed joint. The obtained methods of the constructions turn out to be quite general and are governed by some (linear, in particular) integral equations.

As an example of the construction of processes with continuous time we consider construction of some classes of "bivariate Wiener processes" based on two given Wiener univariates. In applications, the underlying univariate Wiener processes may have econometric meanings. For example, one of the two marginal processes may describe the stock market level time evolution while the other the accompanying gross national product rate at time t. In another application, a process describes the level of inflation at a given time and the other processes the level of employment at the same time. In all cases one is interested in the stochastic dependence between the two marginal processes. That need motivates constructions of the 'joint stochastic processes'. In general, the marginal processes, as applied in econometry, are normal. However, among all the mutual stochastic dependencies at any given time epoch t we choose not to be those described by the classical bivariate normal models, since such dependencies are very well known. Our method provides a large variety of stochastic dependencies, some possibly overlooked in literature.

In another example we consider discrete time bivariate processes associated with first Pareto distributions. All the bivariate stochastic processes considerations are included in section 4.

Random Vectors, A General Approach

2.1 Suppose we are given two arbitrary random variables Y, Z with known marginal survival functions and . Our aim is to provide a general methodology for constructions of various joint survival functions such that all their marginal survival functions are invariant with respect to the performed constructions, always remaining the same as the originally given and . In general, we restrict our considerations to the cases where both the probability distributions, given by the survival functions , , possess probability densities so that the corresponding hazard rates , exist too. In this case, according to the common rule, we can express the considered marginal survival functions as follows:

In the simplest case, i.e., when the random variables Y, Z are independent, their joint survival function satisfies:

To impose a dependence structure for the above "initially" independent random variables Y, Z, we multiply the right hand side of formula (2) by a "dependence factor" . This factor we propose to call the "Aalen factor".

2.2 At this point realize that for any joint survival function having the same fixed marginals , there exist unique Aalen factor . However, in most situations we may encounter, we do not know in advance the joint survival function but rather we have the marginals , . We aim to choose a proper Aalen factor in order to construct the corresponding . To give some more intuitive meaning to the Aalen factor we express it in exponential form:

The function we call the "joiner" of the two marginal survival functions , .

Of course, the joiner is uniquely determined by a given Aalen factor since we have . However, the joiner not always is represented by the integral dudt. The 'integral representation' of a joiner is only valid when such a function , considered as a kind of an "additional hazard rate", exists. Existence of together with the existence of the hazard rates , involves existence of the joint probability density of the considered random vector . Using the function we define the conditional hazard rate of Y | Z = u to be equal to: , and the conditional hazard rate of Z | Y = t as:

For the meaning of the differentials and see section 2.4 and the associated "(microshok □ microdamage) □ microchange in hazard rate" paradigm. (These infinitesimal "microchanges" correspond to the above considered differentials.)

Since we are mostly dealing with the survival functions (not with probability densities), we will rather be interested in the conditional hazard rates of, which involve the form:

and in the conditional hazard rates of using the form

The non-negativity assumption for the function is not mandatory but simplifies many considerations. More general and also necessary assumption for the class of functions are given by the inequalities:

for any nonnegative t, u. As we have already defined it, the left hand sides of the above inequalities represent the (conditional) hazard rates, so they must be non-negative.

Also recall that , are always nonnegative.

Restricting the function to nonnegative values implies considering positive stochastic dependences only between the considered random variables Y, Z.

2.3 Assuming existence of the underlying hazard rates, let us analyze the general form of the joint survival function under the investigation:

with .

First realize that upon setting on the right hand side of (4) z = 0, the expression reduces to the marginal . Likewise, upon the setting y = 0 in (4) turns the expression to the marginal .

Second, upon dividing both sides of (4) by one obtains the conditional survival function

Also dividing (4) by yields the conditional survival function

2.4 Let us explain more the structure of the conditional hazard rate \lambda{1}(y \mid Z \geq z), given "initially" the marginal (baseline) hazard rate \lambda{1}(y). The case of \lambda_{2}(z \mid Y \geq y) can be understood in an analogous way. Example 1. Suppose we have two different objects O1, O2 whose "behaviors" are exhibited by the random quantities Y, Z respectively. Take, at first, object O1 as the object "of main interest", characterized by the quantity Y. Now consider the "activity" of object O2, measured by the random quantity Z. The activity Z is regarded as "stress" that the object O1 is subjected to.

(Such an "activity" may be understood as temperature, for example.) In short, the stress Z is an explanatory variable for the variable of the main interest Y. In the physical absence of the object O2 the variable Y is independent on Z, so the corresponding hazard rate for Y is simply equal to , and the (unconditioned) survival function of Y equals

(The following 'micro-shocks □ micro-damages' mechanism that yields the stochastic dependence is described in more detail in [4].)

When we consider the case where the object O2 is physically accompanying ("connected to") object O1, its activity "produces" at each time epoch 'physical' micro-shock on O1. These micro-shocks result in corresponding micro-damages in the physical structure of O1 which in turn result in micro-changes in the hazard rate of the corresponding random quantity Y. Every such a micro-change (as related to a quantity u) at a given moment t relies on adding to the Y's hazard rate the infinitesimal quantity . This change is, actually, too small to be recognized by any physical measurement, but mathematically we can express it as an infinitesimal differential as proportional to the given "intensity" . The microchanges cumulate as u (for Z = u) "runs" through the interval [o,z].

This phenomenon of the cumulation of the micro-changes (as corresponding to an 'elementary' random events Y = t) results in the conditional hazard rate

which, evidently, is different from the "original" (baseline) .

The foregoing integration relates to the transfer from an elementary event Z = u to the random event [Here realize that since we apply, as the analytic tool, the survival functions, instead of the distribution functions, we eventually will 'calculate' the probability of its complement, i.e., the probability of ].

Integrating (7) with respect to \text{\"time\"} t over the interval [0,y] and applying formula (5), we obtain the conditional survival function P(Y\geq y|Z\geq z). The latter, once multiplied by the marginal distribution in the form P(Z\geq z), results again in the joint survival function of the random vector (Y,Z).

Setting we obtain the formula

which corresponds to the Aalen additive model [1] being the modification of the famous Cox proportional hazards model [2].

Remark 1: In a more general case one may consider the second differential of the joiner (,) as a microchange (corresponding to both elementary events ["actions"] t and u [or u at the time epoch t]) in both hazard rates caused by mutual interaction of the "micro-shocks" described by t and u, caused by activities of the objects O1 and O2. Realize that neither t nor u needs to represent time (they may, for example, be levels of hormones in a human body).

In this framework, that interaction may be better understood if the intensity is chosen to have separated variables, i.e., (see examples in below). Notice also that the secondary differential dtdu is added to both the original hazard rates and , which agrees with the idea of mutual interactions.

On Representation of Bivariate Survival Functions

3.1 The Representation

Write formula (4) in a slightly more general form:

We now depart from the Aalen-like model towards a more general paradigm. We even may drop the assumption on the existence of the hazard rates , , and the function need not to be expressed by the integral dudt. Instead of formulas (4) and (4*) for the joint survival function of (Y,Z) we may use the following, more general, formula:

where is assumed to be a real function, defined for , .

We will require that it satisfies the conditions specified as follows. Thus the function is:

From the last condition we directly obtain the following Property.

Property 1: If the latter condition (3) is satisfied, then both the marginal probability distributions of the joint distribution given by formula (8) are preserved in the sense that they are the same as the baseline distributions , originally given.

If the marginal and the joint probability densities of the considered random variables Y, Z exist, the concern is on non-negativity of them. While the marginal densities, say, , are non-negative due to the common assumption, we need some special condition for to assure non-negativity of the corresponding joint density .

For the joint density to exist, the function must have continuous partial derivatives of first order , and the continuous second order mixed partial derivative .

An additional condition that must be satisfied by has the form of the following inequality:

This follows from the form of the joint density (if it exists):

which must be nonnegative. Here, , and , and . A simpler condition than (9) for the non-negativity of , which is sufficient and necessary too, for the existence of the joint survival function is the following, obtained from (9):

The condition (9) or together with the conditions (1), (2), and (3) are sufficient and necessary for "connecting" the two survival functions , by a given joiner into the bivariate model .

As it was pointed out above, there may be a numerous of such models.

3.2 Particular Cases of the Bivariate Models

A. Obviously, when reduces to zero for all y, z, then model (8) describes independent random variables.

B. If the baseline hazard rates exist and are constant, we call this model "exponential".

In this case we may choose constant to obtain the following special case of model :

where .

This, apparently, represents the first bivariate exponential Gumbel probability distribution (see [5]). C. One also obtains the following 'Weibullian version' of the above bivariate Gumbel:

where and are positive reals.

The latter two models are special cases of representation .

We only need to check when condition is satisfied, i.e.,

To check for that, we properly differentiate the terms of the expression ( over and over and over and ), and set inequality (11*) in the form:

It holds, for every nonnegative y and z, whenever

Thus, if the latter condition is satisfied then model (11) is properly defined.

Remark 2: As for the above given Weibullian version (11) of the Gumbel exponential model (with ), this can be 'partially' generalized by the following model:

Now inequality (12) takes the form:

The necessary condition for (14) to hold for all nonnegative values y, z is that

together with

However, condition (15) makes inequality (14) not true for . The remedy for this is to shift the variables y, z by imposing the additional conditions:

1 = c = y, 1 = c = z for some real c.

Now model (13) is well defined with .

Remark 3: The form of model (8) and Property 1 allows us to construct numerous of bivariate probability distributions. Namely, realize that:

1. For any given fixed 'joiner function' (such that inequality (9) or (9*) is satisfied), may be "connected" by formula (8) to "become" the bivariate survival (distribution) function in which they remain the marginals.

2. For any fixed pair of probability distributions , , there is a wide class of possible joiners [determined by the conditions (1), (2), (3) from section 3.1 together with inequality (9) or ] such that the two distributions can be "connected" into the bivariate model (8) in as many ways as there are possible proper functions .

Finally notice that the above methods for "connecting" any two probability distributions into bivariate distributions, resembles the idea of the copula methodology [6]. It is, however, different. For more remarks on that see [3].

II. ON CONSTRUCTION OF BIVARIATE STOCHASTIC PROCESSES GIVEN ANY TWO UNIVARIATE MARGINAL PROCESSES WHICH SHARE THE SAME TIME

Suppose the (marginal) stochastic processes , are completely defined. So, for every time epoch (T is a nonempty set) two random variables , have known in advance survival functions ,

If the corresponding hazard rates , exist, then, according to what was pointed out in previous sections, for any joiner function [satisfying, for each t, inequality (9) with respect to the , as well as conditions (1)–(3) from section 3.1] there exist a unique joint survival function which is also a function of time t, given by:

Our intention is to consider, for each , the function as the joint survival function (or joint distribution) of the bivariate stochastic process . At this point recall that the marginals , (z) [as functions of time t] are assumed to define completely the processes , respectively.

Suppose time t is continuous. If both survival functions and are continuous functions of the time then we will postulate the joiner to be continuous, as a function of t, as well. Now consider two examples of bivariate stochastic processes whose constructions are based on known initial univariate (marginal) processes.

Example A: Consider two Wiener stochastic processes , with the, given in advance, for each t, survival functions

where , and b, c are the real parameters of the considered Wiener processes.

For all time epochs s, t such that o < s < t, all the distributions of the random vectors and are described by the usual bivariate formulas associated with Wiener processes.

We seek, for every time epoch t, a joiner for the two survival functions (17), which analytically is nice and simple enough, also as a function of t. It needs to satisfy conditions (1) – (3) [section 3.1] as well as inequality (9) or for each t. If found, then, for each t > 0, the joint survival function of the stochastic process would have the form (8) which, in this case, is:

Of course, in the vast majority of cases, (18) is different from classical bivariate normal distributions even if the marginals are normal.

Realize that conditions (1) - (3) from section 3.1 are usually easily checked to be satisfied, so that we really need to check condition (9).

Recall, in the case of stochastic processes, this amounts to checking if:

This inequality needs to be satisfied for each t.

At this point notice that:

These quantities, upon the additional assumption for all nonnegative u, v, are nonnegative and therefore the inequality

always holds.

The latter together with (19), provides solutions of (19) in the form:

where the coefficient function satisfy , and, in particular, may be considered a constant in t, especially one may set for each t.

The model (or rather a candidate for a model in a practical application) we have obtained, has the form of the time dependent joint survival function of the random vectors for . It has the following form:

For the above considered 'bivariate Wiener stochastic process' the marginals and present in the last formula, are given by (17) and (17*).

Also, in the above case, , are the corresponding hazard rates associated with the considered normal distributions. It is clear, however, that the obtained class of models is much wider than the class of the bivariate Wiener. Write inequality (19) in the form:

Another candidate for a set of the models (i.e., set of "bivariate Wiener stochastic processes") will be given by the set of the functions that are solutions of the following integral equation directly derived from inequality :

Where , for .

Realize that the right-hand side of (20) is always less than or equal than the right hand side of . For equation (20) to hold the condition is, in general, essential.

The integral equation (20) is nonlinear. However, it can be simplified to the following:

where is any continuous function.

Equation is linear (nonhomogeneous), but the (in general, known) coefficients are still variable.

The set of solutions of is contained in the set of solutions of equation , so that any solution of is a solution of .

Equation can be simplified more. This will yield a smaller set of solutions being still solutions of . Namely, setting , we obtain the (purely) linear integral equation

So the set of solutions of forms a vector space over the field of real numbers. In the narrowest version of the problem the coefficients can be made constants. Especially the assumptions and comprise the exponential case. But to stick with the "Wiener model" one must keep the hazard rates (coefficients) , in and in as corresponding to the distributions given by (17) and .

Now, having any solution of any integral equation above, we find the joint survival function of the corresponding bivariate stochastic process as expressed by formula (18), where dsdr. The so obtained joiner (or in the case of random vectors only) would be a candidate for the bivariate stochastic model.

Next steps are statistical verifications of an eventual fit of the obtained models to a given data that might yield the eventual choice of the best solution among all the obtained. This subject is, however, out of scope of this paper.

Notice that the above presented method for construction of bivariate Wiener processes is applicable to any two (marginal) stochastic processes , such that each random variable and possesses a hazard rate. Thus, one obtains more models. The method also includes discrete time cases like the case described below.

Example B. In this example we seek an additional method for constructions. From now on we slightly change the notation by replacing the symbols (Y, Z) for the random vectors by the symbols (X,Y). Consider the following first Pareto survival functions for two random variables X, Y:

In this (Pareto) case the variables X, Y usually (but not always) describe income or wealth redistributed within a society. We consider the case where the same society is subdivided in some manner into two groups that differ by professions or some other social indicator (ethnicity, race, gender, age etc...). In the models the differences between the two groups income redistributions are reflected by the parameters , (minimal incomes), and by the shape parameters , . It is reasonable to assume that the probability distributions of the incomes together with the incomes themselves will evolve over time. That is why we consider stochastic (Pareto) processes , description, where the discrete time t is defined as multiplicities t = r, 2r, 3r,... of some time period r such as a year, a quarter, or a month. In our notation we adopt r = 1, and therefore t = 1, 2, 3... We now assume that, for each t, the random variables , are distributed according to the rules:

For simplicity we will assume that , and .

The above univariate "Pareto stochastic processes" are determined by the constants , , and the way we define the functions , for t = 1, 2, 3,.... Our (simplest) proposition is to define

where (as a particular example) , , r = 0.10, s = 0.15.

Now, for each t, we find the joint survival functions of the random vector . Recall the general form of the joint survival functions:

, (here, and the meaning of is different than that of ).

Since we already have both and as given above, we need to find a class of proper joiners to make, for each t, the joint survival function of properly defined. Set .

Recall that we have:

where, for every t, the function must satisfy conditions (1) - (3) from section 3.1. Since the domain of this function is reduced to , with positive and , condition (3) (section 3.1) now takes the form:

Also, for each t, we have and .

Therefore, as pointed out in Property 1 (section 3.1), it follows that, for each t, the functions , are the marginal survival functions for the joint . Moreover, we restrict ourselves to positive stochastic dependences only, so we have, for every t, . [As already mentioned, the last assumption is not necessary (and not necessarily true in this case), but makes the initial considerations clearer.] In order to determine a proper class of the (not necessarily all) joiners [which are equivalent to the corresponding functions ], we need to examine (for each t) inequality ( ) which is

where, in this 'Pareto case', we have and .

Thus inequality (21) takes the form

Like in Example A, the first easiest solution to our problem is

where, for every , .

Recall at this point that like the variables x, y also the variables u, v satisfy , .

The natural assumption in this Pareto case is that both , (i.e., the minimal incomes in both social groups are at least 1 "unit"). Taking the above under the consideration we may choose for (the determinant of the stochastic dependence between and ) the following:

This satisfies (21).

Finally, the full version of the so derived joint survival function of the Pareto random vectors , for every t, is:

Annotation: Just in recent days, when this paper was already finalized, we learned that similar results as for the bivariate distributions (and not bivariate stochastic processes) were considered by Finkelstein [7] in 2003. The author also considered the "Aalen factor" (see formulas (4) and (5) in [7]) but under different names and with no reference to Aalen or Cox. Our impression, however, is that the generality of Finkelstein's approach is somewhat limited. Also, in the references he provides, the results are rather specific and the generality or even universality of this bivariate distributions' representation was probably overlooked. But, anyway, the novum of our results became somehow (but not drastically) limited.