Extreme value modeling with errors-in-variables in detection and attribution of changes in climate extremes

The generalized extreme value (GEV) regression provides a framework for modeling extreme events across various fields by incorporating covariates into the location parameter of GEV distributions. When the covariates are subject to errors-in-variables (EIV) or measurement error, ignoring the EIVs leads to biased estimation and degraded inferences. This problem arises in detection and attribution analyses of changes in climate extremes because the covariates are estimated with uncertainty. It has not been studied even for the case of independent EIVs, let alone the case of dependent EIVs, due to the complex structure of GEV. Here we propose a general Monte Carlo corrected score method and extend it to address temporally correlated EIVs in GEV modeling with application to the detection and attribution analyses for climate extremes. Through extensive simulation studies, the proposed method provides an unbiased estimator and valid inference. In the application to the detection and attribution analyses of temperature extremes in central regions of China, with the proposed method, the combined anthropogenic and natural signal is detected in the change in the annual minimum of daily maximum and the annual minimum of daily minimum.


Introduction
Extreme events arise in many fields and often have catastrophic impacts.Examples include extreme weather (Field et al. 2012), natural hazards (Hamdi et al. 2021), and market crashes (Nolde and Zhou 2021).Statistical methods for prediction and inference in analyzing extreme events need to account for their unique distributional characteristics.The generalized extreme value (GEV) distribution is the natural choice as the limiting distribution of sample maximum (Fisher and Tippett 1928).Covariates can be incorporated into the location parameter to predict extreme events or understand the behavior of extreme observations (Coles et al. 2001).Such GEV regression models have been widely used in an enormous variety of applications (Smith 1989;Wang et al. 2004;Huerta and Sansó 2007;Hyun et al. 2019).
The validity of GEV regression analyses depends on a critical assumption that there is no errors-in-variables (EIV) or measurement error in the data.This assumption may be violated in many applications, such as the detection and attribution analyses in the context of climate extremes (see Sect. 3).One major tool for the detection and attribution analyses is fingerprinting, which was first introduced to study changes in mean climate conditions (e.g., mean temperature) (Hasselmann 1993;Hegerl et al. 1996;Allen and Stott 2003), and later extended to extreme climate conditions (Zwiers et al. 2011;Wang et al. 2017Wang et al. , 2021)).The fingerprinting method fits a GEV distribution to the observed extremes and takes the fingerprints of external forcings of interest as covari-ates in the location parameters.The fingerprint, or signal, of a particular external forcing, is not observed but estimated from climate-model-simulated extremes, which introduces errors (i.e., EIV) and causes bias in the resulting estimator of the coefficients of the fingerprints (Wang et al. 2021).
Despite the voluminous literature on EIV issues in other regression settings, it has been rarely studied in GEV regressions even for independent data, let alone for time series or spatiotemporal data.The main challenge owes to the complex structure of GEV distribution, which prevents a straightforward extension of many methods.In the context of the fingerprinting, for example, though the EIV issue has been solved in the mean climate models under a linear regression setting by total least squares (Allen and Stott 2003), integrated likelihood (Hannart 2016), and Bayesian model averaging (Katzfuss et al. 2017), none of these approaches can be applied to the extremes setting in a relatively straightforward manner.In fact, the bias issue in GEV modeling has been explicitly raised but left unsolved (Wang et al. 2021).Additionally, the correlation in the measurement errors further complicates the problem.Among a wide range of applications of GEV modeling of spatiotemporal data, the EIVs are likely dependent.One example is again the fingerprinting with time series, where the signals are most likely to be temporally dependent as the parameters in the signals are often shared over time.
A promising approach to GEV regression with EIV is the corrected-score method (Nakamura 1990).In principle, this method constructs a corrected score equation whose expectation conditional on the contaminated covariates is the same as the original score equations evaluated at the true covariates.It is attractive because, as a functional approach, it does not require any distributional specification of the covariates as some approaches do.Given the complex expressions of the score functions under the GEV regression, which can be derived from the score functions of GEV distribution (Prescott and Walden 1980) using the chain rule, such construction is challenging because the expectations of the score functions in terms of the contaminated covariates have no closed forms, which prohibits closed-form corrections.The Monte Carlo corrected scores (MCCS) method, however, is well suited for the situation.It approximates the expectations of the score functions in terms of the contaminated covariates by Monte Carlo integration, which facilitates corrected scores (Novick and Stefanski 2002).
Our contributions are two-fold.First, as the first attempt to address the EIV issue in GEV modeling using a functional approach, we propose an MCCS-based approach and extend it to the settings where the EIVs are temporally dependent as in fingerprinting applications.We further develop a small-sample correction to the sandwich variance estimator.This is practically important because the record lengths of extremes are usually insufficient for the variance estimator to work well.Simulation studies show that the proposed method corrects the bias in estimating all the model parameters, especially the regression coefficients in the location parameter.The coverage rates of the confidence intervals based on our corrected sandwich variance estimator are acceptable.Our second contribution is an application of the methods to the detection and attribution analyses of changes in temperature extremes.The method accounts for the fact the unobserved signals are estimated from climate model simulations under the corresponding forcings.In applications we demonstrate the validity of the method via simulation studies that mimic the fingerprinting setting and real data.
The rest of the paper is organized as follows.Section 2 provides the MCCS framework under the GEV setting for both independent and temporally dependent EIVs.Section 3 adapts the method to the specific setting of fingerprinting for climate extremes in detection and attribution analyses.The methods are validated in simulation studies with both independent extremes and fingerprinting settings in Sect. 4. Section 5 reports an application to detection and attribution analyses of changes in temperature extremes in central China.A discussion concludes in Sect.6.

