Estimation and Prediction of Record Values Using Pivotal Quantities and Copulas

Recently, the area of sea ice is rapidly decreasing due to global warming, and since the Arctic sea ice has a great impact on climate change, interest in this is increasing very much all over the world. In fact, the area of sea ice reached a record low in September 2012 after satellite observations began in late 1979. In addition, in early 2018, the glacier on the northern coast of Greenland began to collapse. If we are interested in record values of sea ice area, modeling relationships of these values and predicting future record values can be a very important issue because the record values that consist of larger or smaller values than the preceding observations are very closely related to each other. The relationship between the record values can be modeled based on the pivotal quantity and canonical and drawable vine copulas, and the relationship is called a dependence structure. In addition, predictions for future record values can be solved in a very concise way based on the pivotal quantity. To accomplish that, this article proposes an approach to model the dependence structure between record values based on the canonical and drawable vine. To do this, unknown parameters of a probability distribution need to be estimated first, and the pivotal-based method is provided. In the pivotal-based estimation, a new algorithm to deal with a nuisance parameter is proposed. This method allows one to reduce computational complexity when constructing exact confidence intervals of functions with unknown parameters. This method not only reduces computational complexity when constructing exact confidence intervals of functions with unknown parameters, but is also very useful for obtaining the replicated data needed to model the dependence structure based on canonical and drawable vine. In addition, prediction methods for future record values are proposed with the pivotal quantity, and we compared them with a time series forecasting method in real data analysis. The validity of the proposed methods was examined through Monte Carlo simulations and analysis for Arctic sea ice data.


Introduction
Extreme weather and air pollution have been received steadily increasing attention over the past decade. Examples include extreme temperatures, the exceedances of flood peaks, and pollutant concentrations deviating considerably from expected average levels. In such cases, predicting observations more extreme than the current extreme values is an important issue. The topic of record values was introduced by Chandler [1], and Balakrishnan et al. [2] established some recurrence relationships for single and double moments of lower record values from the Gumbel distribution. Coles and Tawn [3] analyzed a daily rainfall series for modeling the extremes of a rainfall process in a sequence of lower record values is denoted by X L(k) , k = 1, 2, . . . from the original sequence {X 1 , X 2 , . . .} . These record values are heavily related to each other. Imani and Braga-Neto [8] proposed an efficient finite-horizon feedback controller similar to an optimal linear quadratic Gaussian estimator for partially-observed Boolean dynamical systems as a general class of nonlinear state-space model. Imani et al. [9] proposed an optimal Bayesian filter approach to the problem of recursive estimation in partially-observed Boolean dynamical systems. To establish the relationship between the record values, a copulas approach based on the canonical (C) and drawable (D) vine is proposed in this article.
Copulas have recently received much attention as a modeling tool for describing the dependency structure of multivariate data. The notion of a copula function can be found in Sklar [10]. One of the advantages of copula models is to build a variety of dependence structures based on existing parametric or non-parametric models of the marginal distributions. The copulas can be described as follows. Let F be the d-dimensional function of the random vector X = (X 1 , . . . , X d ) T with marginal distributions F X 1 (x), . . . , F X d (x). Then there exists a copula C such that for all x = (x 1 , . . . , x d ) T ∈ [−∞, ∞] d , where the copula C is unique if F X 1 (x), . . . , F X d (x) are continuous by Sklar's theorem (Sklar [10], 1959). In this case, the copula C can be interpreted as the distribution function of a d-dimensional random variable on [0, 1] d with uniform marginal distributions. Tsung et al. [11] conducted a comprehensive literature review of statistical transfer learning methods focusing on statistical models and statistical methodologies, including a Gaussian copula. Rocher et al. [12] proposed a generative copula-based method that can elaborately estimate the likelihood that a particular person will be correctly re-identified, even in a very incomplete dataset. Vine copulas were introduced by Joe [13] to overcome limitations of standard multivariate copulas in higher dimensions, where standard multivariate copulas lack the flexibility of accurately modeling the dependence.
Aas et al. [14] described statistical inference techniques for the C-and D-vine copulas. Berg and Aas [15] and Fischer et al. [16] showed the excellence of the D-vine copula approach, compared to alternative copulas in constructing higher dimensional dependency structures.
To the best of our knowledge, modeling the dependence structure between record values numerically has been little explored. In this paper, we propose an approach with which to model the relationship between the record values based on C-and D-vine copulas and to predict future record values.
The remainder of the article is organized as follows. Section 2 introduces C-and D-vine copulas and provides pivotal-based approaches to estimate the model parameters and to predict future record values by proposing a new algorithm to deal with a nuisance parameter. Section 3 presents simulation studies to validate the proposed approaches. We applied the methods to Arctic sea ice data; see Section 4. Concluding remarks and some discussions are in Section 5.

