\section{Introduction}
\label{sect1}
Late-time cosmic acceleration is an inevitable ingredient of
our Universe directly supported by cosmological observations \cite{HZTEAM,
SCP}. Other observations such as cosmic microwave background (CMB),
baryonic aucostic oscillations (BAO), sloan digital sky survey \cite{cmb1, cmb2,
bao1, bao2, sdss1, sdss2} and many others support this fact indirectly. The
globular cluster reveals that the age of certain objects in the
Universe is larger than the age of Universe estimated in standard model with
normal matter. The only known resolution of the puzzle is provided by
invoking cosmic acceleration at late-times. Although, there are many ways to
explain the accelerating expansion of the Universe (e.g. by adding a source
term in the matter part of Einstein field equations, by modifying the
geometry or by invoking inhomogeneity), the inclusion of a source term with
large negative pressure dubbed ``\textit{dark energy}" (DE) \cite{DE1, DE2, DE3,
DE4, DE5, DE6, DE7, DE8, DE9, DE10} is widely accepted to the theorists.
However, a promising candidate for dark energy (cosmological constant $\Lambda $) is
under scrutiny. A wide varieties of dark energy candidate have been proposed
in the past few years such as cosmological constant \cite{DE3, carollCC,
padmaCC, peeblesCC}, slowly rolling scalar field \cite{quint1, quint2,
quint3, quint4, quint5}, phantom field \cite{phant1, phant2, phant3, phant4,
phant5, phant6}, tachyon field \cite{tachy1, tachy2,alam2017} and chaplygin gas \cite%
{CG1, CG2} etc. (See \cite{DE1, DEREV1, DEREV2, DEREV3} for a detailed list).

The modifications of Einstein's theory of gravity not only account the cosmic acceleration but also resolves many standard problems
such as singularity problem, the hierarchy problem, quantization and
unification with other theories of fundamental interactions. Massive
gravity, Gauss-Bonnet gravity, $f(R)$, $f(T)$, $f(R,T)$ gravities,
Chern-Simon gravity, Galileon gravity are name a few among the various
alternative theories proposed in the past few years. Modified gravity can
also provide unified description of the early-time inflation with that of
late-time cosmic acceleration and dark matter (DM). Traditionally, all these
modifications invoke extra degrees of freedom non-minimally coupled to
matter in the Einstein frame. Generally, it is believed that late-time
acceleration requires the presence of dark energy or the extra degrees of
freedom. Recently, Berezhiani et al. \cite{KHOURY2016} discussed a third
possibility which requires neither any exotic matter nor the large scale
modifications of gravity. They showed that the interaction between the
normal matter components, namely, the dark matter and the baryonic matter (BM)
can also provide the late-time acceleration in Jordan frame. In context of the
coupling, the stability criteria disfavors the conformal coupling while the
maximally disformal coupling can give rise to late-time cosmic acceleration
in Jordan frame but no acceleration in the Einstein frame. Extending the
work of Ref. \cite{KHOURY2016}, Agarwal et al. \cite{ARSSA} have further
investigated the cosmological dynamics of the model obtained by
parameterizing the coupling function. Also, they have shown that the model
exhibits the sudden future singularity which can be resolved by taking a
more generalized parametrization of the coupling function.

In this paper, we shall consider two forms of parametrizations. The first parametrization shows the sudden future singularity, which can be pushed into far future if we consider the second parametrization. We further investigate the models using statefinder and $Om$ diagnostics. The paper is organized as follows. Section \ref{sect:2} is devoted to the basic equations of the models. In section \ref{sec:obs}, we put the observational constraints on the model parameters. The detailed analysis of statefinder and $Om$ diagnostics are presented in sections \ref{sec:state} and \ref{sec:om}, respectively. We conclude our results in section \ref{sec:conc}.
\section{Field equations}
\label{sect:2}
The scenario of the interaction between dark matter and the baryonic
matter as described briefly in Refs. \cite{KHOURY2016} and \cite{ARSSA} in a spatially
flat Freidmann-Lemaitre-Robartson-Walker (FLRW) background
\begin{equation}
ds^{2}=-dt^{2}+a^{2}(t)\left( dx^{2}+dy^{2}+dz^{2}\right) \text{,}  \label{1}
\end{equation}%
yield the field equations,
\begin{equation}
3H^{2}=8\pi G\left( \Lambda _{DM}^{4}\sqrt{\frac{X}{X_{eq}}}\left( \frac{%
a_{eq}}{a}\right) ^{3}-P+QR^{3}\tilde{\rho}_{b}\right) \text{,}  \label{2}
\end{equation}%
and
\begin{equation}
2\frac{\ddot{a}}{a}+H^{2}=-8\pi G(P+P_{b})\text{,}  \label{3}
\end{equation}%
where $Q$ and $R$ are two arbitrary coupling functions.

The Einstein frame metric couples to Jordan frame metric such that $\sqrt{-%
\tilde{g}}=QR^{3}\sqrt{-g}$ ($\tilde{g}$ is the determinant of Jordan frame
metric $\tilde{g}_{\mu \nu }$, constructed from the Einstein frame metric $%
g_{\mu \nu }$) and $X=-g^{\mu \nu }\partial _{\mu }\Theta \partial _{\nu
}\Theta $, $\Theta $ being the dark matter field. The quantities $P$ and $P_{b}$ are
pressures of DM and BM in Einstein frame respectively which are related by $%
P_{b}\equiv QR^{3}\tilde{P}_{b}$. $\tilde{P}_{b}$ and $\tilde{\rho}_{b}$ are
the pressure and density of BM in Jordan frame and are related to BM density
in Einstein frame by the relation $\rho _{b}=QR^{3}\left( \tilde{\rho}%
_{b}\left( 1-2X\frac{Q,_{X}}{Q}\right) +6X\frac{R,_{X}}{R}\tilde{P}%
_{b}\right) $. For the detailed derivation of field equations, see \cite%
{KHOURY2016} and \cite{ARSSA}. Here, we shall note that all the quantity with a
tilde above are in Jordan frame and without tilde are in Einstein frame.

The Jordan frame scale factor dubbed \textit{physical scale factor} is
related a scale factor of Einstein frame  as
\begin{equation}
\tilde{a}=Ra\text{.}  \label{4}
\end{equation}
\qquad Here, we consider the same maximally disformal coupling of BM and DM for which $Q=1$
throughout the evolution and $R=1$ in the early Universe that grows
sufficiently fast such that the physical scale factor $\tilde{a}$ in Jordan
frame experiences acceleration. The conformal coupling is disfavored by the stability criteria \cite{KHOURY2016}. One needs to specify the coupling function $%
R(a)$ to proceed further or equivalently, $a$ can be parametrized in terms
of physical scale factor $\tilde{a}$. The two parametrizations are
\\
(1) Model 1:  $a(\tilde{a})=\tilde{a}+\alpha \tilde{a}^{2}+\beta \tilde{a}^{3}$, where
$\alpha $ and $\beta $ are two model parameters.
\\
(2) Model 2:  $a(\tilde{a})=\tilde{a}e^{\alpha \tilde{a}}$, in this case, only $\alpha $ is a model parameter.

By expanding the functional $\tilde{a}e^{\alpha \tilde{a}}$ in Taylor
series, the first parametrization $\tilde{a}+\alpha \tilde{a}^{2}+\beta
\tilde{a}^{3}$ can be recovered by substituting $\beta=\alpha^2/2$. Agarwal et al. \cite{ARSSA} have studied
various features of the model 1 and constrained the parameters $\alpha $ \& $%
\beta $ by employing the $\chi ^{2}$ analysis using $H(z)+SN+BAO$
datasets. Extending the analysis, we further study some more physical
characteristics of the models 1 and 2 such as the statefinder and $Om$ diagnostics, and also put observational constraints on the parameter $\alpha $ of model 2.

We also need to express the cosmological parameters in terms of redshifts in
both the Einstein frame and Jordan frame which are defined as
\begin{equation}
\tilde{a}=\frac{\tilde{a}_{0}}{1+\tilde{z}}\text{ , \ \ }a=\frac{a_{0}}{1+z}%
\text{,}  \label{5}
\end{equation}%
For both the parametrizations, $\tilde{a}_{0}=1$ but $a_{0}=1+\alpha+\beta\neq 1$ (model 1) and $a_{0}=e^\alpha \neq 1$ (model 2).

For model 1, the explicit expressions for the Hubble and deceleration parameters in Jordan frame are obtained as
\begin{equation}
\tilde{H}(\tilde{a})=\frac{\dot{\tilde{a}}}{\tilde{a}}=\tilde{H}_{0}\frac{%
(1+\alpha +\beta )^{\frac{1}{2}}(1+2\alpha +3\beta )}{\tilde{a}^{\frac{3}{2}}%
\left[ 1+\alpha \tilde{a}+\beta \tilde{a}^{2}\right] ^{\frac{1}{2}}\left[
1+2\alpha \tilde{a}+3\beta \tilde{a}^{2}\right] }\text{,}  \label{eq:H1}
\end{equation}
and%
\begin{equation}
\tilde{q}(\tilde{a})=-\frac{\tilde{a}\ddot{\tilde{a}}}{\dot{\tilde{a}}^{2}}=%
\frac{1}{2}\left( \frac{1+2\alpha \tilde{a}+3\beta \tilde{a}^{2}}{1+\alpha
\tilde{a}+\beta \tilde{a}^{2}}\right) +\left( \frac{2\tilde{a}(\alpha
+3\beta \tilde{a})}{1+2\alpha \tilde{a}+3\beta \tilde{a}^{2}}\right) \text{.}
\label{eq:q1}
\end{equation}%
Using Eq. (\ref{5}), above expressions can be written in terms of redshift $\tilde{z}$ as
\begin{equation}
\tilde{H}(\tilde{z})=\tilde{H}_{0}\frac{(1+\alpha +\beta )^{\frac{1}{2}%
}(1+2\alpha +3\beta )\left( 1+\tilde{z}\right) ^{\frac{9}{2}}}{\left[ \left(
1+\tilde{z}\right) ^{2}+\alpha \left( 1+\tilde{z}\right) +\beta \right] ^{%
\frac{1}{2}}\left[ \left( 1+\tilde{z}\right) ^{2}+2\alpha \left( 1+\tilde{z}%
\right) +3\beta \right] }\text{,}  \label{a3}
\end{equation}%
and%
\begin{equation}
\tilde{q}(\tilde{z})=\frac{\left[ \left( 1+\tilde{z}\right) ^{2}+2\alpha
\left( 1+\tilde{z}\right) +3\beta \right] ^{2}+4\left[ \alpha \left( 1+%
\tilde{z}\right) +3\beta \right] \left[ \left( 1+\tilde{z}\right)
^{2}+\alpha \left( 1+\tilde{z}\right) +\beta \right] }{2\left[ \left( 1+%
\tilde{z}\right) ^{2}+\alpha \left( 1+\tilde{z}\right) +\beta \right] \left[
\left( 1+\tilde{z}\right) ^{2}+2\alpha \left( 1+\tilde{z}\right) +3\beta %
\right] }\text{,}  \label{a4}
\end{equation}%
together with the effective equation of state (EOS) parameter given by
\begin{equation}
\tilde{w}_{eff}(\tilde{z})=\frac{\alpha \left( 5+6\alpha +5\tilde{z}\right)
(1+\tilde{z})^{2}+\beta (14+23\alpha +14\tilde{z})(1+\tilde{z})+18\beta ^{2}%
}{3\{(1+\tilde{z})^{2}+\alpha (1+\tilde{z})+\beta \}\{(1+\tilde{z}%
)^{2}+2\alpha (1+\tilde{z})+3\beta \}}\text{.}  \label{a5}
\end{equation}
Similarly, for model 2, we obtain the expressions for the Hubble and deceleration parameters in Jordan-frame as%
\begin{equation}
\tilde{H}(\tilde{a})=\frac{\tilde{H}_{0}(1+\alpha )e^{\frac{3}{2}\alpha }}{%
\tilde{a}^{\frac{3}{2}}\left[ 1+\alpha \tilde{a}\right] ^{\frac{1}{2}}e^{%
\frac{3}{2}\alpha \tilde{a}}},  \label{eq:H2}
\end{equation}%
and%
\begin{equation}
\tilde{q}(a)=\frac{\left( 1+\alpha \tilde{a}\right) ^{2}+2\tilde{a}(2\alpha
+\alpha ^{2}\tilde{a})}{2\left( 1+\alpha \tilde{a}\right) }\text{.}
\label{eq:q2}
\end{equation}%
with the help of Eq. (\ref{5}), we obtain%
\begin{equation}
\tilde{H}(\tilde{z})=\frac{\tilde{H}_{0}(1+\alpha )e^{\frac{3}{2}\alpha
}\left( 1+\tilde{z}\right) ^{\frac{5}{2}}}{\left[ \left( 1+\tilde{z}\right)
+\alpha \right] e^{\frac{3}{2}\displaystyle\frac{\alpha }{(1+\tilde{z})}}}%
\text{,}  \label{b3}
\end{equation}%
and%
\begin{equation}
\tilde{q}(\tilde{z})=\frac{1+\tilde{z}^{2}+6\alpha +3\alpha ^{2}+\left(
2+6\alpha \right) \tilde{z}}{2\left( 1+\tilde{z}\right) \left( 1+\tilde{z}%
+\alpha \right) }\text{.}  \label{b4}
\end{equation}
The effective equation of state is then given by
\begin{equation}
\tilde{w}_{eff}(\tilde{z})=\frac{5\alpha (1+\tilde{z})+3\alpha ^{2}}{3(1+%
\tilde{z})\left[ (1+\tilde{z})+\alpha \right] }\text{ .}  \label{b5}
\end{equation}
In both the models, another parameter $H_0$ will come in the expressions of $H(z)$, see Eqs. (\ref{a3}) and (\ref{b3}). But, here we focus on the parameters of underlying parametrizations. The first parametrization (model 1) consist of two model parameters (i.e. $\alpha$ and $\beta$) while the second parametrization (model 2) consists of a single model parameter (i.e. $\alpha$). Now we are in position to put the observational constraints on the parameters of model 2 in the following section.

