Efficient Volatility Estimation for Lévy Processes with Jumps of Unbounded Variation

Efficient Volatility Estimation for L\'evy Processes with Jumps of Unbounded Variation

Introduction

Statistical inference for stochastic processes based on high-frequency observations has attracted considerable attention in the literature for more than a decade. Among the many problems studied so far, none has received more attention than that of the estimation of the continuous or predictable quadratic variation of an It\^o semimartingale $X=\{X_{t}\}_{t\geq{}0}$. Specifically, if \[ X_t:=X^{c}_t+J_t:=\int_0^t b_t dt+\int_0^t \sigma_s dW_s+J_t,\qquad t\in[0,T], \] where $W=\{W_t\}_{t\geq{}0}$ is a Wiener process and $J=\{J_{t}\}_{t\geq{}0}$ is a pure-jump It\^o semimartingale, the estimation target is $IV_T=\int_0^T\sigma^2_s ds$. This quantity, also known as the integrated volatility or integrated variance of $X$, has many applications, especially in finance. When $X$ is observed at times $0=t_0<t_1<\ldots<t_n=T$, in the absence of jumps, an efficient estimator of $IV_T$ is given by the realized quadratic variation $\widehat{IV}_T=\sum_{i=1}^{n}(X_{t_i}-X_{t_{i-1}})^2$ in the so-called high-frequency asymptotic regime; i.e., when $\max_{i}(t_i-t_{i-1})\to{}0$ and $T\equiv t_n$ is fixed. In the presence of jumps, $\widehat{IV}_T$ is no longer even consistent for $IV_T$, instead converging to $IV_T+\sum_{s\leq T}(\Delta X_s)^2$, where $\Delta X_s:= X_s-X_{s^-}$ denotes the jump at time $s$. To accommodate jump behavior, several estimators have been proposed, among which the most well-known are the truncated realized quadratic variation and the multipower variations. We focus on the first class, which, unlike the second, is both rate- and variance-efficient, {in the Cramer-Rao lower bound sense,} when jumps are of bounded variation under certain additional conditions.

The truncated realized quadratic variation (TRQV), first introduced by Mancini in mancini2001 and mancini2004, is defined as

where $\varepsilon=\varepsilon_n>0$ is a tuning parameter converging to $0$ at a suitable rate. Above, $\Delta_{i}^{n}X:=X_{t_{i}}-X_{t_{i-1}}$ is the $i$-th increment of $(X_{t})_{t\geq 0}$ based on evenly spaced observations $X_{t_{0}},\ldots,X_{t_{n}}$ over a fixed time interval $[0,T]$ (i.e., $t_{i}=ih_n$ with $h_{n}=T/n$). It is shown in mancini2009non that TRQV is consistent when either the jumps have finite activity or stem from an infinite-activity L\'evy process. For a semimartingale model with L\'evy jumps of bounded variation, cont2011nonparametric showed that the TRQV admits a feasible central limit theorem (CLT), provided that $\varepsilon_n = ch_n^\beta$ with some\footnote{N.b.: in cont2011nonparametric, a different parameterization of the threshold parameter $\beta$ is used.} $\beta \in [\frac{1}{4-Y}, \frac{1}{2})$, where $Y\in [0,1)$ denotes the corresponding Blumenthal-Getoor index. In jacod2008asymptotic, consistency was established for a general It\^o semimartingale $X$, and a corresponding CLT is given when the jumps of $X$ are of bounded variation. In that case, the TRQV attains the optimal rate and asymptotic variance of $\sqrt{ h_n}$ and $2\int_{0}^{T}\sigma_s^4ds$, respectively.