Methods
Let f (x; θ) and F (x; θ) be the marginal density and distribution functions of X, respectively, where θ is an unknown parameter. Then a schematic diagram (Figure 1) of our method is given by The C -and D-vine copulas are first described in the following subsection.

C-and D-Vine Copulas
A vine is a flexible graphical model that decomposes a multivariate probability distribution into bivariate copulas, where each pair-copula can be chosen independently from the others [14]. This article considers C-and D-vine copulas to model the relationship between record values based on C-and D-vine copulas.
The C-vine decomposition is given by Then, we can specify the pairs of the d-dimensional C-vine copula model in the following order: The D-vine decomposition is given by Similarly, the pairs of the d-dimensional D-vine copula model are specified in the following order: To measure the dependence of each pair-copula, we consider tree 1 that can be employed to obtain Kendall's τ (Nelsen [17], 2006) given by In a C-and D-vine, consider the exponentiated Gumbel distribution (EGD) with the CDF where σ and λ are the scale and shape parameters. Then, u i is the value of the marginal distribution of where H(x) = − log F(x) (Ahsanullah [18] and Arnold et al. [19]). That is, for d − 1 pairs of data points {(x L(i) , x L(j) ); i < j}, the corresponding couples {(u i , u j ); i < j} can be computed from the marginal distribution (2). In addition, the marginal density function of X L(i) is given by Note that it is necessary to estimate σ and λ for computing the values of u i and u j .

Pivotal-Based Approach
Here we present a pivotal-based method to estimate the parameters of the CDF (1). First, a lemma is introduced to deal with nuisance parameters in order to establish the relationship between record values. Lemma 1. Let X L(1) , · · · , X L(k) be the lower record values from the CDF (1). Then, (a) T k (σ, λ) = 2λe −X L(k) /σ has a χ 2 distribution with 2k degrees of freedom; has a F distribution with 2k − 2 and 2 degrees of freedom; has a χ 2 distribution with 2k − 2 degrees of freedom.
Proof. Let X L(1) , · · · , X L(k) be the lower record values from the CDF (1). Then, we have that have a standard exponential distribution, and that leads to the following spacings which is independent and identically distributed as the standard exponential distribution with mean 1. From the spacings, the pivotal quantities are derived, which are independent random variables such that there is an χ 2 distribution with 2j(j = 1, . . . , k). From (3), the pivotal quantity (a) is easily proved. In addition, and T 1 (σ, λ) are independent random variables, given the fact that S i (i = 1, . . . , k) are independent and identically random values, as mentioned earlier; both have a χ 2 distribution with 2k − 2 and 2 degrees of freedom, respectively. On top of that, we can derive the following pivotal quantities by using Lemma 2 of Wang et al. [4] in (3) which are independent and identically distributed as the uniform distribution on the interval (0, 1). Then, the pivotal quantity (c) is proved as From Lemma 1(c), the unique solution σ * is given by where Q W follows a χ 2 distribution with 2k − 2 degrees of freedom. Then, for any 0 < α < 1, an exact 100(1 − α)% CI for σ based on W(σ) is given by Note that the exact CI is the equal-tail CI because it splits the probability equally, putting α/2 in each tail of the distribution. For any 0 < α < 1, an exact 100(1 − α)% CI with the shortest-length based on σ * is given by where l * is chosen so that Similarly, the unique solution from V(σ) in Lemma 1 is given by where Q F follows a F distribution with 2k − 2 and 2 degrees of freedom. Then, with the same argument, the exact 100(1 − α)% equal-tailed and shortest CIs for σ based on V(σ) can be constructed. In Section 3, it is found that W(σ) provides a more efficient CI than V(σ) in terms of average lengths (ALs) of the CIs through Monte Carlo simulations, as in the case of Seo and Kim [5].
For λ, we have that by putting T k = T k (σ, λ) in Lemma 1. In addition, let g(W, x L ) be the unique solution of W(σ) = W for W > 0, where x L = x L(1) , · · · , x L(k) . Then, by substituting g(W, x L ) for σ in (4), the following generalized quantity is given by The existing literature (Wang et al. [4] and Wang et al. [20]) supposed that W has a χ 2 distribution with 2k − 2 degrees of freedom, and then obtained the percentiles of ψ generating W and T k independently from the χ 2 distribution with 2k − 2 and 2k degrees of freedom, respectively, althrough W and T k are not independent. As an alternative, the following algorithm is proposed to obtain the percentiles of ψ.
Step 3. Compute W = 2 ∑ k−1 j=1 log T k T j and solve the equation W(σ) = W for σ to obtain g(W, x L ).
From the algorithm, the equal-tailed and shortest CIs for λ based on ψ are given by respectively, where l * is chosen so that