Before proceeding to next section, we consider The DGP model as \cite{dgp}:
\begin{equation}
\label{eq:DGP}
\frac{H(z)}{H_0} = \left[ \left(\frac{1-\Omega_{0 {\rm m}}}{2}\right)+\sqrt{\Omega_{0 {\rm m}} (1+z)^3+
\left(\frac{1-\Omega_{0 {\rm m}}}{2}\right)^2} \right]\,\,
\end{equation}
where $H_0$ and $\Omega_{0 {\rm m}}$ are the present values of Hubble parameter and energy density parameter of matter.
\section{Observational constraints}
\label{sec:obs}
We have already mentioned that the model 1 consists of two parameters, namely, $\alpha$ and $\beta$ which were constrained in Ref. \cite{ARSSA}. In our analysis, we shall use  their best-fit values given as $\alpha =-0.102681$ \& $ \beta =-0.078347$. In this section, we put the constraints on parameters of model 2 by employing the same procedure as in \cite{ARSSA}.

One can use the total likelihood to constrain the parameters $\alpha $ and $H_0$ of model 2. The total likelihood function for a joint analysis can be defined as
\begin{equation}
\mathcal{L}_{tot}(\alpha , H_0 )=e^{-\frac{\chi _{tot}^{2}(\alpha , H_0 )}{%
2}}\text{,~~~~ where }\chi _{\mathrm{tot}}^{2}=\chi _{\mathrm{Hub}}^{2}+\chi _{%
\mathrm{SN}}^{2}+\chi _{\mathrm{BAO}}^{2}\text{.}  \label{o1}
\end{equation}%
Here, $\chi _{\mathrm{Hub}}^{2}$ denotes the chi-square for the Hubble dataset,
$\chi _{\mathrm{SN}}^{2}$ represents the Type Ia supernova and $\chi _{%
\mathrm{BAO}}^{2}$ corresponds to the BAO. By minimizing the $\chi
_{\mathrm{tot}}^{2}$ , we obtain the best-fit value of $\alpha $ and $H_0$. The
likelihood contours are standard i.e. confidence level at 1$\sigma $ and 2$%
\sigma $ are $2.3$ and $6.17$, respectively in the 2D plane.

First, we consider 28 data points of $H(z)$ used by Farooq and Ratra \cite{Farooq:2013hq} in the redshift range $0.07\leq z\leq 2.3$, and use $H_{0}=67.8\pm 0.9~Km/S/Mpc$ \cite{planck2015}. The $\chi ^{2}$, in this case, is defined as
\begin{equation}
\chi _{\mathrm{Hub}}^{2}(\theta )=\sum_{i=1}^{29}\frac{\left[ h_{\mathrm{th}%
}(z_{i},\theta )-h_{\mathrm{obs}}(z_{i})\right] ^{2}}{\sigma _{h}(z_{i})^{2}}%
\,,  \label{o2}
\end{equation}%
where $h=H/H_{0}$ represents the normalized Hubble parameter, $h_{\mathrm{obs}}$ and $h_{\mathrm{th}}$ are the observed and theoretical values of normalized Hubble parameter and $\sigma _{h}=\left(
\frac{\sigma _{H}}{H}+\frac{\sigma _{H_{0}}}{H_{0}}\right) h$. The quantities $\sigma _{H}$ and $\sigma _{H_{0}}$ designate the errors associated with $H$ and ${H_{0}}$,
respectively.

Second, we use $580$ data points from Union2.1 compilation data \cite{Suzuki:2011hu}. The corresponding $\chi ^{2}$ is given as
\begin{equation}
\chi _{\mathrm{SN}}^{2}(\mu _{0},\theta )=\sum_{i=1}^{580}\frac{\left[ \mu
_{th}(z_{i},\mu _{0},\theta )-\mu _{obs}(z_{i})\right] ^{2}}{\sigma _{\mu
}(z_{i})^{2}}\,,  \label{o3}
\end{equation}%
where $\mu _{obs}$, $\mu _{th}$ are the observed, theoretical distance
modulus and $\sigma _{\mu }$ is the uncertainty in the distance modulus, and $\theta $ is an arbitrary parameter. The distance modulus $\mu (z)$ is an observed quantity and related to luminosity distance $D_{L}(z)=(1+z)\int_{0}^{z}\frac{H_{0}dz^{\prime }}{%
H(z^{\prime })}$ as $\mu (z)=m-M=5\log D_{L}(z)+\mu _{0}$; $m$ and $M$ being the apparent and absolute magnitudes of the supernovae, and $\mu _{0}=5\log
\left( \frac{H_{0}^{-1}}{\mathrm{Mpc}}\right) +25$ is a nuisance parameter

Finally, we consider BAO data. The corresponding chi-square ($\chi _{\mathrm{BAO}}^{2}$) is defined by
\cite{Giostri:2012ek}:
\begin{equation}
\chi _{\mathrm{BAO}}^{2}=Y^{T}C^{-1}Y\,,  \label{o4}
\end{equation}%
where
\begin{equation}
Y=\left( \begin{array}{c}
        \frac{d_A(z_\star)}{D_V(0.106)} - 30.95 \\
        \frac{d_A(z_\star)}{D_V(0.2)} - 17.55 \\
        \frac{d_A(z_\star)}{D_V(0.35)} - 10.11 \\
        \frac{d_A(z_\star)}{D_V(0.44)} - 8.44 \\
        \frac{d_A(z_\star)}{D_V(0.6)} - 6.69 \\
        \frac{d_A(z_\star)}{D_V(0.73)} - 5.45
        \end{array} \right)\,,
\end{equation}
and inverse covariance matrix ($C^{-1}$), values $\frac{%
d_{A}(z_{\star })}{D_{V}(Z_{BAO})}$ are taken into account as in \cite%
{Blake:2011en, Percival:2009xn, Beutler:2011hx, Jarosik:2010iu,
Eisenstein:2005su, Giostri:2012ek}, and $z_{\star }\approx 1091$ is the
decoupling time, $d_{A}(z)$ is the co-moving angular-diameter distance and $%
D_{V}(z)=\left( d_{A}(z)^{2}z/H(z)\right) ^{1/3}$ is the dilation scale.

For model 2, we use an integrated datasets of $H(z)+SN+BAO$, and corresponding likelihood contour at 1$\sigma $ and 2$\sigma $ confidence levels are shown in Fig. \ref{fig:cont}. The best-fit values of the model parameters are obtained as $
\protect\alpha =-0.3147$ and $H_{0}=66.84~Km/S/Mpc$.
\begin{figure}[tbph]
\begin{center}
{\includegraphics[width=2.3in,height=2.3in,angle=0]{ah.pdf}}
\end{center}
\caption{This figure shows the 1$\protect\sigma $ (dark shaded) and 2$\protect%
\sigma $ (light shaded) likelihood contours in $\protect\alpha -H_{0}$
plane. The figure corresponds to joint datasets of $H(z)+SN+BAO$. A black dot
represents the best-fit values of the model parameters which are found to be $
\protect\alpha =-0.3147$ and $H_{0}=66.84~Km/S/Mpc$.}
\label{fig:cont}
\end{figure}
The normalized Hubble parameter and the effective EOS ($w_{eff}$) are plotted for both the models with their respective best-fit values of the parameters, and are shown in Fig. \ref{fig:hw}. Both the models represent phantom behavior in future. Model 1 exhibits the sudden singularity in near future while this kind of singularity has been delayed and pushed into far future in case of model 2 that is displayed in the right panel of Fig. \ref{fig:hw}. Fig. \ref{fig:mu} exhibits the error bar plots for models 1 and 2 with $H(z)$ and $SN$ datasets which shows that both the models are consistent with the observations.
\begin{figure}[tbph]
\begin{center}
\begin{tabular}{cc}
{\includegraphics[width=2.3in,height=2.3in,angle=0]{zh.pdf}} & {%
\includegraphics[width=2.3in,height=2.3in,angle=0]{zw.pdf}}%
\end{tabular}%
\end{center}
\caption{The figure shows the evolution of normalized Hubble parameter ($H(z)/H_0$)
and effective EOS ($w_{eff}$) versus redshift ($z$). The
dotted, dashed and dot-dashed lines correspond to DGP, models 1 and 2,
respectively. We use best-fit values of the model parameters. The horizontal line is the phantom divide line.}
\label{fig:hw}
\end{figure}
\begin{figure}[tbph]
\begin{center}
\begin{tabular}{cc}
{\includegraphics[width=2.3in,height=2.3in,angle=0]{Hz-error.pdf}} & {%
\includegraphics[width=2.3in,height=2.3in,angle=0]{sn-error.pdf}}%
\end{tabular}%
\end{center}
\caption{This figure exhibits the error bars of $H(z)$ (left) and $SN$ (right) datasets. In both the panels, the dashed and dot-dashed lines show the best-fitted behavior for models 1 and 2, respectively.}
\label{fig:mu}
\end{figure}