In the presence of jumps of unbounded variation, the situation is notably different, and the problem is comparatively much less studied. In cont2011nonparametric, it is shown that when jumps stem from a L\'evy process with stable-like small-jumps of infinite variation, the TRQV estimator $\widehat{C}_n(\varepsilon)$ converges to $IV_T$ at a rate slower than $\sqrt{h_n}$. Further, in mancini2011speed it is shown that when the jump component $J$ is a $Y$-stable L\'evy process and $\varepsilon_n=h_n^\beta$ with $\beta\in(0,1/2)$, the decomposition $\widehat{C}_{n}( \varepsilon_n)-IV_T=\sqrt{h_n}Z_n+\mathcal{R}_n$ holds, where $Z_{n}$ converges stably in law to $\mathcal{N}(0,2\int_0^T\sigma_s^4 ds)$, while $\mathcal{R}_n$ is precisely of order $\varepsilon_n^{2-Y}$ in the sense that $\mathcal{R}_n=O_P(\varepsilon_n^{2-Y})$ and $\varepsilon_n^{2-Y}=O_P(\mathcal{R}_n)$ (which decays too slowly to allow for efficiency when $Y>1$). In AmorinoGloter20, the smoothed version of the TRQV estimator $\widehat{C}^{Sm}_{n}(\varepsilon)=\sum_{i=1}^{n}(X_{t_i}-X_{t_{i-1}})^2\varphi((X_{t_i}-X_{t_{i-1}})/\varepsilon)$ is considered\footnote{The authors in AmorinoGloter20, in fact, consider the more general estimation of $\int_0^T f(X_s) \sigma_s^2 ds$ for functions $f$ of polynomial growth, for which $IV_T$ is a special case.}, where $\varphi\in C^{\infty}$ vanishes in $\mathbb{R}\backslash{}(-2,2)$ and $\varphi(x)=1$ for $x\in(-1,1)$. In that case, using the truncation level $\varepsilon_n:=h_n^\beta$, it is shown that $\widehat{C}^{Sm}_{n}( \varepsilon_n)-IV_T=\sqrt{h_n}Z_n+\mathcal{R}_n$ with $\mathcal{R}_n$ such that $\varepsilon_n^{-(2-Y)}\mathcal{R}_n\to c_{Y}\int \varphi(u)|u|^{1-Y}du$, for a constant $c_Y$. By taking $\varphi$ such that $\int \varphi(u)|u|^{1-Y}du=0$, a “bias-corrected" estimator was considered under the additional condition that $Y<4/3$. Specifically, the resulting estimator is such that, for any $\tilde{\varepsilon}>0$, $\widehat{C}^{Sm}_{n}( \varepsilon_n)-IV_T=o_{P}(h_n^{1/2-\tilde\varepsilon})$, “nearly" attaining the optimal statistical error $O_{P}(h_n^{1/2})$. Unfortunately, the construction of such an estimator requires knowledge or accurate estimation of the jump intensity index $Y$, and no feasible CLT was proved when jumps are of unbounded variation even assuming $Y$ is known.

Until very recently, in the case of jumps of unbounded variation, the only rate- and variance-efficient estimator of the integrated volatility known in the literature was that proposed by Jacod and Todorov JacodTodorov:2014, under the additional condition that the jump intensity index $Y<3/2$ or the process $X$ has a symmetric jump component\footnote{The authors in JacodTodorov:2014 also constructed an estimator that is rate-efficient even in the presence of asymmetric jumps, but its asymptotic variance is twice as big as the optimal value $2\int_{0}^{T}\sigma_s^4ds$.}. Their estimator is based on locally estimating the volatility from the empirical characteristic function of the process' increments over disjoint time intervals shrinking to $0$, but still containing an increasing number of observations. It requires two debiasing steps, which are simpler to explain for a L\'evy process $X$ with symmetric $Y$-stable jump component $J$. The first debiasing step is meant to reduce the bias introduced when attempting to estimate $\log\mathbb{E}(\cos(\omega X_{h_n}/\sqrt{h_n}))$ with $\log\big\{\frac{1}{n}\sum_{i=1}^{n}\cos\big(\omega \Delta_i^n X/\sqrt{h_n}\big)\big\}$. The second debiasing step is aimed at eliminating the second term in the expansion $-2\log\mathbb{E}(\cos(\omega X_{h_n}/\sqrt{h_n})) /\omega^2=\sigma^2+2|\gamma|^{Y} \omega^{Y-2}h_n^{1-Y/2}+O(h_n)$, which otherwise diverges when multiplied by the optimal scaling $h_{n}^{-1/2}$. Using an extension of this approach, Jacod and Todorov were able to apply these techniques to a more general class of It\^o semimartingales in jacod2016efficient, even allowing any $Y<2$, though only rate-efficient, but not variance-efficient, estimators were ultimately constructed. It is stated therein that both rate- and variance-efficiency can be achieved for symmetric jump components for this more general class of semimartingales.

