Semiparametric Wavelet-based JPEG IV Estimator for endogenously truncated data

Semiparametric Wavelet-based JPEG IV Estimator for endogenously truncated data

\history{This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2019.2929571, IEEE Access. } \doi{10.1109/ACCESS.2019.2929571}

\address[1]{University of Haifa, Haifa, Israel (e-mail: [email removed])} \address[2]{University of Haifa, Haifa, Israel (e-mail: [email removed])}

\markboth

\corresp{Corresponding author: Moshe kim (e-mail: [email removed]).\\ This work was supported by the Research Authority, The University of Haifa.\\ The authors would like to thank seminar participants of the faculty of Mathematics and Computer Science, The Weizmann Institute of Science and Statistic departments of Tel-Aviv and Haifa universities, for very constructive comments.}

\titlepgskip=-15pt

\IEEEpeerreviewmaketitle

Introduction

\linespread{0.8} Scientists routinely try to model and extract causal relations among covariates, rather than merely their correlations.\footnote{For a specific form of causality due to treatment effect see angrist2004semiparametric definition.} In practice however, the presence of endogenous covariates in the model challenges the causal inference due to comovement of the random disturbance with these covariates. We distinguish between population induced comovement and training data (as in machine learning) comovement without having to rely on a covariate shift assumption since the behavioral (causal) model embedded in the training data does not necessarily describing the behavior in the entire population (see discussion in billfeld2019semiparametric).

The common way to overcome the aforementioned, is to generate a variation in the endogenous covariate without introducing variation in the random disturbance. This idea is achieved by employing a proper instrumental variable (IV).\footnote{Note that we deal with endogenously truncated sample selection model to differentiate from censored sample selection models heckman1979sample,newey2009two,powell2001semiparametric, where there exists information pertaining to the non-participants.}

Application of a proper instrumental variable generates variation in the endogenous covariate without introducing variation in the random disturbance and hence is orthogonal to it. Thus, IV should contribute to exogeneity and therefore has been extremely popular in empirical work.

Once we have analytically shown that the IV estimator is no longer valid in an endogenously truncated environment, we offer a truncation-proof estimator, which is a semiparametric wavelet-based JPEG-IV denoising algorithm.\footnote{A wavelet is a bandwidth-free estimator that is based on a multi-scale representation of the data. It is a widely used denoising technique xie2002sar.} This denoising algorithm decomposes the random disturbance into a noise and a systemic bias part, enabling the elimination of the truncation bias. The magnitude of the this bias is captured by the size of the wavelet coefficients which quantify and measure the degree of redundancy hidden in a sequence of data points. Consequently, this algorithm nests the conventional IV estimator as a special case due to the fact that in the absence of systemic endogenous truncation, the wavelet coefficients approach zero except for the intercept which describe the coarse level of the function.\footnote{Unlike Fourier transform, the wavelet estimator preservers not only the data average (coarse) behavior but also its local behavior capturing deviations (details) from the average. This fact renders our denoising suitable for irregular-spaced data which largely depend on local behavior and play an important role in the denoising.}