The above discussions show the validation of our models corresponding to the
observations. In the following sections we shall employ different diagnostics for underlying models.
\begin{figure}[tbph]
\begin{center}
\begin{tabular}{cc}
{\includegraphics[width=2.3in,height=2.3in,angle=0]{sr.pdf}} & {%
\includegraphics[width=2.3in,height=2.3in,angle=0]{qr.pdf}}%
\end{tabular}%
\end{center}
\caption{The figure shows the time evolution of the statefinder pairs $\{r,s\} $ (left) and $\{r,q\}$ (right) for models 1 (dashed), 2 (dot-dashed) and DGP (dotted). In the left panel, the fixed point ($r=1$, $s=0 $) corresponds to $\Lambda$CDM. All the models passes through $\Lambda$CDM. In the right panel, all models diverge from the same point ($r=1$, $q=0.5$) which corresponds to SCDM. The DGP converges to the
point ($r=1,q=-1$) that represents the de-Sitter expansion (dS) whereas models 1 and 2 do not converge to dS due to phantom nature. The dark dots on the curves denote current values $\{r_{0},s_{0}\}$ (left) and $\{r_{0},q_{0}\}$ (right). In all models,
we have taken best-fit values of the model parameters.}
\label{fig:rs}
\end{figure}
\section{Statefinder diagnostic}
\label{sec:state}
The past two decades produced a plethora of theoretical cosmological
models of dark energy with improved quality of observational data. So, there
should be some analysis which can differentiate these models and predict
the deviations from $\Lambda$CDM. Sahni et al. \cite{Sahni1} have pointed out this idea widely known as statefinder diagnostic. The statefinder pairs $\{r,s\}$ and $\{r,q\}$ are the geometrical quantities that are constructed from any space-time metric directly and can successfully differentiate various competing models of dark energy by using the higher order derivatives of scale factor. In the literature, the $\{r,s\}$ and $\{r,q\}$ pairs are defined as \cite{Sahni1}.
\begin{equation}
q=-\frac{\ddot{a}}{aH^{2}}\text{, \ \ } r=\frac{\dddot{a}}{aH^{3}}\text{, \ \ }s=\frac{r-1}{3(q-\frac{1}{2})}\text{.}
\label{d1}
\end{equation}%
The statefinder diagnostic is an useful tool in modern day cosmology and
being used to serve the purpose of distinguishing different dark energy
models \cite{alam1,alam2,alam3}. In this process, different
trajectories in $r-s$ and $r-q$ planes are plotted for various dark energy models
and study their behaviors. In a spatially flat FLRW background, the statefinder pair $\{r,s\}=\{1,0\}$ and $\{1,1\}$ for $\Lambda$CDM and standard cold dark matter (SCDM). In the $r-s$ and $r-q$ planes, the departure of any dark energy model from these fixed points are analyzed. The pairs $\{r,s\}$ and $\{r,q\}$ for models 1 and 2 are calculated as
\begin{equation}
r=\frac{%
\begin{array}{c}
11+101\alpha \tilde{a}+(358\alpha ^{2}+158\beta )\tilde{a}^{2}+(488\alpha
^{3}+1237\alpha \beta )\tilde{a}^{3} \\
+(2440\alpha ^{2}\beta +224\alpha ^{4}+1156\beta ^{2})\tilde{a}%
^{4}+(1424\alpha ^{3}\beta +4087\alpha \beta ^{2})\tilde{a}^{5} \\
+(3346\alpha ^{2}\beta ^{2}+2178\beta ^{3})\tilde{a}^{6}+3375\alpha \beta
^{3}\tilde{a}^{7}+1233\beta ^{4}\tilde{a}^{8}%
\end{array}%
}{2(1+\alpha \tilde{a}+\beta \tilde{a}^{2})^{2}(1+2\alpha \tilde{a}+3\beta
\tilde{a}^{2})^{2}}  \label{d3}
\end{equation}%
\begin{equation}
s=\frac{%
\begin{array}{c}
9+89\alpha \tilde{a}+(332\alpha ^{2}+142\beta )\tilde{a}^{2}+(464\alpha
^{3}+1169\alpha \beta )\tilde{a}^{3} \\
+(216\alpha ^{4}+2348\alpha ^{2}\beta +1112\beta ^{2})\tilde{a}%
^{4}+(1384\alpha ^{3}\beta +3971\alpha \beta ^{2})\tilde{a}^{5} \\
+(3272\alpha ^{2}\beta ^{2}+2130\beta ^{3})\tilde{a}^{6}+3315\alpha \beta
^{3}\tilde{a}^{7}+1215\beta ^{4}\tilde{a}^{8}%
\end{array}%
}{%
\begin{array}{c}
15\alpha \tilde{a}+(63\alpha ^{2}+42\beta )\tilde{a}^{2}+(84\alpha
^{3}+255\alpha \beta )\tilde{a}^{3} \\
+(36\alpha ^{4}+438\alpha ^{2}\beta +222\beta ^{2})\tilde{a}^{4}+(228\alpha
^{3}\beta +693\alpha \beta ^{2})\tilde{a}^{5} \\
+(507\alpha ^{2}\beta ^{2}+342\beta ^{3})\tilde{a}^{6}+477\alpha \beta ^{3}%
\tilde{a}^{7}+162\beta ^{4}\tilde{a}^{8}%
\end{array}%
}  \label{d4}
\end{equation}
and
\begin{equation}
r=\frac{3+(2+18\alpha )\tilde{a}+(4\alpha +41\alpha ^{2})\tilde{a}%
^{2}+(2\alpha ^{2}+33\alpha ^{3})\tilde{a}^{3}+9\alpha ^{4}\tilde{a}^{4}}{2%
\tilde{a}+4\alpha \tilde{a}^{2}+2\alpha ^{2}\tilde{a}^{3}}  \label{d5}
\end{equation}%
\begin{equation}
s=\frac{-3-18\alpha \tilde{a}-41\alpha ^{2}\tilde{a}^{2}-33\alpha ^{3}\tilde{%
a}^{3}-9\alpha ^{4}\tilde{a}^{4}}{-9+(9-30\alpha )\tilde{a}+(18\alpha
-30\alpha ^{2})\tilde{a}^{2}+(9\alpha ^{2}-9\alpha ^{3})\tilde{a}^{3}}
\label{d6}
\end{equation}
The deceleration parameter $q$ is given by Eqs. (\ref{eq:q1}) and (\ref{eq:q2}), respectively. We plot the $r-s$ and $r-q$ diagrams for our models and compare these with the $\Lambda$CDM.

Fig. \ref{fig:rs} shows the time evolution of the statefinder pairs for different DE models. The left panel exhibits the evolution of $\{r,s\}$ while the right one for $\{r,q\}$. In both panels, the models 1 (dashed) and 2 (dot-dashed) are compared with DGP (dotted) and $\Lambda$CDM. In left panel, the fixed point ($r=1$, $s=0$)
corresponds to $\Lambda$CDM, and all models passes through this fixed point. Moreover, one can see that the trajectory of DGP terminates there while the corresponding
trajectories of models 1 and 2 evolve further showing that the phantom behavior in future. In right panel, the fixed point ($r=1$, $q=0.5$) represents SCDM. All the underlying models evolve from this point. The DGP converges to the second point ($r=1,q=-1$) that corresponds to de-Sitter expansion (dS) whereas models 1 and 2 do not
converge to dS due to their phantom behavior. The dark dots
on the curves show present values $\{r_{0},s_{0}\}$ (left) and $%
\{r_{0},q_{0}\}$ (right) for the models under consideration. We chose the best-fit values of the model parameters for these plots.
\section{$Om$ diagnostic}
\label{sec:om}
We shall now diagnose our models with $Om$ analysis which
is also a geometrical diagnostic that explicitly depends on redshift and the Hubble
parameter, and is defined as \cite{Sahni2,Z}:
\begin{equation}
Om\left( z\right) =\frac{\left( \frac{H(z)}{H_{0}}\right)^{2}-1}{\left( 1+z\right)
^{3}-1}\text{}
\label{eq:om}
\end{equation}
\begin{figure}[tbph]
\begin{center}
{\includegraphics[width=2.5in,height=2.5in,angle=0]{om.pdf}}
\end{center}
\caption{This figure shows the evolution of $Om(z)$ versus redshift $z$
for different DE models such as $w=-1$ ($\Lambda$CDM), $-1.4$ (phantom), $-0.6$ (quintessence) and models 1 (dashed) \& 2 (dot-dashed). The horizontal line represents $\Lambda$CDM, and has zero curvature. The DE models with $w>-1$
(quintessence) have negative curvature whereas models with $w<-1$
(phantom) have positive curvature. Models 1 and 2 show the positive curvatures though they have $w_{eff}>-1$ at the present epoch and in future $w_{eff}<-1$
(phantom phase). The vertical solid line represents the present era. The
best-fit values are chosen to plot this figure.}
\label{fig:om}
\end{figure}
The $Om$ diagnostic also differentiates various DE models from $\Lambda$CDM \cite{Amnras}. This is a simpler diagnostic when applied to observations as it depends only on the first derivative of scale factor. The Hubble parameter for a constant EOS is defined as
\begin{equation}
H^2(z) = H_0^2 \left[ \Omega_{0m}(1+z)^3 + (1-\Omega_{0m})(1+z)^{3(1+w)}\right] ,
\label{eq:Hconst}
\end{equation}
The expression for $Om(z)$ corresponding to Eq. (\ref{eq:Hconst}) is written as
\begin{equation}
Om(z) = \Omega_{0m} + (1-\Omega_{0m})\frac{(1+z)^{3(1+w)}-
1}{(1+z)^3-1}
\label{eq:omconst}
\end{equation}
From Eq. (\ref{eq:omconst}), one can notice that, $Om(z) = \Omega_{0m}$ for $w=-1$ ($\Lambda$CDM) whereas $Om(z) > \Omega_{0m}$ for $w> -1$ (quintessence) and $Om(z) < \Omega_{0m}$ for $w < -1$ (phantom). The corresponding evolutions of $Om(z)$ for $w=-0.6$ (quintessence), $-1$ ($\Lambda$CDM) and $-1.4$ (phantom) are shown in Fig. \ref{fig:om}. From Fig. \ref{fig:om}, one can clearly see that the $Om(z)$ has negative, zero and positive curvatures for quintessence, $\Lambda$CDM and phantom, respectively. In contrast, we show the evolution of $Om(z)$ for models 1 and 2 in Fig. \ref{fig:om}. Both models exhibit positive curvatures though thay have $w_{eff}>-1$ (quintessence) at present epoch and in future $w_{eff}<-1$ (phantom phase), see Figs. \ref{fig:hw} and \ref{fig:om}. This is a vital result as in the literature, quintessence does not have positive curvature which is a generic feature \cite{Amnras}.
\section{Conclusion}
\label{sec:conc}
In this work, we have considered the scenario in which cosmic
acceleration might arise due to coupling between known matter components
present in the Universe \cite{KHOURY2016}. To this effect, we further extend
the work of Ref. \cite{ARSSA}. The two models obtained by parameterizing the
coupling function (or correspondingly the Einstein frame scale factor in
terms of physical scale factor) are analyzed by using an integrated observational data.
We used joint data of $H(z)$, Type Ia supernova and BAO, and constrained the model parameters of model 2 (see Fig. \ref{fig:cont}). In this case, the best-fit values of
$\alpha$ and $H_0$ are found to be $\alpha =-0.3147$ and $H_{0}=66.84~Km/S/Mpc$ (model 1 consists of two parameters that were constrained in Ref. \cite{ARSSA} and the best-fit values were obtained as $\alpha =-0.102681$ \& $\beta =-0.078347$). We used best-fit values of the model parameters to carried out the analysis. The time evolutions of Hubble parameter and effective EOS  are shown in Fig. \ref{fig:hw}. From this figure we conclude that the model 1 shows phantom behavior with a pressure singularity in the near future while the model 2 is a generalized case of model 1 that pushes the future sudden singularity to the infinite future. In Fig. \ref{fig:mu}, we have shown error bars of observational data with the models under consideration. One can clearly observe that both the models are compatible with the observations.

In addition, the statefinder diagnostic has been performed for the underlying models. We obtained the expressions for the statefinder pairs and their behaviors have been displayed in $r-s$ and $r-q$ planes as shown in Fig. \ref{fig:rs}. For comparison, we have also shown DGP model in same figure. In the $r-s$ plane, both the models pass
through the fixed point ($r=1, s=0$) and move away from the $\Lambda$CDM while the
DGP terminates at the fixed point. This is due the phantom behavior which is a distinguished characteristic of the underlying scenario. In the $q-r$ plane, it is clearly seen that all the models originated from a
fixed point ($r=1,q=0.5$) that corresponds to SCDM.
The DGP model converges to the second point ($r=1,q=-1$) that represents
the dS while the models under consideration do not
converge to the dS fixed point due to their phantom nature (see right
panel of Fig. \ref{fig:rs}).

The evolution of $Om(z)$ versus redshift $z$ for different DE models are shown in Fig. \ref{fig:om}. We observed that the $\Lambda$CDM, quintessence and phantom have zero, negative and positive curvatures. Models 1 and 2 lie in quintessence regime in the past and remain so till the present epoch, and evolve to phantom in future. Both the models have positive curvature however they lie in the quintessence regime. This is an important result as in the literature, its not possible to have positive curvature for quintessence \cite{Amnras}.

In our opinion, the scenario proposed in \cite{KHOURY2016} and investigated
here is of great interest in view of GW170817. The modification caused by
the interaction between known components of matter does not involve any
extra degree of freedom and falls into the safe category in the light of
recent observations on gravitational waves.

\section*{Acknowledgment}
We are highly indebted to M. Sami for suggesting this problem and constant supervision as well as for providing all necessities to complete this work. SKJP wishes to
thank National Board of Higher Mathematics (NBHM), Department of Atomic
Energy (DAE), Govt. of India for financial support through the post-doctoral
research fellowship.


\section{Introduction}
Let $Y$ be a positive random variable representing the income (or expenditure) distribution in a given population. Suppose that $Y$ is defined on a probability espace $(\Omega,\mathcal{A},\mathbb{P})$, with underlying continuous distribution function $G(y)=\mathbb{P}(Y\leq y),\;\forall y\geq 0.$  Given a   poverty line $z>0$, we say that an individual is poor if his income  is less than $z$. Most of  the poverty measures employed in practice may be represented, in their continuous form, by the following functional :
 \begin{equation} \label{s1}
J(w,f)=\int^{z}_{0}w[G(y),G(z)]f(y,z)dG(y),
\end{equation}
where $w(u,v)$ is a function of $(u,v)\in [0,1]^2$, which is decreasing with respect to its first argument $u$. It is interpreted as a weighting function associated with the kind of poverty measure being considered ; while $f(y,z)$ is called the \textit{poverty deprivation function}, which evaluates the individual contribution of each poor to the global poverty within the population. Note that the function $f(y,z)$ is also decreasing with respect to its first argument $y$.\\

 \indent Formula (\ref{s1}) is quite general and represents a wide class of poverty indices including the most currently used in practice. For instance, the additively decomposable family of poverty measures can be obtained from this formula \eqref{s1} by letting the weighting function $w(u,v)\equiv 1$, for all $(u,v)\in \left[0,1\right]^2$. As well, the non-additively decomposable poverty measures  such as the Sen-like poverty indices may  also be put in the form (\ref{s1}), with specific weighting  and deprivation functions $w(\cdot,\cdot)$ and $f(\cdot,\cdot)$. In the Table \ref{tab0} below, we give some examples of poverty indices with their own weighting and deprivation functions.
 \begin{table}[htbp]
\begin{center}
$$\begin{array}{ccc}
\hline
 \text{Poverty index} &  w(u,v) &  f(y,z)\\
 \hline
\text{FGT}(\alpha)  &         1  &    (1-y/z)^{\alpha}  \\
\text{Sen}  &      2(1- u/v)     &    (1-y/z)   \\
\text{Shorrocks}  & 2(1-u)     & (1-y/z)       \\
\text{Kakwani} (k)     &  (k+1)(1-u/v)^k    &  (1-y/z)^{k}     \\
\text{Watts}    &      1    & \log (z/y)     \\
\hline
\end{array}
$$
\end{center}
\caption{Examples of poverty indices with their weighting and deprivation functions.}
\label{tab0}
\end{table}