Recently, Mies Mies:2020 proposed an efficient estimation method for L\'{e}vy processes based on the generalized method of moments. Specifically, for some suitable functions $f_{1},f_{2},\dots,f_{m}$ and a scaling factor $u_{n}\rightarrow\infty$, Mies:2020 {proposed to search} for the parameter values $\widehat{\boldsymbol{\theta}}=(\widehat{\theta}_{1},\ldots,\widehat{\theta}_{m})$ such that

where $\widetilde{X}$ is the superposition of a Brownian motion and independent stable L\'{e}vy processes closely approximating $X$ in a certain sense. The distribution measure $\mathbb{P}_{\boldsymbol{\theta}}$ of $\widetilde{X}$ depends on some parameters $\boldsymbol{\theta}=(\theta_1,\ldots,\theta_m)$, one of which is the volatility $\sigma$ of $X$, and $\mathbb{E}_{\boldsymbol{\theta}}(\cdot)$ denotes the expectation with respect to $\mathbb{P}_{\boldsymbol{\theta}}$. Though in principle this method is efficient, it suffers from several important drawbacks. First, its finite-sample performance critically depends on the chosen moment functions $f_{1},\ldots,f_{m}$. Secondly, its implementation is computationally expensive and may lead to numerical issues since it involves solving a system of nonlinear equations. Moreover, in {addition} to the required use of a numerical solver to determine the {values} of $\widehat{\pmb{\theta}}$ in (ref), the expectations contained therein need to be numerically evaluated since the moments $\mathbb{E}_{\boldsymbol{\theta}}(f_{j}(u_{n}\Delta^{n}_{i}\widetilde{ X}))$ are typically not explicit. This fact introduces numerical errors that complicates its performance.

In this paper, we consider a new method to estimate the volatility of a L\'evy process under a semiparametric specification of the jump component. To the best of our knowledge, our method, together with of those in JacodTodorov:2014 and Mies:2020, are the only efficient methods to deal with L\'evy processes with jumps of unbounded variation. The idea is natural. We simply apply debiasing steps similar to those of JacodTodorov:2014 to the TRQV of mancini2009non. To give the heuristics as to why this strategy works, let us consider a small-time expansion of the truncated second moment $\mathbb{E}(X_{h_n}^2{\bf 1}_{\{|X_{h_n}|\leq{}\varepsilon_n\}})$ under the asymptotic regime $\varepsilon_n/\sqrt{h_n}\to\infty$. Using a variety of techniques, including a change of probability measure, Fourier-based methods, and small-large jump decompositions, we show the following expansion:

for certain constants $c_1,c_2\neq{}0$. Based on this expansion, it is easy to see that the rescaled bias ${\mathbb{E}[h_n^{-1/2}(\widehat{C}_{n}(\varepsilon)-\sigma^2T)]}$ satisfies

which suggests the necessity of the condition $h_n^{-1/2}\varepsilon_n^{2-Y} =o(1)$ for a feasible CLT for $\widehat C_n(\varepsilon)$ at the rate $\sqrt{h_n}$. However, together with the (necessary) condition $\varepsilon_n/\sqrt{h_n}\to\infty$, this can happen only if $Y<1$, and removal of the first terms in (ref) is necessary for efficient estimation, when jumps are of unbounded variation. To that end, note that for any $\zeta>1$,

