Geyser eruption is one of the signature attractions at the Yellowstone National Park, which is home to two-thirds of the worlds' geysers. Tourists around the world crave to witness this fascinating natural phenomenon. Many researchers are interested in studying geyser eruptions and the underlying mechanisms. Fournier (1969) built a physical model to describe the time interval between eruptions for the Old Faithful geyser, which is one of the most famous geysers in the Yellowstone National Park. Rinehart (1972) showed that the Old Faithful Geyser activities are affected by earth tidal forces, barometric pressure, and tectonic stresses. However, geyser eruptions have not been studied by statistical methods. This paper develops a statistical model for analyzing geyser eruption data from a bivariate geyser system. This work will benefit the geyser study community for understanding and effectively modeling geyser eruption activities.
Typical geyser eruption is a repeating process and hence can be modeled with a recurrent process for events repeatedly occurring over time. Recurrent processes have had broad applications in diverse areas. For example, they have been widely used for studying vehicle failures in warranty studies (Lawless, 1995), relapse biomarkers in cancer research (Schaubel and Cai, 2004), and sports injury analysis (Ullah et al., 2014). Typically the time interval between two consecutive events, which is also referred to as the gap time, is studied to model the event frequency in a recurrent process. The proportional intensity models (Cox, 1972, andAndersen andGill, 1982) are popular for modeling event occurrences of a single type of event. In more sophisticated studies, there are multiple types of recurrent events observed in a single system. The occurrence of any type of event will result in a system event. In addition, in a multi-type recurrent event process, the gap time for different event types could be correlated. For example, the occurrence of one type of event could cause other types of events to occur more frequently. In this case, a multivariate recurrent process should be considered to model the interdependence of multiple event types in the same system.
In many event analyses, covariates are found to be useful for modeling the event occurrence time and frequency. Many recurrent processes are affected by process conditions. For example, some mechanical failures could occur at a higher frequency under a higher temperature, humidity, or pressure. Incorporating covariates into the recurrent process models could improve the model performance and provide a more precise estimation of the event time and frequency. Models of this type are referred to as the covariate-adjusted recurrent event process models.
In this paper, we focus on modeling and analysis of geyser eruptions for a two geyser system in the Yellowstone National Park. In a multi-type recurrent event system, the system events can result from either type of events, and hence any consecutive events could be associated with the same or different event types. In order to describe this bivariate recurrent event process, we need to not only understand the marginal behavior of each type of event, but also understand the interdependence between the two types of events. We consider a bivariate distribution for the gap times between successive events for a bivariate event system.
To improve the estimation, we leverage the covariate information on the eruption duration and develop a covariate-adjusted bivariate recurrent event process model for estimating the eruption gap time of a two geyser system formed by West Triplet Geyser and Grotto Geyser in the Yellowstone National Park.
Recurrent event processes are extensively studied in the areas of reliability, public health, and medical studies. The nonparametric estimation of gap time distribution based on multivariate failure time data was introduced in Schaubel and Cai (2004). Dauxois and Sencey (2009) considered the risks of two nosocomial infections for patients admitted to hospitals. Bouaziz et al. (2013) provided a nonparametric method to estimate the intensity function of a recurrent process. Other than hazard functions and intensity functions, survival status in time is often of interest as well. Huang and Liu (2007) investigated the disease free survival rate in a recurrent heart failure study. Zeng and Lin (2009) and Garre et al. (2008) focused on studying the terminal events in recurrent systems. Meyer and Romeo (2015) presented Bayesian analysis of recurrent event using copulas. An earlier review regarding recurrent events can be found in Lawless (1995) and a review on modeling of repairable systems can be found in Lindqvist (2006). Classical books on recurrent event data analysis include Daley and Vere-Jones (2003), Cook and Lawless (2007), and Duchateau and Janssen (2008).
In medical research, Liu and Huang (2009) presented repeated measurements of biomarker to determine the HIV survival status in a recurrent event system. Sun et al. (2006) applied covariate-adjusted additive hazard model for the data, which is involving recurrent gap times. Prasad and Rao (2002) used a proportional hazard function with covariate adjustment in a repairable system. In another application of recurrent event data, Huzurbazar and Williams (2010) incorporated covariates in a flowgraph model. Yang et al. (2013) introduced multivariate lognormal assumption on event gap times of different event types. Yang et al. (2017) considered a parametric model for the multi-type event recurrent event data without covariates and developed copula function on gap times for the recurrent process in a car body manufacturing process.
Existing research does not consider correlated renewal process with covariate adjustments.
Motivated by the geyser data, we propose the CARP model to study the eruption gap time for a two-geyser system. Application-wise, geyser eruption is rarely studied by statistical models. For the geyser eruption study, the modeling and analysis presented in this paper are new to geyser research.
The rest of this paper is organized as follows. Section 2 provides more details on the geyser eruption data from the Yellowstone National Park. Section 3 discusses the model formulation for the multi-type recurrent event system. Section 4 describes the maximum likelihood approach for estimating the model parameters. A simulation study is described in Section 5 to evaluate the proposed method with model comparisons made under different parameter settings. The modeling and analysis of the Yellowstone geyser data are detailed in Section 6.
Section 7 contains some concluding remarks.
We use the publicly available Yellowstone geyser eruption data, which were collected in 2008 by the Geyser Observation and Study Association (GOSA). By using underground sensors, water levels were measured continuously, and occurrences of geyser eruptions were detected automatically. For each geyser, the GOSA data include the starting time and duration of each eruption. We choose to analyze the data from the West Triplet Geyser and the Grotto Geyser during the study period between June 2008 and November 2008. This particular dataset and study period were chosen to ensure that completely uninterrupted recurrent event data are available over a relatively long time span to allow for the dependence modeling for the two geysers.
We merge their eruption records in a temporal order as illustrated in Table 1. The data include the date and time when an eruption occurred, the eruption duration, and which geyser had the eruption. During the study period, the West Triplet Geyser erupted more frequently than the Grotto Geyser. Also, the eruptions of the West Triplet Geyser also lasted longer than those of the Grotto Geyser on average. More specifically, the average eruption gap time for the West Triplet Geyser is 6.8 hours with a standard deviation of 2.8 hours, while the Grotto Geyser had an average eruption gap time of 9.3 hours with a standard deviation of 8.6 hours. For the eruption duration, the Grotto Geyser eruptions lasted generally longer than the West Triplet Geyser, with the average duration of the Grotto Geyser eruptions being 3.5 hours with a standard deviation of 5.2 hours and the eruption duration of the West Triplet Geyser averaged at 0.7 hours with a 0.5 hours standard deviation.
Figures 1a and1b display the side-by-side boxplots of the time between eruptions and the duration of eruptions, respectively, for the West Triplet and Grotto geysers. We can see compared with the West Triplet Geyser, the Grotto Geyser has a much larger variation of the eruption gap time with a number of extremely long gaps between eruptions and also more variation in the eruption frequency. On the other hand, the West Triplet Geyser has considerably shorter and less variable eruption durations than the Grotto Geyser. In addition, Figure 2 shows the plots of the time between eruptions versus the duration time of the previous eruption for both geysers. We can observe a moderate correlation between these two variables for the West Triplet Geyser and a strong correlation for the Grotto geyser.
Hence, we decided to utilize the duration of the previous eruption for both geysers to help model the eruption gap time.
Figure 3 illustrates the data obtained for the two-geyser system. In this figure, W kj denotes the gap time of the kth eruption (i.e., the time interval between the (k 1)th and kth eruptions) for geyser j where j = 1 for the West Triplet Geyser and j = 2 for the Grotto Geyser. The covariate x i denotes the eruption duration for the ith eruption in the bivariate geyser system, where i = 1, × × ×, n, n = n 1 + n 2 is the total number of eruptions for both geysers, and n j denotes the number of eruptions for the jth geyser. In particular, n 1 = 580, n 2 = 421, and n = 1001 for the dataset we analyzed. Note here the time between eruptions are labeled separately for each individual geyser. To model the bivariate geyser system, we will introduce new notation for jointly describing the event time, event type and covariate information in Section 3.
Suppose that in a bivariate recurrent process with n total events, the systematic event time is described by the variable T i for i = 1, × × ×, n. We use T 0 to denote the starting time or the system installation time, and T n denotes the last event time in the system. In a bivariate system, there are two types of events, and hence we use an indicator variable Δ i }1, 2| to represent the type of event. For an event that occurs at time t i , the covariate vector on the event duration is denoted by X i , which is a vector of the previous duration for both types of events. Therefore, each event can be represented by }T i , Δ i , X i | , where i = 1, × × ×, n. We use }t i , δ i , x i | to denote the observations of the ith event which occurs at time t i and it is from event type δ i with covariates measured as x i . Note in an observed recurrent process with n total events, the last event is observed at time t n , where t n ≥ κ with κ being the pre-defined study termination time.
This sequence of bivariate events can be expressed as a counting process }N (t) : t ∼ 0| , where N (t) denotes the cumulative number of events at time t regardless of the event type.
Similarly, we define the counting process for each individual event type j as }N j (t) : t ∼ 0| , where j = 1 or 2. In addition, we denote the event history of the system up to a time point s ≥ κ by X s = }N (t) : t ≥ s| . Similarly, the covariate history can be expressed as { s = }x t : t ≥ s| . For further discussion, we use H s = }X s , { s | to denote the history including both event and covariate information.
For a bivariate recurrent process system, event time variable T i defined above satisfies
Similarly, for events of type j, the event time variables are defined as T lj , where l = 1, . . . , n j and n j is the number of events from type j. In a bivariate recurrent system, we have the relation n 1 + n 2 = n. The time variable T lj is also presented in temporal order. For the events of type j, j = 1 or 2, we have
Based on the ordered event times, the event gap time for the jth type as the interval between two consecutive events can be calculated as W lj = T l+1,j T lj , where l = 0, 1, . . . , n j 1. Now consider the joint bivariate recurrent system. For the ith event in the system, the two-dimensional event gap time variable is defined as
where l ij = N j (t i ) is the cumulative number of events for jth event type by time t i . For example, at the system installation time where t = 0, the event gap time vector is W 0 = (W 01 , W 02 ) . If the first event occurs at time t 1 and is from event type 1, then the event gap time vector is denoted as W 1 = (W 11 , W 02 ) . In this case, l 11 = 1 since one event of type 1 has occurred as of time t 1 , while l 12 = 0 as no event of type 2 has occurred by that time point.
To link the event gap time W i with the system event time T i , we introduce an age variable at time t which is defined as the time between the time point t and the time of the latest event of each event type. In a bivariate recurrent process, the age variable is a vector of two components denoted as
where A j (t) = t T N j (t), j for j = 1 and 2. For an observed event at time t i , the observed age vector is denoted as
, where a j (t) = t t N j (t), j . For example, at the system installation time at t = 0, the age vector is a 0 = (0, 0) . If the first event occurs from type 1 at time t 1 , then a 1 = (0, t 1 ) . The age variable then connects between the system event time and the event gap time W i through the relation W i ∼ A i . Here, the vector comparison is defined to be an element-wise comparison. In other words, W i ∼ A i suggests W l ij ,j ∼ A j (T i ) for both j = 1 and 2.
In this section, we introduce the proposed covariate-adjusted recurrent process (CARP).
After an event occurred at time t i (or after the system installation at time t 0 = 0), the distribution of the event gap time is given by
(1)
is the joint cumulative distribution function (cdf) of the event gap time conditioned on the age and covariates. At any event time, the information we have about the event type that
has not yet occurred is captured through its age and the conditional probability that is conditioned on the event eruption time is greater than or equal to the age since the last eruption. At each event time, the age is set to be zero for the occurred event type. The event gap time variable W i is adjusted by covariates x i which is further discussed in Section 3.3 in more detail.
Dependence between events from different types in a bivariate system is modeled by implementing distributional assumptions on variable W i = (W i1 , W i2 ) . In this paper, we use a bivariate lognormal distribution and a copula function to model the random vector W i , where in both models, covariates x i are used for adjustment. We refer to the CARP models under these assumptions as the CARP-MLN and CARP-copula models, respectively. Here MLN is short for multivariate lognormal.
For the multivariate lognormal distribution,
where the location parameter in the bivariate lognormal distribution is expressed as a linear form of covariates x i . That is,
In the linear expression above, μ 0 is a vector of baseline location parameters and B is a 2 ≤ 2 coefficient matrix. In the bivariate lognormal assumption, we use a covariance matrix Σ to capture the event dependence between events from two event types. Specifically, the diagonal elements in Σ represent marginal variances while the off diagonal elements stand for the covariances. When using the bivariate lognormal distribution, the covariance matrix is defined as Σ = CC to ensure Σ to be positive definite, where
The correlation ρ is determined by σ 2 and η as
Note the above covariance matrix offers great flexibility to model different correlation relationships of varied size and direction. The sign of the η value determines if the two types of events have a positive or negative correlation. In a bivariate lognormal distribution, the marginal distribution of each dimension also follows a lognormal distribution. Therefore, when the observed marginal distributions do not seem to follow the lognormal distributions, or their dependency cannot be characterized by the covariance matrix Σ, the assumed bivariate lognormal distribution is not appropriate. In this case, the alternative strategy is to define W 1 and W 2 by separate distributions and combine them through a more flexible copula function. This is referred to as the CARP-copula model, which will be introduced in the next section.
Let F 1 and F 2 denote the marginal cdfs of W 1 and W 2 . There always exists a copula function C such that the joint cdf of the two dimensional variable (W 1 , W 2 ) ′ can be written as
where for any unitary uniform variable U j , j = 1, 2, the bivariate copula is defined as
This is also known as the Sklar's Theorem. By using a copula function, we have the flexibility to choose marginal distributions separately for each event type. For instance, a Gamma distribution and a Weibull distribution can be used as marginal distributions for the two types of events, respectively. With selected F 1 and F 2 , one can combine the marginal distributions by using different copula functions.
Note the copula approach allows us to model the marginal distributions separately and then combine them through an appropriate copula function for modeling the dependence structure. For the geyser system, we choose to use parametric distributions for modeling the marginal distributions. However, the method can be easily generalized to using nonparametric methods for modeling the marginal distributions and hence offers great flexibility to be adapted for broad applications.
In the literature, a variety of copula functions has been introduced to capture different dependence patterns among the marginal distributions. A Gaussian copula uses a multivariate normal distribution of transformed marginal distributions based on the inverse cdf of the standard normal distribution. The Archimedean copulas are an associative class of copula functions that model the multivariate dependence through a single parameter. The power variance function copulas including Clayton, Gumbel and Inverse Gaussian are among the most popular ones that are flexible for modeling various dependence structures (e.g., Romeo et al., 2018). In this paper, we use the Gumbel copula function which is popular for modeling stronger dependence in the positive tail. We refer to the CARP model with the use of Gumbel copula as CARP-copula model for the rest of the paper. In fact, the CARP-MLN model is a special case of the CARP-copula model, where in CARP-MLN, the Gaussian copula is applied and lognormal marginal distributions are selected.
The CARP-MLN model characterizes the dependence among event types through the covariance matrix Σ, while the CARP-copula model quantifies dependence using the copula parameter. Particularly, in the Gumbel copula model, the parameter is denoted as α. As a result, in the CARP-MLN and CARP-copula models, we have different parameters to characterize event dependence. In order to compare dependence from different models, we introduce the Kendall's tau.
In the CARP-copula model, the Kendall's tau is expressed as τ = 1 1/α, where α is the copula coefficient in the Gumbel copula. In the CARP-MLN model, the Kendall's tau is calculated as τ = (2/π) arcsin(η/ σ 2 2 + η 2 ), where η and σ 2 can be found in (4). Similar to the CARP-MLN model, for the CARP-copula model, we also use a linear form of the covariates x i to represent location parameters in marginal distributions as in (3). In the CARP-MLN model, we use a two dimensional vector μ(
the location parameter of the lognormal distribution. While in the CARP-copula, μ 1 (x i ) and μ 2 (x i ) stand for the location parameters for the first and second marginal distributions, respectively.
For event gap time variable W i , we define the survival function (sf), cdf and hazard function as follows. We denote the joint sf of W i as,
The joint cdf is
and the corresponding joint probability density function (pdf) is denoted by f W (v). According to (2), the joint pdf of the event gap time variable given all historical events W i H t i is given by
where c j = t l ij ,j + W j , j = 1, 2. The denominator in (7) takes age condition into account, while the numerator builds the relationship among the event gap time W j , the event time t l ij ,j and the age a j (c j ). Covariates x i are used to adjust the gap time W i as discussed in Section 3.3.
For further discussions, let H t -be the event history up to time t. Note that for event type j, T N j (t -),j gives the most recent event time by time t, and the age variable prior to time t is denoted as A - j (t) = t T N j (t -),j , which calculates the cumulative running time upon time t since the last event. We use a -(t) = [a - 1 (t), a - 2 (t)] to denote the age vector for the two event processes prior to time t. Hence, prior to an event time t i , the age vector is denoted as a - i = a -(t i ). For example, at the initial time 0, the vector is a - 0 = (0, 0) . If the event type is δ 1 = 1 at t 1 , then a - 1 = (t 1 , t 1 ) . Note that a 1 updates the age to be zero for the corresponding event type at an event time, while a - 1 does not. This notation is used to derive the likelihood in Section 4, is also used to define the hazard function below.
In the literature, the sub-intensity function (i.e., cause-specific event intensity function) is often used to characterize an event process. In particular, the sub-intensity function for event type j is defined as
where T is the event time, and Δ is the event type. The sub-cumulative intensity function is H j (t) = t 0 h j (s)ds. The hazard function and cumulative hazard function for the system are calculated as the sum of corresponding functions for the two event types,
The sub-intensity function in ( 8) is calculated as
and v = ( v 1 , v 2 ) is a vector with two components. More detail about the calculation of D j ( v) under different models will be discussed in Section 4.
The maximum likelihood (ML) approach is used to estimate model parameters. Parameters in the model include parameters in the joint distribution function, the copula function and the linear covariate transformation function. Given all the event history H τ , the likelihood function is constructed as follows:
where
Let θ denote the vector of all the parameters in the model. The estimated parameters θ are asymptotically normally distributed based on the large sample ML theory (Casella and Berger, 2002). The calculation of the likelihood L i for the proposed CARP models is shown below.
For any observed recurrent process with n total events, the likelihood contribution for i = 1, × × ×, n in (11) is given by
In ( 12), the quantity D δ i (a - i ), which is introduced in (9), is the partial derivative of the bivariate sf S(a - i ). The covariates x i are used to adjust the distribution of W i . Likelihood calculations so far are the same for both the CARP-MLN and CARP-copula models. However, we need different ways to calculate D j (a - i ) for the two CARP models based on how the survival functions are calculated, which are detailed below.
For the CARP-MLN model, the sf is calculated in a closed form as discussed in ( 6), and the partial derivative with regard to the jth event type can be written as
where f j [a - j (t i )] is the jth marginal density function from a bivariate lognormal distribution. In ( 13), the conditional probability can be calculated from the conditional normal distribu-tion with a logarithm transformation. Calculation details can be found in Appendix A.1.
For the CARP-copula model, the likelihood function is calculated based on the relationship between the bivariate sf S(a - i ) and the cdf F (a - i ). In bivariate cases, the sf and cdf have the following relationship:
where the joint cdf can be calculated by the copula as in ( 5). The calculated likelihood will vary for different choices of the marginal distribution and the copula function. In our case, we chose the Gumbel copula as described in Section 3.3. The ML estimates θ are obtained by maximizing the likelihood function in (10).
This section describes the simulation study we conducted for evaluating the proposed method.
We simulated bivariate recurrent system data based on different models and parameter values.
The performance of the proposed method was then evaluated based on the simulated data.
Different parameter values were considered in the simulation for both the CARP-MLN and the CARP-copula models to understand the impact of the model parameters. We evaluated the model goodness of fit by calculating the average AIC and the mean squared error (MSE) of the estimated model parameters. The goal is to demonstrate the performance of the proposed method and the improvement by using the covariate adjustment.
We used both the CARP-MLN and CARP-copula models to simulate the data. For each generated data set, both models were fitted to the data and results were evaluated. The model used to generate data is referred to as the true model, while the models used to fit the data are referred to as the fitted models. We explored changing the sample sizes (n), the Kendall's tau value, the scale parameters (σ j , j = 1, 2), and the linear transformation matrix (B) to understand their impacts on the analysis. First, to understand the impact of sample size, we varied the sample size at n =200, 500, 1000 and 2000 when generating the data using each model. For the CARP-copula model, the lognormal marginal distributions were used for both event types along with the Gumbel copula function. The location and scale parameters for lognormal marginal distributions in the true model were specified at (μ 1 = 1, σ 1 = 0.25) and (μ 2 = 1.5, σ 2 = 0.25). For the linear transformation matrix B, we used the 2 × 2 matrix below
We chose to set the Gumbel copula parameter α at 1.5. When generating the data based on the CARP-MLN model, we used the same location and scale parameters for the marginal lognormal distributions, and adjusted the correlation parameter η to obtain the same Kendall's tau as in the CARP-copula model.
Second, to evaluate the effect of the Gumbel copula coefficient α, we varied its value at 1, 1.12, 1.5 and 2.22 by changing the Kendall's tau parameter in the true models, while keeping all other parameters the same with the sample size fixed at n = 1000. Considering that the CARP-MLN and CARP-copula models quantify the Kendall's tau differently, we used η at 0, 0.0443, 0.1445 and 0.299 for the CARP-MLN model, so that the Kendall's tau from both
CARPs are varied at 0, 0.11, 0.33 and 0.55.
We also explored the impact of the scale parameters by letting σ 1 and σ 2 vary at 0.35, 0.3, 0.25 to 0.20 for the marginal lognormal distributions of the CARP-Copula model while keeping sample size at n = 1000 and the Kendall's tau was kept at 0.33. Similarly for the CARP-MLN model, covariance matrices were adjusted to align with the marginal distributions used in the CARP-Copula model.
Lastly, the effect of covariate adjustment was studied by using different types of B. True models used both non-zero and zero B matrices. In the non-zero B case, we used the coefficient matrix in ( 14) to generate data. In the zero B case, B = 0 was used. The Kendall's tau was set to be 0.33, while other parameters are the same as ones in the sample size case. A sample size of 1000 is used across true models.
The simulation results under different true models and fitted models are summarized in this section. For each true model, the model performance was summarized over 1000 simulated data sets. The simulation size was chosen to ensure a reliable analysis result while balancing the computing time needed to evaluate a broad number of scenarios. Further increasing the simulation size would not result in a meaningful change in the evaluated summary statistics and the drawn conclusions. Under each fitted model, the average AIC was computed. In order to evaluate the performance of parameter estimates, we calculated the MSE of the location parameter µ, scale parameter σ and the linear transformation parameter B for Table 2 shows the average AIC of the fitted models summarized over the 1000 simulated data sets generated using both the CARP-copula and CARP-MLN models at different sample sizes. Both the CARP-MLN and CARP-copula models were applied to each simulated data set. We can see at each fixed sample size, the fitted CARP model that matches the true model used for simulating the data generally outperforms the other model by having a smaller average AIC value, except for the smallest size case at n = 200 where the difference in the average AIC values are extremely small. The advantage of using the true model becomes Table 3 shows the average AIC values of the fitted models when the Kendall's tau for capturing the dependence between two event types varies at the levels of 0, 0.11, 0.33 and 0.55. When the Kendall's tau is not zero which suggests some level of dependence between two marginal variables, using a fitted model that matches the true model for generating the data will result in a smaller average AIC value with a bigger improvement achieved when stronger dependence exists between the two types of events. When the Kendall's tau is zero indicating the marginal distributions are independent, the average AIC values are generally similar regardless which model is used to fit the data.
Table 4 shows the comparison based on using different scale parameters in the true underlying model for both the CARP-copula and CARP-MLN models. Again, across all evaluated scale parameter values, we observe consistently smaller average AIC values when the fitted model matches the true model used for data generation.
Lastly, Table 5 compares the results of using different coefficient matrix B. We simulated data with both non-zero and zero B using both the CARP-MLN and CARP-copula models.
For each simulated data, we fitted both the CARP-MLN and CARP-copula models with non-zero and zero B. Similar patterns can be observed. When the true model was used as the fitted model, it results in the smallest AIC value. When non-zero B was used to generate the data, the use of covariate adjustment led to substantial improvement in the average AIC of the fitted model compared to using zero B. On the other hand, when the data were generated with zero B, using the non-zero B fitted models produced similar AIC values as the zero B models. Therefore, the covariate adjustment is generally recommended due to its potential to substantially improving the model performance by leveraging the additional covariate information.
In this section, we present the analysis for the bivariate geyser system in the Yellowstone National Park using the proposed CARP models. Two adjacent geysers including the West Triplet and the Grotto Geyser are considered for our analysis. Geyser eruptions are highly related to underground water levels, which can be affected by a nearby geyser eruption. Also, it is believed that the eruption duration could affect the gap time until the next eruption. In particular, the longer the current eruption lasts, the longer it will take for the next eruption to occur. This is because a longer eruption usually indicates more water consumption during the eruption and hence a longer water gathering time is expected to reach the next eruption.
We use both the CARP-MLN and the CARP-copula models with lognormal marginal dis- where the correlation estimation is calculated as η/ √ σ 2 2 + η 2 = 0.05. This indicates a small negative correlation between the event gap times of the West Triplet and Grotto Geysers. In other words, a longer eruption gap time for West Triplet Geyser could be associated with a shorter time interval for the next eruption of the Grotto Geyser. The Kendall's tau provides We can see a better agreement between the estimated and the observed cumulative intensity functions based on the CARP-copula model compared with the CARP-MLN model.
This paper introduces two CARP models, CARP-MLN and CARP-Copula, for modeling a bivariate recurrent process. With the covariate adjustment, the CARP models provide improvements over traditional models for capturing the interdependence of the bivariate event system. When the underlying data are close to a bivariate lognormal distribution, both models work similarly well. However, when the real data are inconsistent with a bivariate lognormal distribution, the CARP-copula model is recommended to improve the model performance.
For lognormal cases, the Kendall's tau estimator can be expressed as
where the asymptotic distribution is known from the ML estimator. Using the Delta method,
where
For the Gumbel copula case, τ = 1 1 α .