Pricing perpetual timer options under Heston Model by finite difference method: Theory and implementation

1.   Introduction

1.1.   Joint progressive type-II censoring scheme

Applying censoring techniques to a single population entails specific limitations. While progressive type-II censoring permits the exclusion of specific data, obtaining a sufficient number of observations remains a costly endeavor. Moreover, if our emphasis is on understanding the interactions and interdependencies among populations, experiments conducted solely on a single population may not provide conclusive evidence.

The joint progressive type-II censoring scheme (JPT-II-CS) provides notable benefits for comparing the lifespan distributions of products produced by different units within the same facility. JPT-II-CS has attracted considerable attention in the research community, with numerous authors exploring JPT-II-CS and related inference methods in the literature. For instance, Rasouli and Balakrishnan ^[32], Doostparast et al. ^[18], Balakrishnan et al. ^[15], Mondal and Kundu ^[26], Krishna and Goel ^[24], Goel and Krishna ^[21], and Goel and Krishna ^[22] have contributed to this body of work.

Recently, a multitude of researchers have investigated a range of strategies and different lifetime models in various applications. For more detailed information, please consult the publications by Pandey and Srivastava ^[29], Qiao and Gui ^[31], Abdel-Aty et al. ^[1], Ferreira and Silva ^[20], Celik and Guloksuz ^[16], Chiou and Chen ^[17], Panahi and Lone ^[28], Yan et al. ^[34], Asadi et al. ^{[4,7,8,9,10,11,12]}.

Within the JPT-II-CS framework, two samples, one from Population-A (Pop-A) and the other from Population-B (Pop-B), each comprising $m$ and $n$ units, are amalgamated for a life-testing experiment. Let $k$ be the total number of observed failures in this experiment. Moreover, let $R_1, R_2, ..., R_k$ signify the number of units removed, such that $\sum_{i = 1}^{k} (R_i + 1) = m + n$, where $R_i = S_i + T_i$, and $R_i$ is the sum of units removed from Pop-A $(S_i)$ and Pop-B $(T_i)$ at the i-th stage. Considering the joint sample upon the first failure event denoted as $W_1$, $R_1 = S_1 + T_1$ units are randomly chosen from the remaining pool of $m + n - 1$ surviving units. Here $S_1$ and $T_1$ represent the units removed from Population (A) and Population (B) respectively. In a similar fashion, at the second stage, $R_2 = S_2 + T_2$ units are randomly selected from the remaining pool of $m + n - 2 - R_1$ surviving units and so forth. Finally, at the $k$-th failure, all the remaining $R_k = n + m - k - \sum_{i = 1}^{k} R_i$ surviving units are removed. It is crucial to highlight that when $R_1 = R_2 = ... = R_k = 0$, it implies that $S_i = T_i = 0$ for all $i = 1, 2, ..., k$, simplifying the JPT-II-CS to the scenario of complete samples. Furthermore, if $R_1 = R_2 = ... = R_{k-1} = 0$, resulting in $R_k = n + m - k$, then the censoring scheme transforms back to a conventional Type-II censoring scheme for two distinct samples.

In the context of the JPT-II-CS framework, the observed data comprises $(W, Z, S)$, where $W$ is depicted as $(W_1, W_2, ..., W_k)$, where $k$ is a specified integer such that $1 \leqslant k < m + n.$ $Z$ is denoted as $(Z_1, Z_2, ..., Z_k)$, where

Moreover, $S$ is depicted as $(S_1, S_2, ..., S_k)$. We employ $k_1$ to indicate the number of failed units from Pop-A, and $k_2$ to signify the count of failed units from Pop-B, where $k_1 = \sum_{i = 1}^{k} Z_i$ and $k_2 = k - \sum_{i = 1}^{k} Z_i$, respectively.

1.2.   Nadarajah-Haghighi distribution

