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

Jerzy K. Filus¹, Lidia Z. Filus²

¹Mathematics and Computer Science, Oakton College, Des Plaines, USA.
²Mathematics, Northeastern Illinois University, Chicago, USA.

DOI: 10.4236/jamp.2025.134066                

Abstract

This paper proposes a universal framework for constructing bivariate stochastic processes, going beyond the limitations of copulas and offering a potentially simpler alternative. The achieved generality of the construction methods extends its applicability to diverse stochastic processes also including discrete as well as continuous time cases. The initially given two arbitrary univariate stochastic processes {Y_t}, {Z_t}, are only assumed to share the same time t. When considered as describing (time dependent) random quantities that are physically separated (the baseline case), the processes are independent. From this trivial case we move to the case when physical interactions between the quantities make them stochastically dependent random variables at any moment t. For each time epoch t, we impose stochastic dependence on two “initially independent” random variables Y_t, Z_t by multiplying the product of their survival functions by a proper “dependence factor” ${\varphi }_t\left(y_t,z_t\right)$ , obtaining in this way a universal (“canonical”) form valid for any (!) bivariate distribution. In some known cases, however, this form may become complicated thou it always exists and is unique. The dependence factor, basically, but not always, has the form ${\Phi }_t\left(y,z\right)=\exp \left[-{\displaystyle {\int }_0^y{\displaystyle {\int }_0^z{\Psi }_t}}\left(s;u\right)\text{d}s\text{d}u\right]$ whenever such a continuous function ${\psi }_t\left(s;u\right)$ exists, for each t. That representation of stochastic dependence by the functions ${\psi }_t\left(s;u\right)$ leads, in turn, to the phenomenon of change of the original (baseline) hazard rates of the marginals, similar to those analyzed by Cox and, especially Aalen for single pairs (or sets) of, time independent, random variables. That is why, until Section 4, we consider only single random vectors (Y, Z)' joint survival functions, mostly as a preparation to the theory of bivariate stochastic processes {(Y_t, Z_t)} constructions as initiated in Section 4. The bivariate constructions are illustrated by examples of some applications in biomedical and econometric areas.

Keywords

Bivariate Survival Functions, Bivariate Stochastic Processes’ Constructions, Dependence Functions, Biomedical Applications, Econometrics, Bivariate Wiener and Pareto Stochastic Processes Construction

1. 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 interpretations. 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 z' of Z has an influence on the probabilities P (Y ≥ y₀) (if, for a fixed value z, the random event Z ≥ z occurs and z', such that z' ≥ z, is its elementary realization). In particular, y₀ may be considered a “critical value” of Y (for example, the critical cholesterol level y₀) so that occurrence of the random event Y ≥ y₀ means occurrence of a disease.

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

In the general case of two random variables Y, Z we will be interested in finding analytical formulas for all the conditional probabilities $P\left(Y\ge y|Z\ge z\right)$ for any arbitrary value y of Y and not only for a critical one y₀. 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, [1] and [2]) can be found in [3]. See, also [4].

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 {Y_t}, {Z_t}. Thus, given two marginal (originally independent) univariate stochastic processes, we provide tools for constructing a class of bivariate processes {(Y_t, Z_t)} such that the “initial” processes {Y_t}, {Z_t} 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 process 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 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 already very well known. Our method provides a large variety of stochastic dependencies, some possibly overlooked in literature. Section 2 provides a general approach to the construction of bivariate survival functions. In that we define the method of construction of the conditional survival functions of one random variable given information on the remaining variable. Also, we define there the conditional hazard rates. The section we finalize with analysis of possible reliability application mostly through the described there “micro-shocks → micro-damages” mechanism. In section 3 we continue the problem of bivariate survival functions construction in a more general way, first of all specifying the sufficient and necessary conditions under which various constructions are possible. Examples of some classes of the models are provided. The very topic of construction of bivariate stochastic processes, given two marginal processes, fulfills section 4. After formulating general assumptions for the coming constructions we consider construction of variety of Wiener stochastic processes as Example A. That method of the construction may also comprise other bivariate stochastic processes with continuous time. Next, in Example B, we constructed class of bivariate Pareto stochastic processes with econometric interpretations that may be considered as examples of bivariate processes constructions with discrete time. The paper is finalized with the conclusions.

2. Random Vectors, A General Approach

1) Suppose we are given two arbitrary random variables Y, Z with known marginal survival functions $S_1\left(y\right)=P\left(Y\ge y\right)$ and $S_2\left(z\right)=P\left(Z\ge z\right)$ . Our aim is to provide a general methodology for constructions of various joint survival functions $S\left(y,z\right)=P\left(Y\ge y,Z\ge z\right)$ such that all their marginal survival functions are invariant with respect to the performed constructions, always remaining the same as the originally given S₁(y) and S₂(z). In general, we restrict our considerations to the cases where both the probability distributions, given by the survival functions S₁(y), S₂(z), possess probability densities so that the corresponding hazard rates ${\lambda }_1\left(y\right)$ , ${\lambda }_2\left(z\right)$ exist too. In this case, according to the common rule, we can express the considered marginal survival functions as follows:

$S_1\left(y\right)=\exp \left[-{\displaystyle {\int }_0^y{\lambda }_1}\left(t\right)\text{d}t\right],S_2\left(z\right)=\exp \left[-{\displaystyle {\int }_0^z{\lambda }_2}\left(u\right)\text{d}u\right]$ . (1)

In the simplest case, i.e., when the random variables Y, Z are independent, their joint survival function $S^{*}\left(y,z\right)$ satisfies:

$S^{*}\left(y,z\right)=S_1\left(y\right)S_2\left(z\right)=\exp \left[-{\displaystyle {\int }_0^y{\lambda }_1}\left(t\right)\text{d}t-{\displaystyle {\int }_0^z{\lambda }_2}\left(u\right)\text{d}u\right]$ . (2)

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” $\varphi \left(y,z\right)$ . This factor we propose to call the “Aalen factor”.

2) At this point realize that for any joint survival function $S\left(y,z\right)$ having the same fixed marginals S₁(y), S₂(z) there exist unique Aalen factor $\varphi \left(y,z\right)=S\left(y,z\right)/S_1\left(y\right)S_2\left(z\right)$ . However, in most situations we may encounter, we do not know in advance the joint survival function $S\left(y,z\right)$ but rather we have the marginals S₁(y), S₂(z). We aim to choose a proper Aalen factor $\varphi \left(y,z\right)$ in order to construct the corresponding $S\left(y,z\right)$ . To give some more intuitive meaning to the Aalen factor we express it in exponential form:

$\varphi \left(y,z\right)=\exp \left[-\alpha \left(y,z\right)\right]=\exp \left[-{\displaystyle {\int }_0^y{\displaystyle {\int }_0^z\psi }}\left(t;u\right)\text{d}u\text{d}t\right]$ .

The function $\alpha \left(y,z\right)$ we call the “joiner” of the two marginal survival functions S₁(y), S₂(z).

Of course, the joiner $\alpha \left(y,z\right)$ is uniquely determined by a given Aalen factor $\varphi \left(y,z\right)$ since we have $\alpha \left(y,z\right)=-\log \varphi \left(y,z\right)$ . However, the joiner $\alpha \left(y,z\right)$ not always is represented by the integral

${\displaystyle {\int }_0^y{\displaystyle {\int }_0^z\psi }}\left(t;u\right)\text{d}u\text{d}t$ .

The “integral representation” of a joiner is only valid when such a function $\psi \left(t;u\right)$ , considered as a kind of an “additional hazard rate”, exists. Existence of $\psi \left(t;u\right)$ together with the existence of the hazard rates ${\lambda }_1\left(t\right)$ , ${\lambda }_2\left(u\right)$ involves existence of the joint probability density of the considered random vector (Y, Z). Using the function $\psi \left(t;u\right)$ we define the conditional hazard rate of Y | Z = u to be equal to:

${\lambda }_1\left(t|Z=u\right)={\lambda }_1\left(t\right)+\psi \left(t;u\right)\text{d}u$ ,

and the conditional hazard rate of Z | Y = t as:

${\lambda }_2\left(u|Y=t\right)={\lambda }_2\left(u\right)+\psi \left(t;u\right)\text{d}t$ .

For the meaning of the differentials $\psi \left(t;u\right)\text{d}u$ and $\psi \left(t;u\right)\text{d}t$ 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 Y|Z ≥ z which involve the form:

${\lambda }_1\left(Y=t|Z\ge z\right)={\lambda }_1\left(t\right)+{\displaystyle {\int }_0^z\psi }\left(t;u\right)\text{d}u$ , (3)

and in the conditional hazard rates of Z |Y ≥ y using the form

${\lambda }_2\left(Z=u|Y\ge y\right)={\lambda }_2\left(u\right)+{\displaystyle {\int }_0^y\psi }\left(t;u\right)\text{d}t$ . (3*)

The non-negativity assumption for the function $\psi \left(t;u\right)$ is not mandatory but simplifies many considerations. More general and also necessary assumption for the class of functions $\psi \left(t;u\right)$ are given by the inequalities:

${\lambda }_1\left(t\right)+{\displaystyle {\int }_0^z\psi }\left(t;u\right)\text{d}u\ge 0$ , and

${\lambda }_2\left(u\right)+{\displaystyle {\int }_0^y\psi }\left(t;u\right)\text{d}t\ge 0$ ,

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 ${\lambda }_1\left(t\right)$ , ${\lambda }_2\left(u\right)$ are always nonnegative.

Restricting the function $\psi \left(t;u\right)$ to nonnegative values implies considering positive stochastic dependences only between the considered random variables Y, Z.

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

$\begin{array}{c}S\left(y,z\right)=P\left(Y\ge y,Z\ge z\right)\\ =\exp \left[-{\displaystyle {\int }_0^y{\lambda }_1}\left(t\right)\text{d}t-{\displaystyle {\int }_0^y{\displaystyle {\int }_0^z\psi }}\left(t;u\right)\text{d}u\text{d}t-{\displaystyle {\int }_0^z{\lambda }_2}\left(u\right)\text{d}u\right],\end{array}$ (4)

