Study of the bivariate survival data using frailty models based on Lévy processes

Frailty models allow us to take into account the non-observable inhomogeneity of individual hazard functions. Although models with time-independent frailty have been intensively studied over the last decades and a wide range of applications in survival analysis have been found, the studies based on the models with time-dependent frailty are relatively rare. In this paper, we formulate and prove two propositions related to the identifiability of the bivariate survival models with frailty given by a nonnegative bivariate Lévy process. We discuss parametric and semiparametric procedures for estimating unknown parameters and baseline hazard functions. Numerical experiments with simulated and real data illustrate these procedures. The statements of the propositions can be easily extended to the multivariate case.


Introduction
In survival analysis, frailty is usually defined as a non-observable momentan risk of failure and is included in survival models in the form of an unknown nonnegative random variable or random process characterizing non-homogeneity of population with respect to hazard function. Usually, frailty enters the model multiplicatively to the hazard function and allows us to take into account the correlations between failure times. Observed covariates can also be included in the multivariate survival models in the form of Cox-like regression. Identifiability and other properties of the univariate and multivariate survival models with time-independent frailty have been studied in some depth, and these models now are popular and widely used in survival studies and in the search for genes influencing longevity. Several models with time-dependent frailty have also been suggested. Woodbury and Manton (1997) introduced a stochastic process model of human mortality and ageing. They considered hazard as a quadratic function of stochastically changing unobserved frailty. The model was further extended by Yashin and Manton (1997) to consider a partially observed frailty process. Modifications of this model were successfully applied to the analyses of longitudinal data with informative censoring. Ideas on further development and applications of these results to studying ageing and longevity are summarized by Yashin et al. (2012). Gjessing et al. (2003) considered an approach based on the proportional hazard model with time-dependent frailty given by the formula Here, X (t) is a nonnegative Lévy process (subordinator) with Laplace exponent (c), R(t) = t 0 r (u)du defines the time transformation of subordinator for a nonnegative rate function r (t), and the nonnegative weight function a(u, s) determines the extent to which the previous behavior of transformed subordinator influences the hazard function at time t. The authors derived the expressions for the population survival and hazard functions in a general case: where λ(t) is the baseline hazard function and b(u, t) = t u λ(s)a (u, s−u)ds. Gjessing et al. (2003) showed also that under some conditions, quasi-stationary distributions of survivors can arise. This implies the constant limiting population hazard rate, in spite of the increase of the baseline hazard function. Ata and Özel (2012) considered the proportional hazard model with time-dependent frailty given by the discrete compound Poisson process and applied this model to the earthquake data and to traffic accidents data from Turkey.
The aim of this paper is to study the problem of identifiability for bivariate survival models with/without observed covariates and with time-dependent frailty when this frailty is given by a Lévy process (or Lévy processes). In addition, we demonstrate how these models can be used for longevity datasets based on simulated data.
In Sect. 2, we give the definitions of the univariate survival model under mixed proportional hazard specification and the bivariate correlated frailty model. We then discuss the definitions of the uni-and bivariate survival model with time-dependent frailty given by the nonnegative Lévy processes. At the end of this section, we discuss known findings related to the identifiability of survival models with time-independent frailty, formulate two new propositions about the identifiability of the bivariate survival models with time-dependent frailties, and give the EM algorithm for estimating unknown baseline hazard functions and parameters for the correlated bivariate model with time-dependent frailty. In Sect. 3, we discuss the results of estimations of the baseline hazard functions and unknown parameters (including the Cox-regression coefficients and parameters characterizing the frailty process) in experiments with simulated and real data for parametric and semiparametric approaches. The real-world example was based on the data from the Danish Twin Registry. Conclusions and outlook are presented in Sect. 4.