where “h.o.t." means high order terms. The above formulas motivate the “bias-corrected" estimator

which is the essence of the debiasing procedure of JacodTodorov:2014. As we shall see, the story is more complicated than what the simple heuristics above suggest. Our main result shows that, for a class of L\'evy processes with stable-like small jumps and some additional mild conditions, the estimator (ref) is indeed rate- and variance-efficient provided that $Y<4/3$. Furthermore, if $Y<8/5$, a second bias correcting step will achieve both rate- and variance-efficiency (for the case $ 8/5\leq Y<2$, see remark at the end of Section (ref)). As mentioned above, our estimator provides a simple alternative method to those of JacodTodorov:2014 and Mies:2020. Furthermore, our Monte Carlo experiments indicate improved performance for the important class of CGMY L\'evy processes (cf. CarrGemanMadanYor:2002).

The rest of this paper is organized as follows. Section (ref) introduces the framework and assumptions as well as some known preliminary results from the literature. Section (ref) introduces the debiasing method and main results of the paper. Section (ref) illustrates the performance of our method via Monte Carlo simulations and compares it to the method in JacodTodorov:2014. The proofs are deferred to three appendix sections.

Setting and background

In this section, we introduce the model and main assumptions, some notation, and collect several facts needed from the literature. Our setting is similar to that of FigueroaLopezGongHoudre:2016 and is described here in detail for completeness. We consider a 1-dimensional L\'evy process $X=(X_{t})_{t\in\mathbb{R}_{+}}$ defined on a complete filtered probability space $(\Omega,\mathscr{F},(\mathscr{F}_{t})_{t\in\mathbb{R}_{+}},\mathbb{P})$ that can be decomposed as

where $W:=(W_{t})_{t\in\mathbb{R}_{+}}$ is a Wiener process and $J:=(J_{t})_{t\in\mathbb{R}_{+}}$ is an independent pure-jump L\'{e}vy process with L\'{e}vy triplet $(b,0,\nu)$. The L\'{e}vy measure $\nu$ is assumed to admit a density $s:\mathbb{R}_{0}\rightarrow\mathbb{R}_{+}$ of the form

Here, $C_{\pm}>0$, $Y\in(1,2)$, and $q:\mathbb{R}_{0}\rightarrow\mathbb{R}_{+}$ is a bounded Borel-measurable function satisfying the following assumptions:

These processes are sometimes called “stable-like L\'evy processes" and were studied in FigueroaLopezGongHoudre:2016,FigueroaLopezOlafsson:2019. In simple terms, condition (i) above says that the small jumps of the L\'evy process $X$ behave like those of a $Y$-stable L\'evy process with L\'evy measure

The condition $Y\in(1,2)$ implies that $J$ has unbounded variation in that sense that \[ \sum_{i=1}^{n}|J_{t_i}-J_{t_{i-1}}|\to{}\infty \quad\textnormal{a.s.,} \] as the partition $0=t_0<t_1<\dots<t_n=T$ is such that $\max\{t_i-t_{i-1}\}\to{}0$.

As in FigueroaLopezGongHoudre:2016, it will be important for our analysis to apply a density transformation technique Sato:1999 to “transform" the process $J$ into a stable L\'evy process. Concretely, we can change the probability measure from $\mathbb{P}$ to another locally absolutely continuous measure $\widetilde{\mathbb{P}}$, under which $W$ is still a standard Brownian motion independent of $J$, but, under $\widetilde{\mathbb{P}}$, $J$ has L\'evy triplet $(\tilde{b},0,\tilde{\nu})$, where $\widetilde{\nu}(dx)$ is given as in (ref) and $\widetilde{b}:=b+\int_{0<|x|\leq 1} x(\tilde{\nu}-\nu)(dx)$. For future reference, we recall that (see FigueroaLopezGongHoudre:2016), for any $t\in\mathbb{R}_{+}$,