Our main goal in this paper is to propose an estimator for the theoretical functional $J(\cdot,\cdot)$, and study its asymptotic behavior by considering a relative poverty line. Indeed, empirical studies with fixed poverty lines are plentiful in the literature ( see, e.g.,\cite{kak1}, \cite{bfz},\cite{rg}, \cite{lss}). But most of them do not take account of the sampling error of the poverty line, which may increase or reduce the precision of the poverty estimates (see, e.g. \cite{p}). Investigating relative poverty in US, Zheng \cite{zh2}  proposed an approach which consider the poverty line as a percentage of the mean or a quantile of the distribution function. He dealt with additively separable (or decomposable) poverty measures, and found that the sampling error associated with poverty lines increases the standard error of the poverty estimates.\\

\indent In this paper, we propose an inference procedure which somewhat extends Zheng's \cite{zh2}  approach to  non-additively decomposable poverty measures including  Sen, Shorrocks  and Kakwani indices as well as many other poverty measures available in the literature. Note that the Kakwani's family  is the most interesting family of poverty indices, because it contains the FGT (Foster-Greer-Thorbecke) and Sen measures, and satisfies most of the normative  axioms desirable on on poverty index. Our methodology is inspired by the modern empirical process theory developed in van der Vaart and Wellner \cite{vdv}, which permits us to obtain the uniform weak convergence of a wide class of empirical  poverty estimators.\\

\indent The rest of the paper is organized as follows. In Section 2, we construct an estimator for the functional $J(w,f)$, and establish its asymptotic properties. In Section 3, we derive  testing procedures which allows  to make comparisons based either on one single poverty measure or simultaneously several poverty indices. Section 4 shows, in a simulation study,  the asymptotic normality of the proposed estimator. In Section 5, we give an application to real data sets to estimate the effect of the relative poverty line on the variance of the  some particular poverty estimates. Finally, we prove the main theorem in Appendix.

\section{Construction of the estimator and asymptotic results}
 Let  $Y_1,\cdots, Y_n$ be an independent and identically random sample of the income variable $Y$, with density probability function $g$. Whenever the poverty line $z>0$ is fixed, a direct estimator of \eqref{s1} can be defined as follows :
\begin{equation} \label{s2}
J_n(w ,f) =\frac{1}{n}\sum_{j=1}^{n}
w[G_n(Y_{j,n}),G_n(z)]f(Y_{j,n},z) \mathbb{I}(Y_{j,n}\leq
z),\qquad
\end{equation}
for any specific functions $w$ and  $f$. Here $\mathbb{I}(\cdot)$ designs the indicator function,  $Y_{1,n}\leq\cdots\leq Y_{n,n}$ are  order statistics associated with the sample  $Y_1,\cdots, Y_n$  and $G_n(y)=\frac{1}{n}\sum_{j=1}^{n}\mathbb{I}(Y_{j}\leq y)$ is the  corresponding empirical distribution function.\\
 In contrast,  when we are concerned with relative poverty analysis, the poverty line becomes a functional of the distribution  $G$, say $z=z(G)$, and then must be estimated. Let $\hat{z}$ be a consistent estimator of $z$ such that the following representation $(R)$  (see, e.g. Thuysbaert and Zitikis \cite{tz}) :
$$ (R) \qquad \hat{z}= z + \frac{1}{n}\sum_{j=1}^{n}\zeta(Y_j)+ o_{\mathbb{P}}(n^{-1/2}), $$
where $\zeta(\cdot)$ is a real-valued function such that ${\rm Var}[\zeta(Y)]=\int_0^{\infty} \zeta^2(y)dG(y)<\infty$. \\
The function $\zeta(\cdot)$ may be equal to the constant 0, in which case $\hat{z}=z$ corresponds to an absolute poverty line. If the poverty line is set to a fraction $k$ of the mean of the distribution, i.e. ${z}=k\mu_G$, with $\mu_G=\int ydG(y)$, then $\zeta(y)=ky$. Whenever a fraction $k$ of a $p$-quantile is chosen, i.e. ${z}=kG^{-1}(p)$,  the Bahadur representation yields $\zeta(y)=\frac{k}{g(G^{-1}(p))}\mathbb{I}(y\leq G^{-1}(p) ).$\\
 Substituting $z$ for $\hat{z}$ in \eqref{s2}, we obtain a novel estimator
\begin{equation} \label{s3}
\hat{J}_n(w ,f) =\frac{1}{n}\sum_{j=1}^{n}
w[G_n(Y_{j,n}),G_n(\hat{z})]f(Y_{j,n},\hat{z}) \mathbb{I}(Y_{j,n}\leq
\hat{z}).
\end{equation}
Now, we have to prove that $\hat{J}_n(w,f)$ converges almost surely to $J(w,f)$ for any specific functions $w$ and $f$. Because of certain normative properties desirable on a poverty index, the  functions $w$ and $f$ should satisfy some regularity conditions ; that is they belong to the following classes of functions $\mathcal{W}$ and $\mathcal{F}$, respectively :
\begin{eqnarray*}
\mathcal{W}&=&\{w:[0,1] \times [0,1] \rightarrow \mathbb{R}_{+}, \; w
\;\;\text{  continuous,}\;\;\text{and}\;\;u\mapsto w(u,\cdot)
\;\;\text{is non-increasing}\}\\
\mathcal{F}&=&\{f:\mathbb{R}_{+}\times\mathbb{R}_{+}\rightarrow \mathbb{R}_{+},\;f
\;\;\text{ continuous,}\;\;\text{and}\;\;y\mapsto f(y,\cdot)
\;\;\text{is non-increasing}\}.
\end{eqnarray*}
We also need the following conditions (A.1-2) to establish our asymptotic results :
\begin{itemize}
\item[(A.1)] The functions $(u,v)\mapsto w(u,v)$ and $(y,z)\mapsto f(y,z)$ are differentiable  with continuous first-order partial derivatives ;
\item[(A.2)] $\mathcal{W}$ and $\mathcal{F}$ are both pointwise measurable classes of functions. That is, they contain each one a countable subclass $\mathcal{G}$ such that for all $\phi\in\mathcal{G}$, there exists a sequence $\{\phi_m\}_{m\geq 1} \subset \mathcal{G}$, \noindent with $\phi_m(y)\rightarrow \phi(y)$ for every $y$.\\
\end{itemize}
\indent It is established in \cite{ls} that the estimator $J_n(w,f)$, with $z$ fixed, converges almost surely to $J(w,f)$ for any given functions $w\in \mathcal{W}$ and $f\in\mathcal{F}$.
 In the following proposition, we gives the almost sure consistency of  the estimator $\hat{J}_n(w,f)$ for $J(w,f)$.
\begin{proposition}\label{l1}
 For any couple of functions $(w,f)\in\mathcal{W}\times\mathcal{F}$, one has with probability 1,
 $$ \hat{J}_n(w,f)\longrightarrow J(w,f),\qquad n\longrightarrow\infty.$$
 \end{proposition}
 \textbf{Proof}. It suffices to show that for all $(w,f)\in\mathcal{F\times W}$, $\hat{J}_n(w,f)$ is asymptotically equivalent to ${J}_n(w,f)$ plus an additional quantity of the form $a(\hat{z}-z)$, where the factor $a$ will be specified later on. Let's decompose $\hat{J}_n(w,f)$ as follows :
 \begin{eqnarray*}
 \hat{J}_n(w,f) &= &\frac{1}{n}\sum_{j=1}^{n}w[G_n(Y_{j,n}),G_n({z})]f(Y_{j,n},{z}) \mathbb{I}(Y_{j,n}\leq {z})\\
 &= &\frac{1}{n}\sum_{j=1}^{n} \left\{ w[G_n(Y_{j,n}),G_n(\hat{z})]f(Y_{j,n},\hat{z})- w[G_n(Y_{j,n}),G_n({z})]f(Y_{j,n},{z})\right\} \mathbb{I}(Y_{j,n}\leq {z})\\
 &=& \frac{1}{n}\sum_{j=1}^{n} w[G_n(Y_{j,n}),G_n(\hat{z})]f(Y_{j,n},\hat{z})\left\{\mathbb{I}(Y_{j,n}\leq \hat{z})-\mathbb{I}(Y_{j,n}\leq {z}) \right\} \\
 &=:& I+II+III
 \end{eqnarray*}
One can readily observe that the first term $I$ is exactly $J_n(w,f)$. \\
\indent By applying the mean value theorem, the second term $II$ becomes
\begin{eqnarray*}
II&=&\frac{1}{n}\sum_{j=1}^{n}\left\{ \frac{\partial}{\partial v} w[G_n(Y_{j,n}),G_n(z_0)]f(Y_{j,n},\hat{z})[G_n(\hat{z})-G_n(z)] \right\}\mathbb{I}(Y_{j,n}\leq {z})\\
&+ &  \sum_{j=1}^{n}\left\{\frac{\partial}{\partial z}f(Y_{j,n},z_0)w[G_n(Y_{j,n}),G_n(z)][\hat{z}-z]\right\}\mathbb{I}(Y_{j,n}\leq {z}),
\end{eqnarray*}
where $z_0$ is between $z$ and $\hat{z}$. Recall that  $G_n(y)\rightarrow G(y)$, almost surely for all $y\geq 0$, then we  can write for $n$ large enough $G_n(y)=G(y)+ o(n^{-1/2}), \forall y\geq 0$. Thus, applying again the mean value theorem, we obtain for all large $n$,
\begin{eqnarray*}\label{p1}
 \frac{1}{n}\sum_{j=1}^{n}[\mathbb{I}(Y_{j,n}\leq \hat{z})-\mathbb{I}(Y_{j,n}\leq {z})]= G_n(\hat{z})-G_n(z)&=& G(\hat{z})-G(z)+o(n^{-1/2})\\
& =& g(z_1)[\hat{z}-z]+ o(n^{-1/2}),
\end{eqnarray*}
where $z_1$ is between $z$ and $\hat{z}$. Hence, the second term $II$ can be rewritten into
\begin{eqnarray*}
II&=&[\hat{z}-z]\frac{1}{n}\sum_{j=1}^{n} \frac{\partial}{\partial v} w[G_n(Y_{j,n}),G_n(z_0)]f(Y_{j,n},\hat{z})g(z_1)\mathbb{I}(Y_{j,n}\leq {z})\\
&+ &  [\hat{z}-z]\frac{1}{n}\sum_{j=1}^{n}\frac{\partial}{\partial z}f(Y_{j,n},z_0)w[G_n(Y_{j,n}),G_n(z)]\mathbb{I}(Y_{j,n}\leq {z})+ o(n^{-1/2}).
\end{eqnarray*}

 For the last term $III$, we also make use of  Taylor expansion.  For $Y_{j,n}$ in the vicinity of $z$, we have
\begin{eqnarray*}
 w[G_n(Y_{j,n}),G_n(\hat{z})]f(Y_{j,n},\hat{z})&=& w[G_n(z),G_n(\hat{z})]f(z,\hat{z})\\
 &+ &  \frac{\partial}{\partial u} w[G_n(z),G_n(\hat{z})]f(z,\hat{z})[G_n(Y_{j,n})-G_n(z)] \\
 &+& \frac{\partial}{\partial y}f(z,\hat{z})w[G_n(z),G_n(\hat{z})][Y_{j,n}-z]\\
 &+& o(|G_n(Y_{j,n})-G_n(z)|+|Y_{j,n}-z|).\\
\end{eqnarray*}
Thus $III$ can be transformed into
\begin{eqnarray*}
III&=&\frac{1}{n}\sum_{j=1}^{n}w[G_n(z),G_n(\hat{z})]f(z,\hat{z})[\mathbb{I}(Y_{j,n}\leq \hat{z})-\mathbb{I}(Y_{j,n}\leq {z})]\\
& +&\frac{1}{n}\sum_{j=1}^{n}\frac{\partial}{\partial u} w[G_n(z),G_n(\hat{z})]f(z,\hat{z})[G_n(\hat{z})-G_n(z)][\mathbb{I}(Y_{j,n}\leq \hat{z})-\mathbb{I}(Y_{j,n}\leq {z})] \\
& +& \frac{1}{n}\sum_{j=1}^{n}\frac{\partial}{\partial y}f(z,\hat{z})w[G_n(z),G_n(\hat{z})][Y_{j,n}-z][\mathbb{I}(Y_{j,n}\leq \hat{z})-\mathbb{I}(Y_{j,n}\leq {z})]\\
&+ & \frac{1}{n}\sum_{j=1}^{n}o(|G_n(Y_{j,n})-G_n(z)|+|Y_{j,n}-z|)[\mathbb{I}(Y_{j,n}\leq \hat{z})-\mathbb{I}(Y_{j,n}\leq {z})].
\end{eqnarray*}
Now, we are going show that the last three terms  in the right-hand side of the previous equality are asymptotically negligible.
For the second term we can write, in view of \eqref{p1}, that
 \begin{eqnarray*}