Survival analysis under a frailty setting
Under mixed proportional hazard specification (Abbring and van den Berg 2003), the hazard rates of the failure times for two related individuals depend multiplicatively on the respective baseline hazards λ j (t), regressor functions χ j (u j ) with observed vector of covariates u j , and unobserved nonnegative random variable (frailty) Z j , characterizing the heterogeneity in the population with respect to hazard λ j The function χ j (u j ) is frequently specified as χ j (u j ) = exp(β * j u j ) (the Coxregression term) for some transposed vector parameter β * j , j = 1, 2. The univariate population survivals S j (t|χ j (u j )) = P(T j > t j ) for random times of death T j are with cumulative hazard function j (t) = t 0 λ j (τ )dτ and Laplace transform L(c) = Ee −cZ .
In the bivariate correlated frailty model proposed by Yashin et al. (1995), individual frailties (Z 1 , Z 2 ) were constructed under assumptions for independent nonnegative random variables Y 0 , Y 1 and Y 2 and some nonnegative constants m 0 , m 1 and m 2 . Given frailties (Z 1 , Z 2 ), life spans for both individuals were assumed to be conditionally independent. Since the scale factor common to all subjects in the population can be absorbed into the baseline hazard functions λ j (.), j = 1, 2, we can put EZ 1 = EZ 2 = 1. The heterogeneity in the population can be characterized by the variance of frailties VarZ 1 = σ 2 1 , VarZ 2 = σ 2 2 . The correlation between frailties Corr(Z 1 , Z 2 ) we will denote by ρ.
The assumption that the individual frailty is determined at birth and does not change with age seems to be too strong and unrealistic. To make the approach more flexible, we can weaken this assumption and suppose that the frailty is a random process.
Similarly to Gjessing et al. (2003), we shall consider the frailty process Z = Z (t) : t 0 defined by a nonnegative Lévy process. In accordance with Lévy-Khinchin formula, such a process can be characterized by its Laplace transform with Laplace exponent (c), while c is the argument of Laplace transform. Note that Examples of Lévy processes include the standard compound Poisson process, the compound Poisson process with general jump distribution, gamma processes, stable processes, and PVF (power variance function) processes (Gjessing et al. 2003). In this paper, we shall consider the nonnegative Lévy processes (subordinators) with the Laplace exponent given by nonnegative drift d and the jump measure ν with support (0, ∞) satisfying ∞ 0 min(1, x)ν(dx) < ∞. The Laplace exponent is an increasing and concave function. The Laplace exponent and the jump measure for the gamma process are given by the formulas with the shape ht and the scale parameter γ for the gamma-distributed random variable Z (t). This corresponds to the values of htγ −1 and htγ −2 for mean and variance of Z (t), respectively. To avoid non-identifiability of the model, we shall standardize the frailty distribution and put h = γ . In this case, EZ (t) = t, VarZ (t) = tγ −1 for t 0. In "Appendix 1," we can find the formulas for calculating univariate population survivals in this case for a constant and an exponential (Gompertz) baseline hazard functions. We will denote the Laplace exponent for the univariate frailty processes Z j (t), j = 1, 2 by j (.). Let (Z 1 (t), Z 2 (t)) be a bivariate Lévy process with nonnegative components and the Laplace exponent given by where d, c ∈ R 2 + , < x, y > denotes the dot product of vectors x, y ∈ R 2 , the Lévy measure ν(A) for any Borel set A ∈ R 2 + is the expected number of jumps on the time interval [0,1], whose sizes belong to A, and the following integrability conditions for Lévy measure are satisfied: Note that biv (c 1 , 0) = 1 (c 1 ) and biv (0, c 2 ) = 2 (c 2 ) for marginal Lévy processes.

Model identifiability
The identifiability of the univariate model with unspecified functional form of frailty distribution and baseline hazard has been studied by Elbers and Ridder (1982). This model is identifiable given information on T for finite EZ and is not identifiable when frailty has an infinite mean. Identifiability of the correlated frailty models using data on the pair (T 1 , T 2 ) was proved by Honoré (1993) under the assumption of finite means of Z 1 and Z 2 . Yashin and Iachine (1999) proved the identifiability of the correlated frailty model without observed covariates assuming that Z 1 and Z 2 are gamma distributed. Abbring and van den Berg (2003) studied the identifiability of the mixed proportional hazards competing risks model. We adopt this method to investigate the identifiability of the mixed bivariate survival model for time-dependent correlated frailties.
Then, the mixed bivariate survival model for time-dependent correlated frailties is identified from the bivariate distribution of the failure times (T 1 , T 2 |u 1 , u 2 ).
Differentiating S j (t|χ j ) in t, dividing then by χ j , and setting formally χ j → 0, we get the following equations: Since the expression on the right hand of (9) is observed, we find the hazard function in the form The constant j (0) can be found from the equation j (t * ) = 1, j = 1, 2 (Assumption 2).
If individual frailties (Z 1 , Z 2 ) are constructed under assumptions given by for some positive α, we can weaken the conditions of Proposition 1 and prove the identifiability of the model in the absence of observed covariates.
Proposition 2 Let the following assumptions be satisfied.