where $\varphi(x)=-\ln q(x)$. Under $\widetilde{\mathbb{P}}$, the centered process $S:=(S_{t})_{t\in\mathbb{R}_{+}}$, given by

is a strictly $Y$-stable process with its scale, skewness, and location parameters given by $[(C_{+}+C_{-})\Gamma(-Y)|\cos(\pi Y/2)|]^{1/Y}$, $(C_{+}-C_{-})/(C_{+}+C_{-})$, and $0$, respectively. We recall the following bounds for the tail probabilities and density of $S_1$ under $\widetilde\mathbb{P}$, which hereafter we denote $p_S$:

for a constant $\widetilde{K}<\infty$. The inequalities above can be derived from the asymptotics in Section 14 of Sato:1999 (see FigueroaLopezGongHoudre:2016 for more details).

It will be useful later on to express the processes $S=(S_{t})_{t\in\mathbb{R}_{+}}$ and $U:=(U_{t})_{t\in\mathbb{R}_{+}}$ in terms of the jump measure $N(dt,dx)$ of the process $J$ and its compensator $\widetilde{\nu}(dx)dt$ under $\widetilde{\mathbb{P}}$. Specifically, with $\widetilde{N}(dt,dx):=N(dt,dx)-\widetilde{\nu}(dx)dt$, we have:

where

Under $\widetilde \mathbb{P}$, $S^{+}:=(S^{+}_{t})_{t\in\mathbb{R}_{+}}$ and $-S^{-}:=(-S^{-}_{t})_{t\in\mathbb{R}_{+}}$ are independent one-sided $Y$-stable processes with scale, skewness, and location parameters given by $[C_{\pm}|\Gamma(-Y)\cos(\pi Y/2)|]^{{1/Y}}$, $1$, and $0$, respectively, so that (cf. Zolotarev:1986)

Main results

In this section, we construct an efficient estimator for the volatility parameter $\sigma^2$ based on the well-studied estimator TRQV (ref). Throughout, we assume the process $X=\{X_{t}\}_{t\geq{}0}$ is sampled at $n$ evenly spaced observations, $X_{t_1}, X_{t_2}, \ldots, X_{t_n}$, during a fixed time interval $[0,T]$, where for $i=0,\dots,n$, $t_{i}=t_{i,n}= ih_{n}$ with $h_n=T/n$. As usual, we define the increments of the process as $\Delta_{i}^{n}X:=X_{t_{i}}-X_{t_{i-1}}$, $i=1,\dots,n$ and, for simplicity, assume that $T=1$. In what follows, we sometimes approximate {$X$} by a L\'{e}vy process {$\widetilde{X}:=(\widetilde{X}_{t})_{t\in\mathbb{R}_{+}}$} with characteristic triplet $(0,\sigma^{2},\widetilde{\nu})$, where $\widetilde{\nu}$ is defined in (ref). Note that, under $\widetilde\mathbb{P}$,

where $\stackrel{\mathcal{D}}{=}$ denotes equality in law, and, for notational simplicity, we assume $\widetilde{X}$ has the same distribution under both $\widetilde \mathbb{P}$ and $\mathbb{P}$.