Statistical model
Consider GEV modeling with a regression setting for the location parameter.Let F(• | μ, σ, ξ ) denote a GEV distribution function with location μ, scale σ , and shape ξ , and Y t be observed extreme with p-dimensional error-prone covariates X t and q-dimensional accurate covariates Z t for t = 1, . . ., n. Suppose that Y t follows a GEV distribution with a location parameter that incorporates covariate effects.In particular, where α and β are q-and p-dimensional regression coefficients, respectively.An intercept term is included in α by including a one in Z t .The parameters of interest are θ = (α, β, σ, ξ).
Assumption 1 When X t is observed, there exists an unbiased estimating function ψ(Y t , X t , Z t , θ) with parameters θ such that Then, the parameter θ can be estimated consistently by solving By fulfilling Assumption 1, the score function ψ(Y t , X t , Z t , θ) of GEV loglikelihood is an unbiased estimating function, where and f 1 , f 2 , and f 3 are the partial derivatives of the log density of F(• | μ, σ, ξ ) with respect to μ, σ , and ξ , respectively, with specific expressions given in Appendix A (Prescott and Walden 1980).The maximum likelihood estimator (MLE) for θ can be obtained by solving the score equation ) has expectation zero and finite second moments (Smith 1985;Bücher and Segers 2017).
In situations where instead of X t , a surrogate version W t is observed, substituting X t 's directly with W t in estimating function (2) leads to a biased estimator of θ .The magnitude of the bias depends on the measurement error specification.Assume that where the p-dimensional measurement error e t , independent of X t and Y t , is normally distributed with mean zero and covariance matrix t .
Assumption 2 When X t is not observed, there exists a corrected score function of the observed data such that for Based on Assumptions 1-2, we seek a corrected score function ψ(Y t , W t , Z t , θ), which is conditionally unbiased given the observed quantities and results in estimating equations While the corrected scores for some common families like normal and Poisson are available in a closed-form (Nakamura 1990), it is hard to obtain that of the GEV distribution because of its nonlinear and complicated structure.Therefore, we resort to the MCCS method developed by Novick and Stefanski (2002), where the conditional expectation in (4) is approximated by Monte Carlo integration.This method is relatively robust and can yield satisfying results in certain settings even when the normality assumption of e t is violated.
A further complication in practice, however, is that the measurement errors e t 's in GEV could be correlated in various ways.For instance, in fingerprinting applications, the fingerprint or signal of an external forcing cannot be observed but is generally estimated from climate model simulations.The measurement errors e t 's could be temporally correlated over time and we denote the covariance matrix for e = (e 1 , . . ., e n ) as .This brings an additional complexity compared to the classical MCCS approach in Novick and Stefanski (2002).In the sequel, we address this issue after dealing with independent EIVs.