with $\psi \left(t;u\right)\ge 0$ .

First realize that upon setting on the right hand side of (4) z = 0, the expression reduces to the marginal $\exp \left[-{\displaystyle {\int }_0^y{\lambda }_1}\left(t\right)\text{d}t\right]=S_1\left(y\right)$ .

Likewise, upon the setting y = 0 in (4) turns the expression to the marginal

$\exp \left[-{\displaystyle {\int }_0^z{\lambda }_2}\left(u\right)\text{d}u\right]=S_2\left(z\right)$ .

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

$S_1\left(y|z\right)=P\left(Y\ge y|Z\ge z\right)=\exp \left[-{\displaystyle {\int }_0^y{\lambda }_1}\left(t\right)\text{d}t-{\displaystyle {\int }_0^y{\displaystyle {\int }_0^z\psi }}\left(t;u\right)\text{d}u\text{d}t\right]$ . (5)

Also dividing (4) by S₁(y) yields the conditional survival function

$S_2\left(z|y\right)=P\left(Z\ge z|Y\ge y\right)-\exp \left[-{\displaystyle {\int }_0^z{\lambda }_2}\left(u\right)\text{d}u-{\displaystyle {\int }_0^y{\displaystyle {\int }_0^z\psi }}\left(t;u\right)\text{d}u\text{d}t\right]$ . (6)

Also see the similarity of the general model (4) with the Gumbel bivariate [5] which can be seen as an important specific case of (4).

4) Let us explain more the structure of the conditional hazard rate ${\lambda }_1\left(y|Z\ge z\right)$ , given “initially” the marginal (baseline) hazard rate ${\lambda }_1\left(y\right)$ . The case of ${\lambda }_2\left(z|Y\ge y\right)$ 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 ${\lambda }_1\left(t\right)$ , and the (unconditioned) survival function of Y equals

$S_1\left(y\right)=P\left(Y\ge y\right)=\exp \left[-{\displaystyle {\int }_0^y{\lambda }_1}\left(t\right)\text{d}t\right]$ .

(The following “micro-shocks → micro-damages” mechanism that yields the stochastic dependence is described in more detail in [6].)

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 $\psi \left(t;u\right)\text{d}u$ . This change is, actually, too small to be recognized by any physical measurement, but mathematically we can express it as an infinitesimal differential $\psi \left(t;u\right)\text{d}u$ as proportional to the given “intensity” $\psi \left(t;u\right)$ . The microchanges cumulate as u (for Z = u) “runs” through the interval [0, 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

${\lambda }_1\left(t|Z\ge z\right)={\lambda }_1\left(t\right)+{\displaystyle {\int }_0^z\psi }\left(t;u\right)\text{d}u$ , (7)

which, evidently, is different from the “original” (baseline) ${\lambda }_1\left(t\right)$ .

The foregoing integration relates to the transfer from an elementary event Z = u to the random event 0 ≤ Z < z [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 Z ≥ z].

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

Setting ${\displaystyle {\int }_0^z\psi }\left(t;u\right)\text{d}u=\varphi \left(t,z\right)$ we obtain the formula:

${\lambda }_1\left(t|Z\ge z\right)={\lambda }_1\left(t\right)+\varphi \left(t,z\right)$ (7*)

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

Remark 1. In a more general case one may consider the second differential $\psi \left(t;u\right)\text{d}t\text{d}u$ 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 $\psi \left(t;u\right)$ is chosen to have separated variables, i.e., $\psi \left(t;u\right)={\psi }_1\left(t\right){\psi }_2\left(u\right)$ (see examples in below). Notice also that the secondary differential $\psi \left(t;u\right)\text{d}t\text{d}u$ is added to both the original hazard rates ${\lambda }_1\left(t\right)$ and ${\lambda }_2\left(u\right)$ , which agrees with the idea of mutual interactions.

3. On Representation of the Class of Bivariate Survival Functions

3.1. The Representation

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

$S\left(y,z\right)=P\left(Y\ge y,Z\ge z\right)=\exp \left[-{\displaystyle {\int }_0^y{\lambda }_1}\left(t\right)\text{d}t-\alpha \left(y,z\right)-{\displaystyle {\int }_0^z{\lambda }_2}\left(u\right)\text{d}u\right]$ . (4*)

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 ${\lambda }_1\left(y\right)$ , ${\lambda }_2\left(z\right)$ , and the function $\alpha \left(y,z\right)$ need not to be expressed by the integral ${\displaystyle {\int }_0^y{\displaystyle {\int }_0^z\psi }}\left(t;u\right)\text{d}u\text{d}t$ . Instead of formulas (4) and (4*) for the joint survival function of (Y,Z) we may use the following, more general, formula:

$S\left(y,z\right)=S_1\left(y\right)S_2\left(z\right)\exp \left[-\alpha \left(y,z\right)\right]$ , (8)

where $\alpha \left(y,z\right)$ is assumed to be a real function, defined for y ≥ 0, z ≥ 0.

We will require that it satisfies the conditions specified as follows. Thus the function $\alpha \left(y,z\right)$ is:

1) continuous, for each z, with respect to y and for each y with respect to z,