Assumption 2. Baseline hazard functions
for some real number b > 0 and i = 0, 1, 2. Then, the mixed bivariate survival model for time-dependent correlated frailties is identified from the bivariate distribution of the failure times (T 1 , T 2 ).
The proof of Proposition 2 is given in "Appendix 3."

Model validation
In this section, we will assume that χ(u) = exp(β * u). To validate the regression model, we need to estimate the vectors of Cox-regression parameters β 1 and β 2 , vector parameter defining frailty ζ [in the case of correlated frailty model either (σ 2 1 , σ 2 2 , ρ) for time-independent gamma-frailty or (γ −1 1 , γ −1 2 , ) for gamma-frailty process given by the Laplace exponent (6)], and the baseline hazard function λ 1 (t) and λ 2 (t). If the baseline hazard functions follow a parametric form such as the Gompertz or the Weibull function with vector parameter ξ , we can use the classic maximum likelihood method to estimate all unknown parameters. The log-likelihood in this case is given by Here, θ = (β, ζ, ξ) is the vector of unknown parameters, β = (β 1 , β 2 ) stands for Cox-regression coefficients, ζ for frailty parameters, and ξ for parameters defining the baseline hazard functions λ j (t), j = 1, 2. The data include information about life spans (t i1 , t i2 ), observed covariates (u i1 , u i2 ), and censoring (δ i1 , δ i2 ) for a twin pair i, i = 1, . . . , n. The estimate of the vector parameter θ we can find by maximizing the log-likelihood function (13).
If the form of the baseline hazard functions is not specified, the estimates can be obtained by the EM algorithm. This algorithm combines the maximum partial likelihood estimator of the vector parameter (β, ζ ) with the Breslow estimator (Breslow 1972) of the cumulative baseline hazard function (t) = ( 1 (t), 2 (t)). The EM algorithm is an iterative procedure with two steps-E (expectation) and M (maximization) on each iteration. It works as follows: Let f (z i1 , z i2 |t i1 , t i2 , ζ ) be the density function of the random variable (Z i1 (t i1 ), Z i2 (t i2 )) given parameter vector ζ , i = 1, . . . , n. Denote the estimates of (t), ζ , and β on the l th iteration byˆ l (t), ζ l , andβ l , respectively. Similar to Gorfine and Hsu (2011), we define the failure counting process N i j = δ i j Ind(T i j t) and the at-risk process X i j = Ind(T i j t), where T i j is the random time to the failure of twin j in twin pair i, i = 1, . . . , n, j = 1, 2. Define random processes where the symbol "ˆ" means replacing the unknown frailty Z i j (t) with its conditional expectation given the observed data and the current estimates of (t), ζ and β.
The convergence of the EM algorithm and the properties of parameter estimates have been discussed elsewhere (Zeng and Lin 2007). For the correlated frailty model with time-independent gamma-distributed frailty, the calculation of expression (16) has been discussed in detail by Iachine (1995). The calculation of this expression in the case of the gamma-frailty process for the same model can be found in "Appendix 4." In the parametric approach, the choice of the appropriate baseline hazard function plays an important role. The Gompertz function does not guarantee the good fit of the marginal survival function for real longevity data. The following gamma parameterization of the univariate survival function gives better results by fitting the real data in the model with time-independent frailty is the pseudo-baseline cumulative hazard, t 30, s 2 ,ã,b > 0, and σ 2 > 0 is the variance of the time-independent frailty. Given parameters s 2 , σ 2 ,ã, andb, it is not difficult to find the true baseline cumulative hazard (t) in the form In the case of the time-dependent frailty given by a Lévy process with Laplace exponent (c), we need to consider the following analog of Eq. (17): Unfortunately, Eq. (18) does not have a closed-form solution with respect to (t). In experiments with real data, we will find the approximative solution to (18) such that the function dyn appr (t) is a non-decreasing, nonnegative, piecewise-constant function satisfying the following conditions: dyn appr (0) = 0, for times-to-event t i , i = 1, . . . , 2n, sorted in non-decreasing order, 0 can be calculated recurrently for i = 1, 2, . . . , 2n using a simple bisectional procedure. Note that the function dyn appr (t) converges pointwise to the solution of (19) as n → ∞ and the distance between neighboring moments t i tends to zero.
To compare two approaches, we assume that in the case of the time-independent frailty, the cumulative hazard is also a non-decreasing, nonnegative, piecewiseconstant function stat appr (t) satisfying the following conditions: stat 3 Results