The exponential distribution is renowned for its constant failure rate and memoryless property. However, when exploring phenomena in lifetime and reliability studies, opting for the exponential model may not be conducive, as it lacks the ability to exhibit both monotone (increasing and decreasing) and non-monotone (bathtub and upside-down bathtub) failure rate behaviors. To address this limitation and provide flexibility, alternative generalizations of the exponential distribution have been proposed, such as the Nadarajah–Haghighi distribution (NHD). Introduced by Nadarajah and Haghighi ^[27], this distribution serves as a widespread statistical extension of the conventional exponential distribution and has been recently named the NHD, abbreviating the authors' names. The cumulative distribution function (CDF) of the NHD is expressed as follows (see Figure 1):

Additionally, its corresponding probability density function (PDF) is given by:

In this context, $\gamma$ and $\alpha$ denote the shape and scale parameters, respectively. By setting $\gamma = 1$ for the variable in Eq (1.1), the inverted exponential distribution is presented as a particular case. Nadarajah and Haghighi ^[27] demonstrated that the density of the NHD may exhibit a decreasing trend. Furthermore, its shapes, both unimodal and the hazard rate function, can resemble the increasing, decreasing, or constant patterns observed in gamma, Weibull, and generalized-exponential distributions.

Several researchers have delved into the estimation challenges associated with the NHD. For instance, Selim ^[33] investigated the estimation and prediction of the NHD based on record values. Inferences and ideal censoring techniques for the progressively first-failure censored NHD were covered by Ashour et al. ^[5]. Elshahhat et al. ^[19] used type-II adaptive progressive hybrid censoring with applications to investigate inferences for Nadarajah-Haghighighi parameters. The half logistic inverted NHD was studied by Alotaibi et al. ^[2] under ranked set Sampling with applications. With applications to COVID-19 mortality data and cancer data, Azimi and Esmailian ^[6] examined a novel generalisation of the NHD.

In this publication, we tackle the challenges associated with estimating the NHD. We present both Bayesian and maximum likelihood estimators (MLEs) for the NHD using the JPT-II-CS, accompanied by their corresponding confidence intervals. We derive Bayes estimators for the squared error (SE) and linear exponential (LINEX) loss functions, assuming independent gamma priors. Monte Carlo simulations are conducted to assess the performance of various estimators, evaluating them based on mean squared error (MSE) and average values. Additionally, we scrutinize the average confidence lengths of the 95% two-sided interval estimations to gauge their effectiveness. Finally, we illustrate the practical application of our approach using real-world datasets.

The following outlines the order of the remaining sections in the paper following this introduction: Section 2 examines the concept of sampling significance. Section 3 explains the process of obtaining maximum probability estimates for the unknown parameters. The construction of the Fisher information matrix (FIM) using MLEs is done in Section 4. In Section 5, Bayes estimation is performed with various loss functions using gamma and non-informative priors. Section 6 presents a Monte Carlo simulation analysis along with the results, serving as model validation. Two real data analyses to demonstrate each of the suggested methodologies are carried out in Section 7. Section 8 provides a concise summary of the main points.

2.   Model description

In a particular context, we possess independent and identically distributed (iid) lifetimes, denoted as $X_1, X_2, ..., X_m$, originating from Pop-A. These lifetimes adhere to a NHD characterized by the CDF F(x) and the PDF f(x). In the same vein, we possess $Y_1, Y_2, ..., Y_n$ as iid lifetimes from Pop-B, also adhering to a NHD defined by the CDF G(y) and the PDF g(y). Given a JPT-II-CS ${(m, n); R_1, R_2, ..., R_k}$, we can define $B = (W, Z, S)$ as the JPT-II-CS of size $k$ obtained from Pop-A and Pop-B. It comprises ${(w_1, z_1, s_1), (w_2, z_2, s_2), ..., (w_k, z_k, s_k)}$. The likelihood function based on the observed JPT-II-CS can be expressed as follows:

where $w_1\leq w_2 \leq \ldots \leq w_k$, $\bar{F} = 1-F$, $\bar{G} = 1-G$, $\sum_{i = 1}^k s_i+\sum_{i = 1}^k t_i = \sum_{i = 1}^k R_i$ \text{and} $C = B_1 B_2$ with