MCCS for GEV regression with independent EIV
First, we illustrate the MCCS method for GEV regression under the assumption that measurement error e t 's are mutually independent for t = 1, . . ., n.Let b,t be the bth Monte Carlo copy generated from N (0, t ), where t is the variance matrix for e t , t = 1, . . ., n.The complex variates are constructed as W b,t = W t + ι b,t for b = 1, . . ., B, where B is a large integer and ι = √ −1.Then, the following Lemma was established by Novick and Stefanski (2002).
Lemma 1 (Novick and Stefanski 2002) When the measurement error e t 's are mutually independent, under integrability conditions and the assumption that ψ(Y t , X t , Z t , θ) is an entire function of X t in the complex plane, where (•) is the real part of a complex number.
Consequently, {ψ(Y t , W b,t , Z t , θ)} is a corrected score based on a randomly generated random vector b,t .Then, we use Monte Carlo integration to approximate the conditional expectation as with its specific form in GEV provided in Appendix A. As B → ∞, the MCCS converges to the exact conditional expectation (Novick and Stefanski 2002).The MCCS estimator θn,B solves the estimating equations In practice, t is unknown and needs to be replaced with a consistent estimator in the above derivation.We provide more details on estimating t in Sect.3. Let θ 0 be the true value of θ .When (Y t , W t ) are mutually independent for t = 1, . . ., n, and ψB (Y t , W t , Z t , θ, t ) and ψ(Y t , W t , Z t , θ, t ) are differentiable, there exists a consistent estimator for θ by solving (Nakamura 1990).The solution can be obtained via nleqslv (Hasselman 2018), a nonlinear equation solver package in R.
The proof based on standard estimating equation theory is outlined in Appendix B.1.