Experiments with simulated data
In this subsection, we will discuss the results of the consistency test for the correlated frailty models with time-dependent and time-independent frailties (11). It was assumed that α = 1, VarZ 1 = VarZ 2 = σ 2 , and Corr(Z 1 , Z 2 ) = ρ in the case of the time-independent frailty or VarZ 1 (1) = VarZ 2 (1) = γ −1 and Corr(Z 1 (1), Z 2 (1)) = in the case of the time-dependent frailty, baseline hazard functions λ j (t) followed Gompertz (exponential) form λ j (t) = a exp(bt), and an observed covariate u influenced longevity so that the conditional hazard function was defined by μ j (t j |Z j , u j ) = Z j exp(βu j )λ j (t), j = 1, 2. The covariates were randomly generated from the uniform distribution on the interval [0,1] and were independent for the individuals. The (true) values for data generating are given in Tables 1 and 2 and have been chosen so that the mean and the standard deviation of the generated times-to-event were equal to approximately 75 and 12 years, respectively. The bivariate times-to-event have been generated using formula (4) with σ 1 = σ 2 = σ in the case of the time-independent frailty and using formula (25) given in "Appendix 2" in the case of time-dependent frailty. In both cases, it was assumed that 1 (t) = 2 (t) = (a/b)(exp(bt) − 1) and χ 1 (u) = χ 2 (u) = exp(βu). We have not truncated or censored the generated data. We estimated unknown parameters and cumulative hazard functions using the classic maximum likelihood estimator (parametric method) and the EM algorithm (semiparametric method). In all cases, we simulated 100 datasets for 500 twin pairs. Table 1 shows the results of simulation study for the time-independent gammafrailty model without truncation and censoring. Empirical means and standard deviations of estimates were calculated using the classic maximum likelihood method and the EM algorithm, respectively. Table 2 shows the results of the simulation study for the time-dependent gammafrailty model without truncation and censoring. Empirical means and standard deviations of estimates were also calculated using the classic maximum likelihood method and the EM algorithm, respectively. Analysis of estimates in both tables does not indicate the presence of any bias, and estimates calculated using the classic maximum likelihood estimator are generally more efficient. Furthermore, the estimates of  the Cox-regression coefficients and parameters characterizing the frailty distribution are closer to true vales than those for the EM algorithm. One can see in Fig. 1 that in both cases (the time-dependent and the time-independent frailty), the empirical mean log baseline cumulative hazard trajectory calculated using the EM algorithm fits the true log baseline cumulative hazard trajectory very well.

Experiments with real data
For experiments with real data, we used the datasets from the Danish Twin Registry (DTR). This registry was created in the 1950s. It is one of the oldest populationbased registries in the world and contains information about twins born in Denmark since 1870 and who survived to age 6. Multiple births were manually ascertained in birth registers from all 2200 parishes in Denmark. As soon as a twin was traced, a questionnaire was mailed to the twin, to her/his partner or to their closest relatives if neither of the twin partners were alive. The zygosity of twins was assessed on the basis of questions about phenotypic similarities. The reliability of the zygosity diagnosis was validated by comparing laboratory methods based on the blood, serum, and enzyme group determination. In general, the misclassification rates were less than 5%. Other information includes the data on sex, birth, cause of death, health, and lifestyle. An important feature of the Danish twin survival data is their right censoring and left truncation. In our study, we used the longevity data on the like-sex twins with known zygosity born between 1870 and 1900 and who survived until age 30. This non-censored data include 470 male monozygotic (MZ) twin pairs, 475 female MZ twin pairs, 780 male dizygotic (DZ) twin pairs, and 835 male DZ twin pairs. Further details on the Danish Twin Registry can be found in Hauge (1981).
Since the EM algorithm suffers from its slow convergence, we have estimated unknown parameters for the time-independent and the time-dependent frailty models using the classic maximum likelihood method. Table 3 shows these estimates, the logarithm of the maximum value of the likelihood function, and the value of the AKAIKE Information Criterion (AIC) for the data from the Danish Twin Registry described above. The estimates of parameter s 2 were very close to zero in all experiments with real data. We have put this parameter equal to zero to avoid the efficiency loss of estimates. The AIC values for the model with time-dependent frailty are slightly smaller than the ones for the model with time-independent frailty. That is, the model with time-dependent frailty is slightly more informative than the one with time-independent frailty. Figures 2 and 3 show estimated and empirical bivariate probability density functions for the time-independent and the time-dependent frailty models. Note that the shapes of the estimated bivariate probability density functions are very similar.