3.   Maximum likelihood estimation

MLE stands as a powerful statistical method widely employed in the field of data analysis and parameter estimation. In this manuscript, we employ MLE to estimate the parameters of the NHD. The main principle of MLE is to find the parameter values that maximize the likelihood function, which quantifies the probability of observing the given data under a specific set of parameter values. By maximizing the likelihood, we aim to identify the parameter values that make the observed data most plausible. The MLE estimates are obtained by solving an optimization problem, either analytically or numerically, to find the parameter values that maximize the likelihood function. These estimates are known to possess desirable properties, such as efficiency and consistency, which means that they converge to the true parameter values as the sample size increases. MLE is a versatile and widely applicable method that plays a crucial role in statistical inference and parameter estimation.

The following outcome is obtained by applying the CDF and PDF from Eqs (1.1) and (1.2) to the likelihood equation given in (2.1):

Given that the log-likelihood function exhibits the same monotonic behavior as the likelihood function, it is expressed as follows:

By partially differentiating Eq (3.2) with respect to $\alpha_1, \alpha_2, \gamma_1 \quad \text{and} \quad \gamma_2$ and setting the derivatives equal to zero, the resulting equations are as follows.

It's evident, as indicated in Eqs (3.3)–(3.6), that explicit expressions for the MLEs $\alpha_1, \alpha_2, \gamma_1,$ and $\gamma_2$ are not available. Hence, we recommend using numerical iterative methods, like the Newton-Raphson procedure, to derive the values of $\alpha_1, \alpha_2, \gamma_1,$ and $\gamma_2$.

4.   Fisher information matrix

We investigate approximate confidence intervals in this context for the unknown parameters $(\alpha_1, \alpha_2, \gamma_1, \gamma_2)$, using large sample approximations of the MLEs, which are often known as asymptotic theory. We do this by using the observed FIM to calculate the MLEs' asymptotic variance for the unknown parameters. Next, the observed FIM is represented as:

Derived from the log-likelihood function in (3.2), the following is evident:

Calculating the asymptotic confidence intervals (ACIs) for the parameters $\alpha_1, \alpha_2, \gamma_1,$ and $\gamma_2$ is feasible by leveraging the asymptotic normality of the maximum likelihood estimators (MLEs). This enables us to estimate the $(1-\mu)100\%$ ACIs for $\alpha_1, \alpha_2, \gamma_1,$ and $\gamma_2$ in an approximate manner

Here, $Z_{\mu/2}$ represents the upper $\mu/2$th percentile of the standard normal distribution.

5.   Bayesian inference

In the field of reliability analysis, there are instances where classical estimation using the MLE approach may encounter challenges, particularly when the available data lacks sufficient sampling details. To address this issue, incorporating prior information in conjunction with Bayesian analysis proves beneficial. This section focuses on the Bayesian approach for estimating unknown parameters and deriving the corresponding credible intervals (CRIs).

5.1.   Prior distribution

In Bayesian statistical inference, the role of the prior distribution is pivotal as it represents our existing knowledge or beliefs regarding the parameters, facilitating a more accurate estimation of the posterior distribution. The selection of an appropriate prior distribution is crucial as it can influence the ultimate results of the inference. The gamma distribution, recognized for its flexibility and favorable properties, is a frequently chosen continuous probability distribution for Bayesian parameter priors. The parameters of the gamma distribution can be adjusted to accommodate diverse prior beliefs. Additionally, the gamma distribution exhibits conjugacy, signifying that when employed as a prior distribution, its product with the likelihood function remains a gamma distribution, simplifying posterior distribution calculations.

The choice of hyperparameter values, like any other parameter selection in Bayesian analysis, depends on various factors including prior knowledge, the characteristics of the data, and modeling considerations. Without specific details about the context of your simulation, I can provide some general considerations:

● Prior Knowledge: If you have prior information about the hyperparameters, such as from previous studies or domain expertise, you might set them based on that knowledge. The values could reflect your beliefs about the likely range or distribution of the true hyperparameters.