MCCS for GEV regression with temporally dependent EIV
In many GEV applications, such as the fingerprinting setup, measurement error e t 's could be temporally correlated (Wang et al. 2021).Here we explain how to adapt the aforementioned MCCS approach to time-dependent EIV.To construct Monte Carlo replicates and Z = Z 1 , . . ., Z n to incorporate the temporal correlations.Such practice can be difficult or undesired for GEV modeling as the temporal dependence may be hard to specify.As θ is assumed to be shared across n years, we again consider the corrected score based on each year t as in Eq. ( 5) and obtain the unbiased estimator θn,B .Though the estimating equation is constructed at each time point, the validity of the approach is stated in Proposition 1.
Lemma 2 Assume that ψ(Y t , X t , Z t , θ) is an entire function of X t in the complex plane, and that measurement error e = (e 1 , . . ., e n ) are temporally correlated with covariance matrix .With * b ∼ N (0, ), the following equations still hold: This Lemma is based on representing ψ as a multivariate power series and then applying Lemma 4 from Stefanski and Cook (1995).More details can be found in the Appendix B.3.
Then, we define and the MCCS estimator θ * n,B solves the estimating equations Again, the R package nleqslv can be used to solve Equation (7).
Proposition 2 Under certain regularity conditions, with n → ∞, The detailed proof is provided in Appendix B.2.
Remark 1 Our work is built upon standard estimating equation theory.There are some sufficient but not necessary conditions for our asymptotic theory to hold as below.
1.The parameter space is compact, and there is a unique θ in the parameter space satisfying E{ ψ * B (Y t , W t , Z t , θ, )} = 0 for all t = 1, . . ., n. 2. The estimating equation ( 7) has a unique solution.3. The matrix D * is of full rank within a neighborhood of the true value θ 0 .For sufficiently large n, D n, ψ * is of full rank with eigenvalues bounded away from 0 and ±∞ within a neighborhood of the true value θ 0 .
Sandwich variance estimators are known to under-estimate the variance for small samples (Mancl and   This bootstrap procedure is computationally efficient as it only requires evaluating the outer-products of individual estimating functions instead of solving any estimating equations.Of note, the bootstrap is only used for estimating the middle part of the sandwich estimator.As it does not bootstrap extremes directly, this procedure also avoids some potential issues in sampling extremes (Gilleland 2020), and the confidence intervals of the estimates are still built based on the asymptotic normality.

Application in fingerprinting with time series
Fingerprinting in the context of extremes is a direct application of the proposed method.Fingerprinting aims to quantify the impact of external forcings on the climate variable of interest, in which the expected responses of the climate system to external forcings, also called fingerprints, are treated as predictors (Li et al. 2021(Li et al. , 2023)).These external forcings can be classified into two types: anthropogenic and natural, which are factors that impact the climate system while not being part of the climate system itself, such as greenhouse gas emissions and volcanic eruptions.The fingerprints are not observable but can be inferred from numerical climate model outputs.Climate models typically contain multiple initial condition ensembles, which, most of the time, can be treated as random independent and identically distributed samples from one multivariate spatial-temporal process due to their sensitivity to the initial conditions (Stein 2020).The estimated signals from these outputs are naturally temporally dependent.Here we limit our scope to a setting with time series data instead of spatiotemporal data; handling spatiotemporal EIVs requires a separate investigation.We investigate this application setting as follows.
Suppose the observed climate extremes Y t specified in Sect. 2 are available at a particular grid box for year t = 1, . . ., n.The climate extremes of interest can be station annual maximum of daily maximum (TXx) and minimum (TNx) temperatures or the annual minimum of daily maximum (TXn) and minimum (TNn) temperatures.The errorprone covariates X t in Equation ( 1) are signals of external forcings, and β are the corresponding coefficients of the signals.We also assume a constant shape and scale parameters throughout the year at the site, and an intercept is included in the location parameter, which makes Z t a constant of 1 and α a scalar (Wang et al. 2021;Zwiers et al. 2011).
Since the signal vector X t is unknown, it needs to be estimated from the climate models simulations.Following Zwiers et al. (2011) and Wang et al. (2021), we assume another GEV model for the simulated climate extremes.Let L j denote the total number of available ensemble runs under the jth forcing from the climate models, and U (l) t j extreme observations for the lth ensemble run under the jth forcing in year t, where j = 1, . . ., p, t = 1, . . ., n, and l = 1, . . ., L j .The annual extremes U (l) t j can be modeled as follows, where σ j and ξ j are forcing-specific scale and shape parameters, respectively, and X t j is the fingerprint or the signal of forcing j characterized by function μ with parameter vector γ j .To account for the temporal dynamics, a flexible formulation for μ with B-splines is where S t = (S 1 (t), . . ., S D (t)) , S d (t)'s are a set of B-splines basis with D degree of freedom, and γ j = (γ j,1 , . . ., γ j,D ) is the coefficient vector for the basis.The B-spline representation is assumed to be flexible enough to approximate the temporal pattern in signals.With all the ensemble runs under forcing j as independent replicates, we can obtain the MLE γ j of γ j and its covariance matrix γ j .Therefore, W t j = D d=1 γ j,d S d (t) is an MLE of X t j with covariance matrix S t γ j S t .
Given the B-spline model is correctly specified, using W t = (W t1 , . . ., W t p ) rather than the true X t introduces additional variability to Equation (1), resulting bias when estimating β.It has the same effect as EIV and needs to be recognized.This is exactly the unsolved issue raised by Wang et al. (2021), which results in biased estimation.Here the relationship between X t and W t is the same as in Equation (3), that is, W t = X t + e t .In general, e t and e t for any t = t are not independent as the variance of e t j and e t j both depend on γ j .The covariance of e j = (e 1 j , . . ., e nj ) , can be estimated by S γ j S with S = (S 1 , . . ., S n ) .This fits the situation in Sect.2.3.For each j = 1, . . ., p, we can generate W b, j = ( W b,1 j , . . ., W b,nj ) from a multivariate normal distribution with mean W = (W 1 j , . . ., W nj ) and covariance matrix S γ j S. Since the runs under different forcings are independent, we can generate the Monte Carlo copies separately for each j.Then, W b,t j for any t = 1, . . ., n and j = 1, . . ., p can be used to construct the corrected score in Equation ( 6), yielding the unbiased estimates of θ as explained in Sect. 2.
The coefficients of the fingerprints β, also known as scaling factors, are the primary target of the inference.To account for the time dependency in the measurement error, we rely on the block bootstrap method with 5-year blocks detailed in Sect.2.3 to construct the 90% confidence interval for β, which is frequently used in attribution and detection analyses.If the 90% corresponding confidence interval of a scaling factor is above zero, then the fingerprint is claimed to be detected.Additionally, given that a fingerprint is detected, a 90% confidence interval covering one is the necessary evidence that the observed changes can be attributed to the corresponding external forcing.

Independence setting
We first considered the case with independent measurement errors.For simplicity, X t is a univariate normal random variable with mean 0 and standard deviation 2, and Z t = 1.The observed extremes Y t was subsequently generated from a generalized extreme value distribution with θ = (1, 1, 4, −0.2) by using the R package evd (Stephenson 2002).The observed surrogate W t was then generated from Equation (3) where reduces to σ 2 e I with σ 2 e = 1.To carry out the MCCS method, W b,t was obtained by adding ι b,t to W t , where b,t is also a normal variable with mean 0 and variance σ 2 e = 1, and b = 1, . . ., B.Here we treat σ e as known and focus on estimating θ .In the simulation, we consider the Monte Carlo sampling B ∈ {200, 400} and the sample size n ∈ {100, 200}.For each configuration, 1000 data sets have been generated.Also, for each dataset, we obtained MLE for θ without correcting for W t and θ by implementing MCCS with W b,t .The MLE estimation was carried out by the package ismev (Heffernan and Stephenson 2018), and the MCCS was solved by using nleqslv.
Table 1 summarizes the results in estimating θ .The MLEs of θ are the naive estimates with measurement error ignored.The naive estimates of α, σ , and ξ are unbiased under all configurations.The naive estimate of parameter β, however, is severely biased, and the coverage rate is low across different sample sizes.On the other hand, the proposed method provides an unbiased estimate of β, and the coverage rate is close to the nominal level 90% for all the settings, including the case that B = 200.The other parameters are consistently estimated, and their coverage rates are close to the nominal level as well.In general, as n increases, the convergence rate improves, and the bias, RMSE, ESE, and ASE decreases as n increases.Regarding the performance of the algorithm, the percentage of successfully finding a root using nleqslv is around 95% when n = 100 and improves to 99.5% when n increases to 200 in our simulation.Given a fixed n, however, the performance of MCCS estimators is very similar across different Bs in our setting.Nonetheless, the proposed method works well in correcting bias in β.

Fingerprinting setting
We also carried out simulation studies with fingerprinting setting described in Sect. 3 with one forcing ( p = 1).In order to mimic the real detection and attribution process, we first simulated the climate model extremes U (l) t from a generalized extreme model with σ = 2 and ξ = − 0.2.The location parameter X t is determined by a cubic spline where γ is uniformly distributed in (38,40).The number of ensemble runs L is assumed to be 50, and for each run, we have 100 years of data.Then the observed extreme Y t was generated from a generalized extreme value distribution with {α, β, σ, ξ } = {40, 1, 1, − 0.2}.Again, two levels of the Monte Carlo sampling size, {200, 400} are considered here.
Table 2 summarizes the results in estimating θ in the fingerprinting setup.Both methods provide unbiased estimates to σ and ξ .Unlike the previous results, the bias is observed in both α and β in the naive estimates, resulting in extremely poor coverage rates.This could lead to misleading conclusions in the detection and attribution analysis which mainly relies on the statistical inferences of β.In contrast, MCCS estimates of α and β are nearly unbiased.Even though the bias of α estimated by MCCS is around −1.2 in Table 2, it is relatively unbiased given that the true value of α is 40.
More importantly, the coverage rate of the scaling factor β in MCCS is much closer to 90%, which implies a more reliable decision on the detection and attribution of a signal.Even though the coverage rates of σ and ξ are lower than the nominal level using the MCCS method, they are still around 85%, and their MCCS point estimates are also unbiased.In simulations, the convergence rate of the MCCS method is around 99%, indicating high success in finding a root.

Temperature extremes in central China
The observed annual temperature extremes in the East Asia region over the period of 1951-2010 were extracted from the HadEX2 data (Donat et al. 2013).For each temperature extreme described in Sect.3, grid boxes of resolution 5 Since observations on land are more reliable than those in the ocean, we focused on twelve land grid boxes with latitude ranging from 30N to 40N and longitude covering from 100E to 115E, which covers a large part of inland China (Fig. 1).
The ALL (combination of anthropogenic and natural external) forcing is the single external forcing being considered.To avoid making any assumptions about the exchangeability of different climate models, we only used the secondgeneration Canadian Earth System Model (CanESM2) (Chylek et al. 2011) to estimate signals, which contains 50 runs under the ALL forcing.In order to match the spatial resolution of the HadEx2 data, the CanESM2 data also has been regridded to 5 • × 5 • grid boxes.The 50 runs are independent replicates simulated under both anthropogenic factors (greenhouse gas and aerosol emissions) and natural-only factors (solar radiation and volcano eruptions).The signals under the ALL forcing at each grid box were estimated by fitting (8) to the 50 runs.A cubic spline with a 5-year knots was used to estimate the ALL signal.In regional detection and attribution analyses, it is often assumed that the scaling factors across all the grid boxes are the same (Zwiers et al. 2011;Wang et al. 2017Wang et al. , 2021)), in which case, the spatial dependence among the grid boxes needs to be addressed.Since a spatial extension of the MCCS method needs a full investigation, here we apply the method to each grid box separately.
At each grid box, the fingerprinting method in Sect. 3 was applied to each of the four temperature extremes.The minimum extremes TXn and TNn were negated before the analyses.The size of Monte Carlo sampling was fixed at B = 1, 000 for each analysis.For comparison, we also included the naive estimators, which ignored the measurement errors in the estimated fingerprint and fitted the data using the MLE method.To ensure ξ > −0.5, we use a Beta(9, 6) prior to constraint ξ within a range of (−0.5, 0.5) (Martins and Stedinger 2000) in both methods as done by Wang et al. (2021).
Figure 2 displays the 90% confidence intervals of the scaling factor for each temperature extreme in each of the 12 grid boxes.Additional figures of other parameters are presented in the Appendix 1.The estimates of ξ are negative in most grid boxes, which is consistent with the findings in Huang et al. (2016).The influence of the ALL forcing is detected for all the temperature extremes in most of the 12 grid boxes.In particular, the most detectable influence is on TNx, the annual maximum of daily minimum temperature, and the least detectable influence is on TXx, the annual maximum of daily maximum.In most cases, both naive and proposed methods reach consistent conclusions, although the point estimates of the scaling factor can be quite distinct from both methods in temperature minimums.For the maximums TNx and TXx, for example, the estimates of β are quite close except for the (115E, 30N) case.The analysis on the minimums, especially on the TXn, displays a much larger discrepancy in the β.In the (110E, 40N) case, for instance, the difference in the βs bring about different conclusions in the attribution decision.Moreover, the proposed method tends to produce a larger standard error in most cases, which might lead to different conclusions even though the estimated values of β are close for both approaches.For example, both methods produce similar values in β in the (105E, 40N) grid box, but the lower bound of the 90% confidence interval is much closer to 1 using the proposed method while the counterpart for the naive method is way above 1.Overall, both methods agree on the detection and attribution Fig. 2 The top and bottom graphs present the scaling factor estimation for all the temperature extremes and the corresponding 90% bootstrap confidence intervals for the selected grid boxes results on TNx and TXx for all the grid boxes.ALL forcings are always detected on TNx using both methods, but its evidence of attribution is obvious in the grid boxes with latitude 40N.Different from the results on TNx, the influence of ALL forcings are not detected in grid boxes (110E, 30N), (115E, 30N), (110E, 35N), and (115E, 35N).There is also no evidence that the change in TXx is attributable to ALL forcings.For TNn and TXn, there are a few cases that the ALL forcings are not detectable.For the (115E, 30N) grid box, the influence of the ALL forcing is not detectable using the proposed method, while the naive method reach a different decision.In terms of attribution, only the (115E, 30N) and (115E, 40N) grid boxes show necessary evidence of attributing the ALL influence to the change in TNn.Similar to the results of TNn, the change in TXn can only be attributed to ALL forcings in two grid boxes.Interestingly, both methods agree on the attribution results for the (100E, 30N) grid box, but disagree on the attribution results for (110E, 40N) grid box.
For illustration purposes, we choose the (110E, 40N) grid box and report the point estimates (Est), the corresponding standard errors (SE), and the 90% confidence intervals (CI) of θ for each temperature extreme based on the two methods (Table 3).The estimates of ξ and σ between the two methods are quite similar, which is consistent with the simulation results in Sect. 4. On the other hand, the location parameters α and β can be very different between the naive estimates and the MCCS estimates in some cases.For instance, the estimates of α and β are very different between the two methods on TXn, directly resulting in a different conclusion.The naive method detects the ALL signal at all on TXn, but there is no evidence that the change on TXn is attributable to the influence of ALL forcings.On the opposite side, ALL forcings are not only detected by the proposed method but also demonstrate necessary evidence of attribution on TXn, with the corresponding 90% confidence interval being significantly greater than zero.For extremes other than TXn, both methods detect the ALL signal to the temperature change.