Prediction
Let X L(r) (r = k + 1, k + 2, · · · ) be the future lower record values. Then, for any 0 < α < 1, the conditional quantile is given by where F X L(r) |x L(k) (·) is the conditional distribution of X L(r) given x L(k) . However, the quantile cannot be obtained numerically because it does not have closed forms. Instead, we propose a pivotal approach based on the following lemma.
in the conditional density function of X L(r) given x L(k) be defined by Ahsanullah [18] as Then, Y C has a gamma distribution with the parameters (r − k, 1). (5). Then, the Jacobian of the transformation is

Proof. Let
and the density function of Y C is which is the probability density function of a gamma distribution with parameters (r − k, 1).
With the same argument as Section 2, an algorithm for obtaining the Markov-chain Monte-Carlo (MCMC) samples X i L(r) (i = 1, . . . , N) based on the pivotal quantity Y C is provided as follows.

Simulation Study
A simulation study was performed to examine the validity of the proposed pivotal-based approach in terms of the coverage probabilities (CPs) and ALs of the proposed confidence intervals (CIs). The lower record values with sizes k = 6(2)12 were generated from the standard EGD distribution with λ = 0.5, 0.8, and 1.5. To construct the 95% exact CIs described in Section 2.2, N = 20,000 MCMC samples were generated, and the corresponding CPs and average lengths (ALs) were computed over 10,000 simulations. The results are reported in Table 1 along with those for the classical inference (see Proof) (Appendix A) for comparison. Table 1 shows that the CIs using MCMC samples have nearly same results as those using the classical method, and all considered CIs are well matched to their corresponding nominal levels; however, the CIs based on W(σ) have shorter length than those based on V(σ). In addition, all ALs decrease with an increase in the size of record values k. For ALs, the CIs with the shortest-lengths have shorter lengths than those with equal-tails, as expected.