● Empirical Bayes: If you don't have strong prior knowledge, you might consider using an empirical Bayes approach. This involves estimating the hyperparameters from the data itself. In such cases, the hyperparameters are determined based on the observed data distribution.

● Conjugate Priors: If you choose hyperparameters that make your prior distribution conjugate to the likelihood, it can simplify computations. However, this choice might not always reflect your beliefs about the parameters.

● Non-Informative Priors: A common choice for hyperparameters, especially in the absence of strong prior information, is to set them to values that make the prior distribution non-informative or weakly informative. For instance, setting hyperparameters to achieve a flat or diffuse prior.

In this study, we assert that $c_i$ and $d_i$ are both greater than $0$, where $i$ takes values 1, 2, 3, and 4. This results in the following relationships:

In this scenario, $c_i$ and $d_i$ where $i = 1, 2, 3, 4,$ are introduced to integrate prior information concerning the unidentified parameters. Consequently, the joint prior density function for $\alpha_1$, $\alpha_2$, $\gamma_1$, and $\gamma_2$ is structured as follows:

5.2.   Posterior distribution

The posterior distribution is a fundamental concept in Bayesian statistics, playing a crucial role in the inference process. It represents the updated belief or knowledge about a parameter or set of parameters after taking into account both prior information and observed data. In Bayesian analysis, the posterior distribution is obtained by combining the prior distribution, which encapsulates our initial beliefs, with the likelihood function, which quantifies the probability of observing the data given the parameter values. The posterior distribution provides a comprehensive summary of uncertainty, allowing us to estimate parameters, make predictions, and conduct various analyses. It serves as a bridge between prior knowledge and observed evidence, enabling us to make informed decisions and update our understanding of the underlying phenomenon.

By merging Eqs (3.1) and (5.1), it becomes viable to represent the joint posterior density function for $\alpha_1$, $\alpha_2$, $\gamma_1$ and $\gamma_2$ as following

The joint posterior density function for $\alpha_1$, $\alpha_2$, $\gamma_1$ and $\gamma_2$ can be formulated as:

Evidently, owing to the nonlinear form of (5.3), there is no closed-form solution for the Bayes estimators of $\alpha_1$, $\alpha_2$, $\gamma_1$, or $\gamma_2$ when employing the SE and LINEX loss functions. Consequently, we recommend employing the Markov Chain Monte Carlo (MCMC) approach to acquire the Bayes estimates and establish the corresponding CRIs.

5.3.   MCMC approach

To generate samples using the MCMC approach, we must first determine the conditional posterior distributions of the unknown NHD parameters $\alpha_1$, $\alpha_2$, $\gamma_1$ and $\gamma_2$. The conditions are given by

From Eqs (5.6) and (5.7), it's apparent that the posterior densities of $\gamma_1$ and $\gamma_2$ follow the gamma distribution, given that

and

It's evident that gamma densities can be used to generate random samples of $\gamma_1$ and $\gamma_2$. However, it's not feasible to algebraically simplify the density functions of $\pi_1^*(\alpha_1|\alpha_2, \gamma_1, \gamma_2)$ and $\pi_2^*(\alpha_2|\alpha_1, \gamma_1, \gamma_2)$ to well-known distributions. Therefore, obtaining samples directly through conventional procedures is not practical. Hence, we turn to the Markov Chain Monte Carlo (MCMC) approach and generate a sample using Gibbs sampling with the Metropolis-Hastings (M-H) algorithm (see Metropolis et al. ^[25] and Hastings ^[23]), employing a standard proposal, as described below.

● Step 1. Begin by employing the initial values $\big(\alpha_1^{(0)}, \alpha_2^{(0)}, \gamma_1^{(0)}, \gamma_2^{(0)}\big)$.

● Step 2. Set $j = 1$.

● Step 3. Generate $\gamma_1^{(j)}$ from Gamma $\Big(k_1+c_3, d_3+\sum_{i = 1}^{k} z_i \ln (1+\alpha_1 w_i) \Big)$.