Discussion
Relatively little attention has been paid to understanding and correcting for the EIVs in the location parameters of GEV regressions.Motivated by the optimal fingerprinting method in the detection and attribution analysis of changes in climate extremes, we investigated the impact of measurement error in both the independent and time-correlated measurement error settings through simulation studies.The bias caused by the measurement error can be severe.Our MCCS method takes advantage of the existing score functions of the GEV distribution and the flexibility of Monte Carlo copies to correct for the bias of estimators.For inferences, the confidence intervals based on our block bootstrap approach appear to give acceptable coverage rates.Though our method was applied to single grid boxes, it is able to correct the bias and provide valid inference effectively.Given the fact that the signals at each grid box are much lower than those on a bigger spatial scale, our results are encouraging, indicating that our method can be served as a building block for further research on pooling boxes together to estimate β (Wang et al. 2021).
The corrected score functions also allow further extension of the fingerprinting by incorporating covariates in the scale parameter to account for non-stationarity and other possible influencing factors (Huang et al. 2016).Overall, the MCCS approach can be applied to GEV regressions with independent and temporally dependent EIVs that arise in many fields.
A few future directions are worth pursuing.Computationally, one concern is the numerical instability in estimation when the sample size is small.As many GEV regressions have a limited amount of observations, it would be useful to develop regularized estimation along the line of Firth (1993).For sandwich variance estimation, the "bread" part, or the inverse of the derivative of the corrected scores with respect to the parameters, could be poorly conditioned when the sample size is small.The regularized estimation might be helpful with this issue as well.Due to these computational issues, in our application to central China temperature extremes, we performed only one-signal analysis with the ALL signal.It is desirable to do two-signal analyses with both anthropogenic forcing and natural forcing at the same time.One possibility is to use a longer record of data, such as HadEX3 (Dunn et al. 2020) and CanESM5 (Swart et al. 2019).Methodologically, however, it would be interesting to extend the MCCS method to the spatial setting and restrict the scaling factors across grid boxes to be the same in one regional analysis (Zwiers et al. 2011).Incorporating spatial dependency to improve the effi-ciency of the MCCS method similar to the combined score equations (Wang et al. 2021) merits a full investigation.
Statistics and (2023) 33 :125 The real part of the three functions can be readily obtained by, , and (Z f 1 ) = Zc.The derivation of the case when ξ = 0 can be done in a similar way.