Application: Arctic Sea Ice
Sea ice maintains the Earth's average temperature by reflecting solar energy and keeping the polar regions cool. Currently, the Arctic is warming faster than any other region on earth. The warming of the Arctic Circle leads to a decrease in sea ice, which again causes warming of the Arctic Circle. In addition, it causes global weather changes such as summer heat waves, winter cold waves, and heavy snow. These climate changes are leading to disturbances of ecosystems formed around Arctic sea ice and changes in habitats. For this reason, the importance of sea prediction systems to cope with climate change is increasing. The National Aeronautics and Space Administration (NASA) reported that the area covered by Arctic sea ice has decreased by about ten percent in the last 30 years (Figure 2). This section analyzes the smallest annual Arctic sea ice extent (see Table 2) from October 1978 to October 2018 extracted from the National Snow & Ice Data Center (NSIDC). To measure goodness of fit of the EGD, the replicated data of observed lower record value x L(i) were generated from its marginal density function with σ * and ψ. All results were obtained by generating N = 50,000 MCMC samples. In addition, based on the results in the previous simulation study, σ * from W(σ) was only considered in this data analysis.
The 95% confidence region for the replicated data was plotted in Figure 3. It was found that the confidence regions decreased as the record value of the smallest annual Arctic sea ice decreased. The correlation coefficient between the observed and expected lower record values indicates a strong association.
To examine the relationship between observed record values, Figure 4 plots the first trees of the C-and D-vines with the best copula function in terms of the Akaike information criterion (AIC) and corresponding Kendall τ.  Note that the AIC is defined as −2 ln(L) + 2k, where L is the likelihood function and k is the number of estimated parameters of the model. Therefore, the smaller the AIC, the better. The entire result for the relationships between observed lower record values is reported in Figure 5. It is shown that the observed lower record values have a positive dependence on each other. In addition, the Kendall's τ values increase as the interval between the lower record times decreases. That indicates that x L(i) and x L(j) such that j − i = 1 for j = i have the strongest dependency in terms of the Kendall τ. It is worth noting that the strength of dependency between x L(i) and x L(j) such that j − i = 1 for j = i becomes stronger as the lower record times increase. The 95% exact CIs for σ and λ are reported in Table 3, which shows a similar pattern to the simulation results. For prediction, the last lower record value x L(9) was assumed to be unknown, and a time series analysis was conducted, in which it was expected that differences of the observed lower record values could yield a stationary time series because the observed lower record values had a decreasing pattern. In fact, the ARIMA (0, 1, 0) model was chosen as the best model in terms of the AIC from an ARIMA (p, d, q) model, where p is the autoregressive (AR) model order, d is the difference order, and q is the moving average (MA) model order. Table 4 and Figure 6 present the results for future record values of the least annual Arctic sea ice. Table 4 shows that there is little difference in measures of center such as the mean and median for the predictions of the future lower record values based on the pivotal quantity Y C .  For the last lower record value, the ARIMA (0, 1, 0) model provides a closer predictive value than the mean of X * L(9) |x L(8) to the actual value of 1.29, while the PI from the ARIMA (0, 1, 0) model has a longer length than that for X L(r) |x L(k) based on the pivotal quantity Y C . Finally, Figure 6 shows that as the future record time L(r) increases, the variance of the predicted future record value from the conditional density function increases.

Conclusions
This article proposed a copula approach with which to model the dependence structure between record values from the EGD and to predict future lower record values using a pivotal-based method. In the pivotal-based method, a new algorithm for dealing with a nuisance parameter has been proposed; it not only is very computationally convenient in constructing exact CIs with the shortest lengths, but also provides very satisfactory results in terms of the CPs and ALs, compared with the classical method. In the approach based on the C-and D-vine copulas, we chose the best copula model in terms of the AIC among 40 paircopula families and it showed very intuitive and reasonable results in analysis based on real data. An interesting point is that the strength of the dependency between x L(i) and x L(j) such that j − i = 1 for j = i becomes strong as the lower record times increase in real data analysis. The proposed method is applicable to recording values of other real data that have a probability distribution if the CDF of the probability distribution has a closed form, such as an extreme value distribution. The prediction results of this paper indicate that we should be alert to the decrease in Arctic sea ice extent. In future studies, we envision extending this work to predict the size and decreasing rate of Arctic sea ice extent in real time.  and dθ 2 (α) dθ 1 (α) = g(θ 1 (α)) g(θ 2 (α)) , so that dL dθ 1 (α) = −h(θ 1 (α)) + h(θ 2 (α)) g(θ 1 (α)) g(θ 2 (α)) (X L(1) − X L(k) ), where h(θ 1 (α)) = k − 1 [1 + (k − 1)θ 1 (α)] [log (1 + (k − 1)θ 1 (α))] 2 , h(θ 2 (α)) = k − 1 [1 + (k − 1)θ 2 (α)] [log (1 + (k − 1)θ 2 (α))] 2 .