Discussion
Frailty models are a powerful tool for studying non-observable inhomogeneity in a population related to time-to-failure (e.g., death or disease). Models with timeindependent frailty have been intensively studied over the last decades and have found a wide range of applications in survival analysis and in searching for genes influencing longevity. However, the studies based on the models with time-dependent frailty are scarce. In this paper, we have attempted to improve the knowledge in this area and to study some properties of multivariate survival models with time-dependent frailty components. Proposition 1 we have formulated and proved for the bivariate case. It is not difficult to generalize this result and to prove the identifiability of the frailty model with observed covariates for any number J of related individuals equal or greater than 1 if the time-dependent frailty is a multivariate Lévy process. Similarly, we can generalize Proposition 2 for the case of J 2. However, the number of frailty components in the multivariate analog of the decomposition (11) will be equal to 2 J − 1. The shared frailty model where all individuals in a family or cluster share the same non-observable risk of failure does not meet this problem.
In experiments with simulated data, we tested for consistency and used parametric and semiparametric approaches. In the parametric approach, we assumed that the parametric form of the baseline hazard functions is known and follows the Gompertz form. All unknown parameters characterizing frailty distribution, baseline hazard function, and Cox-regression parameters were estimated directly by maximizing the likelihood function. In the semiparametric approach, we used the EM algorithm and estimated the Cox-regression parameters and the parameters of frailty distribution by maximizing the partial likelihood function. The cumulative baseline hazard function was estimated using the Breslow estimator regarding this function as infinite-dimensional parameters. The EM algorithm suffers from its slow convergence. Moreover, in the semiparametric approach, the number of calculations increases with the number of individuals much more rapidly than in the parametric one. It leads to the drastic slowing of the convergence of the EM algorithm and increases substantially the time of Experiments with real data show that the proposed method and the method with time-independent frailty produce similar shapes of the estimated bivariate probability density functions. The baseline cumulative hazard functions have been chosen so that the estimated marginal survival functions guarantee the best fit to the empirical ones according to . A large degree of similarity of the estimated bivariate density functions for the models with time-dependent and time-independent frailties in the range of ages 30-100 years guarantees the similar bivariate fit. The difference between the two approaches can involve the shape of the baseline hazard functions and the asymptotic behavior of the bivariate probability density functions. The models with time-dependent and time-independent frailties are not nested. Therefore, we cannot compare them using the likelihood ratio test. For this purpose, the AKAIKE Information Criterion can be used. In accordance with this criterion, the model with time-dependent frailty is slightly more informative compared to the one with timeindependent frailty for the data we considered.  Gorfine and Hsu (2011) studied the robustness of the multivariate survival models with frailty components against the violations of the model assumptions. It was found that unnecessary modeling of the dependency between the frailty variates can lead to some efficiency loss of parameter estimates. Misspecification of the frailty distribution can introduce a bias in estimates. Misspecification of the baseline hazard functions can lead to severe bias of all estimates if we use the parametric maximum likelihood estimator, where the baseline hazard functions follow the parametric form. The nonparametric maximum likelihood estimator does not suffer from this drawback. Note that in experiments with real data, we have used a flexible parametrization of the baseline cumulative hazard functions given by formulas (19)-(20). This parameterization does not presume any knowledge about the form of the baseline hazard function. It is sufficient to have a good approximation of the marginal survival function.
An extension of the present study can include the investigation of identifiability of the survival models with competing risks and time-dependent frailty components. The piecewise-constant approximation of the cumulative hazard function has been used in experiments with real data [formulas (19)-(20)]. Other approximative functions such as piecewise linear or piecewise exponential can be used to improve the bivariate goodness-of-fit. Further, numerical experiments with real data are needed to understand whether the proposed method improves the goodness-of fit on the method with time-independent frailty.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix 1: Univariate survival function for frailty model based on a gamma process
0} is a gamma process with Laplace transform given by the Lévy-Khinchin formula where G (c) is the Laplace exponent of the Lévy process Z , the argument for Laplace transform c 0 and Here, ht and γ are the shape and the scale parameters, respectively, for gammadistributed random variable Z (t). Since Z (t) is a.s. increasing in t, we can consider it as subordinator. It is easy to check that where superscript "(l)" stands for lth-order argument derivative. In nonhomogeneous population, the process Z (individual for each member of population) can characterize the individual risk of mortality (frailty). In accordance with the definition of the frailty model under mixed proportional hazard specification, we have for covariate vector u, parameter vector β, and baseline hazard λ(t). In Gjessing et al. (2003) (Theorem 1), it was shown that the population survival in this case can be calculated For constant λ(t) = λ 0 , it holds that (τ, t; u) = λ 0 e β * u (t − τ ), Using the following formula for indefinite integral we get for the baseline hazard given by the Gompertz function λ(t) = a exp(bt) Here, f (t, u) = (aχ(u)/(bγ (1 + aχ(u)e bt /(bγ )) and Li 2 (z) is the dilogarithm function defined by This function is an analytical extension of the infinite series
If the baseline hazard function is given by the Gompertz function, we can use formula (22) to calculate the bivariate population survival.

Appendix 3: The proof of Proposition 2
To prove Proposition 2, we need to check several auxiliary statements.
Lemma 1 Let c 1 , c 2 0. Then, Proof Assume that c 1 c 2 0 and set c = c 1 − c 2 . Then, It is easy to check that Therefore, Last result holds also in the case of c 2 > c 1 0.
Lemma 2 Let (c) be a Laplace exponent given by (5), and (t) is an increasing continuous function defined in R + . Then, Proof The statement of lemma follows directly from formula (5) and monotone increase of the function (.).
The bivariate survival function for correlated frailty processes given by (11) is where i (.), i = 0, 1, 2, is the Laplace exponent for the frailty process Y i (t) and Using formula (28) and Assumptions 2-3 of Proposition 2, it can be shown that Note that expressions on the left hand of (29) are observable, g 1 (t) = g 0 (t, 0) and g 2 (t) = g 0 (0, t).
Let L 1 be a set of all continuous increasing functions (t) on R + satisfying (0) = 0 and (t) → ∞ as t → ∞. Given Laplace exponent 0 , consider the map A 1 defined by Lemma 3 A 1 is a bijective map from L 1 to image A 1 (L 1 ) of L 1 under A 1 .
Proof It is easy to see that (A 1 )(t) is a continuous monotone increasing function with (A 1 )(0) = 0 and (A 1 )(t) → ∞ as t → ∞ for ∈ L 1 . That is, A 1 (L 1 ) ⊂ L 1 . For any 0 < T < ∞, the set L T 1 of all continuous monotone increasing functions on [0, T ] satisfying (0) = 0 is a complete metric space with usual supremum distance Taking into account that (t) is differentiable in t a.e. and 0 (0) = 0, we get that and, therefore, In accordance with the mean value theorem applied to the integral on the right hand of (31), there exists τ ∈ [0, τ ] such that Since ∞ 0 xν(dx) is finite, the function on the left hand of (33) is integrable on [0, t], t < ∞. Assume that two solutions a , b of (32) satisfy a (t) = b (t) on the interval [0, t 0 ], t 0 0. We shall show now that these solutions satisfy this property on the interval [t 0 , t 0 + t ε ] for some t ε > 0 too. Let the map B 1 be given by Taking into account g (τ ) 0 and inequalities (26) and (27), we get that for all t ∈ [0, t 0 + t ε ] and some positive q < 1. It means that Since this map is also surjective, it is bijective Define the set L 2 of all functions 12 (t 1 , t 2 , α) such that for 1 , 2 ∈ L 1 , and α > 0. Given Laplace exponent 0 , consider the map A 2 defined by Proof Since (t) = 12 (t, t) ∈ L 1 , 1 (t) = 12 (t, 0) ∈ L 1 , and α 2 (t) = 12 (0, t) ∈ L 1 for t ∈ [0, ∞), in accordance with Lemma 3, these functions can be uniquely defined by their images g 0 (t, t) = (A 2 1,2 )(t, t) = (A 1 )(t), g 0 (t, 0) = (A 2 12 )(t, 0) = (A 1 1 )(t), and g 0 (0, t) = (A 2 (α 12 ))(0, t) = (A 1 (α 2 ))(t), respectively. Assume now that t 2 > t 1 and the solution 12 (t 1 , t) to (34) is uniquely defined for 0 t t 0 with t 0 t 1 . Define the map B 2 for t > t 0 by Similarly to the proof of Lemma 3, it can be shown that the solution to is also uniquely defined on the interval [t 0 , t 0 + t ε ] for some t ε > 0 and, therefore, is uniquely defined for all t t 1 . The same result can be checked in the case of t 1 > t 2 .
For the sake of contradiction, assume now that there are functions a 0 (.
, and real positive numbers α a and α b such that holds for all t 1 , t 2 0. In accordance with Lemmas 3-4, the maps A a i and A b i are invertible on their images of L i , i = 1, 2. Moreover, Note that i 1 (t 1 ) = i 12 (t 1 , 0) and α i i we get from (36) Now we will prove that the action of the map T ba 2 (resp. T ab 2 ) on the function a 12 (resp. b 12 ) calculated in the point (t 1 , t 2 ) depends only on the value of b 12 (t 1 , t 2 ) (resp. a 12 (t 1 , t 2 )).
Lemma 5 Given the maps a 0 (.), b 0 (.), a 1 (.), b 1 (.), a 2 (.), b 2 (.), and real positive numbers α a and α b , there exist uniquely defined, monotone increasing, continuous functions f ab : R 1 + → R 1 + and f ba : Proof Put α a 1 = α b 1 = 1. From (33) and Lemma 3, it follows that given the functions are strictly monotone increasing, continuous, and uniquely defined. It holds that . Therefore, the maps f ba , are also strictly monotone increasing and continuous. Since for some continuous monotone increasing functions f ba : R 1 + → R 1 + and f ab : R 1 + → R 1 + . Substituting in (39) t 2 = 0 or t 1 = 0, we get that f ab 1 = f ab 2 = f ab and f ba Proof of Proposition 2 Note firstly that in accordance with Assumption 4, the functions i (c), i = 0, 1, 2, are analytic ones on Re(c) > −b. Since f ba (.) and f ab (.) are defined on R 1 + , additive and continuous, they are linear and it holds that f ab (x) = C ab 1 x + C ab 2 and f ba (x) = C ba 1 x + C ba 2 for some constant C ab 1 , C ab 2 , C ba 1 , and C ba 2 . As f ab (0) = f ba (0) = 0, we have that C ab 2 = C ba 2 = 0. Moreover, C ab 1 = C ba 1 = 1 and α a = α b because a i (t * ) = b i (t * ) = 1 for some t * > 0 and i = 1, 2. From here, it follows that f ab (.) and f ba (.) are identity transformations and a i (t) = b i (t) = i (t) for all t 0, i = 1, 2. It holds also that (A a 1 i )(t) = (A b 1 i )(t) for all t 0, i = 1, 2. To complete the proof of Proposition 2, we need to show that a 0 (x) = b 0 (x) for any x ∈ R 1 + . Indeed, the maps a 0 and b 0 are real analytic functions on R 1 + and, therefore, ( 1 (t) − 1 (τ )) n dτ for small values of 1 (t). This holds iff ( a 0 ) (n) (0) = ( b 0 ) (n) (0) for all nonnegative integer n and means that a 0 (x) = b 0 (x) = 0 (x) in some neighborhood of 0. Since 0 (x) is a real analytic function on R 1 + , this function is uniquely defined in this area. Similarly, it can be proved that the functions i (x), i = 1, 2, are also uniquely defined in R 1 + .