2) nondecreasing in y and z,

3) $\alpha \left(0,z\right)=\alpha \left(y,0\right)=0$ .

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 S₁(y), S₂(z) 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, f(y), g(z) are non-negative due to the common assumption, we need some special condition for $\alpha \left(y,z\right)$ to assure non-negativity of the corresponding joint density $h\left(y,z\right)$ .

For the joint density $h\left(y,z\right)$ to exist, the function $\alpha \left(y,z\right)$ must have continuous partial derivatives of first order ${\alpha }_{y,z}\left(y,z\right)$ , ${\alpha }_z\left(y,z\right)$ and the continuous second order mixed partial derivative ${\alpha }_{y,z}\left(y,z\right)=\psi \left(y,z\right)$ . An additional condition that must be satisfied by $\alpha \left(y,z\right)$ has the form of the following inequality:

$\left[{\lambda }_1\left(y\right)+{\alpha }_y\left(y,z\right)\right]\cdot \left[{\lambda }_2\left(z\right)+{\alpha }_z\left(y,z\right)\right]\ge {\alpha }_{yz}\left(y,z\right)=\psi \left(y,z\right)$ . (9)

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

$\begin{array}{c}h\left(y,z\right)={\partial }^2/\partial y\partial zS\left(y,z\right)\\ = \left\{\left[{\lambda }_1\left(y\right)+{\alpha }_y\left(y,z\right)\right]\cdot \left[{\lambda }_2\left(z\right)+{\alpha }_z\left(y,z\right)\right]-{\alpha }_{yz}\left(y,z\right)\right\}\\ \times \exp \left[-{\Lambda }_1\left(y\right)-{\Lambda }_2\left(z\right)-\alpha \left(y,z\right)\right]\end{array}$

which must be nonnegative.

Here, $\text{d}/\text{d}y{\Lambda }_1\left(y\right)={\lambda }_1\left(y\right)$ , and $\text{d}/\text{d}z{\Lambda }_2\left(z\right)={\lambda }_2\left(z\right)$ , and ${\alpha }_{y,z}\left(y,z\right)=\psi \left(y,z\right)$ .

A simpler condition than (9) for the non-negativity of $h\left(y,z\right)$ , which is sufficient and necessary too, for the existence of the joint survival function (4*) is the following, obtained from (9):

${\lambda }_1\left(y\right){\lambda }_2\left(z\right)\ge {\alpha }_{yz}\left(y,z\right)=\psi \left(y,z\right)$ . (9*)

The condition (9) or (9*) together with the conditions (1), (2), and (3) are sufficient and necessary for “connecting” the two survival functions S₁(y), S₂(z) by a given joiner $\alpha \left(y,z\right)$ into the bivariate model $S\left(y,z\right)$ .

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 $\alpha \left(y,z\right)$ 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 $\psi \left(y,z\right)=a=\text{constant}$ to obtain the following special case of model (4*):

$S\left(y,z\right)=\exp \left[-{\lambda }_{1}y-ayz-{\lambda }_{2}z\right]$ , (10)

where $0\le a\le {\lambda }_1{\lambda }_2$ .

This, apparently, represents the first bivariate exponential Gumbel probability distribution

C) One also obtains the following “Weibullian version” of the above bivariate Gumbel:

$S\left(y,z\right)=\exp \left[-{\lambda }_1{y}^{{\gamma }_1}-a{y}^{{\gamma }_1}{z}^{{\gamma }_2}-{\lambda }_2{z}^{{\gamma }_2}\right]$ , (11)

where ${\gamma }_1$ and ${\gamma }_2$ are positive reals.

The latter two models are special cases of representation (4*).

We only need to check when condition (9*) is satisfied, i.e.,

${\lambda }_1\left(y\right){\lambda }_2\left(z\right)\ge \psi \left(y,z\right)$ . (11*)

To check for that, we properly differentiate the terms of the expression $-{\lambda }_1{y}^{{\gamma }_1}-a{y}^{{\gamma }_1}{z}^{{\gamma }_2}-{\lambda }_2{z}^{{\gamma }_2}$ (${\lambda }_1{y}^{{\gamma }_1}$ over y and ${\lambda }_2{z}^{{\gamma }_2}$ over z and $a{y}^{{\gamma }_1}{z}^{{\gamma }_2}$ over y and z), and set inequality (11*) in the form:

${\lambda }_1{\lambda }_2{\gamma }_1{\gamma }_2{y}^{{\gamma }_1-1}{z}^{{\gamma }_2-1}\ge a{\gamma }_1{\gamma }_2{y}^{{\gamma }_1-1}{z}^{{\gamma }_2-1}$ . (12)

It holds, for every nonnegative y and z, whenever $a\le {\lambda }_1{\lambda }_2$ . 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 $a\le {\lambda }_1{\lambda }_2$ ), this can be ‘partially’ generalized by the following model:

$S\left(y,z\right)=\exp \left[-{\lambda }_1{y}^{{\gamma }_1}-a{y}^{{\delta }_1}{z}^{{\delta }_2}-{\lambda }_2{z}^{{\gamma }_2}\right]$ . (13)