& \left|\frac{1}{n}\sum_{j=1}^{n}\frac{\partial}{\partial u} w[G_n(z),G_n(\hat{z})]f(Y_{j,n},\hat{z})[G_n(\hat{z})-G_n(z)][\mathbb{I}(Y_{j,n}\leq \hat{z})-\mathbb{I}(Y_{j,n}\leq {z})] \right|& \\
& \leq w[G_n(z),G_n(\hat{z})]f(z,\hat{z})|G_n(\hat{z})-G_n(z)|^2 &\\
& \leq w[G_n(z),G_n(\hat{z})]f(z,\hat{z})[g^2(z_1)|\hat{z}-z|^2+ o(n^{-1/2})].
 \end{eqnarray*}
 Thus, the second term tends to 0, because $\hat{z}\rightarrow z, \; n\rightarrow\infty$.
Analogously,  by using \eqref{p1} and assuming that the $Y_{j,n}$ are lying between $z$ and $\hat{z}$, we can major the third term by a quantity converging to 0.
\begin{eqnarray*}
& \left|\frac{1}{n}\sum_{j=1}^{n}\frac{\partial}{\partial y}f(z,\hat{z})w[G_n(Y_{j,n}),G_n(\hat{z})][Y_{j,n}-z][\mathbb{I}(Y_{j,n}\leq \hat{z})-\mathbb{I}(Y_{j,n}\leq {z})]\right| &\\
& \leq \frac{\partial}{\partial y}f(z,\hat{z})w[G_n(z),G_n(\hat{z})]|\hat{z}-z|[G_n(\hat{z})-G_n(z)] &\\
& \leq \frac{\partial}{\partial y}f(z,\hat{z})w[G_n(z),G_n(\hat{z})]g(z_1)|\hat{z}-z|^2 + o(n^{-1/2}).&\\
\end{eqnarray*}
For the last and fourth term, one has
\begin{eqnarray*}
&\left|\frac{1}{n}\sum_{j=1}^{n}o(|G_n(Y_{j,n})-G_n(z)|+|Y_{j,n}-z|)[\mathbb{I}(Y_{j,n}\leq \hat{z})-\mathbb{I}(Y_{j,n}\leq {z})]\right|& \\
& \leq g(z_1)o(|\hat{z}-z|^2)+ o(n^{-1/2}).&
\end{eqnarray*}
Thus the fourth term also tends to 0, as $n\rightarrow\infty$. It follows from this that  $III$ is asymptotically equivalent to
\begin{equation}
III=w[G_n(z),G_n(\hat{z})]f(z,\hat{z})g(z_1)[\hat{z}-z)]+ o(n^{-1/2}).
\end{equation}
Finally, we obtain, for $n$ large enough, the following decomposition for  $\hat{J}_n(w,f)$ :
\begin{eqnarray*}
 \hat{J}_n(w,f) &= &\frac{1}{n}\sum_{j=1}^{n}w[G_n(Y_{j,n}),G_n({z})]f(Y_{j,n},{z}) \mathbb{I}(Y_{j,n}\leq {z})\\
&+ &[\hat{z}-z]\frac{1}{n}\sum_{j=1}^{n}\frac{\partial}{\partial v} w[G_n(Y_{j,n}),G_n(z_0)]f(Y_{j,n},\hat{z})g(z_0)\mathbb{I}(Y_{j,n}\leq {z})\\
&+ &[\hat{z}-z]\frac{1}{n}\sum_{j=1}^{n}\frac{\partial}{\partial z}f(Y_{j,n},z_0)w[G_n(Y_{j,n}),G_n(z)]\mathbb{I}(Y_{j,n}\leq {z})\\
 &+& w[G_n(z),G_n(\hat{z})]f(z,\hat{z})g(z_1)[\hat{z}-z)]+ o(n^{-1/2}).
\end{eqnarray*}
Combining the continuity of the partial derivatives of $w$ and $f$ on the compact interval $[0,z]$ and the fact that $z_0$ converges to $z$, we can approximate the second and the third summation terms in the second member of the above equality  respectively by the integrals $$\int_{0}^{z}\frac{\partial}{\partial v} w[G(y),G(z)]f(y,z)g(z)dG(y)$$ and $$\int_{0}^{z}\frac{\partial}{\partial z}f(y,z) w[G(y),G(z)]dG(y).$$
Besides, since $\hat{z}$ converges almost surely to $z$ and $f$ continuous,  $f(z,\hat{z})$ converges almost surely to $f(z,z)$, which is equal to 0, in virtue of the normative focalisation axiom on poverty indices.  Thus, the fourth  term converges to 0, and $\hat{J}_n(w,f)$ becomes
\begin{equation}\label{r1}
\hat{J}_n(w,f)={J}_n(w,f)+a(\hat{z}-z)+ o(n^{-1/2}),
\end{equation}
where
$$a=\int_{0}^{z}\left(\frac{\partial}{\partial v} w[G(y),G(z)]f(y,z)g(z)+ \frac{\partial}{\partial z}f(y,z) w[G(y),G(z)]\right) dG(y).$$
From this, we can conclude that $\hat{J}_n(w,f)$ converges almost surely to ${J}(w,f)$, as $n\rightarrow\infty$.$\square$\\

 Now, we state our main result  which is the weak convergence of the normalized and centered process $\{\sqrt{n}[\hat{J}_{n}(w,f)-J(w,f)]: w\in \mathcal{W}, f\in\mathcal{F}\}$ in $l^{\infty }(\mathcal{W}\times
\mathcal{F})$, the set of all real-valued and bounded functions defined on $\mathcal{W}\times \mathcal{F}.$
\begin{theorem}\label{t1}
Let $G(y)$ be a continuous distribution function with probability density  $g(y)$. If assumptions (A.1-2)  hold, then the process
$\{\sqrt{n}[\hat{J}_{n}(w,f)-J(w,f)]: w\in \mathcal{W}, f\in\mathcal{F}\}$ converges weakly in $l^{\infty }(\mathcal{W}\times\mathcal{F})$ to a zero-mean Gaussian process with covariance function defined, for any $(w,f),(\widetilde w,\widetilde f)\in\mathcal{W\times F}$, as
\begin{eqnarray*}\label{st}
\Gamma \left[(w,f); (\widetilde w,\widetilde f)\right] &=& \Sigma \left[(w,f); (\widetilde w,\widetilde f)\right] +a\left( \int_{0}^{\infty} h(y)\zeta(y)dG(y)-\mathbb{E}[\zeta(Y)]J(w,f)\right)  \\ \nonumber
& +  & a\left( \int_{0}^{\infty} \widetilde{h}(y)\zeta(y)dG(y)-\mathbb{E}[\zeta(Y)]J(w,f)\right)+ a^2{\rm Var}[\zeta(Y)],
\end{eqnarray*}
with
$$ h(y)=w[G(y),G(z)]f(y,z)\mathbb{I}(y\leq z) \quad; \quad\widetilde{h}(y)=\widetilde w[G(y),G(z)]\widetilde f(y,z)\mathbb{I}(y\leq z)$$
and
\begin{eqnarray*}
\Sigma \left[(w,f); (\widetilde w,\widetilde f)\right] &=&
\int^{z}_{0}w[G(y),G(z)]f(y,z)\widetilde w[G(y),G(z)]\widetilde f(y,z)dG(y)\\
 &- &\int^{z}_{0}w[G(y),G(z)]f(y,z)dG(y)\int^{z}_{0}\widetilde w[G(y),G(z)] \widetilde f(y,z)dG(y)\\
 &+&\int_0^{z}\int_0^{z} a_1(x,y)[G(x)\wedge G(y) - G(x)G(y)]dG(x)dG(y), \nonumber  \\
&+& [1-G(z)]\int_0^{z}\int_0^{z} a_2(x,y)G(y)dG(x)dG(y) \nonumber  \\
&+&[1- G(z)] \int_0^{z}\int_0^{z} a_3(x,y)G(x)dG(x)dG(y)  \nonumber  \\
&+&  G(z)[1-G(z)] \int_0^{z}\int_0^{z} a_4(x,y)dG(x)dG(y), \nonumber
\end{eqnarray*}
where
\begin{eqnarray*}
a_1(x,y) & =& \frac{\partial w}{\partial u}[G(x),G(z)]f(x,z)\frac{\partial \widetilde w}{\partial u}[G(y),G(z)] \widetilde f(y,z),\\
a_2(x,y) & =& \frac{\partial w}{\partial
u}[G(x),G(z)]f(x,z)\frac{\partial \widetilde w}{\partial v}[G(y),G(z)] \widetilde f(y,z),\\
a_3(x,y)& =&\frac{\partial w}{\partial v}[G(x),G(z)]f(x,z)\frac{\partial \widetilde w}{\partial u}[G(y),G(z)] \widetilde f(y,z),\\
a_4(x,y)& =& \frac{\partial w}{\partial v}[G(x),G(z)]f(x,z)\frac{\partial \widetilde w}{\partial v}[G(y),G(z)] \widetilde f(y,z).\\
\end{eqnarray*}
\end{theorem}
\textbf{Remark 1.} For any given functions $w$ and $f$, Theorem \ref{t1} gives the variance of the estimator $\hat{J}_n(w,f)$ which is equal to
\begin{eqnarray}\label{stt}
\Gamma \left[(w,f); (w, f)\right] &=& \Sigma \left[(w,f); ( w,f)\right] +2a\left(\int_{0}^{\infty} h(y)\zeta(y)dG(y) -\mathbb{E}[\zeta(Y)]J(w,f)\right)  \nonumber\\
 &&+ a^2{\rm Var}[\zeta(Y)].
\end{eqnarray}
This means that the variance of the poverty estimate $\hat{J}_n(w,f)$ is increased by a term $$\Delta=2a\left(\int_{0}^{\infty} h(y)\zeta(y)dG(y)-\mathbb{E}[\zeta(Y)]J(w,f)\right)+ a^2{\rm Var}[\zeta(Y)]$$ whenever the poverty line $z$ is estimated from the sample.\\

\textbf{Remark 2.} If the function $w\equiv 1$, which corresponds to the additively decomposable measures, one can observe that all the integrals  with $a_i$-term $i=1,2,3,4$  vanish.Then the remainding terms in \eqref{stt} are exactly the expressions in Equations (21) and (22) in \cite{zh2} for suitable functions $\zeta$. The quantity corresponding to sum of all  integrals  with $a_i$-term is due to the weight function $w$,
when it is considered, as for example, in the case of Sen poverty index.\\

\textbf{Proof.} It relies on the following decomposition which is obtained from \eqref{r1}
\begin{equation}\label{r2}
\sqrt{n}[\hat{J}_n(w,f)-J(w,f)]=\sqrt{n}[{J}_n(w,f)-J(w,f)]+ \frac{a}{\sqrt{n}}\sum_{j=1}^{n}\zeta(Y_j) +o_{\mathbb{P}}(1).
\end{equation}
Observe that the  term $\frac{a}{\sqrt{n}}\sum_{j=1}^{n}\zeta(Y_j)$ in the right-hand side of \eqref{r2} is a sum of independent and identically random variables with mean $\mathbb{E}[\zeta(Y)]$ and finite variance ${\rm Var}[\zeta(Y)]=\int \zeta^2(y)dG(y)$. Then, by applying the central limit theorem, it converges in law to a Gaussian random variable, with  variance $a^2 {\rm Var}[\zeta(Y)]$.\\
\indent Next, by using the modern theory of empirical processes indexed by functions, we prove in Appendix (see, also \cite{ls}) that the centered and normalized process $\{\sqrt{n}[J_{n}(w,f)-J(w,f)]: w\in \mathcal{W}, f\in\mathcal{F}\}$ converges weakly in $l^{\infty }(\mathcal{W}\times\mathcal{F})$ to a tight Gaussian process with zero-mean and covariance function given by the kernel $\Sigma(\cdot,\cdot)$. This entails that the process $\{\sqrt{n}[J_{n}(w,f)-J(w,f)]: w\in \mathcal{W}, f\in\mathcal{F}\}$ is asymptotically tight. Since the second term $\frac{a}{\sqrt{n}}\sum_{j=1}^{n}\zeta(Y_j)$ does not depend on the indexing parameter $(w,f)$, we can infer that the sum process $\{\sqrt{n}[{J}_{n}(w,f)-J(w,f)]+\frac{a}{\sqrt{n}}\sum_{j=1}^{n}\zeta(Y_j): w\in \mathcal{W}, f\in\mathcal{F}\}$ is asymptotically tight.
 Moreover,  the $o_{\mathbb{P}}(1)$-term converges uniformly to 0 in $(w,f)$, as it does  not depend on $(w,f)$.  Thus, since the finite margins of this process are Gaussian (by applying the multivariate central limit theorem), we can conclude that $\{\sqrt{n}[\hat{J}_{n}(w,f)-J(w,f)]: w\in \mathcal{W}, f\in\mathcal{F}\}$ in distribution to a limit Gaussian process. By independence of the $Y_{j}$'s, the cross covariance of the two terms in the right hand side of \eqref{r2} is given for all $(w,f)$ by