● Step 4. Generate $\gamma_2^{(j)}$ from Gamma $\Big(k_2+c_4, d_4+\sum_{i = 1}^{k} (1-z_i) \ln (1+\alpha_2 w_i) \Big)$.

● Step 5. Equations (5.4) and (5.5) can be employed to generate $\alpha_1^{(j)}$ and $\alpha_2^{(j)}$ using the Metropolis-Hastings (M-H) algorithm. The recommended normal distributions to use are $N(\alpha_1^{(j-1)}, \text{var}(\alpha_1))$ and $N(\alpha_2^{(j-1)}, \text{var}(\alpha_2))$. In this instance, the major diagonal of the inverted FIM can be utilized to calculate $\text{var}(\alpha_1)$ and $\text{var}(\alpha_2)$.

(I) Produce proposed values $\alpha_1^*$ and $\alpha_2^*$ from the corresponding normal distributions.\\ (II) Determine the acceptance probabilities using the following procedure:

(III) Produce a random value $u$ from a uniform distribution within the range of $(0, 1)$.

(IV) If $u \leq r_1$, approve the proposal and assign $\alpha_1^{(j)}$ as $\alpha_1^*$; otherwise, maintain $\alpha_1^{(j)}$ as $\alpha_1^{(j-1)}$.

(V) If $u \leq r_2$, approve the proposal and assign $\alpha_2^{(j)}$ as $\alpha_2^*$; otherwise, maintain $\alpha_2^{(j)}$ as $\alpha_2^{(j-1)}$.

● Step 7. Set j = j + 1.

● Step 8. Repeat 2–7, for V times. Therefore, the estimated posterior means of $(\alpha_1, \alpha_2, \gamma_1, \gamma_2)$ represented by $\lambda$ under the squared error loss function, can be determined as:

Lastly, calculate the Bayesian estimates of $\lambda$ utilizing the LINEX loss function:

Here, $N$ denotes the burn-in period.

6.   Numerical study

This section's goal is to assess the potency of the various estimation techniques covered in earlier sections. We demonstrate this by examining an actual dataset and conducting a simulated experiment to evaluate the statistical performance of the estimators under the JPT-II-CS. The computations were carried out using Mathematica ver. 10 software.

6.1.   Simulation study

In this section, we conduct simulation studies to assess the performance of the estimation methods developed in the preceding sections. We explore various sample sizes for the two populations, including $(m, n)$ = $(10, 20), (20, 30), (40, 50)$, and different numbers of failures for each sample size, such as $(15, 20, 30), (35, 40, 50), (60, 70, 90)$, respectively. The parameter values for the two populations are set as $(\alpha_1, \alpha_2, \gamma_1, \gamma_2)$ = $(1.5, 1.3, 0.4, 0.3)$.