B.1 Argument for Proposition 1
Denote ψB,t (θ ) = ψB (Y t , W t , Z t , θ, t ).By Taylor series,  As the measurement error e t and t Zb,t are independent and identically distributed random vectors generated from N (0, ), the last equality is a result of Lemma 4 in Stefanski and Cook (1995).Hence, we finish the proof for Proposition 2.
Fig. 3 The top and bottom graphs present the location parameter α estimation for all the temperature extremes and the corresponding 90% bootstrap confidence intervals for the selected grid boxes Fig. 4 The top and bottom graphs present the scale parameter σ estimation for all the temperature extremes and the corresponding 90% bootstrap confidence intervals for the selected grid boxes Fig. 5 The top and bottom graphs present the scale parameter ξ estimation for all the temperature extremes and the corresponding 90% bootstrap confidence intervals for the selected grid boxes θn,B )] for m = 1, . . ., M, where W * m,b,t is generated for the mth bootstrap sample.5. Determine the covariance matrix of ψ m,B using M samples and multiplying the covariance matrix with n to get C n, ψ * .
DeRouen 2001; Westgate 2012).A major reason is that the middle matrix C * is under-estimated.To mitigate the potential underestimation of C n, ψ * when the sample size is small, and deal with the time dependency in extreme observations if that exists in the data, we propose a block bootstrap to account for the time dependency in ψ * B , rendering an estimate for the variance matrix of ψ * B .To create M bootstrap samples, we use the following bootstrap procedure for each m = 1, . . ., M and t = 1, . . ., n: 1. Calculate residuals by subtracting Z t α + W t β from Y t .2. Group residuals into non-overlapping blocks of samples in specified length with replacement.3. Obtain bootstrapped data Y m,t by adding the grouped residuals to the Z t α + W t β. 4. Calculate the ψ m,B

Table 1
Summary of simulation results based on 1000 replications Bias is the bias of point estimates; RMSE is the root mean square error; ESE is the empirical standard error; ASE is the average of the standard errors of the estimator; and CR is the coverage rate of 90% confidence interval

Table 2
Summary of simulation results based on 1000 replicationsBias is the bias of point estimates; RMSE is the root mean square error; ESE is the empirical standard error; ASE is the average of the standard errors of the estimator; and CR is the coverage rate of 90% confidence interval Fig. Twelve grid boxes in central China were selected for analysis