\begin{eqnarray*}
\sigma_{w,f,\zeta} &=&\frac{a}{n} \sum_{i=1}^{n}\sum_{j=1}^{n}{\rm cov}(w[G_n(Y_{i,n}),G_n({z})]f(Y_{i,n},{z}), \zeta(Y_j))\\
&=& a.{\rm cov}(w[G_n(Y_{i,n}),G_n({z})]f(Y_{i,n},{z}), \zeta(Y_j))\\
&=& a\left(\int_{0}^{z}w[G_n(y),G_n(z)]f(y,{z}) \zeta(y)dG(y) \right.\\
 &  & \left. -\int_{0}^{z}w[G_n(y),G_n(z)]f(y,{z})dG(y)\int_{0}^{\infty}\zeta(y)dG(y)\right)\\
&\longrightarrow & a\left(\int_{0}^{z}w[G(y),G(z)]f(y,{z}) \zeta(y)dG(y)-\mathbb{E}[\zeta(Y)]J(w,f)\right),\; n \rightarrow\infty.
\end{eqnarray*}

 \section{ Testing procedures }
 Inference procedures for  testing poverty  usually allow one to say that there is less or more poverty in a given population than in another, but  do not permit to answer questions of type : \textit{How much poverty has been decreased or increased ?} Therefore, it is not possible to use  these procedures in order to quantify the poverty variation (or change) between two populations . In this section, we propose a testing procedure which allows us to evaluate the poverty change between two populations, by checking  for whether there exists a proportionality relation between their associated poverty indices. That is, we aim to test the following hypotheses :
   $$H_0: J_F=\alpha J_G,\qquad \text{versus}\qquad  H_1: J_F\neq \alpha J_G,$$
where $\alpha$ is a positive real number, and  $J_F$,  $J_G$ are aggregated scalar poverty  indices computed from distributions $F$ and $G$. Note that $J_F$ and $J_G$ must be defined with the same specific functions $w$ and $f$ indicating the type of poverty measure being considered
The acception of the null hypothesis $H_0$ leads to the estimation of the relative poverty variation between the two distributions $F$ and $G$. For example, if $\alpha=1/2$, we can say that poverty has decreased by an half, if the reference distribution is $G$.  These kind of conclusions are important for policy makers, as they enable to show the effect of poverty reducing strategies. Indeed, our approach  may be applied to check for the poverty Millennium Development Goals (MDG) which consisted of halving poverty in the world by 2015.\\

\indent Consider now two independent populations with cumulative distribution functions $F$ and $G$, and relative poverty lines $z_F$ and $z_G$, respectively. Assume that two independent and identically distributed samples of sizes $n_F$ and $n_G$ are respectively drawn from thereof. Denote by $\widehat{J}_F=\widehat{J}_F(w,f)$ and $\widehat{J}_G=\widehat{J}_G(w,f)$ the estimators of $J_F$ and $J_G$, respectively. By
{Theorem \ref{t1}}, $\widehat{J}_F$ and $\widehat{J}_G$ are
asymptotically normally distributed, with  variances $\sigma^2_F$ and $\sigma^2_G$ that can be readily determined from \eqref{st} by computing $\Gamma[(w,f);(w,f)]$ with the right distribution $F$ or $G$.\\

Let $\Delta \widehat{J}=\widehat{J}_F- \alpha\widehat{J}_G$. Then, under $H_0$, $\Delta \widehat{J}$ is
asymptotically normally distributed with zero mean and variance $\sigma^2$ which,
by independence of the two samples, is equal to
$$\sigma^2=\frac{1}{n_F}\sigma^2_F +\frac{\alpha^2}{n_G}\sigma^2_G.$$
A consistent estimator for $\sigma^2$, may be defined as
$$
\widehat{\sigma}^2=\frac{1}{n_F}\widehat{\sigma}^2_F+\frac{\alpha^2}{n_G}\widehat{\sigma}^2_G,
$$
where $\widehat{\sigma}^2_F$ and $\widehat{\sigma}^2_G$ are consistent estimators for ${\sigma}^2_F$ and ${\sigma}^2_G$, respectively. Thus, for checking the null hypothesis $H_0$, we may use the following test statistic :
$$ \widehat{T}=\frac{\Delta \widehat{J}}{\widehat{\sigma}},$$
which, by Slutsky's Theorem, converges in law to the standard normal distribution, under $H_0$.\\

\indent The previous test is distribution-free and may be extended to a vector of several particular poverty indices. To perform a joint test using simultaneously several poverty indices, we consider two finite $d$-dimensional vectors ($d$ is a positive integer) of poverty indices denoted by $I_F$ and $I_G$ and associated with distributions $F$ and $G$, respectively. The hypotheses we wish to test are  of the form :
$$H_0: I_F= M I_G,\qquad \text{versus}\qquad  H_1: I_F\neq M I_G,$$
where $M={\rm diag}(\alpha_1,\cdots,\alpha_d)$ is a diagonal matrix of positive coefficients $\alpha_i$. When the matrix $M$ coincides with the identity matrix $\textbf{I}_d$, hypothesis $H_0$ corresponds to the equivalence of the two distributions in terms of poverty. The test can performed by making use of the Wald test statistic which is defined as
$$ \widehat{W}= (\widehat{I}_F-M \widehat{I}_G)'\widehat{\Gamma}^{-1}_{FG}(\widehat{I}_F-M \widehat{I}_G), $$
where $x'$ designs the transpose of a vector $x\in\mathbb{R}^{d}$,$\; \widehat{I}_F, \widehat{I}_G$ are consistent estimators of ${I}_F$ and ${I}_G$, respectively.  $\widehat{\Gamma}_{FG}$ is a consistent estimator of the covariance matrix of the vector $I_F-M I_G$.  It is clear that, under $H_0$,  the statistic $\widehat{W}$  converges weakly  to a chi-square distribution with $d$ degrees of freedom, $\chi^{2}(d)$. Thus, at level of significance $\alpha$, the critical region is of the form $\{\widehat{W} > c\},$ where $c$ is the $(1-\alpha)$-quantile of $\chi^{2}(d)$.\\

 \indent To estimate the covariance matrix $\Gamma$ in Theorem \ref{t1}, denote $\widehat{\Gamma}_{k,l}$ the estimators of the entries ${\Gamma}_{kl}:=\Gamma[(w_l,f_l),(w_k,f_k)],\, k,l=1,\cdots,d$. Take $z$ equal to the quantile of order $q/n$ of a the considered distribution $F$ or $G$, represented by the sample  $Y_1,\cdots,Y_n$, with $q$ a positive  integer less than $n$. Then, the number of poor in the sample is equal to $q$, and a consistent estimators of the elements $\Sigma_{k,l}:=\Sigma\left[(w_{k},f_{k}); (w_{l},f_{l})\right]$,  $1\leq k,l\leq d$, may be defined as
 \begin{eqnarray} \label{eq:om}
\widehat{\Sigma}_{k,l} &= &
\frac{1}{n}\sum_{j=1}^{q}w_k\left(\frac{j}{n},\frac{q}{n}\right)f_k(Y_{j,n},z)w_l\left(\frac{j}{n},\frac{q}{n}\right)f_l(Y_{j,n},z) \nonumber \\
 &- & \frac{1}{n^2}\sum_{i=1}^{q}\sum_{j=1}^{q}w_k\left(\frac{i}{n},\frac{q}{n}\right)f_k(Y_{i,n},z)w_l\left(\frac{j}{n},\frac{q}{n}\right)f_l(Y_{j,n},z) \nonumber \\
& + & \frac{1}{n^2}\sum_{i=1}^{q}\sum_{j=1}^{q}a_1(Y_{i,n},Y_{j,n})\left(\frac{i}{n}\wedge\frac{j}{n}-\frac{i}{n}\frac{j}{n}\right)\nonumber \\
& +& \left(1-\frac{q}{n}\right) \frac{1}{n^2}\sum_{i=1}^{q}\sum_{j=1}^{q}a_2(Y_{i,n},Y_{j,n})\frac{j}{n}  \nonumber \\
& + & \left(1-\frac{q}{n}\right)\frac{1}{n^2}\sum_{i=1}^{q}\sum_{j=1}^{q}a_3(Y_{i,n},Y_{j,n})\frac{i}{n} \nonumber \\
& + &\frac{q}{n}\left(1-\frac{q}{n}\right)\frac{1}{n^2}\sum_{i=1}^{q}\sum_{j=1}^{q}a_4(Y_{i,n},Y_{j,n}),
\end{eqnarray}

where for all $i$, $w_i$ and $f_i$ represent respectively the weighting and deprivation functions of a particular poverty index, $a_r, r=1,2,3,4$ are real-valued functions given in Theorem \ref{t1},  and $Y_{1,n}\leq\cdots\leq Y_{n,n}$
are the order statistics associated with the sample $Y_{1},\cdots,Y_{n}$. It follows from this, that consistent estimators for the entries $\Gamma_{k,l},\, 1\leq k,l\leq d$, are
\begin{eqnarray}\label{vaR}
\widehat{\Gamma}_{k,l}&=&\widehat{\Sigma}_{k,l}+a\left(\frac{1}{n}\sum_{j=1}^{n}h_k(Y_{j,n})\zeta(Y_{j,n}) -\frac{1}{n^2}\sum_{j=1}^{n}h_k(Y_{j,n})\sum_{j=1}^{n}\zeta(Y_{j,n})\right) \nonumber \\
&+&a\left( \frac{1}{n}\sum_{j=1}^{n}h_l(Y_{j,n})\zeta(Y_{j,n})-\frac{1}{n^2}\sum_{j=1}^{n}h_l(Y_{j,n})\sum_{j=1}^{n}\zeta(Y_{j,n})\right) \nonumber \\
&+& \frac{a^2}{n}\sum_{j=1}^{n}\left(\zeta(Y_{j,n})-\frac{1}{n}\sum_{j=1}^{n}\zeta(Y_{j,n})\right)^2.
\end{eqnarray}
which lead to a consistent and non-parametric estimator for the covariance matrix $\Gamma$.

\section{Simulation study}
Here, we make some experiments for showing the asymptotic normality of our estimator $\hat{J}_n$ in relatively small samples of sizes $n=50,100,150.$
We essentially deal with two simple cases ; that is the case where the relative poverty line $z$ is taken equal to the mean of the distribution and the case where $z$ is set to the median of the distribution. The simulation procedure is the following :
\begin{itemize}
\item Generate data from a known distribution with positives values ;
\item Calculate the estimator $\hat{J}_n(w,f)$ and the theoretical indice $J(w,f)$ for specific functions $w,f$ ;
\item Compute the variance, say $\sigma^2_n(w,f)$, by using \eqref{vaR} ;
\item Determine the statistic test $T_n=\left|\frac{\hat{J}_n(w,f)- J(w,f)}{\sigma_n(w,f) }\right|$ ;
\item compute the p-value $p=2*(1-\phi(T_n))$, where $\phi$ is the standard Gaussian distribution;
\item Repeat all these steps $B$ times.
\end{itemize}
We generate  data from two distributions : Exponential(1/2) and Lognormal(0,1). In case of the  poverty line being equal to the mean of the distribution, the function $\zeta(y)=y$. While the median case corresponds to $\zeta(y)=\frac{1}{g(G^{-1}(1/2))}\mathbb{I}(y\leq G^{-1}(1/2))$. The tables below present the p-values of the normality test for different indices : FGT(1), FGT(2) and Sen, by considering a number of replications $B=1000$.

\begin{table}[htbp]
$$\begin{array}{|c|c|c|c||c|c|c|}
\hline
       &\multicolumn{3}{c||}{\text{Exponential}(1/2)} &  \multicolumn{3}{c|}{\text{Lognormal}(0,1)} \\
\hline
\text{Size} & \text{FGT(1)} &  \text{FGT(2)} & \text{Sen} &  \text{FGT(1)} &  \text{FGT(2)} & \text{Sen}\\
 \hline
n=50& 0.87   &0.86    &0.92   &0.89   &0.93    &0.72\\
 \hline
n=100&  0.93  &0.92    & 0.94   &0.91   &0.96    &0.81\\
 \hline
n=150& 0.95   & 0.93 & 0.94  &0.99   & 0.99   &0.89\\
\hline
\end{array}
$$
\caption{$p$-values of test in case the poverty line $z$ is equal to the mean of the distribution.}
\label{tab1}
\end{table}

\begin{table}[htbp]
$$\begin{array}{|c|c|c|c||c|c|c|}
\hline
       &\multicolumn{3}{c||}{\text{Exponential}(1/2)} &  \multicolumn{3}{c|}{\text{Lognormal}(0,1)} \\
\hline
\text{Size} & \text{FGT(1)} &  \text{FGT(2)} & \text{Sen} &  \text{FGT(1)} &  \text{FGT(2)} & \text{Sen}\\
 \hline
n=50& 0.24 &0.33    &0.37    &0.03   &0.21   &0.06\\
 \hline
n=100&  0.62  &0.74   & 0.80   &0.08   &0.40    &0.20\\
 \hline
n=150& 0.86  & 0.88   & 0.95  &0.17  & 0.62   &0.36\\
\hline
\end{array}
$$
\caption{$p$-values of test in case the poverty line $z$  is equal to the median of the distribution.}
\label{tab2}
\end{table}