Our main contribution to the biorthogonal wavelet estimator is by formulating analytically (instead of approximating) the inverse of the transpose of wavelet transform without involving matrices which are computationally cumbersome voronin2015compression\footnote{“The implemented routine for the inverse transpose transform is approximate." voronin2015compression, p.285.} as well as the management of irregular-spaced data.\footnote{In the orthogonal wavelets design various interpolation methods are used to alleviate these irregularities hall1997interpolation and specific methodologies can be used to extend the Haar wavelet transform to the unequally spaced case sardy1999wavelet.} Additionally, our proposed methodology enables the combination of several penalty functions in the estimation procedure which are resolution-dependent.\footnote{It is known that soft thresholding provides smoother results relative to the hard thresholding because it is continuous. The latter, however, provides better edge preservation in comparison with the former.}

Wavelets are useful in denoising data. Several image quality assessment (IQA) measures have been introduced to choose the optimal level of denoising. These approaches can be classified to full-reference (FR) in cases the original image (noise-free) is observed; reduced-reference (RR) in cases where there is a partial information about the reference; reference-free (RF) in cases where the original image is not accessible carnec2003full,farah2014full. As we deal with truncated distributions, the source (the complete non-truncated distribution) is intrinsically unobservable and thus, we cannot assess the success of the denoising by comparing it to the original non-truncated distribution. Therefore, we select both the thresholding (tuning) parameter as well as the penalty function using a reference-free criterion function.

The proposed JPEG IV is biorthogonal, thus preserving both the symmetry (the original shape of the data) and compact support (small number of coefficients) properties of the data.\footnote{Our JPEG IV is a biorthogonal wavelet as it requires two sets of vectors, which are the dual basis and the series expansion sets, to obtain a denoised representation of the data. The elements in the former set are orthogonal to the corresponding elements in the latter set. See cohen1992biorthogonal for a formal definition of biorthogonality.} Importantly, the proposed methodology is easy to compute by precluding the need to find an optimal bandwidth as conventionally done.\footnote{Kernel estimation involves computational burden due to the necessity of finding the optimal bandwidth ichimura1993semiparametric. Unlike the nonparametric case, in the semiparametric context there is no “protocol" for finding the optimal bandwidth, as the traditional bandwidth choice methods might lead to bias estimates due to improper bandwidth choice lewbel2007nonparametric.} These properties make it suitable for denoising by alleviating both the problem of coefficient expansion as well as border discontinuities usevitch2001tutorial. The proposed algorithm corrects for both sources of bias: the endogeneity of covariates as well as the endogenous self-selection biases.

We run Monte Carlo simulations to measure the magnitude of the potential bias in the parameters' estimates under endogenous truncation, obtained by employing a conventional IV to eliminate the endogeneity bias. Our empirical implementation shows that even under mild correlation between the random disturbances, the resulting bias in the estimated parameter of the endogenous covariate in the substantive equation can amount to almost tenfold the true parameter value. Further, for sake of generality of the offered estimator, we subject it to various distributions in which the disturbances are neither jointly nor marginally normally distributed. These disturbances are constructed as realizations of non-symmetric and non-unimodal distribution functions.\footnote{Unlike the practice in some other studies applying only normally distributed disturbances.}

The rest of this paper is organized as follows. The methodology is presented is section (ref). Section (ref) prepares the ground for the biorthogonal wavelet. Section (ref) presents our proposed JPEG algorithm. In section (ref) we employ Monte Carlo simulations to validate our estimator performance. Section (ref) concludes.

Methodology

As discussed above, the IV is based on the following basic requirements: it is correlated with the endogenous covariate, as well as orthogonal to the random disturbance. Additionally, it must satisfy the exclusion restriction, such that in the presence of the endogenous covariate, the IV must be excluded from the regression. The IV is allowed to affect the dependent variable only through its effect on the endogenous covariate. However, the orthogonality condition is rarely satisfied in the presence of endogenous truncation, which is very frequently the nature of data used in empirical research, and therefore the IV will not provide a solution for the endogeneity problem. In what follows, we demonstrate the shortcoming of the conventional IV estimator, as well as potential bias generated in an environment of endogenously truncated data.

Suppose that there is a population random variable $\mathbf{\mathlcal{w}}=(\mathrm{z};\mathrm{x_{1}},\boldsymbol{{\mathrm{x_{-1}}}};\boldsymbol{{\mathrm{w}}})$ and that there is an independent and identically distributed sample $\left \{{z_i,x_{1i},\boldsymbol{{x_{-1_i}}}, \boldsymbol{{w_i}}} \right\}_{i=1}^N$ drawn from this population, referred to as the complete data set consisting of $N$ observations.\footnote{Capital letters indicate random variables; lower case letters indicate realizations of these random variables.} The instrumental variable is $\mathrm{z}$, the endogenous variable is $\mathrm{x_{1}}$ and the exogenous random variables are $(\boldsymbol{{\mathrm{x_{-1}};\mathrm{w}}})$, and where $\boldsymbol{{\mathrm{w}}}\in \mathbb{R}^l$ is a covariate vector.

Let $\xi_{1i}$, $\xi_{2i}$ and $\mathrm{v}_i$ be jointly dependent random disturbances with the respective marginal distribution functions $F_{\xi_1}$, $F_{\xi_2}$ and $F_{\mathrm{v}}$. Their joint distribution function is $F_{\xi_1, \xi_2,\mathrm{v}}$. The model is semiparametric, as neither the marginals nor the joint distribution function are required to be specified by the researcher.

The underlying model is composed of two parts. The first part consists of a selection equation, while the second part consists of the substantive (of interest) equation.

The population (non-truncated) selection equation is defined as:

where $\boldsymbol{{\gamma}}\in\mathbb{R}^l$ and $\boldsymbol{{w_i}}\in\mathbb{R}^l$ are the selection equation's coefficients and covariates vector, respectively. The selection equation's random disturbance is denoted by $\xi_{2i}$.

The substantive equation and the endogenous variable equation are defined as a system of equations:

where $\boldsymbol{{\beta}}\in\mathbb{R}^{p_1}$ and $\boldsymbol{{\delta}}\in\mathbb{R}^{p_2}$ are covariates vectors, $x_{1i}$ is an endogenous variable included in vector $\boldsymbol{{x}}_i\in\mathbb{R}^{p_1}$, and the exogenous variables are denoted by $\boldsymbol{{x}}_{-1_i}^T$. The substantive equation's random disturbances are denoted by $\xi_{1i}$ and $\mathrm{v}_{1i}$.

However the variables $y_{1i}^*, y_{2i}^*, x_{1i}^*$ are latent in the truncated environment and their respective observed realizations are denoted by $y_{1i}, y_{2i}, x_{1i}$, defined in (ref) and (ref) to follow.

The variable $y_{2i}^*$ is latent, while $y_{2i}$ is observed and defined as:

In the next section we reformulate the substantive equation as a partially linear single index model.

Semiparametric selectivity bias correction

The key difference between censored and truncated sample selection models is that in the former the entire covariate set (including the non-participants) and the selection variable are fully observed. In the latter, the entire data are truncated. Nevertheless, in both cases, the substantive equation can be represented as a partially linear regression, in which the dependent variable is observed only for the participants, as we are about to show. Following robinson1988root, the conditional expectation of the substantive equation in semiparametric (censored)\footnote{His approach is a generalization of the well-known inverse-mills ratio estimator introduced by heckman1979sample for the substantive equation's bias term $\mathbb{E}\left[ {\xi_{1i}|\xi_{2i} > -\boldsymbol{{w}}_i^T\boldsymbol{{\gamma}}} \right]$ in the case of a censored sample selection model. Note the difference between censored data and truncated data, which is the case we deal with.} sample selection models is some generally unknown function $\mathcal{M}_1(.)$ (to be estimated) of the selection equation's covariates variables $\boldsymbol{{w_i}}$:

such that $\boldsymbol{{\gamma}}$ is the selection equation's coefficient vector. Since $y_{1i}$ is observed only if $i$ is a participant, the substantive equation's dependent variable obtains the following functional form:

The regression equation in (ref) is referred to as a semiparametric partially linear regression (SP-NLS), in which the non-linear part is the bias term function. This regression can be estimated semiparametrically in cases of a truncated sample selection model using a non-linear least squares procedure as suggested by ichimura1993semiparametric.

Both ichimura1993semiparametric and robinson1988root models involve a kernel function estimation. However, kernel estimates' accuracy is sensitive to the bandwidth selected. This entails a potential problem of finding the optimal bandwidth resulting in computational complexity.\footnote{There is an open question whether there is a way to choose a bandwidth sequence that is optimal for the estimation of the parameters ichimura1991semiparametric.} Due to the lack of applicability of the traditional bandwidth selection methods in the semiparametric context, informal methods are being used, that may lead to a non-ignorable bias in the estimates lewbel2007simple.\footnote{“The well known bandwidth selection rules used in non-parametric estimation, such as cross validation, are not generally applicable to semiparametric settings.” lewbel2007simple} In order to avoid the problems involved with kernel estimation, our methodology relies on a (thresholding-propagated) nonlinear wavelet-based JPEG IV estimator to approximate the bias term (in (ref)).

The substantive equation depicted in (ref) deals with endogenous truncation bias, assuming that the random disturbance and the covariates are not jointly dependent. However, in cases where this random disturbance is jointly dependent with one (or more) of the covariates there will emerge two bias terms: the first one is propagated by the endogenous truncation and the second one is propagated by the endogenous covariate. Next we present a decomposition Theorem (ref), which enables reformulating the substantive equations as a partially linear single index model in the presence of an endogenous covariate.

Decomposition of the substantive equations

Next we formulate the relationship between the covariates and dependent variables in the equations to be estimated, in the presence of an endogenous covariate in the substantive equation under truncation.

Truncated sample selection model with an endogenous covariate

In cases where the substantive equation's dependent variable is a function of an endogenous covariate $x_{1i}$, both $x_{1i}$ as well as $y_{1i}$ (as in (ref)) are truncated, we face a truncated sample selection model with an endogenous covariate.

Thus, the semiparametric partially linear index model in a truncated environment consists of the following system of equations:

where $\epsilon_{1i}^{**}$ and $\epsilon_{2i}$ are two jointly dependent random disturbances,\footnote{There is dependence of these two random disturbances due to the dependence between $\mathrm{v_i}$ and $\mathrm{\xi_{1i}}$ (as in (ref)) in the complete (non-truncated) data.} and by construction are independent of the random variables vector $\boldsymbol{{\mathrm{w}}}$.\footnote{Not to be confused with its realization $\boldsymbol{{w_i}}$.} The intrinsic endogeneity in the model is captured by the joint dependence of $\epsilon_{1i}^{**}$ and the covariates.\footnote{The intrinsic model's endogeneity is related to the joint dependence of the random disturbance and the covariates in the population, unlike a conditional joint dependence of the random disturbance and the covariates given participation in the sample.} The presence of the function $\mathcal{M}_2(.)$ implies that we allow for a dependence between $\mathrm{v_i}$ (the endogenous part of $x_i$) and the selection equation's random disturbance $\xi_{2i}$ (in (ref), the complete, non-truncated, sample selection equation).

Our primary interest is to show that the instrumental variable and the random disturbance might be correlated in a truncated environment as will be depicted in Theorem (ref) to follow. By doing this, we denote a truncated environment using the indicator (selection variable) $\mathrm{s}=I(\xi_{2i}>-\boldsymbol{{\mathrm{w}}}^T\boldsymbol{{\mathrm{\gamma}}})$ and postulate the following assumptions:

These two assumptions implies that the conditional expectation of the instrumental variable, given the selection variable, is a function of the random variable vector $\boldsymbol{{\mathrm{w}}}$, as the following proposition argues:

In Theorem (ref) to follow we use proposition (ref) and present our primary argument: in truncated sample selection models, the orthogonality condition of the instrumental variable with respect to the random disturbance might be violated. This violation stems from a dependency between the instrumental variables and the selection equation's covariates.

However, the orthogonality condition can be satisfied by removing the contamination factor, which is the covariate generating the comovement between the random disturbance and the instrumental variable, as shown in the following Theorem (ref).

Therefore, a valid instrumental variable $\mathrm{z}$ is orthogonal to the truncated distribution (non-contaminated) disturbance $\epsilon_{1i}^{**}$ in (ref), even though $\mathrm{z}$ and $\boldsymbol{{\mathrm{w}}}$ are dependent.

The joint dependence of $(\xi_1,\xi_2,\mathrm{v})$ implies the violation of zero mean expectation (under truncation) in the $x_{1i}$ regression equation (ref), such that $\mathbb{E}[\mathrm{v}|\xi_{2}>-\boldsymbol{{w'\gamma}}]=\mathcal{M}_2(\boldsymbol{{w}}^T\boldsymbol{{\gamma}}) \ne \mathbb{E}[\mathrm{v}]=0$. That is, the conditional expectation of $\mathrm{v}$ given an endogenous truncation is a function of the covariate vector $\boldsymbol{{w}}$, while in the population it does not depend on $\boldsymbol{{w}}$ and has a zero mean expectation. This violation is a precondition for the endogeneity of $\left \{{\boldsymbol{{x_{-1i},z_i}}} \right\}$ with respect to $\mathrm{v_i}$ given participation in the regression of $x_{1i}$.\footnote{As been discussed in heckman1979sample, the fact that the conditional disturbance (given participation) in the substantive equation of $x_{1i}$ is a function of the selection equation's covariates, leads to a potential correlation between the disturbance and the substantive equation's covariates. This correlation implies the endogeneity of the substantive equation's covariates $\left \{{\boldsymbol{{x_{-1i},z_i}}} \right\}$ with respect to its random disturbance $\mathrm{v_i}$ given participation.} The following theorem indicates that such violation is also obtained in cases where the comovement of $\mathrm{v}$ and $\xi_2$ is entirely related to a variation in $\xi_1$.

Next we show that the conventional IV estimator is inconsistent in the presence of a truncated environment in which the expectation of the instrumental variable and the random disturbance are functions of the selection equation's covariates vector $\boldsymbol{{\mathrm{w}}}$. The proof in section (ref) to follow, relies on a linear dependence assumption between these two functions of $\boldsymbol{{\mathrm{w}}}$. The rationale for the linear dependence is due to the fact that the random disturbance's ($\xi_1$) conditional expectation generally satisfies monotonicity with respect to the index variable $\boldsymbol{{\mathrm{w}'\gamma}}$. Therefore, it is enough to assume that, on average, $\mathrm{z}$ is affected monotonically by the index variable $\boldsymbol{{\mathrm{w}'\gamma}}$ to generate a linear dependence between $\mathrm{z}$ and the conditional expectation of $\xi_1$ given participation.\footnote{Both functions are dependent through $\boldsymbol{{\mathrm{w}}}$ by construction, generally leading to some degree of linear dependence.}

The conventional IV estimator's asymptotic bias

The IV estimator's asymptotic bias is:

Given any correlation between $\boldsymbol{{\mathrm{z}}}$ and $\mathcal{M}_1(\boldsymbol{{\mathrm{w}}}^T\boldsymbol{{\mathrm{\gamma}}})$, $\underset{N\to \infty}{\mathrm{plim}}\left[ {\boldsymbol{{\mathrm{z}}}^T\mathcal{M}_1(\boldsymbol{{\mathrm{w}}}^T\boldsymbol{{\mathrm{\gamma}}})} \right]\not\to 0$. Thus, the $\boldsymbol{{\mathrm{\widehat{\beta}_{iv}}}}$ estimator is an inconsistent estimator for $\boldsymbol{{\mathrm{\beta}}}$.

Next we discuss the two types of joint dependence which are present in our model. This is done is order to facilitate the understanding of our proposed procedure, which is intended to correct for the bias propagated by each type of joint dependence.

Preliminaries

The objective is to eliminate the selection bias term captured by $\mathcal{M}_1(\cdot)$ in (ref). As we don't want to impose a specific distribution function on the random disturbances, the aforementioned elimination should be performed in a nonparametric manner. This can be achieved using a semiparametric estimation method, which is distribution-free. However, the bias term might be a discontinuous function with different levels of smoothness that must be considered. These issues can be alleviated using multi-resolution analysis by employing the wavelet estimator haar1910theorie. Wavelet is a bandwidth-free estimator, that is based on the idea of multi-scale representation of the data delouille2006second\footnote{Due to its multi-scale property, we can distinguish between the important information, the function's average behavior, from the noise. The coarse scales (lower resolution-levels) usually convey important information, while at fine scales there is usually more noise.} and is used as a denoising technique by simple thresholding, which is based on the concept of sparsity.\footnote{Sparsity implies that the majority of wavelet coefficients are small, and can be replaced by zero vanraes2002stabilised.}

The applicability of the classical wavelet estimator is problematic in several important aspects. First, it limits the sample size to be represented as $2^J$, with $J$ a non-negative integer, and the observations to be equispaced, which challenges the estimation in case of irregular-spaced data.\footnote{The observations location in space or time must be of equal distance.} Second, the classical wavelet estimator imposes the parametric assumption that the disturbances are independent identically distributed normal variables silverman1999wavelets. Lastly, there are the problems of coefficient expansion and border discontinuities.\footnote{The standard orthogonal wavelet transform has the shortcoming in that it requires a large number of coefficients (coefficient expansion) to represent the original data usevitch2001tutorial.}

In order to overcome these limitations, second generation wavelets have been introduced daubechies1998factoring which define wavelets in terms of lifting-steps instead of matrices to reduce computational complexity.\footnote{The lifting-steps are consecutive operations of prediction (scaled-moving average) and update (scaled-first difference) to obtain the wavelet coefficients.} An earlier attempt to deal with irregular-spaced data using second generation wavelets is presented in delouille2006second by postulating a prior distribution function for the wavelet coefficients.\footnote{delouille2006second adopt the parametric Bayesian denoising approach introduced by johnstone2005empirical, johnstone2004needles to obtain the wavelet coefficients assuming the coefficients are distributed according to a continuous mixture of a normal by a Beta density.} Alternative approaches extend Haar wavelet transform to accommodate for irregular data delouille2004smooth.

Both first as well as second generation wavelet estimation methods involve three steps: coefficient estimation (forward transform); (ii) denoising by using element-wise thresholding (coefficients selection) and (iii) reconstruction of the data without the noise (inverse transform). It is important to notice that the sequential nature of the estimation that relies on element-wise thresholding is applicable for limited types of wavelets, referred to as orthogonal wavelets which consist of the above described limitations. The main shortcoming of orthogonal wavelets is that the compact support and the symmetry properties which are useful in denoising are conflicting.\footnote{Unlike biorthogonality, orthogonality and symmetry are conflicting properties for design of compactly supported nontrivial wavelets (see Theorem 8.1.4 in daubechies1992ten).} To preserve both these properties, the biorthogonal wavelet-based JPEG is used chern1999interpolating, wei1998new.\footnote{The JPEG algorithm used here is termed `wavelet CDF 9/7'.}

In what follows we briefly explain the concept of biorthogonality. Denote a set of functions $\left \{{\varphi_{_{k}}(t)} \right\}$ which spans a vector space $\mathcal{F}$, referred to as the expansion set. By construction, any function $g(t)\in \mathcal{F}$ can be expressed by using a series expansion, such that $g(t)=\sum_{k}\eta_{_{k}}\varphi_{_{k}}(t)$,

where $\eta_{_{k}}$ and $\varphi_{_{k}}$ are the expansion coefficients and expansion functions, respectively. The set $\left \{{\varphi_{_{k}}(t)} \right\}$ is biorthogonal to the set $\left \{{\tilde{\varphi}_{_{k}}(t)} \right\}$ if $\langle \varphi_{_{k}},\tilde{\varphi}_{_{k'}}\rangle = \mathfrak{d}(k-k')$ $\forall k$ and $k'$, with $\langle \cdot \rangle$ being the $L_2$ inner product and the function $\mathfrak{d}(\cdot)$ is the Kronecker delta.\footnote{$\mathfrak{d}(k)=

$. The set $\varphi_{_{k}}(t)$ is orthogonal if $\langle \varphi_{_{k}},\varphi_{_{k'}}\rangle = 0$ $\forall k\ne k'$.} These two sets form a biorthogonal system, in which $\left \{{\tilde{\varphi}_{_{k}}(t)} \right\}$ is referred to as the dual basis of $\left \{{\varphi_{_{k}}(t)} \right\}$. Thus, we get the following unique representation:

Substituting each $\eta_k$ coefficient with its analytic expression in (ref), to obtain:

Obviously, in the present case of biorthogonality, the coefficients in (ref) are obtained by using the dual basis and the function is reconstructed in (ref) by using another basis which is the expansion set. In cases where $\left \{{\tilde{\varphi}_{_{k}}(t)} \right\}=\left \{{\varphi_{_{k}}(t)} \right\}$ we have an orthogonal basis $\left \{{\varphi_{_{k}}(t)} \right\}$, which is referred to as self-dual. Therefore, biorthogonality is a generalization of orthogonality that allows for a larger class of expansions.

Recall that our objective is to estimate the bias term for an unknown functional form, captured by $\mathcal{M}_1(\cdot)$ in (ref), the conditional expectation of $\varepsilon_{i}$, given participation defined as $\mathbb{E}[\varepsilon_{i}|y_{2i}=1]=\mathcal{M}_1(\boldsymbol{{w}}_i^T\boldsymbol{{\gamma}})$.\footnote{For brevity, we present $\mathcal{M}_1(\cdot)$ only. Identical treatment is applied to $\mathcal{M}_2(\cdot)$. } In what follows, we attend to the estimation of $\mathcal{M}_1(\cdot)$ using the wavelet estimator.

We use the concept of a frame in (ref) to define Riezs basis in (ref). Riezs basis is a building block in the definition of biorthogonal wavelets in (ref) to follow.

Let $\mathbb{H}$ be a separable Hilbert space with inner product $\langle\cdot,\cdot\rangle$ and a norm $\Big\lVert{\cdot}\Big\rVert_{2}^{2}$. We denote a sequence $\mathcal{F}=\left \{{f_k, k\in \Lambda} \right\}\subset\mathbb{H}$, in which $\Lambda \subset \mathbb{Z}$.

We use the following frame and Riesz basis definitions zalik1999riesz:

Let $L_2(\mathbb{R})$ be the space of square integrable and real-valued functions on $\mathbb{R}$. We use the following biorthogonal wavelet definition dicken1996wavelet:

It is worth noticing that both existence as well as uniqueness of the series representation are satisfied in definition (ref). However, our proposed nonparametric estimator might be unstable, rendering the estimation problem ill-posed, which is one of the challenges in nonparametric estimation of unknown functions.\footnote{An estimator violating at least one of the requirements: existence, uniqueness and stability is referred to as ill-posed.} This ill-posed problem can be alleviated by employing regularization on the wavelet series expansion coefficients abramovich1998wavelet,horowitz2014adaptive.

Let $u_i=y_{1i}-\boldsymbol{{x}}_i^T\boldsymbol{{\beta}}$ be the $i$'th element in vector $\boldsymbol{{u}}$ of size $n\times 1$, which satisfies:

where $\left \{{t_i} \right\}_{i=1}^{n}$ is a sequence in which the $i$'th element satisfies $t_i=\boldsymbol{{w}}_i^T\boldsymbol{{\gamma}}$ and $\epsilon_{1i}$ is the white noise described in (ref).\footnote{A more general formulation solves the ill-posed problem by employing regularization in cases where a linear transform of the unknown function replaces the original function abramovich1998wavelet.}

We use $\boldsymbol{{\Phi}}_{_I}$ and $\boldsymbol{{\Phi}}_{_F}\equiv \boldsymbol{{\Phi}}_{_I}^{-1}$ to denote the inverse and forward transformation matrices, respectively, each of size $n\times n$ he2011computer in an orthogonal wavelet.\footnote{The forward transform is referred to as the Discrete Wavelet Transform (DWT).} We note that using orthogonal wavelets one obtains the close-form solution to the wavelet coefficients as follows:

where $\widehat{\mathcal{M}_1}(\cdot)$ is the estimate of the unknown function $\mathcal{M}_1(\cdot)$, and $\rho_{_{\lambda}}(\cdot)$ represents the element-wise thresholding generated by some penalty function, in which the tuning parameter is represented by $\lambda$. The procedure in (ref) to obtain $\widehat{\mathcal{M}_1}(\cdot)$ by employing a thresholding operator is referred to as denoising.

We depart from the denoising procedure in (ref) by employing biorthogonal wavelets, as we are interested in the applicability of the general case where the penelization is not an element-wise due to correlations among wavelet regressors.\footnote{It has been shown that the performance of wavelet estimator can be improved when the dependencies among coefficients were taken into account tomassi2015wavelet.} In such cases, there is no such a close-form solution, which necessitates the regularized least squares optimization method to follow.

Let $\boldsymbol{{\Psi}}_{_I}$ and $\boldsymbol{{\Psi}}_{_F}\equiv(\boldsymbol{{\Psi}}_{_I}^T\boldsymbol{{\Psi}}_{_I})^{-1}\boldsymbol{{\Psi}}_{_I}^T$ denote the inverse and forward transformation matrices, respectively, each of size $n\times n$ of the wavelet-based JPEG, which is a biorthogonal wavelet.\footnote{For a definition of biorthogonal wavelets see daubechies1998factoring.}$^{,}$\footnote{ Unlike biorthogonality, orthogonality implies that the wavelet regressors are mutually uncorrelated and that the inverse transform is the transpose of the forward transform. This simplifies the computation as the wavelet coefficients are obtained analytically (closed-form) using element-wise thresholding operators (e.g., hard and soft thresholding operators). However, we opted for the biorthogonality wavelet to exploit the correlation structure of the regressors. Biorthogonal wavelets preserve the perfect reconstruction property (by employing dual-filters) as well, but is more flexible in that the inverse of $\boldsymbol{{X}}$ is not required to be its transpose. Consequently, the thresholding is applied to the entire coefficient vector. }

Let $\rho_{\lambda,\gamma} (\cdot)$ be the minimax concave penalty (MCP) function zhang2010nearly, defined as breheny2011coordinate:\footnote{The penalty function in (ref) represents a family of penalty functions as a generalization of the soft thresholding (if $\gamma\to\infty$) and hard thresholding (if $\gamma\to 1^{+}$) breheny2011coordinate.}

where $\theta\in(-\infty,\infty)$ is the parameter to be penalized, $\lambda>0$ and $\gamma\in (1,\infty)$.

We define resolution-dependent regularized least squares at resolution levels $1,...,J$:

where $\boldsymbol{{\delta}}=[\boldsymbol{{\delta}}_{1}^T,\boldsymbol{{\delta}}_{2}^T,...,\boldsymbol{{\delta}}_{J}^T]^T$ is the wavelet coefficient vector of size $n\times 1$ and $\boldsymbol{{\delta}}_j$ is of size $n_j\times 1$. $\left\lVert\boldsymbol{{\cdotp}}\right\rVert_{2}$ is the usual $\ell_2$ (Euclidean) norm, defined as $\left\lVert\boldsymbol{{b}}\right\rVert_{2}=\left( {\sum_{i=1}^{n}\left| {b_{i}} \right|^2} \right)^{1/2}$. The penalty function is $P_{\lambda_j,\gamma_j}({\boldsymbol{{\delta}}_j})=\sum_{k=1}^{n}\rho_{\lambda_j,\gamma_j}(\left| {{\delta_{j,k}}} \right|)$.

It is evident that when the $\boldsymbol{{\delta}}_j\to 0$, the bias propagated by the endogenous truncation approaches zero and thus, our algorithm is reduced to the conventional IV estimator.

The univariate solution of a regularized least squares problem using the penalty function in (ref) is denoted by $S_{\alpha}(\cdot)$ and defined as:\footnote{In order to utilize the min-max concave (MCP) penalty function in (ref), we depart from the regularized least squares algorithm in yang2015sparse, as it is limited to its special case of the LASSO penalty function. We introduce $\alpha$ as an approximation to the Hessian of the least squares problem in order to obtain an element-wise thresholding. This amounts to a dimensional reduction technique for reducing computational complexity. For the special case of $\alpha=1$, see breheny2011coordinate. }

where $\alpha\in(0,\infty)$. It is worth noting that if $\gamma\to \infty$ the solution is soft-thresholding introduced by donoho1994ideal; in case that $\alpha\gamma\to 1^{+}$ the solution is hard-thresholding (see proof in the Appendix (ref)).

To reduce computational complexity, the optimization problem in (ref) is reformulated as:

where $\Psi_{_I}^T$ is the transpose of matrix $\Psi_{_I}$, $\alpha \boldsymbol{{I}}$ is an approximation of the Hessian and $\boldsymbol{{I}}$ is the identity matrix of size $n\times n$. The number of iterations is denoted by the integer $iter$.

For brevity, we divide the argument to be minimized by $\alpha$ and complete the squares using the expressions in $\Big\lVert{\cdot}\Big\rVert_{2}^{2}$ yang2015sparse to get:

The iterative procedure performs MCP-thresholding on a proximal gradient-descent update for $k=1,...,n$ (see Algorithm (ref) in the Appendix):

where $\delta_{j,k}^{(iter)}$ is the $k$'th coefficient in vector $\boldsymbol{{\delta}}_{j}^{(iter)}$ and $\boldsymbol{{\psi}}_{_{I_k}}$ is $k$'th column in $\boldsymbol{{\Psi}}_{_I}$ (the inverse wavelet transform). The notation $\boldsymbol{{\psi}}_{_{I_k}}^T$ implies the transpose of $\boldsymbol{{\psi}}_{_{I_k}}$. We use (ref) to update the wavelet coefficients iteratively until the update is negligible, such that the following convergence criterion is satisfied:

where $\tau$ is the tolerance which is a positive real number that we arbitrarily set to $10^{-16}$.

The optimization method in (ref) involves matrices multiplication which is computationally infeasible for large data sets. To alleviate this computational complexity we develop a lifting scheme to be employed in order to perform simultaneously the transposed-inverse of the wavelet transform, consisting of lifting steps (see Algorithms (ref)-(ref) to follow). Conventionally, a lifting step can be either a prediction, that is a procedure generating a smoothed version of the data (the scaled coefficients), or an update that is the procedure to generate the remainder (the detail coefficients) between the data and its smoothed version. For the present case we define a new operator because the existing lifting steps do not provide an analytic representation of the transposed-inverse, as discussed in voronin2015compression.

In the next section we discuss the main idea behind lifting steps, in order to obtain analytically the transposed-inverse transform. First we describe the lifting steps in a regular-spaced data given a sample size of $2^{J}$ for a non-negative integer $J$. Then in equations (ref)-(ref) to follow, we alleviate these two restrictions by formulating our proposed algorithm.

Lifting steps to obtain the wavelet coefficients

Let $w=\left( {w_1,...w_n} \right)$ be a discrete sequence of data consisting of $n$ real numbers, such that the sequence is referred to as dyatic iff $n=2^J$ for some integer $J\ge 0$. The sequence can be expressed uniquely in terms of detail (difference) and summation coefficients denoted by $\left \{{d_{_{J-1},k}} \right\}_{k=1}^{n/2}$ and $\left \{{c_{_{J-1},k}} \right\}_{k=1}^{n/2}$, respectively. The former capture the variation in the sequence at different scales and locations and the latter are a smooth representation of the original sequence.

The multi-scale representation of a function $g\in L_2(\mathbb{R})$ is obtained as follows:

The first set of terms, $\phi_{_{0,k}}$, represents the average level of function $g$ and the second set of terms $\varphi_{_{j,k}}$ represents its details by accumulating information at a set of scales $j\in\mathbb{Z}$.

Let $\left \{{0,...,J-1} \right\}$ denotes a set of scales (resolution levels). We define $d_{_{J-1},k}$ and $c_{_{J-1},k}$ as follows nason2010wavelet:

The key idea is that a lower detail coefficient $d_{_{J-1},k}$ implies that $w_{2k}$ is very close to $w_{2k-1}$ and visa versa, as such a smoother function is represented by a small sequence of detail coefficients.

In order to represent the sequence in a coarser-scale (using a lower resolution), we define the coefficients:

By repeating the procedure in (ref) we obtain detailed and smoothed coefficients for lower resolutions. The multiscale algorithm stops when the $c_{_{0},1}$ coefficient is produced.

Next we discuss how to select optimally the thresholding (tuning) parameter in (ref) for each resolution-level.

Optimal thresholding by a reference-free criterion function

Since we deal with truncated distributions, the source (the complete non-truncated distribution) is intrinsically unobservable and thus, we cannot assess the success of denoising by comparing it to the original non-truncated distribution. Therefore, we utilize the “two-fold cross-validation" in nason1996wavelet which is a reference-free criterion function assesing the quality of the function estimated by denoising:\footnote{The methodology implemented in nason1996wavelet chooses one threshold that is applicable to all resolution levels in the wavelet transform. In the present case, however, we select a threshold for each level in order to implement multi-resolution analysis increasing our proposed estimator's accuracy. }

where $\lambda_{j,\gamma_{_j}}$ is the tuning (thresholding) parameter being used in (ref), in which $\gamma_{_j}$ is a specific penalty function. The odd sample and even sample are denoted by $\boldsymbol{{u}}_j^o$ and $\boldsymbol{{u}}_j^e$, respectively, and their corresponding estimates are $\widehat{f}_{\lambda_{j,m}}^o$ and $\widehat{f}_{\lambda_{j,m}}^e$. These estimates are obtained by employing the iterative procedure in (ref).

\color{black} As previously discussed, our proposed truncation-proof IV estimator requires controlling for the bias terms $\mathcal{M}_1(\cdot)$ and $\mathcal{M}_2(\cdot)$. For generality and applicability purposes of the proposed estimator, we adopt a semiparametric approach which is not subjected to distributional assumptions and consequently, does not require specifying the functional form of these unknown functions.

The wavelet-based JPEG semiparametric estimator is chosen for its many advantages. It enables a multi-resolution representation of the noisy data points, implying that the data points are characterized both globally as well as locally.\footnote{Global representation is a weighed average (smoothing) of the data, while local representation consists of more detailed information regarding first differences between neighboring data points. } Such a multi-resolution decomposition facilitates distinguishing between the noise and the systemic part. The systemic part is the functional relationship between the covariates and the dependent variables in the regression equations in (ref).

An additional advantage of incorporating the aforementioned newly introduced JPEG estimator is that it accommodates for various data set forms of different types of irregularities, such as non-equispaced design that will be described in section (ref) to follow. These irregularities are alleviated by introducing the locations in space of the various data points as an additional covariate that is unique for each resolution-level. An additional merit of our approach is enabling a group-wise denoising rather than the traditional element-wise JPEG denoising in the cases of image processing. Group-wise denoising plays an important role in data denoising, as it takes into account potential dependencies among the various data points. Thus, our contribution to the JPEG algorithm are controlling for irregularities, group-wise thresholding on the entire data and utilizing a reference-free criterion function to choose the optimal thresholding. \color{black}

In next section we describe the JPEG algorithm which is introduced to estimate each of the bias terms $\mathcal{M}_1(\cdot)$ and $\mathcal{M}_2(\cdot)$. Although our proposed denoising procedure can be applicable to both even as well as odd sample sizes (as will be demonstrated in Algorithm (ref) to follow), for ease of presentation and without loss of generality, the denoising procedure is formulated as a function of a data set consisting of $2n$ observations.

The wavelet-based JPEG denoising

Let $\left \{{(u_i, t_i)} \right\}_{i=1}^{2n}$ be a pairwise sequence of $2n$ data points as described in (ref), such that $t_i<t_j$ $\forall i<j$. The sequence $\left \{{u_i} \right\}_{i=1}^{2n}$ indicates the noisy data points (or colors of pixels in image processing) and their respective locations in space are represented by $\left \{{t_i} \right\}_{i=1}^{2n}$. The JPEG algorithm is a procedure generating a multi-resolution denoised representation of the sequence $\left \{{u_i} \right\}_{i=1}^{2n}$, which is denoted by the sequence $\left \{{\widehat{u}_i} \right\}_{i=1}^{2n}$. The purpose of the present section is three-fold: first, to describe the JPEG algorithm to be employed in order to obtain a noise-free representation of the noisy data; second, to extend the JPEG algorithm to be compatible with irregularities in the data;\footnote{We define $\Delta_i\equiv t_{i}-t_{i-1}$ $\forall$ $2\le i\le n$, such that equispaced (regular) data is a sequence of data points satisfying $\Delta_i=\Delta_j$ $\forall i$ and $j$. Other cases are referred to as non-equispaced (irregular) spaced data.} and thirdly, to incorporate a reference-free criterion to evaluate the denoising procedure accuracy.

Applying the conventional JPEG algorithm on a vector of data points is equivalent to employing three different procedures on the noisy data: (i) the JPEG forward transform $\mathcal{T}_{_{F}}:\mathbb{R}^{2n\times 2}\to \mathbb{R}^{2n\times 1}$ to obtain the wavelet-based JPEG coefficients (as will be shown in (ref) to follow); (ii) coefficients selection $\mathcal{T}_{_{S}}:\mathbb{R}^{2n\times 2}\to \mathbb{R}^{2n\times 1}$ by applying a thresholding procedure (as will be shown in (ref)) and (iii) the JPEG inverse transform $\mathcal{T}_{_{I}}:\mathbb{R}^{2n\times 2}\to \mathbb{R}^{2n\times 1}$, which recovers the noise-free data by utilizing the selected coefficients (as will be shown in (ref) to follow).

\color{black} In the ensuing section \color {black} we introduce auxiliary matrices to be used in each of the JPEG transforms, which are essential to construct the covariate matrix in the wavelet-based JPEG regression (in (ref) to follow).

Auxiliary matrices for the JPEG alogrithm

The implementation of the JPEG algorithm necessitates the construction of $\mathcal{T}_{_{F}}$ and $\mathcal{T}_{_{I}}$ operators. For this purpose, we construct auxiliary matrices $\mathcal{A}_{_{2n}}$, $\mathcal{S}_{_{2n}}$, $\left \{{\mathcal{H}^{\boldsymbol{{(t)}}}_{_{2n,\ell}}} \right\}_{\ell=1}^{4}$ which are the shifting, rescaling and smoothing operator matrices, respectively, each of size $2n\times 2n$.

Let $\mathcal{S}_{_{2n}}$ and $\mathcal{S}_{_{2n}}^{-1}$ be the rescaling and inverse-rescaling matrices, respectively each of size $2n\times 2n$. Its elements are defined for $m=0,...,n$ as:

The rescaling operator $\boldsymbol{{\tilde{v}}}=\mathcal{S}_{_{2n}}\boldsymbol{{v}}$ takes a vector $\boldsymbol{{v}}$ of size $2n\times 1$ and return a rescaled vector $\boldsymbol{{\tilde{v}}}$ of the same size, such that even and odd elements of the original vector are multiplied by the scalars $1/\varphi$ and $\varphi$, respectively.

Let $\mathcal{A}_{_{2n}}$ be a shifting operator matrix of size $2n\times 2n$, its elements are defined for $m=0,...,n$ as:

The $\boldsymbol{{\tilde{v}}}=\mathcal{A}_{_{2n}}\boldsymbol{{v}}$ operator takes a vector $\boldsymbol{{v}}=[v_{_{1}},...,v_{_{2n}}]^T$ of size $2n\times 1$ and return the vector $\boldsymbol{{\tilde{v}}}=[\boldsymbol{{v}}_{\text{odd}}^T, \boldsymbol{{v}}_{\text{even}}^T]^T$. The vectors $\boldsymbol{{v}}_{\text{odd}}=[v_{_{1}},...,v_{_{2n-3}},v_{_{2n-1}}]^T$ and $\boldsymbol{{v}}_{\text{even}}=[v_{_{2}},...,v_{_{2n-2}},v_{_{2n}}]^T$ consist of the odd and even elements of $\boldsymbol{{v}}$, respectively.

Unlike the conventional JPEG, we allow for data irregularities by controlling for the data set location in space. For doing so, we denote a sequence of matrices $\left \{{\mathcal{H}^{\boldsymbol{{(t)}}}_{_{2n,\ell}}} \right\}_{\ell=1}^4$, such that the elements of matrix $\mathcal{H}^{\boldsymbol{{(t)}}}_{_{2n,\ell}}$ $\forall \ell\in\left \{{1,2,3,4} \right\}$ of size $2n\times 2n$ are defined for $m = 1,...,n$ as:

where each of sequences $\left \{{\omega^{\boldsymbol{{(t)}}}_{_{\ell,l}}} \right\} \hspace{0.5em}\forall\ell \in\left \{{1,2,3,4} \right\}$ are the interpolation weights (defined in (ref) to follow) to control for the location in space of data points (enabling irregular non-equispaced data to be used) and $\pi_1,\pi_2,\pi_3,\pi_4$ are scalar constants described in schelkens2009jpeg, which are referred to as the filter coefficients of the wavelet-based JPEG. In the special case in which $\omega^{\boldsymbol{{(t)}}}_{_{\ell,l}}=0.5$ $\forall l$ and $\ell \in\left \{{1,2,3,4} \right\}$ the algorithm is reduced to the regular-spaced wavelet-based JPEG.

We define the linear interpolation weights:

For tractability, we formulate the JPEG coefficients estimation problem as a linear regression estimation, which necessitates obtaining a closed-form expressions of the JPEG forward and inverse transforms. These closed-form expression are required to characterize the JPEG covariate matrix to be used in the wavelet-based JPEG regression. \color{black} In the following section \color {black} we express analytically each of the forward and inverse transforms using matrix notation as a function of the auxiliary matrices presented above.

The various JPEG transforms in matrix notation

\color{black}Employing our proposed JPEG algorithm on a data set involves representation of data set in multiple resolution levels, a property which referred to as a multi-resolution analysis. Let $J$ be the highest resolution level, which requires the same number of data points as in the noisy data set. The data set representation in $j$'th resolution-level $\forall j<J$ is a transformation of the data set representation in the finer (higher) resolution-level $j+1$. Consequently, the JPEG noise-free representation $\forall j<J$ can be formulated recursively. However, the implication of this formulation is that the location in space of the data points in any given resolution-level is also determined recursively. This fact stems from depicting the noisy data set $\boldsymbol{{u}}=[u_{_{1}},...,u_{_{2n}}]^T$ and its location in space $\boldsymbol{{t}}=[t_{_{1}},...,t_{_{2n}}]^T$ as a pairwise sequence. For ease of notation we construct the adjusted space location operator $\forall j<J$: \color{black}

where $\boldsymbol{{\mathfrak{{m}}}}(j)\equiv \lceil2n/2^{^{J-j}}\rceil$, $J=\lceil\log_2\left( {2n} \right)\rceil$ is the number of resolution-levels and $j$ is a specific resolution level.\footnote{The operator's notation $\lceil\cdot \rceil$ represents the ceiling of a real number.}

This recursive formulation takes the the noisy data points locations in space as control variables, which are essential for alleviating irregularities in the noisy data.

Using the adjusted space location sequence $\left \{{\boldsymbol{{t}}_{_j}} \right\}$ in (ref), we define matrix $\boldsymbol{{\Psi_{_F}^{\boldsymbol{{(t)}}}}}ol{{(t)}}}}$ (to be used in (ref) to follow) for $J\in\left \{{1,...,\lceil\log_2\left( {2n} \right)\rceil} \right\}$ resolution levels as:

where $\boldsymbol{{\Phi}}^{(\boldsymbol{{t}})}_{m\times m}\equiv\mathcal{A}_{{_m}}\mathcal{S}_{{_m}} \mathcal{H}^{\boldsymbol{{(t)}}}_{_{m,4}} \mathcal{H}^{\boldsymbol{{(t)}}}_{_{m,3}} \mathcal{H}^{\boldsymbol{{(t)}}}_{_{m,2}} \mathcal{H}^{\boldsymbol{{(t)}}}_{_{m,1}}$ and $\boldsymbol{{I}}_{m\times m}$ is the identity matrix of size ${m\times m}$. It worth noticing that $\boldsymbol{{\Phi}}^{\boldsymbol{{(t)}}}_{m\times m}$ is a product of invertible matrices and consequently, its inverse is characterized as:

The multi-resolution JPEG forward transform operator $\mathcal{T}_{_{F}}(\boldsymbol{{u}},\boldsymbol{{t}})$ is the linear transform:

where $\boldsymbol{{u}}$, $\boldsymbol{{t}}$ and $\boldsymbol{{\delta}}$ are the noisy data, the location in space and the JPEG coefficient vectors, respectively, each of size $2n\times 1$.

Similarly, the multi-resolution JPEG inverse transform operator $\mathcal{T}_{_{I}}(\boldsymbol{{\delta}},\boldsymbol{{t}})$ is the linear transform:

where $\boldsymbol{{u}}$, $\boldsymbol{{t}}$ and $\boldsymbol{{\delta}}$ are the noisy data, the location in space and the JPEG coefficient vectors, respectively, each of size $2n\times 1$.

Matrix $\boldsymbol{{\Psi_{_I}^{\boldsymbol{{(t)}}}}}ol{{(t)}}}}$ in (ref) is defined for $J\in\left \{{1,...,\lceil\log_2\left( {2n} \right)\rceil} \right\}$ resolution levels as:

which $\boldsymbol{{\Psi_{_I}^{\boldsymbol{{(t)}}}}}ol{{(t)}}}}$ is the analytic inverse transform operator.

We introduce the wavelet-based JPEG nonparametric regression given a non-equispaced irregular data set:

where $\boldsymbol{{\delta}}$ consists of the sequences of coefficients $\left \{{c_{_{0,k}}} \right\}$ and $\left \{{d_{_{j,k}}} \right\}$, which capture the function's average behavior and details, respectively. The matrix $\boldsymbol{{\Psi_{_I}^{\boldsymbol{{(t)}}}}}ol{{(t)}}}}$ consists of the sequences of covariates $\left \{{\phi_{_{0,k}}(t)} \right\}$ and $\left \{{\phi_{_{j,k}}(t)} \right\}$. Unlike the present case which employs a group-wise denoising on the entire data, in the conventional JPEG the coefficients selection operator, $\mathcal{T}_{_{S}}( \boldsymbol{{\delta}},\boldsymbol{{t}})$, is constructed to be used element-wise (for each resolution-level $j$), e.g, $\mathcal{T}_{_{S}}( \delta_{ij},\boldsymbol{{t}})=S_{\alpha}(\delta_{ij}; \lambda_{j,\gamma_{_j}}, \gamma_{_j})$ using $S_{\alpha}(\cdot)$ in (ref), where $\lambda_{j,\gamma_{_j}}$ and $\gamma_{_j}$ are defined in (ref) using $\alpha=1$, such that $\left \{{\delta_{ij}} \right\}_{i=1}^{\boldsymbol{{\mathfrak{{m}}}}(j)}$ is a subset of vector $\boldsymbol{{\delta}}$ consisting of the $j$'th resolution-level coefficients.

In the following section, we discuss about the JPEG group-wise coefficients selection to perform denoising.

JPEG group-wise coefficients selection for irregular data denoising

Lastly, we formulate the transpose of the inverse wavelet transform, in order to employ the group-wise denoising procedure depicted in (ref):

where $\scaleobj{0.8}{\left[ {\left( {\boldsymbol{{\Phi}}^{\boldsymbol{{(t)}}}_{m\times m}} \right)^{-1}} \right]^T}$ is constructed as:

The JPEG estimator, denoted by $\widehat{\boldsymbol{{\delta}}}$, is obtained by minimizing the objective function in (ref) using the iterative procedure in (ref) given the chosen thresholding level. The latter is determined by minimizing the reference-free criterion function depicted in section (ref).

The construction of the wavelet-based JPEG transforms matrix involves computational complexity, a problem which can be alleviated by employing a faster algorithm, referred to as `a lifting step'. \color{black} In the succeeding section \color {black} we discuss the algorithm to obtain a denoised representation of non equispaced data design by using a procedure that does not necessitate matrix operation to reduce computational complexity.

The irregular forward transform

The irregular forward transform in (ref) is the procedure to obtain the wavelet coefficients, $\boldsymbol{{\delta}}$, as follows:

The Filter function in algorithm (ref) defined as follows:

The irregular inverse transform

The irregular inverse transform is the procedure to reconstruct the vector $\boldsymbol{{u}}$ in (ref), as follows:

The irregular transpose of the inverse transform

The irregular transpose of the inverse transform in (ref) enables to obtain the vector $\boldsymbol{{\tilde{\boldsymbol{{u}}}}}mbol{{u}}}}$, as follows (see transposed-inverse filter function, TransInvFilter, in algorithm (ref)):

In the ensuing section \color {black} we describe the estimation procedure of the parameter vector of interest, $\boldsymbol{{\beta}}$, using the wavelet-based JPEG estimate of $\widehat{\mathcal{M}}_1(\cdot)$ obtained from (ref).

The instrument variable estimator: truncated sample

Denote the truncated data by a sequence of observations $\left \{{y_{1i},\boldsymbol{{x}}_i,\boldsymbol{{w}}_i,\boldsymbol{{z}}_i} \right\}_{i=1}^n$, such that each observation is an independent realization of the conditional joint distribution function of the random variables $\left \{{\mathrm{y}_1,\boldsymbol{{\mathrm{x}}},\boldsymbol{{\mathrm{w}}},\boldsymbol{{\mathrm{z}}}} \right\}$ given that they are selected into the sample ($\mathrm{y_{2}}=1)$. The endogenous variable is denoted by $\mathrm{x}_1$ and is included in vector $\boldsymbol{{\mathrm{x}}}$. There are two types of joint dependence between the covariate vector and the substantive equation's random disturbance. The first type is intrinsic in the model and is generated by a variation in $\mathrm{v}$ (the endogenous part of $\mathrm{x_1}$) leading to a comovement between $\mathrm{x_1}$ and $\xi_{1}$. The second type is related to the sample selection and is generated by a variation in $\boldsymbol{{\mathrm{w}}}$ leading to a comovement between the covariate vector $\boldsymbol{{\mathrm{x}}}$ and $\xi_{1}$. This implies that there are two sources of endogeneity to be taken into consideration: the first source is related to the endogenous covariate, while the second source is due to the truncation environment of the data.

Next we discuss the two-step estimation procedure to be employed for the correction of both endogeneity and truncation bias propagated by truncation.

The estimation procedure

In this section we introduce a two-step estimation procedure to eliminate the two sources of bias discussed. To eliminate the endogeneity bias term we adapt a similar approach to the two step procedure in zhou2016estimation for a partially linear single index model estimation, in which the first stage is a regression of the endogenous covariate on all the exogenous covariates and the instrumental variable. In the second stage, the endogenous covariate is substituted with the fitted values obtained from the first stage. However, the estimation approach in zhou2016estimation cannot be implemented in a truncated environment, because it treats the first stage regression as a linear population regression (as if the entire covariates distribution function is observed). We alleviate this by modeling both the first as well as the second stage equations as endogenously truncated equations. In order to eliminate the endogenous truncation bias, we control for this source of bias by including the truncation bias term as an additional covariate in the substantive equations, as depicted in (ref). Thus, the partial linearity is applied to both the first as well as the second stage equations.

In the first stage, we regress the endogenous covariate on the instrumental and exogenous variables, by minimizing the partially linear index model:

In the second stage, the endogenous variable is replaced by its predicted value obtained from the first stage in (ref), and we minimize the following function:

As can be seen in (ref) the two sources of endogeneity bias we deal with are: (i) the bias propagated by the endogenous covariate is alleviated by utilizing the covariate set $\left[ {\widehat{x_{1i}},\boldsymbol{{x}}_{-1_i}^T} \right]$ consisting entirely of exogenous covariates and (ii) the bias propagated by the endogenous truncation is alleviated by controlling for the selection bias term $\widehat{\mathcal{M}}_1(\cdot)$.

Next we present Monte Carlo simulation to examine our semiparametric IV estimator's performance in a truncated environment.

Simulation

In this section, we generate multiple random data sets to be used for the examination of our model's performance, using different sample sizes.

First, we discuss the procedure for the data generation process (DGP).

Data generation process

Denote the sample size by $N\in\bigl\{500$, $2000$, $3000$, $5000$, $8000$, $10000\bigr\}$. In order to not restrict the data generation process to the family of symmetric unimodal distribution functions, a mixture of distribution functions is utilized to generate each of the selection model's disturbances that are jointly dependent (as will be discussed in section (ref) to follow). In order to verify that our proposed model performs well under different data generating processes (DGP), we construct a data set consisting of 2,000,000 distribution functions,\footnote{The estimates obtained given the various data distribution functions will be supplied upon request.} practically generating 100 millions realizations which are not i.i.d. By construction, each observation is randomly drawn from a unique mixture of distribution functions.

The disturbances' joint distribution function

Each triple of disturbances $\left \{{\xi_{1i},\xi_{2i},\mathrm{v}_i} \right\}$ is randomly and independently drawn from $F_{\xi_1,\xi_2,\mathrm{v}}$, which is the substantive and participation equations' disturbances joint distribution function. The aforementioned joint density function consists of two components: a Copula function,\footnote{Any continuous joint distribution function can be characterized by a set of marginal distribution functions and a joint distribution function determining the dependence structure which is referred to as a Copula function (Sklar's Theorem sklar1959fonctions). } which characterizes the disturbances' dependence structure, and three marginal distribution functions $F_{\xi_1}$, $F_{\xi_2}$ and $F_{\mathrm{v}}$. In order to verify our model's performance in the presence of random disturbances' distribution functions that are not restricted to the family of symmetric and unimodal distribution functions, each one of the sample selection model's disturbances $\xi_1$ and $\xi_2$ is marginally-distributed according to a mixture of three different distribution functions: (i) a normal distribution function with expectation and standard deviation parameters $(\mu, \sigma_a)$ denoted by $\mathcal{N}(\mu,\sigma_{a}^2)$; (ii) a normal distribution function with expectation and standard deviation parameters $(-\mu, \sigma_b)$ denoted by $\mathcal{N}(-\mu,\sigma_{b}^2)$; (iii) a gamma distribution function with scale and shape parameters $(\mu\varphi,\varphi)$ denoted by $\Gamma_{\scaleto{\mathcal{\text{Gamma}}}{4pt}}\left( {\mu\varphi, \varphi} \right)$\footnote{The scale and shape parameters imply that the expectation and standard deviation parameters are $(\mu ,\sqrt{\mu/\varphi})$, respectively.}. This mixture distribution function is defined as:

where $\mathbb{E}\left[ {\mathrm{\xi_{j}}} \right]=0$ and $\mathbb{E}\left[ {\mathrm{v}} \right]=0$ .

The parameters set $(\mu,\sigma_a,\sigma_b,\varphi,\sigma_{\mathrm{v}})=(4,2.5,1.5,2,1)$ is arbitrarily chosen. Due to its simplicity, the Clayton Copula (as will be discussed in section (ref) to follow) with a degree of dependence parameter is set to equal $1$, assuring a mild correlation between the disturbances, is used for controlling the dependence structure. Choosing a mild correlation, is important in order to be conservative by examining the potential bias in the parameter estimates under conditions which are not extreme.

Next we employ a function characterizing the dependence properties of the Copula mcneil2009multivariate, referred to as a generator function to construct the joint dependence of the random disturbances in (ref).

Archimedean Copula function

An Archimedean Copula is a Copula characterized by a non-increasing, continuous generator function $\psi$: $[0,\infty] \to [0, 1]$, which satisfies $\psi(0) = 1$, $\psi(\infty) = 0$ and is strictly decreasing on $[0, \inf\left \{{t : \psi(t) = 0} \right\}]$. In particular, we are interested in the $d$ dimensional Archemdean Copula family ($3$ in the present case\footnote{$d=3$ representing the three-dimensional vector of random disturbances $(\mathrm{v_i},\xi_{1i},\xi_{2i})$.}) which has the simple algebraic form mcneil2009multivariate:\footnote{Knowing the distribution corresponding to a generator $\psi$, marshall1988families presented a sampling algorithm for exchangeable Archimedean copulas which does not require the knowledge of the copula density. This algorithm is therefore applicable to large dimensions hofert2008sampling.}

where $\psi$ is a specific function known as the generator of $\mathcal{C}$. To generate the disturbances, the Clayton Copula's generator $\psi(t)=(1+t)^{-1/\theta}$ is chosen.

The covariates vector of $i$'th observation is a realization of the random variables $\left[ {\mathrm{z},\mathrm{x_2},\mathrm{w_1},\mathrm{w_2}} \right]$ which are jointly distributed with their corresponding marginal distribution functions $\left \{{\mathlcal{F}_{\mathrm{z}}^i, \mathlcal{F}_{\mathrm{x_2}}^i, \mathlcal{F}_{\mathrm{w_1}}^i,\mathlcal{F}_{\mathrm{w_2}}^i} \right\}$. The dependence structure is modeled by utilizing a Gaussian Copula which is a convenient way to generate high dimensional data. By construction, each datum is generated by utilizing a different sequence of marginal distribution functions (constructed as a finite mixture drawn from a menu of $2,000,000$ continuous distribution functions). These random variables expectation vector $\boldsymbol{{\mathrm{\mu}}}=[0,0,0,0]^T$ and a covariance matrix $\boldsymbol{{\mathrm{\Sigma_{4\times 4}}}}$. The arbitrarily chosen covariance matrix is:

We generate the data $y_{1i}, y_{2i}, x_{1i}$ according to the following data generation process (DGP) escanciano2017simple:

where each $i$ element in the sequence $\left \{{x_{2i},z_i,w_{1i},w_{2i}} \right\}_{i=1}^N$ is an independent realization of the random variables $(\mathrm{x_2,z,w_1,w_2})$. We choose the parameter setting $[\alpha_1,\alpha_2,\beta_1,\beta_2,\delta_1,\delta_2,\gamma_1,\gamma_2]=[2, 0.5, 1, 1.25, 0.5, 1 ,2, -1]$.

The truncated data set is characterized by the following equations:

where $x_{1i}$ is an endogenous variable included in vector $\boldsymbol{{x}}_i\in\mathbb{R}^p$, in which all the elements (except for $x_i$) are exogenous variables and $\boldsymbol{{\beta}}\in\mathbb{R}^p$ is a covariates vector. The substantive equation's random disturbance is denoted by $\xi_{1i}$.

Simulations result

We have randomly generated for each sample size $N\in\left \{{500,2000,3000,5000,8000,10000} \right\}$, a total of $10,000$ data sets using the data generation process elaborated on in (ref). For a given number of observations $N$, different models are estimated: (i) an $OLS$ estimator utilizing a sample consisting of random realizations from the complete distribution function, without correcting for the endogeneity of $x_{1i}$ covariate; (ii) a conventional IV estimator, correcting for the endogeneity of $x_{1i}$ covariate using the aforementioned entire distribution function; (iii) both $OLS$ as well as a conventional IV estimators are applied to a truncated portion of the data distribution function consisting of participants only (without correcting for the self-selection bias); (iv) truncated sample model's estimates using the developed wavelet-based JPEG IV estimator, correcting for both truncation as well as endogeneity biases.\footnote{For sake of brevity, we delegate results of the first stage to Appendix (ref).}

Table (ref) presents summary statistics of estimates for models (i) and (ii), while Table (ref) presents summary statistics of estimates for models (iii) and (iv). In Table (ref), different convergence measures of these estimates are presented.

\linespread{0.1}

Entries in Table (ref) indicate that regardless of sample size, the means of the $OLS$ estimates are biased, such that e.g., for a sample size of $3000$ observations $\beta_1=2.6807$ and $\beta_2=-0.6988$, while the mean of the full sample IV's estimates are $\beta_1=1.0025$ and $\beta_2=1.2457$ for the same sample size. For a sample size of $500$ observations, the standard deviation obtained for $\beta_1$ (the endogenous covariate's coefficient) using the IV estimator is $2.73$ times larger than in the $OLS$ estimator and decreases from $0.4397$ to $0.1006$, when the sample size increases from $500$ to $10,000$ observations.

In Table (ref) to follow, the estimates obtained by using the truncated data are presented. It is evident that the mean of the truncated sample $OLS$ estimates are $\beta_1=2.26$ and $\beta_2=-0.3084$ for a sample size of $500$ observations, where as the true parameter values are $\beta_1=1$ and $\beta_2=1.25$, respectively. This is a huge bias which is hardly improved as the sample size increases. Further, applying a conventional IV produces estimates which still represent a huge bias, particularly $\beta_1=0.1084$ and $\beta_2=2.0544$ for the same sample size (500 observations). \linespread{0.1}

Entries in Table (ref) indicate that regardless of sample size, the means of the truncated sample IV's estimates are biased (ranges from one-tenth to one-fifth of the estimate that would have emerged by employing the conventional IV method in the absence of truncation).\footnote{For sake of brevity, we have omitted the estimates of the nuisance parameters which can be furnished upon request as well as the parameter estimates of the first stage, which are delegated to Table (ref) Appendix (ref).} Note that estimates' accuracy hardly improved as sample size increases. This is due to the presence of two sources of bias. The mean estimate of $\beta_1$ (the endogenous covariate's parameter) which is obtained from implementing our proposed methodology, basically mimics the results obtained using a random sample from the entire data distribution function for sample sizes, above $2,000$ observations. The standard deviations of this estimate for sample sizes of $500$ and $10,000$ observations are $0.7325$ and $0.1577$, respectively. For a sample size of $5,000$ observations (or above), the mean estimate of $\beta_2$ (the exogenous covariate's parameter) approximates the estimate obtained by employing the conventional IV, using a random sample from the entire data distribution function. However, the estimate of $\beta_2$ obtained by employing the conventional IV in a truncated sample is biased even for $10,000$ observations.

We conduct sensitivity test to measure the influence of an increase in number of observations on the accuracy of the truncated sample's parameter estimates.

The first accuracy measure we use is the standardized root mean square error, $RMSE_j$, measuring the bias in the truncated regression estimate relative to the true parameter value that would have been obtained in an non truncated distribution, defined as:

where $\hat\beta_{i,j}^s$ and $\beta_{j}^s$ stand for the substantive (s) equation's $j$'th coefficient estimated in the $i$'th sample and the coefficient in the theoretical model that would have been obtained in the entire population, respectively. $\Omega$ is the number of data sets generated for the Monte Carlo simulations, which is $5000$ data sets (each one consists of $N$ observations).

Another measure is based on a formula similar to the one described in (ref), and is intended to find the relative accuracy of the truncated sample's estimates, in comparison to full sample estimates, defined as:

where $\hat\beta_{i,j}^{ts}$ and $\hat\beta_{i,j}^s$ stand for the substantive (s) equation's $j$'th coefficient estimated using the truncated (t) sample and the full sample, respectively. This measure evaluates the relative model's performance in the truncated sample, with respect to the conventional IV using the full sample.

The last estimates' accuracy measure is the $\delta$ coefficient used for the calculation of the estimators' standard deviations convergence rate $n^\delta$ with respect to the sample size. It depicts the speed of standard deviation's shrinkage resulting from increasing the sample size. This coefficient is calculated based on the following ratio:

where $\sigma_1$ and $\sigma_2$ are the estimate's standard deviations that are calculated for data sets with $n_1$ and $n_2$ number observations, respectively (calculated for a given estimate).

\linespread{0.1}

Table (ref) entries indicate that the root mean squares error (RMSE) measure of the estimates obtained by employing the conventional IV estimator, using a random sample from the entire data distribution function, gets smaller as the sample size increases, as can be expected. However, applying the same procedure to the truncated data set leads to RMSE measures, which are in the range of 2 to 8-fold larger, given a sample size of $2,000$ to $10,0000$ observations, respectively. This is indeed a huge bias generated by the conventional IV, which is not immune to truncation bias. Additionally, the RMSE measures show negligible improvements as a function of number of observations for the conventional IV, whereas there is a huge improvement of the RMSE, as a function of the number of observations for the JPEG IV estimator provided by our model. The proximity between the JPEG IV and the full sample IV estimates increases with the sample size, as reflected by the $R_j$ proximity measure. Using the same measure, we find that there is a much smaller improvement in the proximity between the truncated sample conventional IV and the full sample IV, relative to the improvement in the proximity between the JPEG IV and the full sample IV estimates.

It is evident that JPEG IV is a $\sqrt{n}$ consistent estimator, as depicted by the $\delta$ consistency measure, which is about $0.5$, implying that multiplying the sample size by 2 shrinks the estimators' standard deviations by $2^{\delta}=\sqrt{2}$. It is also evident that the truncated data conventional IV is poorly functioning in terms of consistency, as is shown by entries in Table (ref).

Conclusion

We provide an analytical proof showing that in an endogenously truncated data the conventional IV estimator does not perform the task it was intended to, but rather introduces an additional unintended bias into the parameter estimates of the substantive equation. The instrumental variable is endogenous by itself in the context of endogenously truncated data due to a comovement between the instrumental variable and the substantive equation's random disturbance, generated by mediating covariates. We offer a truncation-proof JPEG IV, shown to be a proper estimator under endogenous truncation. Employing Monte Carlo simulations attests to the JPEG IV estimator's high accuracy and its $\sqrt{n}$ consistency. These results have been verified by utilizing 2,000,000 different distribution functions (not restricted to the unimodal symmetric family), generating 100 million realizations to construct the covariates' data sets which are not imposed to be i.i.d. The various distribution functions attest to a very high accuracy of the model as depicted by the parameter estimates that closely mimic the true parameters.

\parskip = 0pt