As mentioned in the introduction, the TRQV estimator $\widehat{C}_n(\varepsilon)$ is not efficient since it possesses a bias that vanishes at a rate slower than $n^{-1/2}$, the rate at which the “centered” TRQV, \[ \overline{C}_{n}(\varepsilon)=\sum_{i=1}^{n}\mathopen{}\mathclose\bgroup\originalleft\{\big(\Delta_{i}^{n}X\big)^{2}{\bf 1}_{\{|\Delta_{i}^{n}X|\leq\varepsilon\}}- \mathbb{E}\big(X^2_{h_n}{\bf 1}_{\{|X_{h_n}|\leq\varepsilon\}}\big)\aftergroup\egroup\originalright\}, \] converges to Gaussianity. To overcome this, our idea is to apply the debiasing procedure of JacodTodorov:2014 to the TRQV. Unlike the referred work, our procedure is simpler, as it does not require an extra debiasing step to correct the nonlinear nature of the logarithmic transformation employed therein nor does it require a symmetrization step to deal with asymmetric L\'evy measures.

Before constructing our estimator, we first establish an infeasible CLT for an approximately centered estimator, where the truncated moment $\mathbb{E}\big(X^2_{h_n}{\bf 1}_{\{|X_{h_n}|\leq\varepsilon\}}\big)$ is replaced by the quantity $\widetilde\mathbb{E}\big(\widetilde{X}^2_{h_n}{\bf 1}_{\{|\widetilde{X}_{h_n}|\leq\varepsilon\}}\big)$. Below and throughout the rest of the paper, we use the usual notation $a_{n}\ll b_n$, whenever $a_n/b_n\to0$ as $n\to\infty$.

The role of the next result is twofold. First, it establishes the asymptotic normality of $\widehat C_n(\varepsilon)$ using a fully specified (unobservable) centering quantity $A(\varepsilon,h)$ rather than the inexplicit centering term $\mathbb{E}\big(\widetilde{X}^2_{h_n}{\bf 1}_{\{|\widetilde{X}_{h_n}|\leq \varepsilon_n\}}\big)$ of Proposition (ref) (the quantity $A(\varepsilon,h)$ will subsequently be estimated and removed through our proposed debiasing method). {Secondly}, subject to this centering, it provides the joint asymptotic behavior of $\widehat C_n(\varepsilon)$ and the difference $\widehat C_n(\zeta\varepsilon) -\widehat C_n(\varepsilon)$, for some $\zeta>1$; this is the main technical result from which we deduce the efficiency of our debiased estimator.

We are now in a position to introduce our proposed estimator. Let $\zeta_1,\zeta_2>1$, and set

The next theorem is the main result of the paper. It establishes the rate- and variance-efficiency of the two-step debiased estimator $ \widetilde C_n ''(\varepsilon, \zeta_2,\zeta_1)$ provided $Y<8/5$.

It is customary to use power thresholds of the form $\varepsilon_n = c_0 h_n^\beta$, where $c_0>0$ and $\beta>0$ are some constants\footnote{ Though, it is shown in gong2021 that under the setting of Section (ref), the optimal threshold $\varepsilon_n^{*}$ is such that $\varepsilon_n^{*}\sim\sqrt{(2-Y)\sigma^{2}h_n\ln(1/h_n)}$.}. In that case, the assumption $h_n^{\frac{4}{8+Y}}\ll \varepsilon_n\ll h_n^{\frac{1}{4-Y}\vee \frac{1}{2+Y/2}}$ in Theorem (ref) becomes

Note that the value of $\beta=5/12$ satisfies the above constraint for any value {$Y$} of the possible range considered by Theorem (ref).

Monte Carlo performance for CGMY L\'evy processes