We computed MLEs along with 95% CIs for the parameters $(\alpha_1, \alpha_2, \gamma_1, \gamma_2)$ across all specified scenarios. This process was repeated $1000$ times, and we calculated the mean values of MLEs and their respective lengths. The results are presented in Tables 1–4. Additionally, we employed informative gamma priors for for $\alpha_1$, $\alpha_2$, $\gamma_1$, and $\gamma_2$ in the context of Bayesian estimation under SE and LINEX loss functions. The hyperparameters were set as $c_i = 0.8$ and $d_i = 2.5$, where $i = 1, 2, 3, 4$, with $a = 3$ denoting overestimation and $a = -3$ representing underestimation. We employed the(MCMC approach with $11, 000$ samples to derive Bayesian estimates for $\alpha_1$, $\alpha_2$, $\gamma_1$, and $\gamma_2$, along with 95% CRIs, through $1000$ simulations. The initial $1000$ values were excluded due to "burn-in". MSE was utilized to evaluate the performance of the generated estimators for $\alpha_1$, $\alpha_2$, $\gamma_1$, and $\gamma_2$. Following $1000$ repetitions of this process, we computed the mean values of MLEs and their respective lengths. Tables 1–4 present the results.

Based on the aforementioned data, several conclusions can be drawn:

(1) Tables 1 to 4 reveal that, in most cases, Bayesian estimates outperform MLEs in terms of MSEs.

(2) Examining Tables 1 to 4, it becomes apparent that CRIs exhibit the shortest average length of intervals compared to the average length of CIs, indicating the superiority of Bayesian estimators over MLEs.

(3) Notably, Bayesian estimation using the LINEX loss function at $a = 3$ surpasses Bayesian estimators using the LINEX loss function at $a = -3$ and the SE loss function in terms of MSEs.

(4) The results indicate that both MLEs and Bayesian estimators yield favorable outcomes for two sample lines, whether they have the same or different population numbers, suggesting the model's suitability for various sample situations.

6.2.   Data analysis

The data represent the air-conditioning system failure times (in hours) for planes 7913 and 7914, originally sourced from Proschan ^[30]. The assumption is made that the two datasets are independent, and within each dataset, the failure times are also considered independent. The data is provided below.

Data 1. (Plane 7913): 1, 4, 11, 16, 18, 18, 18, 24, 31, 39, 46, 51, 54, 63, 68, 77, 80, 82, 97, 106, 111, 141, 142, 163, 191, 206, 216.

Data 2. (Plane 7914): 3, 5, 5, 13, 14, 15, 22, 22, 23, 30, 36, 39, 44, 46, 50, 72, 79, 88, 97, 102, 139, 188, 197, 210.

Table 5 displays the outcomes of the Kolmogorov-Smirnov (K-S) test, employed to assess the data's adherence to the NHD.

Taking into account the details in Table 5, it is apparent that the computed K-S values for the data are lower than the corresponding values expected at a significance level of $5\%$. Additionally, we have noted notably high P-values. Consequently, we can reasonably conclude that the NHD serves as a well-fitted model for the data. Moreover, we have generated empirical $S(x)$ and fitted $S(x)$ for each dataset, as depicted in Figures 2 and 3, respectively. These plots provide additional confirmation that the NHD model offers a superior fit to the data.

We generated a JPT-IISC sample from the aforementioned datasets using the following censoring scheme. For the first sample, set $m = 27$, and for the second sample, set $n = 24$, implementing JPT-IISC with $r = 20$. The censoring vectors are defined as:

Here are the generated datasets:

We obtained estimates for $\alpha_1, \alpha_2, \gamma_1$, and $\gamma_2$ using the MLE method, depending on the data type used in this study. The corresponding results are displayed in Table 6, whereas Tables 7 and 8 present the $95\%$ ACIs for $\alpha_1, \alpha_2, \gamma_1$, and $\gamma_2$. We utilized the MCMC method for Bayesian estimation conducting $20 000$ iterations to ensure convergence we excluded the initial $5000$ iterations as 'burn in' We have selected hyperparameters $c_i$ and $d_i$ as $0.0001$, approaching values near zero for the prior distributions. Bayesian estimates for $\alpha_1, \alpha_2, \gamma_1$, and $\gamma_2$ were derived using both the SE loss and LINEX loss functions, with the corresponding results detailed in Table 6. Additionally, the 95% CRIs for $\alpha_1, \alpha_2, \gamma_1$, and $\gamma_2$ are provided in Tables 7 and 8.

7.   Conclusions

In this paper, we explored statistical inference for two populations characterized by NHD. These distributions feature distinct shape and scale parameters, and our analysis focused on a JPT-II-CS. Assuming life distribution for both populations, we derived maximum likelihood estimates for unknown parameters and performed Bayesian estimation under gamma and non-information priors. Two loss functions, namely the SE and LINEX loss functions, were employed. To assess the effectiveness of the proposed estimates, we conducted Monte Carlo simulation experiments, revealing that Bayes estimates, along with their associated CRIs, outperform other estimators. Finally, we presented a numerical example to illustrate the inferential results established in this study.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research project was supported by the Researchers Supporting Project Number (RSP2024R488), King Saud University, Riyadh, Saudi Arabia.

Conflict of interest

The authors declare no conflict of interest.