Tables \ref{tab1} and \ref{tab2} allow us to accept the asymptotic normality of  our estimator for sample sizes greater than or equal to $n=100$ at a nominal level of 5\%.

\section{Application to genuine data}
Here, we use real data to estimate the quantity $\Delta$, representing the added term to the variance of the poverty estimate when we deal with  a relative poverty line $z$.
We employ Senegalese households expenditure data, which consist of a sample of size $n=3163$ provided  by a national survey entitled \textit{ESAM 2}, that was conducted in 2001 by the National Statistic Agence (ANSD). We consider two relative poverty lines : the mean and the median of the observed distribution. Denote the data by
$y_1\leq y_2\leq \cdots \leq y_n$. Then, an estimation of the quantity $\Delta$ is given by
\begin{eqnarray*}
\hat{\Delta}&=& 2\hat{a}\left[\frac{1}{n}\sum_{i=1}^{q} w\left(\frac{i}{n},\frac{q}{n}\right)f(y_i,z)\zeta(y_i) - \overline{\zeta(Y)}. \frac{1}{n}\sum_{i=1}^{q} w\left(\frac{i}{n},\frac{q}{n}\right)f(y_i,z)\right]\\
&& + \hat{a}^2\left[\frac{1}{n}\sum_{i=1}^{n}\left(\zeta(y_i)-\overline{\zeta(Y)}\right)^2\right],
\end{eqnarray*}
where
$$ \hat{a}= \frac{1}{n}\sum_{i=1}^{q}\left[ \frac{\partial}{\partial v} w\left(\frac{i}{n},\frac{q}{n}\right)f(y_i,z) +\frac{\partial}{\partial z} f(y_i,z)w\left(\frac{i}{n},\frac{q}{n}\right)\right],$$
$\overline{\zeta(Y)}=\frac{1}{n}\sum_{i=1}^{n}\zeta(y_i)$ and $q$ is the number of poor in the sample and satisfies : $G_n(z)=q/n$.\\
The following table \ref{tab3} gives estimations for the quantity $\Delta$ and the variance of the poverty estimate when the poverty line is fixed. The results concern the Sen index and the FGT indices of parameter $\beta=1,2.$\\

\begin{table}[htbp]
$$
\begin{array}{|c|c|c|c|}
\hline
    & \text{FGT(1)} &  \text{FGT(2)} &  \text{Sen}  \\
 \hline
\multirow{2}{*}{z=\text{mean}}& \hat{\Delta}=0.004 & \hat{\Delta}=0.01  & \hat{\Delta}=-0.05       \\
                                               & \hat{var}=0.07 &  \hat{var}=0.03   & \hat{var}=0.23 \\
\hline
\multirow{2}{*}{z=\text{median}}& \hat{\Delta}=1.44 & \hat{\Delta}=0.38  &\hat{\Delta}=1.63     \\
                                                   & \hat{var}=0.05  &     \hat{var}=0.01  &  \hat{var}=0.15 \\
\hline
\end{array}
$$
\caption{Estimation of $\Delta$ and  variances of poverty estimates when the poverty line is fixed.}
\label{tab3}
\end{table}

We observe that when the poverty line $z$ is set to the median of the distribution, the sampling error due the estimation of $z$ i increases the variance of the poverty estimate  for all the considered indices. In contrast, when the poverty line is set to the mean of the distribution, the variance of the FTG indices increases while the variance of the Sen measure decreases. This may be due to the fact that, the Sen measure affects a weight which depends on the poverty line $z$ and the individual's ranks.