{In this section, we study the performance of the two-step debiasing procedure introduced in the previous section to the case of a L\'evy process with a CGMY jump component $J$ (cf. CarrGemanMadanYor:2002). Specifically, we work with simulated data from the model (ref) where $\{J_t\}_{t\geq{}0}$ is a CGMY process, independent of the Brownian motion $\{W_{t}\}_{t\geq{}0}$, with L\'evy measure having a $q$-function, in the notation of (ref), of the form:

and $C_{+}=C_{-}=C$. Thus, the conditions of Assumption (ref) are satisfied with $\alpha_{+}=-M$ and $\alpha_{-}=G$. We adopt the parameter setting

which are similar to those used in FigueroaLopezOlafsson:2019\footnote{FigueroaLopezOlafsson:2019 considers the asymmetric case $\nu(dx)=C_{{\rm sgn} (x)}\bar{q}(x)|x|^{-1-Y}\,dx$ with $C_+=0.015$ and $C_-=0.041$. Here, we take $C=(C_++C_-))/2$ in order to simplify the simulation of the model. The parameter values of $C_+$, $C_-$, $G$, and $M$ used in FigueroaLopezOlafsson:2019 were taken from Kawai, who calibrated the tempered stable model using market option prices.}, and take $\sigma=0.2, 0.4$ and $Y=1.25, 1.35, 1.5, 1.7$, respectively. We take $T=1$ year and $n=252(6.5)(12)$, which corresponds to a frequency of $5$ minutes (assuming $252$ trading days and a $6.5$-hour trading period each day).

In a fashion similar to JacodTodorov:2014, for the threshold $\varepsilon = c_0h^\beta$, we take {$c_0 = \sigma_{BV}$, where \[ \sigma_{BV}^2= \frac{\pi}{2} \sum_{i=2}^n |\Delta_{i-1}^nX||\Delta_i^nX|, \] which} is the standard Bipower variation estimator of $\sigma^2$ introduced by Barndorff. {For the value of $\beta$ we take $\beta = \frac{5}{12}$, which, as mentioned above, satisfies the condition (ref) for any $Y\in(1,8/5)$.} We compare the performance of the following estimators:

We further compare the simulated performance of the above estimators to the estimator implemented in the Monte Carlo study in JacodTodorov:2014. Specifically, we use the equation (5.3) in the paper JacodTodorov:2014:

where $\widehat C_{\text{JT}}$ denotes their “nonsymmetrized" two-step debiased estimator. For the parameter settings, we take $$ \zeta=1.5, \quad u_{n}=(\ln(1/h_{n}))^{-1/30}/\sigma_{BV}, \quad p_0=0.5, $$ where the values of $\zeta$ and $u_n$ were those suggested by JacodTodorov:2014 and the value of $p_0$ was chosen for favorable estimation performance. Note that, since a L\'evy model has constant volatility, it is not necessary to localize the estimator and, hence, we treat the 1-year data as one block, which corresponds to $k_n=252(6.5)(12)$ in the notation of JacodTodorov:2014.

The simulation results are summarized in Tables (ref)-(ref). We report the sample means, standard deviations {(SDs)}, the average and {SD} of relative errors, the MSEs, and median of absolute deviations (MADs) of each of the four estimators described above, based on 1000 simulations.

When $\sigma=0.2$ and $Y=1.25$ or $1.35$, or when $\sigma=0.4$ and $Y=1.25$, $1.35$, or $1.5$, $\widetilde C_{n,nb} ''$ significantly outperforms $\widetilde C_{n,pb}'$. Also, $\widetilde C_{n,nb} ''$ has superior performance compared to $\widehat C_{\text{JT}, 53}$ as measured by MSE and MAD. Table (ref) also shows that, though when $\sigma=0.2$ and $Y=1.5$, $\widetilde C_{n,nb} ''$ has larger MSE and MAD than $\widehat{C}_n$ and $\widetilde C_{n,pb}'$, it still performs better than $\widehat C_{\text{JT}, 53}$: the MSE and MAD of $\widetilde C_{n,nb} ''$ in this case are approximately {53% and 73%} of those of $\widehat C_{\text{JT}, 53}$, respectively. Tables (ref) and (ref) show that, when $Y=1.7$, which is not covered by our theoretical framework (c.f. Remark (ref)), $\widetilde C_{n,nb} ''$ has larger MSE and MAD than $\widehat C_{\text{JT}, 53}$. Overall, we conclude that our debiasing procedure outperforms $\widehat C_{\text{JT}, 53}$ when $Y \in (1, 1.5]$}.