Now inequality (12) takes the form:

${\lambda }_1{\lambda }_2{\gamma }_1{\gamma }_2{y}^{{\gamma }_1-1}{z}^{{\gamma }_2-1}\ge a{\delta }_1{\delta }_2{y}^{{\delta }_1-1}{z}^{{\delta }_2-1}$ . (14)

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

$1\le {\delta }_i\le {\gamma }_i$ , for $i=1,2$ , (15)

together with

$a\le {\lambda }_1{\lambda }_2$ . (16)

However, condition (15) makes inequality (14) not true for 0 ≤ y, z < 1. The remedy for this is to shift the variables y, z by imposing the additional conditions: $1\le c\le y$ , $1\le c\le z$ for some real $c$ .

Now model (13) is well defined with $1\le {\delta }_i\le {\gamma }_i$ . 

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” $\alpha \left(y,z\right)$ , ‘any’ pair of two probability distributions [not necessarily both from the same class of the distributions], determined by S₁(y), S₂(z) (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 S₁(y), S₂(z), there is a wide class of possible joiners $\alpha \left(y,z\right)$ [determined by the conditions (1), (2), (3) from section 3.1 together with inequality (9) or (9*)] such that the two distributions can be “connected” into the bivariate model (8) in as many ways as there are possible proper functions $\alpha \left(y,z\right)$ .

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

4. On Construction of Bivariate Stochastic Processes Given Any Two Univariate Marginal Processes Which Share the Same Time

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

$S_t\left(y\right)=P\left(Y_t\ge y\right)$ ,

$R_t\left(z\right)=P\left(Z_t\ge z\right)$ .

If the corresponding hazard rates ${\lambda }_1\left(y,t\right)$ , ${\lambda }_2\left(z,t\right)$ exist, then, according to what was pointed out in previous sections, for any joiner function ${\alpha }_t\left(y,z\right)=\alpha \left(y,z,t\right)$ [satisfying, for each t, inequality (9) with respect to the ${\lambda }_1\left(y,t\right)$ , ${\lambda }_2\left(z,t\right)$ 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:

$S_t\left(y,z\right)=P\left(Y_t\ge y,Z_t\ge z\right)=S_t\left(y\right)R_t\left(z\right)\exp \left[-\alpha \left(y,z,t\right)\right]$ .

Our intention is to consider, for each $t\in T$ , the function S_t(y, z) as the joint survival function (or joint distribution) of the bivariate stochastic process ${\left\{\left(Y_t,Z_t\right)\right\}}_{t\in T}$ . At this point recall that the marginals S_t(y), R_t(z) [as functions of time t] are assumed to define completely the processes {Y_t}, {Z_t} respectively.

Suppose time t is continuous. If both survival functions S_t(y) and R_t(z) are continuous functions of the time then we will postulate the joiner $\alpha \left(y,z,t\right)$ 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 {Y_t}, {Z_t} with the, given in advance, for each t, survival functions

$P\left(Y_t\ge y\right)=1-\varphi \left(\left(y-bt\right)/{\sigma }_1\sqrt{t}\right)$ and (17)

$P\left(Z_t\ge z\right)=1-\varphi \left(\left(z-ct\right)/{\sigma }_2\sqrt{t}\right)$ , (17*)

where

$\varphi \left(x\right)=\left(1/\sqrt{2\pi }\right){\displaystyle {\int }_{-\infty }^x\exp \left[-w^2/2\right]\text{d}w}$ ,

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

For all time epochs s, t such that 0 < s < t, all the distributions of the random vectors (Y_s, Y_t) and (Z_s, Z_t) are described by the usual bivariate formulas associated with Wiener processes.

We seek, for every time epoch t, a joiner $\alpha \left(y,z,t\right)$ for the two survival functions (17), (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 (9*) for each t. If found, then, for each t > 0, the joint survival function of the stochastic process {(Y_t, Z_t)} would have the form (8) which, in this case, is:

$S_t\left(y,z\right)=S\left(y,z,t\right)=P\left(Y_t\ge y\right)P\left(Z_t\ge z\right)\exp \left[-\alpha \left(y,z,t\right)\right]$ . (18)

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 whether:

$\begin{array}{l}{\partial }^2/\partial y\partial z\alpha \left(y,z,t\right)={\alpha }_{yz}\left(y,z,t\right)=\Psi \left(y,z,t\right)\\ \le \left[{\lambda }_1\left(y,t\right)+{\alpha }_y\left(y,z,t\right)\right]\times \left[{\lambda }_2\left(z,t\right)+{\alpha }_z\left(y,z,t\right)\right]\end{array}$ (19)

For safety reasons we additionally propose adopting the assumption: $\Psi \left(y,z,t\right)\le \Psi \left(y,z\right)$ , where the function $\Psi \left(y,z\right)$ is constant over the time.

The inequality (19), however, needs to be satisfied for each t.

At this point notice that:

${\alpha }_y\left(y,z,t\right)={\displaystyle {\int }_{-\infty }^z\Psi }\left(y,v,t\right)\text{d}v$ and

${\alpha }_z\left(y,z,t\right)={\displaystyle {\int }_{-\infty }^y\Psi }\left(u,z,y\right)\text{d}u$ .

These quantities, upon the additional assumption $\Psi \left(u,v,t\right)\ge 0$ for all nonnegative $u,v$ , are nonnegative and therefore the inequality

${\lambda }_1\left(y,t\right){\lambda }_2\left(z,t\right)\le \left[{\lambda }_1\left(y,t\right)+{\alpha }_y\left(y,z,t\right)\right]\times \left[{\lambda }_2\left(z,t\right)+{\alpha }_z\left(y,z,t\right)\right]$

always holds.

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

$\Psi \left(y,z,t\right)=\omega \left(t\right){\lambda }_1\left(y,t\right){\lambda }_2\left(z,t\right)$ ,

where the coefficient function ω(t) satisfy 0 ≤ ω(t) ≤ 1, and, in particular, ω(t) may be considered a constant in t, especially one may set ω(t) = 1 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 (Y_t,Z_t) for t ≥ 0. It has the following form:

$\begin{array}{c}{S}_t\left(y,z\right)=P\left(Y_t\ge y\right)P\left(Z_t\ge z\right)\exp \left[-\omega \left(t\right){\displaystyle {\int }_{-\infty }^y{\displaystyle {\int }_{-\infty }^z{\lambda }_1}}\left(u,t\right){\lambda }_2\left(v,t\right)\text{d}u\text{d}v\right]\\ =\exp \left[-{\displaystyle {\int }_{-\infty }^y{\lambda }_1}\left(u,t\right)\text{d}u-\omega \left(t\right){\displaystyle {\int }_{-\infty }^y{\displaystyle {\int }_{-\infty }^z{\lambda }_1}}\left(u,t\right){\lambda }_2\left(v,t\right)\text{d}u\text{d}v-{\displaystyle {\int }_{-\infty }^z{\lambda }_2}\left(v,t\right)\text{d}v\right]\end{array}$

For the above considered ‘bivariate Wiener stochastic process’ the marginals P(Y_t ≥ y) and P(Z_t ≥ z) present in the last formula, are given by (17) and (17*).

Also, in the above case, ${\lambda }_1\left(u,t\right)$ , ${\lambda }_2\left(v,t\right)$ 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:

$\Psi \left(y,z,t\right)\le \left[{\lambda }_1\left(y,t\right)+{\displaystyle {\int }_{-\infty }^z\Psi }\left(y,v,t\right)\text{d}v\right]\times \left[{\lambda }_2\left(z,t\right)+{\displaystyle {\int }_{-\infty }^y\Psi }\left(u,z,t\right)\text{d}u\right]$ (19*)

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 $\Psi \left(y,z,t\right)$ that are solutions of the following integral equation directly derived from inequality (19*):

$\begin{array}{c}\Psi \left(y,z,t\right)=\left[{\lambda }_1\left(y,t\right)+{\omega }_1\left(t\right){\displaystyle {\int }_{-\infty }^z\Psi }\left(y,v,t\right)\text{d}v\right]\\ \cdot \left[{\lambda }_2\left(z,t\right)+{\omega }_2\left(t\right){\displaystyle {\int }_{-\infty }^y\Psi }\left(u,z,t\right)\text{d}u\right]\end{array}$ (20)

where $0\le {\omega }_i\left(t\right)\le 1$ , for $i=1,2$ .

Realize that the right-hand side of (20) is always less than or equal than the right hand side of (19*).

For equation (20) to hold the condition $\Psi \left(y,z,t\right)\ge 0$ is, in general, essential.

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

$\begin{array}{c}\Psi \left(y,z,t\right)=g\left(y,z,t\right)+{\omega }_1\left(t\right){\lambda }_2\left(z,t\right){\displaystyle {\int }_{-\infty }^z\Psi }\left(y,v,t\right)\text{d}v\\ +{\omega }_2\left(t\right){\lambda }_1\left(y,t\right){\displaystyle {\int }_{-\infty }^y\Psi }\left(u,z,t\right)\text{d}u,\end{array}$ (20*)

where $0\le g\left(y,z,t\right)\le {\lambda }_1\left(y,t\right){\lambda }_2\left(z,t\right)$ is any continuous function.

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

The set of solutions of (20*) is contained in the set of solutions of equation (19*), so that any solution of (20*) is a solution of (19*).

Equation (20*) can be simplified more. This will yield a smaller set of solutions being still solutions of (20*). Namely, setting $g\left(y,z,t\right)=0$ , we obtain the (purely) linear integral equation

${\omega }_1\left(t\right){\lambda }_2\left(z,t\right){\displaystyle {\int }_{-\infty }^z\Psi }\left(y,v,t\right)\text{d}v+{\omega }_2\left(t\right){\lambda }_1\left(y,t\right){\displaystyle {\int }_{-\infty }^y\Psi }\left(u,z,t\right)\text{d}u=\Psi \left(y,z,t\right)$ . (20**)

So the set of solutions of (20**) 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 ${\lambda }_1\left(y,t\right)={\lambda }_1$ and ${\lambda }_2\left(z,t\right)={\lambda }_2$ comprise the exponential case. But to stick with the “Wiener model” one must keep the hazard rates (coefficients) ${\lambda }_1\left(y,t\right)$ , ${\lambda }_2\left(z,t\right)$ in (20*) and in (20**) as corresponding to the distributions given by (17) and (17*).

Now, having any solution $\Psi \left(y,z,t\right)$ of any integral equation above, we find the joint survival function of the corresponding bivariate stochastic process as expressed by formula (18), where

$\alpha \left(y,z,t\right)={\displaystyle {\int }_{-\infty }^y{\displaystyle {\int }_{-\infty }^z\Psi }}\left(s,r,t\right)\text{d}s\text{d}r$ .

The so obtained joiner $\alpha \left(y,z,t\right)$ (or $\alpha \left(y,z\right)$ 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 {Y_t}, {Z_t } such that each random variable Y_t and Z_t 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:

$S\left(x\right)=P\left(X\ge x\right)={\left(x_0/x\right)}^{\alpha }$ ,

$R\left(y\right)=P\left(Y\ge y\right)={\left(y_0/y\right)}^{\beta }$ ,

where x ≥ x₀ > 0, y ≥ y₀ > 0, α > 0, β > 0.

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 x_0, y₀ (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 {X_t}, {Y_t} 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 X_t, Y_t are distributed according to the rules:

$S_t\left(x\right)=P\left(X_t\ge x\right)={\left(x_0\left(t\right)/x\right)}^{\alpha \left(t\right)}$ ,

$R_t\left(y\right)=P\left(Y_t\ge y\right)={\left(y_0\left(t\right)/y\right)}^{\beta \left(t\right)}$ .

For simplicity we will assume that x₀(t) = x₀ = constant, and y₀(t) = y₀ = constant.

The above univariate “Pareto stochastic processes” are determined by the constants x_0, y₀, and the way we define the functions α(t), β(t) for t = 1, 2, 3, …. Our (simplest) proposition is to define

$\alpha \left(t\right)=A+\left(t-1\right)r,\beta \left(t\right)=B+\left(t-1\right)s$ ,

where (as a particular example) 1 ≤ A ≤ 2, 1 ≤ B ≤ 2, r = 0.10, s = 0.15.

Now, for each t, we find the joint survival functions $U_t\left(x,y\right)$ of the random vector $\left(X_t,Y_t\right)$ . Recall the general form of the joint survival functions:

$U_t\left(x,y\right)=S_t\left(x\right)R_t\left(y\right)\exp \left[-{\alpha }^{*}\left(x,y,t\right)\right]$ ,

(here, ${\alpha }^{*}\ne \alpha$ and the meaning of ${\alpha }^{*}$ is different than that of $\alpha$ ).

Since we already have both S_t(x) and R_t(y) as given above, we need to find a class of proper joiners ${\alpha }^{*}\left(x,y,t\right)$ to make, for each t, the joint survival function $U_t\left(x,y\right)$ of $\left(X_t,Y_t\right)$ properly defined. Set ${\alpha }^{*}\left(x,y,t\right)={\alpha }_t^{*}\left(x,y\right)$ .

Recall that we have:

${\alpha }_t^{*}\left(x,y\right)={\displaystyle {\int }_{x_0}^x{\displaystyle {\int }_{y_0}^y{\Psi }_t}}\left(u,v\right)\text{d}u\text{d}v$ ,

where, for every t, the function ${\alpha }_t^{*}\left(x,y\right)$ must satisfy conditions (1) - (3) from section 3.1. Since the domain of this function is reduced to $\left[x_0,+\infty \right)\times \left[y_0,+\infty \right)$ , with positive x0 and y0, condition (3) (section 3.1) now takes the form:

${\alpha }_t^{*}\left(x_0,y\right)={\alpha }_t^{*}\left(x,y_0\right)=0$ , for each $t$ .

Also, for each $t$ , we have $S_t\left(x_0\right)=1$ and $R_t\left(y_0\right)=1$ .

Therefore, as pointed out in Property 1 (section 3.1), it follows that, for each t, the functions S_t(x), R_t(y) are the marginal survival functions for the joint $U_t\left(x,y\right)$ . Moreover, we restrict ourselves to positive stochastic dependences only, so we have, for every $t$ , ${\Psi }_t\left(u,v\right)\ge 0$ . [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 ${\alpha }_t^{*}\left(x,y\right)$ [which are equivalent to the corresponding functions ${\Psi }_t\left(u,v\right)$ ], we need to examine (for each t) inequality (9*) which is

${\Psi }_t\left(u,v\right)\le {\lambda }_1\left(u,t\right){\lambda }_2\left(v,t\right)$ , (21)

where, in this ‘Pareto case’, we have ${\lambda }_1\left(u,t\right)=\alpha \left(t\right)/u$ and ${\lambda }_2\left(v,t\right)=\beta \left(t\right)/v$ .

Thus inequality (21) takes the form ${\Psi }_t\left(u,v\right)\le \alpha \left(t\right)\beta \left(t\right)/uv$ , for every $t$ .

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

${\Psi }_t\left(u,v\right)=\omega \left(t\right)\alpha \left(t\right)\beta \left(t\right)/uv$ ,

where, for every $t$ , $0\le \omega \left(t\right)\le 1$ .

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

The natural assumption in this Pareto case is that both x₀ ≥ 1, y₀ ≥ 1 (i.e., the minimal incomes in both social groups are at least 1 “unit”). Taking the above under the consideration we may choose for ${\Psi }_t\left(u,v\right)$ (the determinant of the stochastic dependence between X_t and Y_t) the following:

${\Psi }_t\left(u,v\right)=\alpha \left(t\right)\beta \left(t\right)/u^{{\gamma }_1}{v}^{{\gamma }_2}$ , for ${\gamma }_1\ge 1$ , ${\gamma }_2\ge 1$ .

This ${\Psi }_t\left(u,v\right)$ satisfies (21).

Finally, the full version of the so derived joint survival function of the Pareto random vectors $\left(X_t,Y_t\right)$ , for every t, is:

$U_t\left(x,y\right)=S_t\left(x\right)R_t\left(y\right)\exp \left[-\alpha \left(t\right)\beta \left(t\right){\displaystyle {\int }_{x_0}^x{\displaystyle {\int }_{y_0}^y{u}^{-{\gamma }_1}{v}^{-{\gamma }_2}\text{d}u\text{d}v}}\right]$ , ${\gamma }_1\ge 1$ , ${\gamma }_2\ge 1$ .

Annotation. Just in recent days, when this paper was already finalized, we learned that similar results as for the bivariate distributions (but not bivariate stochastic processes) were considered by Finkelstein [4] in 2003. The author also considered the “Aalen factor” (see formulas (4) and (5) in [4]) but under different name 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 apparently was overlooked. But, anyway, the novum of our results became somehow (although not drastically) limited. What makes our own contribution to the theory of bivariate probability distributions (and not only to the theory of bivariate stochastic processes), besides the indication on its universality, is some more in detail performed analysis of the conditions under which continuous functions of two real variables may serve as proper joiners for any two, given in advance, fixed univariate distributions or univariate stochastic processes. Such conditions often rely on being solutions of nontrivial integral equations as, for example, the equations given by the formulas (20), (20*), (20**). Their eventual solutions are, at the same time, solutions of the most basic inequality (19*) that ensures legitimacy of the so obtained bivariate distributions. An underlying theory of those integral equations seems to be a nontrivial one. There mentioned, integral equations can actually be applied to both: bivariate distributions as well as to the considered bivariate stochastic processes.

Also, on our side lies general indication on the joiners’ theory as the one competitive to the copula methodology.

5. Conclusions

1) Applying the joiners as measures of stochastic dependencies for both, bivariate random vectors and bivariate stochastic processes enable an easy method for underlying constructions as well as provides a general theory on the topic, commonly known under the name of “multivariate stochastic analyses”. To a measure, this general theory may, possibly, be competitive to the theory associated with the notion of copula.

2) This can easily be seen, that the performed in this paper bivariate constructions can be extended to k-dimensional (k > 2) random vectors and processes. However, when the dimension was increased, the phenomena of r-dependence (k ≥ r ≥ 3) should be taken into account (for that see [8]).

3) Using the joiners as functions of time, new multivariate time dependent stochastic models can be constructed. The range of applications of those models comprises almost all areas of real life wherever multivariate data are to be analyzed.

4) Statistical analysis must follow the construction of new models. New statistics for estimating or testifying various types of the joiners and their parameters should be invented and analyzed. This rich subject was, however, out of scope of this paper which has a rather purely probabilistic character.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Cox, D.R. (1972) Regression Models and Life-Tables. Journal of the Royal Statistical Society Series B: Statistical Methodology, 34, 187-202.[CrossRef]
[2]	Aalen, O.O. (1989) A Linear Regression Model for the Analysis of Life Times. Statistics in Medicine, 8, 907-925.[CrossRef] [PubMed]
[3]	Filus, J.K. and Filus, L.Z. (2017) The Cox-Aalen Models as Framework for Construction of Bivariate Probability Distributions, Universal Representation. Journal of Statistical Science and Application, 5, 56-63.[CrossRef]
[4]	Finkelstein, M.S. (2003) On One Class of Bivariate Distributions. Statistics & Probability Letters, 65, 1-6.[CrossRef]
[5]	Gumbel, E.J. (1960) Bivariate Exponential Distributions. Journal of the American Statistical Association, 55, 698-707.[CrossRef]
[6]	Filus, J.K. and Filus, L.Z. (2008) On Multi-Component System Reliability with ‘Microshocks–Microdamages’ Type of Components Interaction. Proceedings of the International Multi Conference of Engineers and Computer Scientists, Hong Kong, 19-21 March 2008, 1946-1951.
[7]	Sklar, A. (1959) Fonctions de Répartition à n Dimensions et Leurs Marges. Publications de l’Institut Statistique de l’Université de Paris, 8, 229-231.
[8]	Filus, J.K. and Filus, L.Z. (2023) Stochastic “Three-Dependence” in Application to System Reliability. Research Gate Blog, Jerzy Kazimierz Filus Account.