\section*{Appendix}
\subsection{Proof of Theorem \ref{t1} }
First, recall the definition of the classes of functions $\mathcal{W}$ and $\mathcal{F}$
$$\mathcal{W}=\{w:[0,1] \times [0,1] \rightarrow \mathbb{R}_{+}, \; w
\;\;\text{  continuous,}\;\;\text{and}\;\;u\mapsto w(u,\cdot)
\;\;\text{is non-increasing}\}$$
$$\mathcal{F}=\{f:\mathbb{R}_{+}\times\mathbb{R}_{+}\rightarrow \mathbb{R}_{+},\;f
\;\;\text{ continuous,}\;\;\text{and}\;\;y\mapsto f(y,\cdot)
\;\;\text{is non-increasing}\}$$
Next, introduce the  class of functions
$$\mathcal{K}=\{k:\mathbb{R} \rightarrow [0,1] \;  \text{increasing} \}.$$
For  $z>0$ fixed, $ w\in\mathcal{W},f\in\mathcal{F}, k \in\mathcal{K}$  define the real-valued function
$$ h_{w,f,k}(y) = w[k(y),k(z)]f(y,z) \mathbb{I}(y <z),\quad\text{for all}\; y\in\mathbb{R}_{+}$$
and let $\mathcal{H}_z $ be the class of functions defined as
$$
\mathcal{H}_z  =  \{y\mapsto h_{w,f,k}(y) : w\in\mathcal{W},f\in\mathcal{F}, k \in\mathcal{K}\}.
$$
According to the sketch given at the end of the statement of the Theorem \ref{t1}, we split the proof  into four parts. In the first, we establish the Donsker property for the class $\mathcal{H}_z$, and derive from this, that the empirical process $\{\mathbb{G}_n(h_{w,f,G}) : w\in\mathcal{W},f\in\mathcal{F}\}$ converges weakly
to a limit Gaussian process $\mathbb{G}(h_{w,f,G})$.
In the second, we show that
\begin{equation}\label{84'}
\sup_{(w,f) \in \mathcal{W}\times\mathcal{F}}|\mathbb{G}_{n}(h_{w,f,G_n}-h_{w,f,G})|\longrightarrow_{p} 0,\quad n\rightarrow\infty,
\end{equation}
where "$\longrightarrow_{p}$" denotes the convergence in probability.
In the third part, we prove the weak convergence of the process $\mathbb{W}_n(w,f)$ to a zero-mean Gaussian process $\mathbb{W}(w,f)$ in $l^{\infty}(\mathcal{W}\times\mathcal{F})$. Finally in the last part, we prove that the joint process $(\mathbb{G}_n,\mathbb{W}_n)$ converges weakly to $(\mathbb{G},\mathbb{W})$ which is a zero-mean Gaussian process.

 \subsection{Part I}
Recall that $P$ is the common probability law of the $Y_j's$ and $G$ stands for its cumulative distribution function.
We have to prove that the class of functions $\mathcal{H}_z$ is $P$-Donsker. This will be done if we prove that  the bracketing integral
 $$
 J_{[]}(\infty,\mathcal{H}_z,L_2(P))=\int_0^{\infty}\sqrt{\log
N_{[]}(\epsilon,\mathcal{H}_z, L_{2}(P))}d\epsilon
 $$
 is finite, where $N_{[]}(\cdot)$ denotes the bracketing number. Before proving this, observe that the elements of
$\mathcal{H}_z$ are continuous and increasing functions, bounded on $\mathbb R_+$ by $w[k(0), k(z)]f(0,z)$, for every $(w,f,k) \in \mathcal {W \times F \times K}$. By assumption (A), the classes of functions $\mathcal{W}$ and $\mathcal{F}$ are pointwise measurable.
Further, Lemma 2.2 of \cite{vdg} entails that the $\delta$-entropy, relatively to the supremum norm,  of the class of increasing functions $\mathcal{K}$ is finite for any $\delta>0$.
 That is, the class $\mathcal K$ is totally bounded relatively to the supremum norm, and hence is pointwise measurable.
This enables us to take the supremum over the set $\mathcal{W} \times \mathcal{F}\times \mathcal{K}$ as equal to the supremum over a countable subset $\mathcal{G}_0\subset\mathcal{W} \times \mathcal{F}\times \mathcal{K}$. Since for  $z>0$ fixed, the quantity $w[k(0), k(z)]f(0,z)$ is finite for any $(w,f,k) \in \mathcal {W \times F \times K}$, we may define the constant function
$$
H(y)=\sup_{(w,f,k) \in \mathcal{W} \times \mathcal{F}\times \mathcal{K}}w[k(0), k(z)]f(0,z),\;\forall y\in\mathbb R_+,
$$
as an envelope function for the class $\mathcal{H}_z$. Then $\mathcal{H}_z$ is uniformly bounded by $H(y)$, and we may assume without loss of generality that
$H(y)\equiv 1$. Thus, $\mathcal{H}_z$ is a subset of the class of monotone functions defined on $\mathbb R$ with values in $[0,1]$. It follows from Theorem 2.7.5, page 159 of \cite{vdv} that for all $\epsilon >0$,
 \begin{equation}\label{sw00}
\log N_{[]}(\epsilon,\mathcal H_z, L_{2}(P))< C\epsilon^{-1},
\end{equation}
where $C$ is a positive constant.\\
From the fact that the elements of $\mathcal{H}_z$ take their values in $[0,1]$, for $\epsilon > 1 $ the number of $\epsilon$-brackets needed to cover $\mathcal{H}_z$ is just 1. Then $J_{[]}(\infty,\mathcal{H}_z,L_2(P))$ would be finite if $$ \int_0^{1}\sqrt{\log
N_{[]}(\epsilon,\mathcal{H}_z, L_{2}(P))}d\epsilon < \infty.$$
 Now, integrating both sides of (\ref{sw00}), one obtains
$$
\int_0^{1}\sqrt{\log
N_{[]}(\epsilon,\mathcal{H}_z, L_{2}(P))}d\epsilon< \sqrt{C}\int_0^1 \epsilon^{-1/2}d\epsilon=2\sqrt{C}<\infty.
$$
 That is, $J_{[]}(\infty,\mathcal{H}_z,L_2(P))$ is finite and the class $\mathcal{H}_z$ is $P$-Donsker.
In particular for $k=G$ (the distribution function  associated with the probability law $P$), the class  $\mathcal{H}_{z}$ restricts to
$$
\mathcal{H}_{z,G}=\{h_{w,f,G}:w\in\mathcal{W},f\in\mathcal{F}\},
$$
which may be identified to $\mathcal{W}\times\mathcal{F}$. Since $\mathcal{H}_{z,G}\subset\mathcal{H}_z$ is $P$-Donsker, so is the class $\mathcal{W}\times\mathcal{F}$. Then it follows that the empirical process $\{\mathbb{G}_n(h_{w,f,G}) : w\in\mathcal{W},f\in\mathcal{F}\}$
converges weakly in $l^{\infty}(\mathcal{H}_{z,G})$ to a tight limit process $\mathbb{G}$ , which is a zero- mean Gaussian process with covariance function defined, for all $(w,f)$ and $(\widetilde{w},\widetilde{f})$, by
\begin{eqnarray*}
 {\rm cov}(\mathbb{G}(h_{w,f,G}),\mathbb{G}(h_{\widetilde{w},\widetilde{f},G}))
&= & Ph_{w,f,G}h_{\widetilde{w},\widetilde{f},G}-Ph_{w,f,G}Ph_{\widetilde{w},\widetilde{f},G}\\
& = & \int^{z}_{0}w[G(y),G(z)]f(y,z)\widetilde{w}[G(y),G(z)]\widetilde{f}(y,z)dG(y) \\
 &-&\int^{z}_{0}w[G(y),G(z)]f(y,z)dG(y) \int^{z}_{0}\widetilde{w}[G(y),G(z)]\widetilde{f}(y,z)dG(y).
\end{eqnarray*}

 \subsection{Part II}
For establishing (\ref{84'}), we first remark that
for any $(w,f) \in \mathcal W \times \mathcal F$, the functions $h_{w,f,G}$  and $h_{w,f,G_n}$ are elements of $\mathcal H_z$, which is shown to be $P$-Donsker according to the preview part. Since $h_{w,f,G}$  and $h_{w,f,G_n}$ are bounded, they are in $L_2(P)=L_2(G)$.  Now, one has
\begin{eqnarray*}
\int_0^{\infty} \left[h_{w,f,G_n}(y)- h_{w,f,G}(y) \right]^2 dG(y) & \le & \sup_{y \le z} \left[h_{w,f,G_n}(y)- h_{w,f,G}(y) \right]^2\\
&\leq &\sup_{y \le z} |w(G_n(y),G_n(z))-w(G(y),G(z))|^2f^2(y,z)\\
&\le & f^2(0,z)\sup_{y \le z} |w(G_n(y),G_n(z))-w(G(y),G(z))|^2,
\end{eqnarray*}
which tends almost surely to 0, as $n\rightarrow\infty$, by  continuity of the function $w$ and the fact that
the empirical distribution function $G_n(y)$ converges almost surely to $G(y)$ for all $y\in\mathbb{R}$. Thus, as $n$ tends to infinity,
$\int_0^{\infty} \left[h_{w,f,G_n}(y)- h_{w,f,G}(y) \right]^2 dG(y)$
 converges almost surely and hence in probabilty to zero. It follows from Lemma 19.24 of \cite{vdv1} that $\mathbb{G}_n(h_{w,f,G_n}- h_{w,f,G}){\longrightarrow}_p 0$, $n \rightarrow \infty$ which, by the continuous mapping theorem, implies that
$$ \sup_{(w,f) \in \mathcal W \times \mathcal F}| \mathbb{G}_n(h_{w,f,G_n}- h_{w,f,G})|{\longrightarrow}_p 0, \ n \rightarrow \infty.$$
This establishes the second part of our proof.

\subsection{Part III}
 For any given functions $(w,f)\in\mathcal{W\times F}$, we define on the class $\mathcal K=\{k:\mathbb{R} \rightarrow [0,1] ,\;  \text{increasing} \}$ the following operator
 $$\phi_{w,f} : k \mapsto \phi_{w,f}(k)=\int_0^z w[k(y),k(z)]f(y,z)dP(y)=P h_{w,f,k},$$

Recall that
${\partial \over \partial u} \xi(a,b)$ and ${\partial \over \partial v}\xi(a,b)$ are the partial derivatives of a differentiable function $\xi(u,v)$ with respect to its first and second arguments, taken at $(u,v)=(a,b)$.
Let $$\mathcal{K}'=\{k\in\mathcal{K},\;k \;\text{ continuous} \}.$$
For all $k \in \mathcal K$, and $s_t \in \mathcal K$ such that $k+ts_t \in \mathcal K$ and $s_t \rightarrow s \in \mathcal K',$ as $ t \rightarrow 0$, one has by a first-order Taylor expansion of $w$, for some functions $\zeta$ and $\pi$ defined on $\mathbb R$,  with values in $ (0,1)$ :
\begin{eqnarray*}
 { \phi_{w,f}(k+ts_t) -  \phi_{w,f}(k) \over t}
&=&
\int_0^z  s_t(y) {\partial \over \partial u}w \left[k(y)+t\pi(t)s_t(y),k(z)+t\zeta(t)s_t(z)\right] f(y,z)dG(y)\\
&&+\int_0^z   s_t(z) {\partial \over \partial v}w\left[k(y)+t\pi(t)s_t(y),k(z)+t\zeta(t)s_t(z)\right] f(y,z)dG(y) \\
&=:& I_t+II_t.
\end{eqnarray*}
Now, we have to show that as $t \rightarrow 0$,
\begin{eqnarray*}
I_t & \longrightarrow & I=\int_0^z  s(y) {\partial \over \partial u}w \left[k(y),k(z)\right] f(y,z)dG(y) \\
II_t & \longrightarrow & II=\int_0^z   s(z) {\partial \over \partial v}w\left[k(y),k(z)\right] f(y,z)dG(y) .
\end{eqnarray*}
We only establish the first result as the other can be handled with the same techniques. By assumption (A.1) the function $w$ and its first-order partial derivatives are bounded on $(0,z]$ and one has :
\begin{eqnarray*}
|I_t-I| &\le&\sup_{y \le z}\Big|  s_t(y)
{\partial \over \partial u}w \left[k(y)+t\pi(t)s_t(y),k(z)+t\zeta(t)s_t(z)\right]\\&&- s(y) {\partial \over \partial u}w \left[k(y),k(z)\right]  \Big | f(y,z)\int_0^z dG(y).
\end{eqnarray*}
Adding and substracting appropriate terms and observing that both $k$ and
$s$ are bounded by 1, one has :
\begin{eqnarray*}
|I_t-I|
&\le&
\sup_{y \le z}\Big| s_t(y)-s(y) \Big |\times  \sup_{y \le z} \Big \{
\Big|{\partial \over \partial u}w \left[k(y)+t\pi(t)s_t(y),k(z)+t\zeta(t)s_t(z)\right]\Big|f(y,z) \Big \}
\\&&+ \sup_{y \le z} \Big \{\Big |
{\partial \over \partial u}w \left[k(y)+t\pi(t)s_t(y),k(z)+t\zeta(t)s_t(z)\right]-{\partial \over \partial u}w \left[k(y),k(z)\right] \Big |f(y,z) \Big \}.
\end{eqnarray*}
The fact that $s_t \longrightarrow s$, as $t \rightarrow 0$ entails that $|s_t(y)-s(y) | \longrightarrow 0$, as $t \rightarrow 0$. Consequently, the first term in the right-hand side of the above inequality tends to 0, as $t$ tends to 0. The second term also tends to 0, as $t$ goes to 0. This is due to the continuity of $w$ and its first-order partial derivatives.
It results from above that, $\phi$ is Hadamard-differentiable at $k \in \mathcal K$, tangentially to $\mathcal K'$, with      derivative $ \phi'_{w,f}[k]$, given for all $s \in \mathcal K$ by
\begin{eqnarray*}
 \phi'_{w,f}[k](s)&=&
\int_0^z \left\{ s(y) {\partial \over \partial u}w[k(y),k(z)] + s(z) {\partial \over \partial v}w[k(y),k(z)] \right \}f(y,z)dG(y).
\end{eqnarray*}
Since $\sqrt{n}[G_n - G]$ converge weakly to $ \mathbb B \circ G$, where $ \mathbb B$ stands for the standard Brownian bridge, it follows from the  functional delta method (see, e.g., \cite{vdv}) that $\sqrt{n}[\phi_{w,f}( G_n)-\phi_{w,f}(G)]=\sqrt{n}[Ph_{w,f,G_n}-Ph_{w,f,G}]=\mathbb{W}_n(w,f)$ converges in distribution to the Gaussian variable
\begin{eqnarray*}
 \phi'_{w,f}[G](\mathbb B \circ G)&=&
\int_0^z \left\{  \mathbb B \circ G (y) {\partial \over \partial u}w[G(y),G(z)] + \mathbb B \circ G(z) {\partial \over \partial v}w[G(y),G(z)] \right \}f(y,z)dG(y)\\
 &=: &\mathbb{W}(w,f).
\end{eqnarray*}
Since the class of functions $\mathcal{W\times\mathcal{F}}$ is shown to be Donsker according to Part $I$, we can infer that the process $\{\mathbb{W}_n(w,f): w\in\mathcal{W}, f\in\mathcal{F}\}$ converge in distribution to $\mathbb{W}(w,f)$ which is a zero-mean Gaussian process, with covariance kernel given, for all $(w,f)$ and $(\widetilde w, \widetilde f)$, by
\begin{eqnarray*}
{\rm cov}(\mathbb{W}(w,f),\mathbb{W}(\widetilde w, \widetilde f))& = &
\int_0^{z}\int_0^{z} a_1(x,y)[G(x)\wedge G(y) - G(x)G(y)]dG(x)dG(y), \nonumber  \\
&& + [1-G(z)]\int_0^{z}\int_0^{z} a_2(x,y)G(y)dG(x)dG(y) \nonumber  \\
&& +[1- G(z)] \int_0^{z}\int_0^{z} a_3(x,y)G(x)dG(x)dG(y)  \nonumber  \\
&& + G(z)[1-G(z)] \int_0^{z}\int_0^{z} a_4(x,y)dG(x)dG(y), \nonumber
\end{eqnarray*}
where
\begin{eqnarray*}
a_1(x,y) & =& \frac{\partial w}{\partial u}[G(x),G(z)]f(x,z)\frac{\partial \widetilde w}{\partial u}[G(y),G(z)] \widetilde f(y,z)\\
a_2(x,y) & =& \frac{\partial w}{\partial
u}[G(x),G(z)]f(x,z)\frac{\partial \widetilde w}{\partial v}[G(y),G(z)] \widetilde f(y,z)\\
a_3(x,y)& =&\frac{\partial w}{\partial v}[G(x),G(z)]f(x,z)\frac{\partial \widetilde w}{\partial u}[G(y),G(z)] \widetilde f(y,z)\\
a_4(x,y)& =& \frac{\partial w}{\partial v}[G(x),G(z)]f(x,z)\frac{\partial \widetilde w}{\partial v}[G(y),G(z)] \widetilde f(y,z).
\end{eqnarray*}

\subsection{Part IV}
Here we show that the couple of processes $(\mathbb{G}_n, \mathbb{W}_n)$  converges weakly to joint process $(\mathbb{G},\mathbb{W})$ which is a zero-mean Gaussian process. To this end, we show that it is tight and that its finite marginal distributions converge to those of a Gaussian process.

The tightness follows immediately  from Parts $I$ and $III$ where, it is proved that $\mathbb{G}_n$ converges weakly to a tight Gaussian process $\mathbb{G} \in l^{\infty}(\mathcal{H}_{z,G})$,  and $\mathbb{W}_n $ converges weakly to a tight Gaussian process $\mathbb{W} \in l^{\infty}(\mathcal{W\times F})$.

For the study of the finite dimensional distributions, we have to show that for all $\alpha_1, \ldots, \alpha_m,$   $\beta_1, \ldots, \beta_{\ell} \in \mathbb R$ and $(w_l,f_l), (\widetilde w_i, \widetilde f_i) \in \mathcal {W \times F}$, $l=1, \ldots,m$, $i=1, \ldots, \ell$, the linear combination
\begin{equation}\label{cl}
\sum_{l=1}^m \alpha_l\mathbb{G}_n(h_{w_l,f_l,G})+\sum_{i=1}^{\ell} \beta_i
\mathbb{W}_n(\widetilde w_i, \widetilde f_i)
\end{equation}
is asymptotically Gaussian. For this, we make use of the asymptotic linearity of the two processes $\mathbb{G}_n$ and $\mathbb{W}_n$. For larger values of $n$, the latter can be expressed in terms of the former. Indeed for all $(w,f) \in \mathcal {W \times F}$ denote by  $L_{w,f}$ the Hadamard derivative of $\phi_{w,f}$ at $G$ ; that is $L_{w,f}=\phi'_{w,f}[G]$. Then for larger values of $n$ one has
$$
\mathbb{W}_n(w,f)= \sqrt{n}(Ph_{w,f,G_n}-Ph_{w,f,G})=
L_{w,f}(\sqrt{n}[G_n -G])+o_P(1).
$$
Since $G_n (\cdot)=n^{-1}\sum_{j=1}^n \mathbb{I}(Y_j\leq \cdot)=n^{-1}\sum_{j=1}^n \mathbb{I}_{[Y_j, \infty)}(\cdot)$, using the linearity of $L_{w,f}$, we obtain for $n$ large enough that
\begin{eqnarray*}
\mathbb{W}_n (w,f
& = & \frac{1}{\sqrt{n}}\sum_{j=1}^n \left[ L_{w,f}( \mathbb{I}_{[Y_j, \infty)}(\cdot) - L_{w,f}(G)\right] +o_P(1)\\
&= & \mathbb{G}_n\left( L_{w,f}(\mathbb{I}_{[Y_j, \infty)}(\cdot)\right) +o_P(1).
\end{eqnarray*}
Combining this with the linearity of $\mathbb{G}_n$ we obtain, for all $\alpha_1, \ldots, \alpha_m,$ $ \beta_1, \ldots, \beta_{\ell} \in \mathbb R$ and $(w_l,f_l)$, $(\widetilde w_i, \widetilde f_i)$ $\in$ $ \mathcal {W \times F}$, $l=1, \ldots,m$, $i=1, \ldots, \ell$, that
\begin{eqnarray*}
&&\sum_{l=1}^m \alpha_l \mathbb{G}_n(h_{w_l,f_l,G})+\sum_{i=1}^{\ell} \beta_i
\mathbb{W}_n(\widetilde w_i, \widetilde f_i)\\
&=&
\sum_{l=1}^m \alpha_l \mathbb{G}_n(h_{w_l,f_l,G}) +
\sum_{i=1}^{\ell} \beta_i
\mathbb{G}_n \left(L_{\widetilde w_i, \widetilde f_i}(\mathbb{I}_{[Y_j, \infty)}(\cdot)\right)+o_P(1)\\
&=&
\mathbb{G}_n \left(\sum_{l=1}^m \alpha_l h_{w_l,f_l,G} +
\sum_{i=1}^{\ell} \beta_i L_{\widetilde w_i, \widetilde f_i}(\mathbb{I}_{[Y_j, \infty)}(\cdot)\right)+o_P(1).
\end{eqnarray*}
Recall that $\mathbb{G}_n$ is the empirical process and that the function $\sum_{l=1}^m \alpha_l h_{w_l,f_l,G} +
\sum_{i=1}^{\ell} \beta_i L_{\widetilde w_i, \widetilde f_i}(\mathbb{I}_{[Y_j, \infty)}(\cdot)$ belongs to $L^2(P)$.
Then it follows that the random variable defined in (\ref{cl})
is asymptotically Gaussian, and hence the finite marginal distributions of the process $(\mathbb{G}_n, \mathbb{W}_n)$,
$$\left( \mathbb{G}_n(h_{w_1,f_1,G}),\cdots,\mathbb{G}_n(h_{w_m,f_m,G}),\mathbb{W}_n(\widetilde w_1, \widetilde f_1),\cdots,\mathbb{W}_n(\widetilde w_{\ell}, \widetilde f_{\ell}) \right)  $$  are asymptotically Gaussian too. Combining this with the tightness argument enable us to conclude that the joint process $(\mathbb{G}_n, \mathbb{W}_n)$ converges weakly to the process $(\mathbb{G}, \mathbb{W})$ which  is Gaussian and centered.