The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
121,100 characters
Semiparametric inference for partially linear regressions with Box-Cox transformation
\let\WriteBookmarks\relax
\title{\bf
Semiparametric inference for partially linear regressions with Box-Cox transformation}
\shorttitle{SmoothMD for semiparametric partially linear regressions with Box-Cox transformation}
\author[1]{Daniel Becker}[]
\address[1]{PhD-Student in Economics, University of Bonn, Germany, Adenauerallee 24-26, 53113 Bonn(e-mail: [email removed])}
\author[2]{Alois Kneip}[]
\address[2]{Statistic Professor, University of Bonn, Germany, Adenauerallee 24-26, 53113 Bonn(e-mail: [email removed])}
\author[3]{Valentin Patilea}[]
\address[3]{Statistic Professor, CREST(Ensai), France, Campus de Ker-Lann, Rue Blaise Pascal - BP 37203 (e-mail: [email removed])}
\shortauthors{Becker et~al.}
\begin{abstract}
In this paper, a semiparametric partially linear model in the spirit of Robinson (1988) with Box-Cox transformed dependent variable is studied.
Transformation regression models are widely used in applied econometrics to avoid misspecification.
In addition, a partially linear semiparametric model is an intermediate strategy that tries to balance advantages and disadvantages
of a fully parametric model and nonparametric models.
A combination of transformation and partially linear semiparametric model is, thus, a natural strategy.
The model parameters are estimated by a semiparametric extension of the so called smooth minimum distance (SmoothMD) approach
proposed by \citet{lavergne2013smooth}. SmoothMD is suitable for models defined by conditional moment conditions
and allows the variance of the error terms to depend on the covariates. In addition, here we allow for infinite-dimension nuisance parameters.
The asymptotic behavior of the new SmoothMD estimator is studied under general conditions and new inference methods are proposed.
A simulation experiment illustrates the performance of the methods for finite samples.
\end{abstract}
\begin{keywords}
Semiparametric partially linear model, Nonparametric kernel estimators, Root N-consistent estimation, Conditional estimating equations, Hypothesis testing
\end{keywords}
\maketitle
\section{Introduction} \label{intro}
The data consists of independent copies of a response variable $Y$ and a random covariate vector
$\left(\boldsymbol X^T,\boldsymbol Z^T \right)^T \in \mathbb{R}^p \times \mathbb{R}^q$.\footnote{Herein, vectors are column matrices and for any matrix $\boldsymbol A$, $\boldsymbol A^T$ denotes its transpose.} We adopt a semiparametric approach to model the dependence of $Y$ on $X$ and $Z$.
Starting with \citet{robinson1988root}, the use of partial linear models of the form $Y= \boldsymbol X^T\boldsymbol \beta + m(\boldsymbol Z) + \varepsilon$ for some unknown, nonparametric function $m$ has become popular in this context. For an overview consider \citet{hardle2012partially} and \citet{li2007nonparametric}.
In many important economic applications, the response variable $Y$ is positive, i.e. $P(Y>0)=1$. It is then often questionable to assume (partial) linear regression models. A commonly used remedy is to apply a partial linear model to a suitable transformation of the response variable.
In econometric practice the dependent variable is then frequently $\log$-transformed (see for example \citet{acemoglu2001colonial} and \citet{autor2013putting}). But usually no substantial knowledge exists ensuring that this specific transformation leads to the correct model. The transformation proposed by \citet{box1964analysis}
\begin{align*} \label{Box_Cox}
T(Y,\lambda) =
\begin{cases}
\frac{Y^\lambda -1}{\lambda}&, \lambda \neq 0 \\
\log(Y) &, \lambda = 0
\end{cases}
\end{align*}
offers much more flexibility. Specifying the transformation up to a parameter and estimating the parameter together with the regression coefficients, leads to more reliable results at the cost of having to estimate only one additional parameter. In addition, the common $\log$-transformation is nested in the Box-Cox transformation and, thus, can be confirmed by a statistical test.
This motivates the approach adopted in this paper. To model the relationship between a positive response variable and the covariate vectors, we consider a \emph{transformation partially linear} mean regression model given by
\begin{equation} \label{main_model}
T(Y,\lambda) = \boldsymbol X^T\boldsymbol \beta + m(\boldsymbol Z) + \varepsilon,
\end{equation}
where $m(\cdot)$ is an unknown function and
\begin{align} \label{mom_con}
E[\varepsilon\mid\boldsymbol X,\boldsymbol Z] = 0.
\end{align}
The true values $\lambda_0$ and $\boldsymbol \beta_0$ of the parameters $\lambda>0$ and $\boldsymbol \beta\in \mathbb{R}^p$ are unknown and have to be estimated from an i.i.d. sample $(Y_i,\boldsymbol X_i^T,\boldsymbol Z_i^T)$, $i=1,\dots,n$, of $(Y,\boldsymbol X^T,\boldsymbol Z^T)$.
We impose no further assumption on the conditional distribution of $\varepsilon$. In particular, we allow for heteroscedasticity of unknown form. The vector $\boldsymbol Z$ contains only continuous variables, but the components of $\boldsymbol X$ need not be continuous.
We wish to note that the Box-Cox transformation is widely used in applications, and is discussed in various textbooks, e.g. \citet{amemiya1985advanced, greene2003econometric,horowitz2012semiparametric,showalter1994monte,wooldridge1992some}. Furthermore, there exist several empirical studies that employ the Box-Cox transformation. See, for instance, \citet{berndt1993empirical,heckman1974empirical} or \citet{keane1988real}. For an overview of the Box-Cox transformation, consider \citet{horowitz2012semiparametric} and \citet{sakia1992box}. The semiparametric partially linear specification of the conditional mean of the response seems to be quite appealing as it allows a linear dependency on a subvector $\boldsymbol X$ of covariates, which could include discrete variables, and meanwhile allows a nonparametric additive effect of the covariates $\boldsymbol Z$. These features could help practitioners faced with a large cross-sectional data set with independent observations including many candidate explanatory variables, who, on the basis of economic theory or past experience with similar data, feel able to parameterize only some of them.
Similar to other semiparametric models, the major challenge is to develop $\sqrt{n}$-consistent estimators and corresponding inference procedures for the parameters $(\lambda_0,\boldsymbol \beta_0^T)$.
This requires a careful
methodological development, since to our knowledge, there is no established procedure which can readily be applied under our general setup.
Despite its popularity, even in a purely parametric framework, estimation and inference in a Box-Cox transformation model is a difficult statistical problem and is usually based on quite restrictive assumptions on the conditional law of the response. See for instance chapter 5 of \citet{horowitz2012semiparametric} for an illuminating discussion. In particular, ordinary least squares estimation of
$(\lambda_0,\boldsymbol\beta_0^T)$ may lead to inconsistent results, and more sophisticated procedures, such as the nonlinear two-stage least squares (NL2SLS) estimator introduced by \citet{amemiya1981comparison}, have to be applied. The problem becomes even more complex in the case of the semiparametric regression (\ref{main_model}) where one only assumes the minimal identification condition (\ref{mom_con}).
In order to motivate our procedure, let us first consider the special case that the true value $\lambda_0$ is known a priori.
With $Y_0=T(Y,\lambda_0)$, we then arrive at a standard partial linear model
$Y_0=\boldsymbol X^T\boldsymbol \beta + m(\boldsymbol Z) + \varepsilon$. Consequently, $m(Z)=E(Y_0|Z)-E[\boldsymbol X\mid \boldsymbol Z])^T\boldsymbol \beta_0$ and
\begin{align*}
0= E[\varepsilon|\boldsymbol X,\boldsymbol Z] =E\bigg[ Y_0 - E[Y_0\mid \boldsymbol Z]- (\boldsymbol X-E[\boldsymbol X\mid \boldsymbol Z])^T\boldsymbol \beta_0\bigg| X,Z \bigg] .
\end{align*}
The basic idea of \citet{robinson1988root} now consists in disentangling the nonparametric estimation of unknown functions and the parametric estimation of the coefficient vector $\boldsymbol\beta$. In a first step Nadaraya-Watson kernel estimators are used for nonparametric estimation
of the functions $E[Y_0\mid \boldsymbol Z]$ and $E[\boldsymbol X\mid \boldsymbol Z]$. Plugging in these nonparametric function estimates,
an OLS regression of $Y_0- E[Y_0\mid \boldsymbol Z]$ on $\boldsymbol X-E[\boldsymbol X\mid \boldsymbol Z]$ then leads to an estimator
$\widehat{\boldsymbol\beta}$. If $q<4$ and standard nonparametric estimators based on second order kernels are applied, then, assuming homoscedasticity and suitable bandwidth sequences, Robinson showed that $\widehat{\boldsymbol\beta}$ is a $\sqrt{n}$-consistent, asymptotically normally distributed and efficient estimator of $\boldsymbol\beta_0$. If the dimensionality of $Z$ is larger, i.e. $q\ge 4$, then $\sqrt{n}$-consistent estimators can only be achieved by using higher order kernels.
The quite straightforward way to build efficient estimators made the partially linear model quite a popular. Versions of this model have also been studied by \citet{engle1986semiparametric,heckman1986spline,shiller1984smoothness} and \citet{wahba1984partial}. In order to avoid the trimming introduced by \citet{robinson1988root} to ensure that the estimate of the density of $\boldsymbol Z$, $f_z(\boldsymbol Z)$, stays away from zero, \citet{li1996root} considered as starting point the unfeasible OLS regression of $(Y_0- E[Y_0\mid \boldsymbol Z])f_z(\boldsymbol Z) $ on $(\boldsymbol X-E[\boldsymbol X\mid \boldsymbol Z]) f_z(\boldsymbol Z) $. Premultiplying by the density of $\boldsymbol Z$ does not break the consistency of the unfeasible OLS estimator since $E[f_z(\boldsymbol Z) \varepsilon\mid \boldsymbol X,\boldsymbol Z] = f_z(\boldsymbol Z) E[\varepsilon\mid \boldsymbol X,\boldsymbol Z] = 0.$ Next, \citet{li1996root} proposed to build OLS estimates using standard kernel estimators instead of the unfeasible response and covariates. This new estimator is still $\sqrt{n}$-consistent and asymptotically normally distributed. Moreover, \citet{li1996root} relaxed the condition on the bandwidth with the consequence that the smoothing requires higher order kernels only if the dimension of $ \boldsymbol Z$ is larger than 5, instead of larger than 3 as required in \citet{robinson1988root}.
Let us now return to the general model \eqref{main_model} with unknown parameter $\lambda$. Adopting Li's idea of premultiplying with the density $f_z$ of $\boldsymbol Z$, the conditional moment condition $E[\varepsilon|\boldsymbol X,\boldsymbol Z] = 0$ leads to
\begin{equation} \label{CMC1}
E\left(\bigg(T(Y,\lambda) - E[T(Y,\lambda)\mid \boldsymbol Z]- (\boldsymbol X-E[\boldsymbol X\mid \boldsymbol Z])^T\boldsymbol \beta\bigg)f_z(\boldsymbol Z)\bigg| X,Z\right)=0\quad \Longleftrightarrow
\quad \lambda=\lambda_0, \boldsymbol\beta=\boldsymbol\beta_0.
\end{equation}
It is now immediately seen that the unknown, additional parameter $\lambda$ introduces a major complication. Unlike the standard partial linear model, there is no way to disentangle nonparametric estimation of unknown functions and parametric estimation of the coefficient vector $(\lambda_0,\boldsymbol\beta_0^T)$. The reason is that $ E[T(Y,\lambda)\mid \boldsymbol Z]$ depends on $\lambda$. Indeed, there does not seem to exist a straightforward generalization of Robinson's approach which is able to cope with model \eqref{main_model}.
On the other hand, by \eqref{CMC1}, our model belongs to the large class of models identified by conditional moment restrictions. Methodologically however, we have to deal with
the obvious facts that
a) the model is highly nonlinear in $\lambda$ and b) \eqref{CMC1} incorporates an infinite dimensional nuisance parameter $\boldsymbol \eta_\lambda$ consisting of the functions
$f_z(z)$, $E[T(Y,\lambda)\mid \boldsymbol Z=z]$ and $E[\boldsymbol X\mid \boldsymbol Z=z]$.
Even if $\boldsymbol\eta_\lambda$ were known a priori, any use of the generalized method of moments (GMM) runs into the problem that the conditional moment restrictions identifying our model imply an infinite number of unconditional moment restrictions, since the conditioning variables have a support with infinite cardinality. But GMM relies only on a finite number of instruments and, thus, in general, consistency of GMM requires regularization and additional assumptions. See \citet{dominguez2004consistent}. This problem has already been pointed out for the Box-Cox transformation by \citet{foster2001estimation} and \citet{shin2008semiparametric} in the linear case. See also \citet{horowitz2012semiparametric}.
More recent work explicitly focuses on regularization techniques in order to account for complex conditional moment conditions.
Some methods rely on increasing the number of considered unconditional estimating equations (or instruments) with the sample size, such as the sieve minimum distance (SMD) approach of \citet{ai2003efficient}, or generalizations of GMM and empirical likelihood (EL) by \citet{donald2003empirical} and \citet{hjort2009extending}. \citet{carrasco2000generalization} use a regularization approach to generalize the GMM approach to a continuum of estimating equations. Other EL-type estimators use nonparametric smoothing to estimate conditional equations, such as \citet{antoine2007efficient}, \citet{kitamura2004empirical}, and \citet{smith2007efficient, smith2007local}. All these approaches share one common feature. The estimators’ sensitivity to the user-chosen parameter (number of estimating equations, regularization parameter, or smoothing parameter) remains largely unknown.
In this paper, we rely upon the SmoothMD approach proposed by \citet{lavergne2013smooth} to estimate $(\lambda_0,\boldsymbol\beta_0^T)$. Roughly speaking, SmoothMD can be seen as a new technique to translate conditional moment conditions into unconditional ones which can be approximated by sample averages. Although the method involves some tuning parameters, an attractive feature consists of the fact that a practical choice is quite uncritical, since asymptotic results can be established for a wide range of possible values of these tuning parameters (including values independent of the sample size). SmoothMD thus bridges a gap between Dominguez and Lobato’s method, which does not require a user-chosen parameter, and the competing SMD estimator and EL and GMM-type methods that rely on smoothing with restrictive conditions on the choice of smoothing parameters. Furthermore, although \citet{lavergne2013smooth} rely on a more standard setup. We will show that this technique can be well adapted to deal with complex functional nuisance parameters.
In order to explain the methodology, we introduce some abbreviations in order to simplify the lengthy expressions in \eqref{CMC1}. Let $\boldsymbol W=\left(\boldsymbol X^T,\boldsymbol Z^T \right)^T \in \mathbb{R}^p \times \mathbb{R}^q$ and $\boldsymbol U=\left(Y,\boldsymbol W^T \right)^T$.
Moreover, set $\boldsymbol\theta = \left(\lambda, \boldsymbol{\beta}^T\right)^T$, and for a real value $\gamma$ define
\begin{equation}\label{def1}
g(\boldsymbol U; \boldsymbol \theta,
\boldsymbol \eta_\lambda) = \left(T(Y,\lambda) - E[T(Y,\lambda)\mid \boldsymbol Z]- (\boldsymbol X-E[\boldsymbol X\mid \boldsymbol Z])^T\boldsymbol \beta\right)f_z(\boldsymbol Z) .
\end{equation}
Recall that $\boldsymbol \eta_\lambda$ is an infinite-dimensional nuisance parameter defined by
$$\boldsymbol \eta_\lambda=\boldsymbol\eta_\lambda(z)=(f_z(z), E[T(Y,\lambda)\mid \boldsymbol Z=z],E[\boldsymbol X\mid \boldsymbol Z=z]^T).$$
Condition \eqref{CMC1} is then equivalent to requiring
\begin{equation} \label{CMC2}
E\left(g(\boldsymbol U; \boldsymbol \theta,
\boldsymbol \eta_\lambda)\big| W\right)=0\quad \Longleftrightarrow \quad
\boldsymbol\theta = \left(\lambda, \boldsymbol{\beta}^T\right)^T=
\left(\lambda_0, \boldsymbol{\beta}_0^T\right)^T =:\boldsymbol\theta_0 .
\end{equation}
SmoothMD is based on the following insight: Let $\boldsymbol U_1$ and $\boldsymbol U_2$ be two independent copies of $\boldsymbol U$ with corresponding subvectors $\boldsymbol W_1$ and $\boldsymbol W_2$. For any symmetric function $\omega(\cdot)$ of $\boldsymbol W$ with positive Fourier transform, we than have that
\begin{equation*}
E\left(g(\boldsymbol U; \boldsymbol \theta,
\boldsymbol \eta_\lambda)\big| W\right)=0\quad \text{if and only if}\quad Q(\boldsymbol \theta) =
E[g(\boldsymbol U_1; \boldsymbol \theta,
\boldsymbol \eta_{\lambda,1} )g(\boldsymbol U_2; \boldsymbol \theta,
\boldsymbol \eta_{\lambda,2})
\omega ( \boldsymbol W_1 - \boldsymbol W_2) ] = 0.
\end{equation*}
Here $\boldsymbol \eta_{\lambda,1} $ and $\boldsymbol \eta_{\lambda,2} $ are the vector of nuisance functions corresponding to $\boldsymbol Z_1$ and $\boldsymbol Z_2$.
Indeed, in Lemma \ref{lem_ident} of Section \ref{identification}, it will be shown that
\begin{align} \label{MGA}
Q(\boldsymbol \theta) =
\begin{cases}
\hskip 0.35cm 0 &, \text{ if } \boldsymbol\theta = \left(\lambda, \boldsymbol{\beta}^T\right)^T=
\left(\lambda_0, \boldsymbol{\beta_0}^T\right)^T , \\
>0 &, \text{ else}.
\end{cases}
\end{align}
\citet{lavergne2013smooth} list several possible ways to define $\omega(\cdot)$, but throughout this paper we will rely on the simple choice $\omega (\boldsymbol W) := \exp\left\{-\boldsymbol W^T \boldsymbol D\boldsymbol W\right\}$, where $\boldsymbol D$ is a diagonal matrix whose positive diagonal elements $d_1,\dots,d_{p+q}$ represent user selected tuning parameters. A sensible choice consists of using the
standard deviations of the components of the vectors $(\boldsymbol X_i^T, \boldsymbol Z_i^T)^T$.
If the functions which define $\boldsymbol \eta_\lambda$ were known, then a data-based estimator
of the unconditional moment $Q(\boldsymbol \theta)$ could be obtained by the sample averages
\begin{equation} \label{MGA1}
Q(\boldsymbol \theta) = \frac{1}{n^2} \sum_{1\leq i,j \leq n}g(\boldsymbol U_i; \boldsymbol \theta, \boldsymbol \eta_{\lambda,i} )
g(\boldsymbol U_j; \boldsymbol \theta, \boldsymbol \eta_{\lambda,j} ) \omega ( \boldsymbol W_i - \boldsymbol W_j),
\end{equation}
where $\boldsymbol \eta_{\lambda,i} $ and $\boldsymbol \eta_{\lambda,j} $ are the nuisance parameter values corresponding to $\boldsymbol Z_i$ and $\boldsymbol Z_j$, respectively.
The average scheme proposed by \citet{lavergne2013smooth} relies on leaving out diagonal elements with $i=j$. In our setup, inclusion of these diagonal terms provides more stable and reliable estimators. We will show that the resulting bias is asymptotically negligible. Note that by definition of $\omega(\cdot)$
$$
E\left[ Q_n\left(\boldsymbol \theta\right)\right]
=
Q(\boldsymbol \theta ) + \frac{1}{n} E\left( g(\boldsymbol U; \boldsymbol \theta,\boldsymbol \eta_\lambda)^2\right)
=
Q(\boldsymbol \theta) + O(n^{-1}).
$$
Under model \eqref{main_model}, the nuisance parameter $\boldsymbol \eta_\lambda$ is unknown, and $Q_n(\cdot)$ cannot be directly computed. We therefore use kernel estimation to determine nonparametric estimators $\widehat{\boldsymbol \eta}_\lambda$. This then leads to a feasible version $\widehat{Q}_n(\cdot)$. More precisely, our estimation procedure can be described as follows. For each $\lambda$, we define the map
$$
\boldsymbol \beta \mapsto \widehat Q_n \left(\left( \lambda, \boldsymbol{\beta}^T\right)^T
\right)
= \frac{1}{n^2} \sum_{1\leq i,j \leq n}g(\boldsymbol U_i; \boldsymbol \theta,
\widehat{\boldsymbol \eta}_{\lambda,i})g(\boldsymbol U_j; \boldsymbol \theta,
\widehat{\boldsymbol \eta}_{\lambda,j}) \omega ( \boldsymbol W_i - \boldsymbol W_j),
$$
which is quadratic with an explicit unique minimum $\widehat {\boldsymbol \beta}(\lambda)$. Thus, we define a profile SmoothMD estimator of $\lambda_0$ as
$$
\widehat \lambda = \arg\min_{\lambda } \widehat Q_n\left(\left(\lambda,\widehat {\boldsymbol \beta}(\lambda)^T\right)^T \right),
$$
and, with at hand the estimate $\widehat{\lambda}$, we eventually calculate $\widehat{\boldsymbol \beta} (\widehat \lambda)$, the semiparametric SmoothMD estimate of $\boldsymbol \beta_0$. In a final step, the function $m(\cdot)$ in model \eqref{main_model} can be estimated from the residuals $\hat \epsilon_i=T(Y_i,\widehat{\lambda})-\boldsymbol X_i^T\widehat{\boldsymbol \beta}(\widehat{\lambda})$ by using any established smoothing procedure.
Details of the method are described in Section \ref{prelim}. Under mild regularity conditions it is then shown that our procedure leads to consistent estimators. Using standard kernel estimators based on second order kernels and suitable bandwidth sequences for nonparametric function estimation, we then establish $\sqrt{n}$-consistency and asymptotic normality, provided that the dimension of $\boldsymbol Z$ is $q<4$. Corresponding test procedures are described in Section 4. All theoretical results are derived uniformly for all possible choices in a compact set for the $d_1,\dots,d_{p+q}$ used in \eqref{MGA}. This provides theoretical grounds for a sample-based choice of $d_1,\dots,d_{p+q}$, such as the sample standard deviations of the components of the vectors $(\boldsymbol X_i^T, \boldsymbol Z_i^T)^T$.
We wish to note that the generality of our approach implies that the method may be used as a powerful tool to check parametric models. For example, in addition to verifying a log-transformation to the response variable, one may check linearity assumptions. The latter may be done by comparing the outputs of the parametric model with the results of a semiparametric analysis, where some of the regressors enter the model nonparametrically and define a corresponding vector $Z$. This is exemplified by our real data application in Section 4.
The remainder of the paper is organized as follows. In Section \ref{prelim}, we present our new estimation method and establish identification of the model parameters. Theoretical properties of the estimators are derived in Section \ref{con_asy_norm}, while in Section \ref{Test}, we investigate a distance-metric procedure for testing restrictions on parameters. In Section \ref{small_sample_study}, we study the finite sample behavior by a simulation study and apply the estimator to a real data sample. Our estimator performs well in our experiments and our tests yield accurate levels and good power in moderate samples. Finally, in Section \ref{discussion} we formulate few conclusions and discuss the extension of our approach to higher-dimension vectors $\boldsymbol Z$ using higher-order kernels, as well as efficiency aspects. The proofs are left to the Appendix.
\section{The semiparametric SmoothMD approach} \label{prelim}
In this section we formally define our semiparametric estimator. First, we investigate two issues. On the one hand, we prove identification of the true value $\boldsymbol \theta_0=(\lambda_0,\boldsymbol \beta_0^T)^T$ of the parameter of interest. Next, we discuss the recommendation appearing in the literature for normalizing the response variable. This issue is specific to the Box-Cox transformation, though similar problems occur with other families of transformations. Finally, we define our semiparametric SmoothMD estimator.
Before proceeding with this plan, let us slightly modify the definition of the conditional moment equation. The function $g(\boldsymbol U; \boldsymbol \theta, \boldsymbol \eta_\lambda)$ defined in \eqref{def1} has zero-mean if $\boldsymbol \theta=\boldsymbol \theta_0$, but there is no reason to expect its sample version to be centered, as is the case when $\lambda_0$ is given and one uses least squares; see \citet{li1996root}. We therefore propose to introduce an intercept and hereafter replace $g(\boldsymbol U; \boldsymbol \theta, \boldsymbol \eta_\lambda)$ with
\begin{equation}\label{def1_b}
g(\boldsymbol U; \boldsymbol \theta,
\gamma,
\boldsymbol \eta_\lambda) = \left(T(Y,\lambda) - E[T(Y,\lambda)\mid \boldsymbol Z]- (\boldsymbol X-E[\boldsymbol X\mid \boldsymbol Z])^T\boldsymbol \beta\right) f_z(\boldsymbol Z)
- \gamma.
\end{equation}
The true value of $\gamma$ is known to be $\gamma_0=0$, but this intercept slightly improves the results with finite samples, while it does not introduce any additional theoretical or computational complexity.
We use the following notation throughout the remainder of the paper. For $d_l, d_c\geq 1,$ let $\mathbb{R}^{d_l\times d_c}$ denote the set of $d_l\times d_c-$ matrices with real elements. Let $\boldsymbol{1}_{d_l}$ (resp. $\boldsymbol{0}_{d_l}$) denote the vector with all components equal to 1 (resp. 0), $\boldsymbol{0}_{d_l \times d_c}$ the $d_l\times d_c-$null matrix and $\boldsymbol{I}_{d_l \times d_c}$ the identity matrix with dimension $d_l \times d_c$. For a matrix $\boldsymbol A$, $\lVert \boldsymbol A\rVert$ is the Frobenius norm and $\lVert \boldsymbol A\rVert_{\rm{Sp}}$ the spectral norm. Below, $\boldsymbol D={\rm diag}(\boldsymbol d)$ is some positive definite diagonal matrix with $\boldsymbol d\in \mathcal{D} \subset \mathbb{R}^{p+q}_+$ being a diagonal vector with strictly positive components. Herein, $\mathcal{D}$ is a compact set and our asymptotic results are derived uniformly with respect to $\boldsymbol d\in \mathcal{D}$.
\subsection{Identification} \label{identification}
Let $-\infty <\lambda_{\rm min}< \lambda_0 < \lambda_{\rm max} <\infty$, with $\lambda_{\min}< 0$ and $\lambda_{\max}> 0$. For any $\lambda\in
[\lambda_{\rm min},\lambda_{\rm max}]$, let
\begin{equation}\label{beta_l}
(\gamma(\lambda), \boldsymbol \beta (\lambda)^T)^T = \arg\min_{\gamma\in\mathbb{R},\boldsymbol\beta\in \mathbb{R} ^p} E[g(\boldsymbol U_1; \boldsymbol \theta,\gamma, \boldsymbol \eta_{\lambda,1})g(\boldsymbol U_2; \boldsymbol \theta,\gamma, \boldsymbol \eta_{\lambda,2})\omega ( \boldsymbol W_1 - \boldsymbol W_2) ],
\end{equation}
with $ g(\boldsymbol U_1; \boldsymbol \theta,\gamma, \boldsymbol \eta_{\lambda,1}) $ and $ g(\boldsymbol U_2; \boldsymbol \theta, \gamma , \boldsymbol \eta_{\lambda,2})$ being independent copies of $g(\boldsymbol U; \boldsymbol \theta,\gamma, \boldsymbol \eta_{\lambda}) $ defined in equation \eqref{def1_b} with $\boldsymbol U=\left(Y,\boldsymbol W^T \right)^T$, $\boldsymbol W=\left(\boldsymbol X^T,\boldsymbol Z^T \right)^T $, $\boldsymbol \theta=(\lambda,\boldsymbol \beta^T)^T$ and $\boldsymbol \eta_{\lambda,k} =\boldsymbol \eta_{\lambda} (\boldsymbol Z_k)$, $k=1,2,$ where
\begin{equation*}\label{eq_eta}
\boldsymbol \eta_{\lambda} (\boldsymbol z)= (f_z(\boldsymbol z), E[T(Y,\lambda)\mid \boldsymbol Z= \boldsymbol z\;], E[\boldsymbol X\mid \boldsymbol Z=\boldsymbol z\;]^T)^T.
\end{equation*}
\begin{assumption}\emph{Data Generating Process}
\begin{enumerate}
\item The observations $\left(Y_i, \boldsymbol X_i^T, \boldsymbol Z_i^T \right)^T$ , $1 \leq i \leq n$, are i.i.d. copies of $\left(Y, \boldsymbol X^T,\boldsymbol Z^T \right)^T \in \mathbb{R} \times \mathbb{R}^p \times \mathbb{R}^q$. Moreover, there exists a constant $c>0$ such that $\mathbb{P}(Y > c)=1$.
\item The covariate vector $\boldsymbol Z$ admits a bounded density in $\mathbb{R}^q$. The covariate vector $\boldsymbol X$ is split into two subvectors $\boldsymbol X_c \in \mathbb{R}^{p_c}$ and $\boldsymbol X_d \in \mathbb{R}^{p_d} $ with $0\leq p_c, p_d\leq p$ and $p_c+p_d=p$. The subvector $\boldsymbol X_c$ admits a bounded density in $\mathbb{R}^{p_c}$. The subvector $\boldsymbol X_d$ takes values in a finite set.
\item The diagonal of the matrix $\boldsymbol D$ belongs to the $(p+q)-$dimension cube $\mathcal{D}= [d_L,d_U]^{p+q}$, with some fixed $0<d_L<d_U<\infty$.
\end{enumerate}
\label{ass_dgp}
\end{assumption}
The assumption that the discrete components of $\boldsymbol X$ take values in a finite set, is a technical condition that simplifies the proofs without significant restriction of the generality of the applications.
\begin{assumption}\emph{Identification}
\begin{enumerate}
\item $E\left[\| \boldsymbol X \|^2\right]<\infty$, $E\left[\| \boldsymbol Z \|^2\right]<\infty$, and $Var\left[\boldsymbol X - E[\boldsymbol X|\boldsymbol Z]\right]$ has full rank.
\item The true value $\boldsymbol\beta_{0,c}\in \mathbb{R}^{p_c}$ of the subvector of coefficients corresponding to $\boldsymbol X_c$ is not equal to $\boldsymbol 0_{p_c}$.
\item The continuous random subvector $\boldsymbol X_c$ is such that, for any $ \boldsymbol b \in \mathbb{R}^{p_c}$, $\boldsymbol b \neq \boldsymbol 0_{p_c}$, the variable
$\boldsymbol X^T_c \boldsymbol b $ is continuous with the support equal to the whole real line.
\item Whenever $\lambda\neq \lambda_0$, for any $\boldsymbol z$ in the support of $\boldsymbol Z$ and $\boldsymbol x_d$ in the support of the discrete subvector $\boldsymbol X_d$, the set of values of the map
$
\boldsymbol x_c \mapsto E\left[ T(Y,\lambda) - T(Y,\lambda_0)\mid \boldsymbol X_c=\boldsymbol x_c, \boldsymbol X_d=\boldsymbol x_d,
\boldsymbol Z=\boldsymbol z
\right] $,
$\boldsymbol x_c \in \mathbb{R}^{p_c},
$
is unbounded.
\item $E\left[Y^{2C_\lambda}\right]<\infty$, where $C_\lambda=\max(|\lambda_{min}|,\lambda_{max})<\infty$.
\end{enumerate}
\label{ass_ident}
\end{assumption}
Note that $Var\left[(\boldsymbol{X}^T, \boldsymbol{Z}^T)^T\right]$ necessarily has full rank, by Assumption \ref{ass_ident}.1 and the fact that $\boldsymbol Z$ admits a density.
The complete justification of this statement is given in the Appendix.
With all this in hand, we can now state the following identification result.
\begin{lem} Suppose that Assumptions \ref{ass_dgp} and \ref{ass_ident} hold true.
Let
$\gamma(\lambda)$ and $\boldsymbol \beta (\lambda)$ be defined as in
\eqref{beta_l}. Then, $\gamma(\lambda_0)=0$ and $\boldsymbol {\beta} (\lambda_0) = \boldsymbol\beta_0$ and
\begin{align*}
\mathbb{P} \left(E\left[ (T(Y,\lambda) - E[T(Y,\lambda)\mid \boldsymbol Z] ) f_z(\boldsymbol Z) - \gamma -
(\boldsymbol X - E[\boldsymbol X\mid\boldsymbol Z])^T\boldsymbol\beta f_z(\boldsymbol Z) \mid \boldsymbol X,\boldsymbol Z \right] = 0 \right) < 1,
\end{align*}
for all $\gamma\in \mathbb{R}$ and $\boldsymbol \theta = (\lambda,\boldsymbol \beta^T)^T \in [\lambda_{\rm min},\lambda_{\rm max}]\times \mathbb{R}^p$ such that $(\gamma,\boldsymbol \theta^T)^T\neq (0,\boldsymbol {\theta}_0^T)^T$. Moreover, for any $\varepsilon >0$,
\begin{multline}
\inf_{
|\lambda-\lambda_0| \geq \varepsilon} \; \inf_{\boldsymbol d \in\mathcal{D} } E\left[g\left(\boldsymbol U_1; (\lambda, \boldsymbol \beta (\lambda)^T)^T , \gamma (\lambda), \boldsymbol \eta_{\lambda,1}\right)g\left(\boldsymbol U_2; (\lambda, \boldsymbol \beta (\lambda)^T)^T , \gamma (\lambda), \boldsymbol \eta_{\lambda,2}\right)\right. \\ \times \left. \exp \left\{- (\boldsymbol W_1 - \boldsymbol W_2)^T \boldsymbol D (\boldsymbol W_1 - \boldsymbol W_2)\right\} \right] > 0.\label{well_sep_0}
\end{multline}
\label{lem_ident}
\end{lem}
\subsection{Box-Cox transformation and standardized responses} \label{stand_Y}
Let us note that
$$
\lim_{\lambda \uparrow \infty }\frac{y^\lambda -1}{\lambda} = 0 \quad \text{if $0< y<1$} \qquad and \qquad \lim_{\lambda \downarrow -\infty }\frac{y^\lambda -1}{\lambda} = 0 \quad \text{if $y>1$}.
$$
In classical estimation approaches for parametric regression models with Box-Cox transformed response, this is likely to induce instability for the estimation of the parameter $\lambda$. See, e.g., \citet{khazzoom1989note}, \citet{powell1996rescaled} and \citet{showalter1994monte} for a discussion of this well-known issue. In order to avoid such problems, the common recommendation is to standardize the response by some constant, say $s$, such that
$$
\mathbb{P}\left( Y/s < 1 \right) >0 \qquad \text{and} \qquad \mathbb{P}\left( Y/s > 1 \right) >0.
$$
The constant $s$ could be, for instance, the mean of $Y$ or the geometric mean of $Y$.\footnote{The geometric mean is defined as $G(Y)=\exp\{E(\log(Y)\}$ and the sample counterpart is $G_n = \prod\limits_{i=1}^{n} Y_i^{1/n}$.}
With finite samples, the practitioner would first estimate such a constant using the sample, and next would normalize the responses. The same type of problems might occur in our semiparametric extension of the Box-Cox transformation model. For this reason, we will replace our function $g(\boldsymbol U; \boldsymbol \theta,\gamma,\boldsymbol \eta_\lambda) $ by a family of functions $s^{-\lambda} g(\boldsymbol U; \boldsymbol \theta,\gamma, \boldsymbol \eta_\lambda)$ also indexed by $s$ that we shall let depend on the sample. This change of the family of functions is equivalent to changing $Y$ to $Y/s$ in the definition \eqref{def1_b}, and a rescaling of the parameters $\boldsymbol \beta$ and $\gamma$. By the profiling-based construction of our SmoothMD estimator, the replacement of the response $Y$ by $Y/s$ matters only for computing $\widehat \lambda$. Clearly, the identifiability property established in Lemma \ref{lem_ident} is preserved. In the remainder of the paper, we provide asymptotic results that are uniform with respect to $s$ in order to allow for a data-driven choice of $s$, such as for instance, the sample geometric mean of the response.
\subsection{The estimator} \label{estimator}
Given an independent sample $\left(Y_1,\boldsymbol X^T_1,\boldsymbol Z^T_1 \right)^T,\ldots, \left(Y_n,\boldsymbol X^T_n,\boldsymbol Z^T_n \right)^T$ from $\left(Y,\boldsymbol X^T,\boldsymbol Z^T \right)^T \in \mathbb{R} \times \mathbb{R}^{p+q}$, let us define
$$
\widehat{\mathbb{Y}}_n(\lambda) = \left( (T(Y_1,\lambda) - \widehat{E}[T(Y_1,\lambda)\mid \boldsymbol Z_1] ) \widehat{ f} _z(\boldsymbol Z_1),\ldots, (T(Y_n,\lambda) - \widehat{E} [T(Y_n,\lambda)\mid \boldsymbol Z_n]) \widehat{f}_z(\boldsymbol Z_n) \right)^T\in\mathbb{R}^n,
$$
and
$$
\widehat{\mathbb{X}}_n = \left( (\boldsymbol X_1-\widehat{E}[\boldsymbol X_1\mid \boldsymbol Z_1]) \widehat{f}_z(\boldsymbol Z_1),\ldots, (\boldsymbol X_n-\widehat{E}[\boldsymbol X_n\mid \boldsymbol Z_n]) \widehat{f}_z(\boldsymbol Z_n)\right)^T\in \mathbb{R}^{n\times p}.
$$
For $1\leq i \leq n$, $\widehat{\boldsymbol \eta}_{\lambda,i} = (\widehat{f}_z(\boldsymbol Z_i) , \widehat{E}[T(Y_i,\lambda)\mid \boldsymbol Z_i], \widehat{E}[\boldsymbol X_i\mid \boldsymbol Z_i]^T)^T$ are nonparametric kernel estimates of $ \boldsymbol \eta_{\lambda,i} =(f_z(\boldsymbol Z_i), E[T(Y_i,\lambda)\mid \boldsymbol Z_i], E[\boldsymbol X_i\mid\boldsymbol Z_i]^T)^T$. More precisely,
$$
\widehat{f}_z(\boldsymbol Z_i) = \frac{1}{nh^q}\sum_{j=1}^n K\left( \frac{\boldsymbol Z_i-\boldsymbol Z_j}{h} \right),\quad
\widehat{E}[T(Y_i,\lambda)\mid \boldsymbol Z_i]\widehat{f}_z(\boldsymbol Z_i)= \frac{1}{nh^q}\sum_{j=1}^nT(Y_j,\lambda) K\left( \frac{\boldsymbol Z_i-\boldsymbol Z_j}{h} \right),
$$
and
$$
\widehat{E}[\boldsymbol X_i\mid \boldsymbol Z_i]\widehat{f}_z(\boldsymbol Z_i)= \frac{1}{nh^q}\sum_{j=1}^n\boldsymbol X_j K\left( \frac{\boldsymbol Z_i-\boldsymbol Z_j}{h} \right).
$$
Here $K(\cdot)$ is a multivariate kernel function and $h$ is the bandwidth. Let $\boldsymbol \Omega_n$ be the $n\times n-$ symmetric matrix with elements
$$
\boldsymbol{\Omega}_{n,ij} = \exp \{ - (\boldsymbol X_i^T-\boldsymbol X_j^T, \boldsymbol Z_i^T-\boldsymbol Z_j^T)\boldsymbol D (\boldsymbol X_i-\boldsymbol X_j, \boldsymbol Z_i-\boldsymbol Z_j) \}, \qquad 1\leq i, j \leq n.
$$
Typically, the components of the vector $ \boldsymbol d$ defining the diagonal matrix $\boldsymbol D$, could be taken as proportional to the standard deviation of the components of the vectors $(\boldsymbol X_i^T, \boldsymbol Z_i^T)^T$. The definition of $\boldsymbol \Omega_{n,i,j}$ allows also to take into account, discrete components of $\boldsymbol X$. For finite support discrete covariates, one could set some large value for the corresponding diagonal element of $\boldsymbol D$, which in practice would be equivalent to an indicator of the event that the observations $i$ and $j$ have the same value for that covariate.
We can now define, for any $\lambda$, the estimates of $(\gamma(\lambda), \boldsymbol \beta (\lambda)^T)^T\in \mathbb{R}^{1+p}$ introduced in
\eqref{beta_l}. For any $s>0$, let
$$
\widehat Q_n\left(\left(\lambda, {\boldsymbol \beta}^T\right)^T, {\gamma} ;s \right) = n^{-2} s^{-2\lambda}\left( \widehat{\mathbb{Y}}_n(\lambda) - \gamma\boldsymbol{1}_n - \widehat{\mathbb{X}}_n \boldsymbol \beta \right)^T\boldsymbol \Omega_n \left( \widehat{\mathbb{Y}}_n(\lambda) - \gamma\boldsymbol{1}_n - \widehat{\mathbb{X}}_n \boldsymbol \beta\right).
$$
For fixed $s$ and $\lambda$, consider the generalized least-squares problem
\begin{equation}\label{gls}
\min_{\gamma,\boldsymbol \beta} \widehat Q_n\left(\left(\lambda, {\boldsymbol \beta}^T\right)^T, {\gamma} ;s \right).
\end{equation}
The solution of this problem does not depend on $ s^{-\lambda}$ and has the form of standard generalized least-squares estimators:
\begin{equation*}\label{gamma_l_hat}
\widehat \gamma(\lambda,{\boldsymbol \beta} (\lambda))
= \frac{1}{\boldsymbol{1}_n^T \boldsymbol\Omega_n \boldsymbol{1}_n} \boldsymbol{1}_n^T \boldsymbol\Omega_n \left(\widehat{\mathbb{Y}}_n(\lambda) - \widehat{\mathbb{X}}_n {\boldsymbol \beta} (\lambda) \right),
\end{equation*}
and
\begin{equation*}\label{beta_l_hat}
\widehat{\boldsymbol \beta} (\lambda) = \left(\widehat{\mathbb{X}}_n ^T{\mathbb{D}} _n\widehat{\mathbb{X}}_n \right)^{-1} \widehat{\mathbb{X}}_n ^T {\mathbb{D}} _n \widehat{\mathbb{Y}}_n(\lambda),
\end{equation*}
with
\begin{equation} \label{d_n}
{\mathbb{D}} _n = \boldsymbol\Omega_n - \frac{1}{\boldsymbol{1}_n^T \boldsymbol\Omega_n \boldsymbol{1}_n} \boldsymbol\Omega_n \boldsymbol{1}_n \boldsymbol{1}_n^T \boldsymbol\Omega_n \in\mathbb{R}^{n\times n}.
\end{equation}
Next, plugging $(\widehat{\gamma}(\lambda,\widehat{\boldsymbol \beta} (\lambda)), \widehat{\boldsymbol \beta} (\lambda)^T)^T$ into the problem \eqref{gls}, for a given $s$, we define the SmoothMD estimator of $\lambda_0$ as
\begin{equation}\label{lambda_l_hat}
\widehat \lambda = \widehat \lambda (s) = \arg\min_{\lambda\in\Lambda } s^{-\lambda} \widehat{\mathbb{Y}}_n(\lambda)^T \; \widehat {\mathbb{B}} _n \; s^{-\lambda} \widehat{\mathbb{Y}}_n(\lambda),
\end{equation}
with
$$
\widehat {\mathbb{B}}_n = {\mathbb{D}} _n - {\mathbb{D}} _n \widehat{\mathbb{X}}_n \left(\widehat{\mathbb{X}}_n ^T{\mathbb{D}} _n \widehat{\mathbb{X}}_n \right)^{-1} \widehat{\mathbb{X}}_n ^T {\mathbb{D}} _n \in\mathbb{R}^{n\times n}.
$$
Note that, by construction,
\begin{equation*}
{\mathbb{D}} _n \boldsymbol{1}_n = \widehat {\mathbb{B}} _n \boldsymbol{1}_n = \boldsymbol 0_n \quad \text{ and } \quad
\widehat {\mathbb{B}}_n \widehat{\mathbb{X}}_n = \boldsymbol {0}_{n\times p}.
\end{equation*}
Finally, the SmoothMD estimator of $\boldsymbol \beta_0$ is $\widehat{\boldsymbol \beta} (\widehat \lambda)$. We close this section by showing that our estimator is well-defined.
\begin{lem}\label{omega_mat}
If Assumptions \ref{ass_dgp}.3 and \ref{ass_ident} hold true, then, for each $n\geq 1$,
\begin{enumerate}
\item the matrices $\boldsymbol \Omega_n$ and $\widehat{\mathbb{X}}_n ^T\mathbb D_n \widehat{\mathbb{X}}_n $ are positive definite with probability 1. In particular, $\boldsymbol{1}_n^T \boldsymbol \Omega_n \boldsymbol{1}_n>0$ and $\widehat{\mathbb{X}}_n ^T\mathbb D_n \widehat{\mathbb{X}}_n $ is invertible with probability 1.
\item the matrix $\widehat {\mathbb{B}} _n $ is positive semi-definite with probability 1.
\end{enumerate}
\end{lem}
\begin{rmk}
The matrix ${\mathbb{D}}_n$ is defined in equation \eqref{d_n} and has dimension $n\times n$. Therefore, one might imagine that it becomes difficult to work with this matrix when the sample size is large. This is not the case because it is not necessary to estimate the matrix ${\mathbb{D}}_n$ itself. It suffices to compute $\widehat{\mathbb{X}}_n ^T{\mathbb{D}} _n$ and $\widehat{\mathbb{Y}}_n(\lambda)^T{\mathbb{D}} _n$ to be able to calculate $\widehat{\boldsymbol \beta} (\lambda)$ and $\widehat \lambda$. In chapter \ref{small_sample_study} we shall show that the estimator can be easily applied even for $n > 100,000$.
\end{rmk}
\section{Consistency and asymptotic normality} \label{con_asy_norm}
The estimator $\widehat \lambda(s)$ depends on $s$, a value that in practice could be calculated from the sample, such as for instance, the sample geometric mean. For this reason, our asymptotic results are stated uniformly with respect to $s$.
Our asymptotic results are also stated uniformly with respect to the diagonal of the matrix $\boldsymbol D$. This ensures that we can use a data driven estimate of $\boldsymbol D$ proportional to the empirical standard deviation of $\boldsymbol X$ or $\boldsymbol Z$.
Let's introduce some more notation: for each $\lambda \in\Lambda $, let
$$
{\mathbb{Y}}_n(\lambda) = \left( (T(Y_1,\lambda) - {E}[T(Y_1,\lambda)\mid \boldsymbol Z_1] ) { f} _z(\boldsymbol Z_1),\ldots, (T(Y_n,\lambda) - {E} [T(Y_n,\lambda)\mid \boldsymbol Z_n]) {f}_z(\boldsymbol Z_n) \right)^T\in\mathbb{R}^n,
$$
and
$$
{\mathbb{X}}_n = \left( (\boldsymbol X_1- {E}[\boldsymbol X_1\mid \boldsymbol Z_1]) {f}_z(\boldsymbol Z_1),\ldots, (\boldsymbol X_n- {E}[\boldsymbol X_n\mid \boldsymbol Z_n]) {f}_z(\boldsymbol Z_n)\right)^T\in \mathbb{R}^{n\times p}.
$$
Moreover,
$$
{\mathbb{B}} _n = {\mathbb{D}} _n - {\mathbb{D}} _n {\mathbb{X}}_n \left( {\mathbb{X}}_n ^T{\mathbb{D}} _n {\mathbb{X}}_n \right)^{-1} {\mathbb{X}}_n ^T {\mathbb{D}} _n \in\mathbb{R}^{n\times n},
$$
with $ {\mathbb{D}}_n $ defined in equation \eqref{d_n}. Again, by construction ${\mathbb{B}} _n \boldsymbol{1}_n = \boldsymbol 0_n$ and
${\mathbb{B}}_n {\mathbb{X}}_n = \boldsymbol {0}_{n\times p}$.
\begin{assumption}\label{ass_con} \emph{Consistency}
\begin{enumerate}
\item The functions $f_z(\cdot)$, $(mf_z)(\cdot)$, $E[\| \boldsymbol X \|^2\mid \boldsymbol Z=\cdot\;]f_z(\cdot)$ and $\sup_{\lambda \in\Lambda}(\partial^2 /\partial \lambda^2) E[T(Y, \lambda)\mid \boldsymbol Z=\cdot\;]f_z(\cdot)$ have H\"older continuous partial derivatives of order four.
\item The kernel $K(\cdot)$ is the product of $q$ univariate kernel functions $\widetilde K$ of bounded variation. Moreover, $\widetilde K$ is a symmetric function with integral equal to one and $\int t^2 \widetilde K(t) dt <\infty$.
\item The bandwidth $h$ belongs to a range $\mathcal{H}_{c,n}=[c_{min}n^{-\alpha}, c_{max}n^{-\alpha}],$ with $0 < \alpha < 1/q$ and $c_{min}$, $c_{max}$ positive constants.
\end{enumerate}
\end{assumption}
With all this in hand, we can now state the consistency of our estimator.
\begin{thm}[Consistency]\label{consist}
Assume that Assumptions \ref{ass_dgp}, \ref{ass_ident} and \ref{ass_con} hold true. Let $s_0$ be some normalizing value such that
$\mathbb{P}\left( Y/s_0 < 1 \right) >0$ and $\mathbb{P}\left( Y/s_0 > 1 \right) >0$ and let $S_n$ be an arbitrary $o_{\mathbb{P}}(1)$ neighborhood of $s_0$. Then
$$
\sup_{h\in\mathcal{H}_{c,n}} \sup_{s\in S_n} \sup_{\boldsymbol d \in \mathcal{D} }\left|\widehat \lambda - \lambda_0 \right|= o_{\mathbb{P}}(1) \quad \text{ and } \quad \sup_{h\in\mathcal{H}_{c,n}} \sup_{s\in S_n} \sup_{\boldsymbol d \in \mathcal{D}} \left\| \widehat{\boldsymbol{\beta}} (\widehat \lambda) - \boldsymbol{\beta}_0 \right\| = o_{\mathbb{P}}(1).
$$
\end{thm}
In Theorem \ref{consist} we require that $h\in\mathcal{H}_{c,n}$.
This implies that $nh^q \rightarrow \infty$ and $h \rightarrow 0$ for $n \rightarrow \infty$ which is in line with \citet{robinson1988root} and \citet{li1996root}.
Next, we prove asymptotic normality for our estimator. For this purpose, we first derive the asymptotic linear representation of $\widehat \lambda$ and and $\widehat {\boldsymbol \beta}(\widehat \lambda )$ from which the $\sqrt{n}-$asymptotic normality follows. In the following result, we show that $\widehat \lambda$ and $\widehat {\boldsymbol \beta}(\widehat \lambda )$ are asymptotically not equivalent to the infeasible estimators of $\lambda_0$ and $\boldsymbol \beta_0$, one would obtain when the infinite-dimensional parameter $\boldsymbol \eta_\lambda$ is given and the intercept $\gamma$ is equal to 0. This is in contrast to the results of \citet{li1996root} and \citet{robinson1988root}. The reason is that they can use the fact that $E[{\mathbb{X}}_{n,i}\mid \boldsymbol{Z}_i] = 0$ when controlling higher order terms. In our case, we weight the observations by $\boldsymbol{\Omega}_{n,ij}$ and, in general, $E[{\mathbb{X}}_{n,i}\boldsymbol{\Omega}_{n,ij}\mid \boldsymbol{Z}_i] \neq 0$ for $i \neq j$. This is also the reason why we need to ask for $q < 4$ instead of $q < 6$ as in \citet{li1996root}. Therefore, we require that $h\in\mathcal{H}_{sc,n}$, where $\mathcal{H}_{sc,n}=[c_{min}n^{-\alpha}, c_{max}n^{-\alpha}]$, with $\alpha \in (1/4,1/q)$. This is the small price we have to pay for getting what, to the best of our knowledge, is the first consistent estimation procedure for the semiparametric transformation model we investigate. However, as discussed in section \ref{discussion}, one could still apply our approach with $q\geq 4$, provided one uses higher-order kernels.
\begin{assumption} \label{ass_asy_norm} \emph{Asymptotic Normality}
\begin{enumerate}
\item $E\left[m(\boldsymbol{Z})^4\right] < \infty$ and $Var\left[\frac{\partial }{\partial \lambda} {T(Y,\lambda_0)}\right] > 0$.
\item $E\left[\varepsilon^4\right] < \infty$ and $E\left[\varepsilon^2\mid\boldsymbol{X},\boldsymbol{Z} \right] = \sigma^2(\boldsymbol{X},\boldsymbol{Z})$ is in $L^1 \cap L^2$.
\item The bandwidth $h$ belongs to $\mathcal{H}_{sc,n}=[c_{min}n^{-\alpha}, c_{max}n^{-\alpha}],$ with $\alpha \in (1/4, 1/q)$ and $c_{min}$, $c_{max}>0$.
\end{enumerate}
\end{assumption}
The results are again obtained uniformly with respect to the elements on the diagonal of the matrix $\boldsymbol D$ that determines $\boldsymbol{\Omega}_n$ and with respect to the scaling factor $s$ that could be used for numerical stability, as mentioned in Section \ref{stand_Y}. In addition, let $K_h (\cdot) = h^{-q}K(\cdot/h)$ and, for any $1\leq i, j \leq n$, let
$$
K_{h,ij}=K_h (\boldsymbol Z_i-\boldsymbol Z_j).
$$
\begin{prop}[Asymptotic representation]\label{AN_prop}
Assume that the conditions of Theorem \ref{consist} and Assumption \ref{ass_asy_norm} hold true. Then, uniformly with respect to $h\in\mathcal{H}_{sc,n}$, $\boldsymbol d\in \mathcal{D}$ and $s\in S_n$,
\begin{equation*}\label{lin_lambda}
\widehat \lambda - \lambda_0 = - \left[ \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0)^T \; {\mathbb{B}} _n \; \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) \right]^{-1} \; \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n( \lambda_0)^T \;{\mathbb{B}} _n \left[(\boldsymbol{\varepsilon f_z})_n - \left(\boldsymbol{\widehat{\varepsilon}_{|z}\widehat{f}_z}\right)_n \right]+ o_{\mathbb{P}}(n^{-1/2}) = O_{\mathbb{P}}(n^{-1/2}),
\end{equation*}
and
\begin{equation*}\label{beta_l_hat_0}
\widehat {\boldsymbol \beta}(\widehat \lambda ) - \boldsymbol {\beta}_0 =
\left( {\mathbb{X}}_n ^T {\mathbb{D}}_n {\mathbb{X}}_n \right)^{-1} {\mathbb{X}}_n ^T {\mathbb{D}}_n \left[ (\boldsymbol{\varepsilon f_z})_n -\left(\boldsymbol{\widehat{\varepsilon}_{|z}\widehat{f}_z}\right)_n
+ \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) \left( \widehat \lambda - \lambda_0 \right)\right] + o_{\mathbb{P}}(n^{-1/2})
= O_{\mathbb{P}}(n^{-1/2}),
\end{equation*}
where $(\boldsymbol{\varepsilon f_z})_n = (\varepsilon_1 {f}_z(\boldsymbol Z_1),\ldots,\varepsilon_n {f}_z(\boldsymbol Z_n))^T$
and $\left(\boldsymbol{\widehat{\varepsilon}_{|z}\widehat{f}_z}\right)_n = \left(\frac{1}{n}\sum\limits_{k=1, k\neq 1}^n \varepsilon_k K_{h,1k}, \ldots, \frac{1}{n}\sum\limits_{k=1, k\neq n}^n \varepsilon_k K_{h,nk} \right)^T$.
\end{prop}
Note that the asymptotic representation of $\widehat \lambda$ does not depend on $s_0$, i.e. the choice of $s_0$ does not influence the asymptotic behavior of $\widehat \lambda$. This result is in line with the result of \citet{powell1996rescaled}.
Next, we state the asymptotic normality of our estimator. We use the notation $\boldsymbol\Omega_{n,i,j}(\boldsymbol{d})= \boldsymbol\Omega_{n,i,j}$ and
\begin{equation*}
{\mathbb{D}}_n(\boldsymbol{d}) = \boldsymbol\Omega_n(\boldsymbol{d}) - \frac{1}{\boldsymbol{1}_n^T \boldsymbol\Omega_n(\boldsymbol{d}) \boldsymbol{1}_n} \boldsymbol\Omega_n(\boldsymbol{d}) \boldsymbol{1}_n \boldsymbol{1}_n^T \boldsymbol\Omega_n(\boldsymbol{d}),
\end{equation*}
to make the dependence of $\boldsymbol\Omega_{n}$ on $\boldsymbol{d}$ explicit. Note that with
\begin{align*}
\boldsymbol \Omega^X_{n,ij} (\boldsymbol d) = \boldsymbol \Omega^X_{n,ij} &= \exp \{ - (\boldsymbol X_i-\boldsymbol X_j)^T{\rm diag}(d_1,\ldots,d_p) (\boldsymbol X_i-\boldsymbol X_j) \} \quad\qquad \text{and}\\
\boldsymbol \Omega^Z_{n,ij}(\boldsymbol d) = \boldsymbol \Omega^Z_{n,ij} &= \exp \{ - (\boldsymbol Z_i-\boldsymbol Z_j)^T{\rm diag}(d_{p+1},\ldots,d_{p+q}) (\boldsymbol Z_i-\boldsymbol Z_j) \}, \quad 1\leq i, j \leq n,
\end{align*}
$\boldsymbol \Omega_{n,ij}(\boldsymbol d) = \boldsymbol \Omega^X_{n,ij}(\boldsymbol d) \boldsymbol \Omega^Z_{n,ij}(\boldsymbol d)$. Furthermore, we define, for $1 \leq i \leq n$,
\begin{align*}
\!\!\!\!\!\! \boldsymbol{\tau}_i(\boldsymbol{d})
:=
\Bigg(\!
\left(\!\frac{\partial}{\partial \lambda}\mathbb{Y}_{n,i} -\frac{1}{E\left[\boldsymbol{1}_n^T \boldsymbol\Omega_n(\boldsymbol{d}) \boldsymbol{1}_n\right]}E\left[\frac{\partial}{\partial \lambda}\mathbb{Y}_n ^T\boldsymbol\Omega_n(\boldsymbol{d}) \boldsymbol{1}_n \right]\right) ,
-\left(\!\mathbb{X}_{n,i}^T -\frac{1}{E\left[\boldsymbol{1}_n^T \boldsymbol\Omega_n(\boldsymbol{d}) \boldsymbol{1}_n\right]}E\left[ \boldsymbol{1}_n^T\boldsymbol\Omega_n(\boldsymbol{d}) {\mathbb{X}}_n \right]\right)
\!\Bigg)^{\!T}\!,
\end{align*}
where $\frac{\partial}{\partial \lambda}\mathbb{Y}_{n,i}(\lambda) = \left(\frac{\partial}{\partial \lambda}T(Y_i,\lambda) - {E}[\frac{\partial}{\partial \lambda}T(Y_i,\lambda)\mid \boldsymbol Z_i]\right) { f} _z(\boldsymbol Z_i)$ and $\mathbb{X}_{n,i} = (\boldsymbol X_i- {E}[\boldsymbol X_i\mid \boldsymbol Z_i]) {f}_z(\boldsymbol Z_i)$. In addition, let
\begin{align*}
\boldsymbol \Phi^X_{n,ij}(\boldsymbol d) = \boldsymbol \Omega^X_{n,ij}(\boldsymbol d) - E\left[\boldsymbol \Omega^X_{n,ik}(\boldsymbol d)\mid \boldsymbol X_i\right].
\end{align*}
Finally, for any vector $\boldsymbol a$, we denote $\boldsymbol a^{\otimes 2} = \boldsymbol a \boldsymbol a^T$. With all this in hand, we can state the following result.
\begin{thm}[Asymptotic normality]\label{AN}
Assume that the conditions of Proposition \ref{AN_prop} hold true. Then, uniformly with respect to $h\in\mathcal{H}_{sc,n}$, $\boldsymbol d
\in \mathcal{D}$ and $s\in S_n$,
\begin{align*}
\sqrt{n}\left((\widehat \lambda, \widehat {\boldsymbol \beta}(\widehat \lambda )^T)^T - (\lambda_0, \boldsymbol {\beta}_0 ^T)^T\right) + \boldsymbol{V}(\boldsymbol{d})^{-1}\left(\frac{1}{\sqrt{n}}\sum_{j=1}^{n}\varepsilon_j{f}_z(\boldsymbol Z_j)E\left[\boldsymbol{\tau}_i(\boldsymbol{d})\,\boldsymbol\Omega_{n,ij}^Z(\boldsymbol{d}) \boldsymbol \Phi^X_{n,ij}(\boldsymbol d)\mid \boldsymbol X_j, \boldsymbol Z_j\right]\right) =
o_{\mathbb{P} }\left(1\right),
\end{align*}
where
$$
\boldsymbol V(\boldsymbol{d}) = \underset{n \rightarrow \infty}{\lim}
\begin{pmatrix}
E\left[n^{-2}\frac{\partial}{\partial \lambda}\mathbb{Y}_n(\lambda_0)^T{\mathbb{D}}_n(\boldsymbol{d})\frac{\partial}{\partial \lambda}\mathbb{Y}_n(\lambda_0) \right]&
- E\left[n^{-2}\frac{\partial}{\partial \lambda}\mathbb{Y}_n(\lambda_0)^T{\mathbb{D}}_n(\boldsymbol{d}) \mathbb{X}_n\right]\\
- E\left[n^{-2}\mathbb{X}_n^T{\mathbb{D}}_n(\boldsymbol{d})\frac{\partial}{\partial \lambda}\mathbb{Y}_n(\lambda_0)\right] &E\left[n^{-2}\mathbb{X}_n^T{\mathbb{D}}_n(\boldsymbol{d})\mathbb{X}_n\right]
\end{pmatrix}.
$$
As a consequence, $\sqrt{n}\left((\widehat \lambda, \widehat {\boldsymbol \beta}(\widehat \lambda )^T)^T - (\lambda_0, \boldsymbol {\beta}_0 ^T)^T\right)$ converges in distribution to a $(p+1)-$dimension centered Gaussian vector with variance $\boldsymbol V(\boldsymbol{d})^{-1}\boldsymbol\Delta(\boldsymbol{d})\boldsymbol V(\boldsymbol{d})^{-1}$
where
$$
\boldsymbol\Delta(\boldsymbol{d}) = E\left\{Var\left[\varepsilon_j\mid \boldsymbol X_j, \boldsymbol Z_j\right]{f}^2_z(\boldsymbol Z_j)
\left( E \left[\boldsymbol{\tau}_i(\boldsymbol{d})\,\boldsymbol\Omega_{n,ij}^Z(\boldsymbol{d}) \boldsymbol \Phi^X_{n,ij}(\boldsymbol d)\mid \boldsymbol X_j, \boldsymbol Z_j\right]\right)^{\otimes 2}
\right\}.
$$
\end{thm}
If $\boldsymbol{\eta}_\lambda$ were known, which corresponds to the case studied by \citet{lavergne2013smooth}, $\boldsymbol \Phi^X_{n,ij}(\boldsymbol d)$ should be replaced by $\boldsymbol \Omega^X_{n,ij}(\boldsymbol d)$ in the expression of $\boldsymbol\Delta(\boldsymbol{d})$.
We can estimate the covariance matrix by $\widehat{\boldsymbol V}(\boldsymbol{d})^{-1} \widehat{\boldsymbol\Delta}(\boldsymbol{d}) \widehat{\boldsymbol V}(\boldsymbol{d})^{-1}$, where
\begin{equation}\label{variance_est}
\begin{aligned}
\widehat{\boldsymbol V}(\boldsymbol{d}) &=
\begin{pmatrix}
n^{-2}\frac{\partial}{\partial \lambda}\widehat{\mathbb{Y}}_n(\widehat \lambda)^T{\mathbb{D}}_n(\boldsymbol{d})\frac{\partial}{\partial \lambda}\widehat{\mathbb{Y}}_n(\widehat\lambda) & - n^{-2}\frac{\partial}{\partial \lambda}\widehat{\mathbb{Y}}_n(\widehat\lambda)^T{\mathbb{D}}_n(\boldsymbol{d}) \widehat{\mathbb{X}}_n\\
- n^{-2}\widehat{\mathbb{X}}_n^T{\mathbb{D}}_n(\boldsymbol{d})\frac{\partial}{\partial \lambda}\widehat{\mathbb{Y}}_n(\widehat\lambda) &n^{-2} \widehat{\mathbb{X}}_n^T{\mathbb{D}}_n(\boldsymbol{d})\widehat{\mathbb{X}}_n
\end{pmatrix}
\\ \text{and} \quad \\
\widehat{\boldsymbol\Delta}(\boldsymbol{d})
& = n^{-3}\left(\frac{\partial}{\partial \lambda}\widehat{\mathbb{Y}}_n(\widehat \lambda), -\widehat{\mathbb{X}}_n\right)^T{\mathbb{D}}_{n,inf}(\boldsymbol{d})\widehat{\boldsymbol \Phi}_{n}(\boldsymbol{d})\widehat{\boldsymbol{\Sigma}}_n \widehat{\boldsymbol \Phi}^T_{n}(\boldsymbol{d}){\mathbb{D}}^T_{n,inf}(\boldsymbol{d})\left(\frac{\partial}{\partial \lambda}\widehat{\mathbb{Y}}_n(\widehat \lambda), -\widehat{\mathbb{X}}_n\right).
\end{aligned}
\end{equation}
Here, $\widehat{\boldsymbol \Phi}^X_{n}$ and $\widehat{\boldsymbol \Phi}_{n}$ are the $n\times n-$ symmetric matrices with elements
\begin{align*}
\widehat{\boldsymbol \Phi}^X_{n,ij}(\boldsymbol d) &= \boldsymbol \Omega^X_{n,ij}(\boldsymbol d) - \frac{1}{n}\sum\limits_{k=1}^{n} \boldsymbol \Omega^X_{n,ik}(\boldsymbol d),
\quad 1\leq i, j \leq n \\
\widehat{\boldsymbol \Phi}_{n,ij}(\boldsymbol d) &= \widehat{\boldsymbol \Phi}^X_{n,ij}(\boldsymbol d)\boldsymbol \Omega^Z_{n,ij}(\boldsymbol d),
\hskip 1.74cm 1\leq i, j \leq n \\
and \quad
{\mathbb{D}}_{n,inf}(\boldsymbol{d}) &= \boldsymbol{I}_{n\times n} - \frac{1}{\boldsymbol{1}_n^T \boldsymbol\Omega_n(\boldsymbol{d}) \boldsymbol{1}_n} \boldsymbol\Omega_n(\boldsymbol{d}) \boldsymbol{1}_n \boldsymbol{1}_n^T.
\end{align*}
$\widehat{\boldsymbol{\Sigma}}_n = $ \rm{diag}$\left(\widehat{Var}\left[\varepsilon_1 f_z(\boldsymbol Z_1)\mid \boldsymbol X_1, \boldsymbol Z_1\right], \ldots, \widehat{Var}\left[\varepsilon_n{f}_z(\boldsymbol Z_n)\mid \boldsymbol X_n, \boldsymbol Z_n\right]\right)$ is an estimator of
$\rm{diag}\Big(Var\left[\varepsilon_1 f_z(\boldsymbol Z_1)\mid \boldsymbol X_1, \boldsymbol Z_1\right],\\ \ldots, Var\left[\varepsilon_n{f}_z(\boldsymbol Z_n)\mid \boldsymbol X_n, \boldsymbol Z_n\right]\Big)$.
One can use a nonparametric estimator for the conditional variance or alternatively, use an estimate of the error terms to approximate the conditional variance in the spirit of the Eiker-White variance estimator. Consistency of the above estimators is straightforward to establish.
\section{Testing based on SmoothMD for parameter restrictions} \label{Test}
In Section \ref{con_asy_norm} we established consistency and asymptotic normality of our estimator.
The asymptotic behavior of our estimator is not influenced by the standardization with $s$ but the asymptotic variance is affected by the estimation of $\boldsymbol{\eta}_\lambda$.
In addition, the behavior of our estimator is, even asymptotically, influenced by the vectors $\boldsymbol{d}_1$ and $\boldsymbol{d}_2$. When developing a test theory, we should
take that influence into account in order to get reliable results. That's what we do in the following.
\subsection{Testing the transformation parameter}
When it comes to testing parameter restrictions in the semiparametric partially linear regression model with Box-Cox transformation we might be mainly interested in testing if $\lambda$ is zero or not and if the components of $\boldsymbol{\beta}$ are zero. However, we shall consider here a more general approach to allow for more complex hypotheses as well. We separate the discussion into two parts. In the first part, we consider only restrictions for $\lambda$ and in the second part we consider restrictions for $\boldsymbol\beta$ with and without restricting $\lambda$.
Suppose we want to test the restriction for $\lambda$ given by
\begin{align}\label{test_lambda}
H_0: \lambda_0 = \lambda_R.
\end{align}
In order to test this restriction, we can use the distance metric statistic proposed by \citet{lavergne2013smooth}. Adapted to our case and for testing \eqref{test_lambda}, we consider the distance
\begin{align*}
DM_{\lambda} =
\frac{1}{n}\widehat{\mathbb{Y}}_n(\lambda_R)^T \; \widehat {\mathbb{B}} _n \widehat{\mathbb{Y}}_n(\lambda_R)
-
\frac{1}{n}\widehat{\mathbb{Y}}_n(\widehat\lambda)^T \; \widehat {\mathbb{B}} _n \widehat{\mathbb{Y}}_n(\widehat\lambda).
\end{align*}
The distance metric is based on the object that is minimized to get the estimate for $\lambda$, see equation \eqref{lambda_l_hat}. However, the test statistic is not standardized by $s^{-\lambda}$ as we only need this for the estimation of $\lambda$.
Let
\begin{align*}
\boldsymbol{A}_n =
\begin{pmatrix}
\frac{\partial}{\partial \lambda}\mathbb{Y}_n(\lambda_0)^T\\
- \mathbb{X}_n^T
\end{pmatrix}
{\mathbb{D}}_n\left( (\boldsymbol{\varepsilon f_z})_n - \left(\boldsymbol{\widehat{\varepsilon}_{|z}\widehat{f}_z}\right)_n\right).
\end{align*}
We can therefore now state the following Proposition.
\begin{prop}\label{prop_test_lambda}
Assume that the conditions of Proposition \ref{AN_prop} hold true. Then, uniformly with respect to $h\in\mathcal{H}_{sc,n}$, $\boldsymbol d\in \mathcal{D}$ and $s\in S_n$,
$$
DM_{\lambda} -
(1, \boldsymbol{0}_p^T)\boldsymbol V(\boldsymbol{d})^{-1} n^{-3/2}\boldsymbol{A}_n n^{-3/2}\boldsymbol{A}_n^T \boldsymbol V(\boldsymbol{d})^{-1}(1, \boldsymbol{0}_p^T)^T
E\left[ \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) ^T {\mathbb{B}} _n \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) \right]
= o_{\mathbb{P}}(1),
$$
under $H_0$ and $\mathbb{P}(n^{-1}DM_{\lambda} > c) \rightarrow 1$ for any $c>0$ if $H_0$ does not hold.
\end{prop}
The process $(1, \boldsymbol{0}_p^T)\boldsymbol V(\boldsymbol{d})^{-1} n^{-3/2}\boldsymbol{A}_n n^{-3/2}\boldsymbol{A}_n^T \boldsymbol V(\boldsymbol{d})^{-1}(1, \boldsymbol{0}_p^T)^T E\left[ \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) ^T {\mathbb{B}} _n \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) \right]$ is asymptotically tight and for each $\boldsymbol{d}$ behaves asymptotically as a chi-square times \\
$(1, \boldsymbol{0}_p^T)\boldsymbol V(\boldsymbol{d})^{-1} \boldsymbol\Delta(\boldsymbol{d})
\boldsymbol V(\boldsymbol{d})^{-1}(1, \boldsymbol{0}_p^T)^T E\left[ \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) ^T {\mathbb{B}} _n \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) \right]$, see \citet{johnson1995}. The distribution of the distance metric statistic therefore is, in general, non-pivotal. Determining critical values requires the estimation of
$(1, \boldsymbol{0}_p^T)\boldsymbol V(\boldsymbol{d})^{-1}\boldsymbol \Delta(\boldsymbol{d})
\boldsymbol V(\boldsymbol{d})^{-1}(1, \boldsymbol{0}_p^T)^T E\left[ \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) ^T {\mathbb{B}} _n \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) \right]$, which can rely on the estimators stated in \eqref{variance_est}.
\subsection{Testing the slope coefficients}
In the next part, we consider restrictions for $\boldsymbol\beta$. Suppose we want to test $r$ linear restrictions for $\boldsymbol\beta$ given by
\begin{align}\label{test_beta}
H_0: \boldsymbol{R}\boldsymbol\beta_0 = \boldsymbol c,
\end{align}
where $\boldsymbol{R}$ is a $r \times p-$ matrix of full rank and $\boldsymbol c \in \mathbb{R}^r$. In order to test the restrictions, we need to find the restricted estimators for $\boldsymbol{\beta}_0$, $\widehat{\boldsymbol{\beta}}_R(\lambda)$, and $\lambda_0$, $\widehat \lambda_R$. We minimize
\begin{align*}
n^{-2} s^{-2\lambda}\left( \widehat{\mathbb{Y}}_n(\lambda) - \widehat{\mathbb{X}}_n \boldsymbol \beta \right)^T\boldsymbol
{\mathbb{D}}_n \left( \widehat{\mathbb{Y}}_n(\lambda) - \widehat{\mathbb{X}}_n \boldsymbol \beta\right)
\quad s.t. \quad
\boldsymbol{R}\boldsymbol{\beta} = \boldsymbol c,
\end{align*}
with respect to $\boldsymbol \beta$ and get that
\begin{align*}
\widehat{\boldsymbol{\beta}}_R(\lambda) = \widehat{\boldsymbol{\beta}}(\lambda) - \left(\widehat{\mathbb{X}}_n^T{\mathbb{D}}_n\widehat{\mathbb{X}}_n\right)^{-1} \boldsymbol{R}^T
\left(\boldsymbol{R} \left(\widehat{\mathbb{X}}_n^T{\mathbb{D}}_n\widehat{\mathbb{X}}_n\right)^{-1} \boldsymbol{R}^T \right)^{-1}
\left(\boldsymbol{R}\widehat{\boldsymbol{\beta}}(\lambda) - \boldsymbol c\right).
\end{align*}
The restricted estimator for $\lambda_0$ is then given by
\begin{equation}\label{lambda_l_hat_R}
\widehat \lambda_R = \widehat \lambda_R (s) = \arg\min_{\lambda\in\Lambda }
s^{-\lambda}\left( \widehat{\mathbb{Y}}_n(\lambda) - \widehat{\mathbb{X}}_n \widehat{\boldsymbol{\beta}}_R(\lambda)\right)^T {\mathbb{D}}_n \,
s^{-\lambda}\left( \widehat{\mathbb{Y}}_n(\lambda) - \widehat{\mathbb{X}}_n \widehat{\boldsymbol{\beta}}_R(\lambda)\right).
\end{equation}
With all the estimators in hand, we can now define our distance metric statistic for testing \eqref{test_beta}.
\begin{align*}
DM_{\boldsymbol\beta} = \frac{1}{n}
\left( \widehat{\mathbb{Y}}_n(\widehat \lambda_R) - \widehat{\mathbb{X}}_n \widehat{\boldsymbol{\beta}}_R(\widehat \lambda_R)\right)^T
{\mathbb{D}}_n \,
\left( \widehat{\mathbb{Y}}_n(\widehat \lambda_R) - \widehat{\mathbb{X}}_n \widehat{\boldsymbol{\beta}}_R(\widehat \lambda_R)\right)
-
\frac{1}{n}\widehat{\mathbb{Y}}_n(\widehat\lambda)^T \; \widehat {\mathbb{B}} _n \widehat{\mathbb{Y}}_n(\widehat\lambda).
\end{align*}
The distance metric is based on the object that is minimized to get the restricted estimate for $\lambda$, see equation \eqref{lambda_l_hat_R}. Once again the test statistic is not standardized by $s^{-\lambda}$ as we only need this for the estimation of $\lambda$.
Let,
$$
{\mathbb{B}}_{n,R} = {\mathbb{B}}_n +
{\mathbb{D}}_n {\mathbb{X}}_n\left({\mathbb{X}}_n^T{\mathbb{D}}_n{\mathbb{X}}_n\right)^{-1} \boldsymbol{R}^T
\left(\boldsymbol{R} \left({\mathbb{X}}_n^T{\mathbb{D}}_n{\mathbb{X}}_n\right)^{-1} \boldsymbol{R}^T \right)^{-1}
\boldsymbol{R} \left({\mathbb{X}}_n^T{\mathbb{D}}_n{\mathbb{X}}_n\right)^{-1}{\mathbb{X}}_n^T{\mathbb{D}}_n,
$$
and
$$
\boldsymbol{V}_R(\boldsymbol{d}) =
\left(
E\left[\frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0)^T{\mathbb{B}}_{n,R}\frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0)\right]^{-1},
E\left[\frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0)^T{\mathbb{B}}_{n,R}\frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0)\right]^{-1}
E\left[\frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0)^T {\mathbb{D}}_n{\mathbb{X}}_n{\mathbb{B}}_n^+\right]
\right),
$$
where
$$
{\mathbb{B}}_n^+ = \left({\mathbb{X}}_n^T{\mathbb{D}}_n{\mathbb{X}}_n\right)^{-1}
- \left({\mathbb{X}}_n^T{\mathbb{D}}_n{\mathbb{X}}_n\right)^{-1} \boldsymbol{R}^T
\left(\boldsymbol{R} \left({\mathbb{X}}_n^T{\mathbb{D}}_n{\mathbb{X}}_n\right)^{-1} \boldsymbol{R}^T \right)^{-1}
\boldsymbol{R} \left({\mathbb{X}}_n^T{\mathbb{D}}_n{\mathbb{X}}_n\right)^{-1}.
$$
We can therefore now state the following proposition.
\begin{prop} \label{prop_test_beta}
Assume that the conditions of Proposition \ref{AN_prop} hold true. Then, uniformly with respect to $h\in\mathcal{H}_{sc,n}$, $\boldsymbol d\in \mathcal{D}$ and $s\in S_n$,
\begin{align*}
DM_{\boldsymbol\beta}&\\
-
&n^{-3/2}\boldsymbol{A}_n^T\Bigg(
\left(\boldsymbol{0}_{p\times 1}, \boldsymbol{I}_{p\times p}\right)^T E\left[{\mathbb{X}}_n^T{\mathbb{D}}_n{\mathbb{X}}_n\right]^{-1} \boldsymbol{R}^T
\left(\boldsymbol{R} E\left[{\mathbb{X}}_n^T{\mathbb{D}}_n{\mathbb{X}}_n\right]^{-1} \boldsymbol{R}^T \right)^{-1}
\boldsymbol{R} E\left[{\mathbb{X}}_n^T{\mathbb{D}}_n{\mathbb{X}}_n\right]^{-1}\left(\boldsymbol{0}_{p\times 1}, \boldsymbol{I}_{p\times p}\right)\\
&\quad -
\boldsymbol{V}_R(\boldsymbol{d})^T\boldsymbol{V}_R(\boldsymbol{d})
E\left[\frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) ^T {\mathbb{B}}_{n,R} \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) \right] \\
&\quad + \boldsymbol V(\boldsymbol{d})^{-1}(1, \boldsymbol{0}_p^T)^T (1, \boldsymbol{0}_p^T)\boldsymbol V(\boldsymbol{d})^{-1}
E\left[ \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) ^T {\mathbb{B}} _n \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) \right]
\Bigg)\boldsymbol{A}_nn^{-3/2} = o_{\mathbb{P}}(1),
\end{align*}
under $H_0$ and $\mathbb{P}(n^{-1}DM_{\boldsymbol\beta} > c) \rightarrow 1$ for any $c>0$ if $H_0$ does not hold.
\end{prop}
The process in Proposition \ref{prop_test_beta} is asymptotically tight and for each $\boldsymbol{d}$ behaves asymptotically as a weighted sum of $p +1 - r$ independent chi-squares, where the weights are the positive eigenvalues of
\begin{equation*}\label{dist_beta}
\begin{aligned}
&\left(\boldsymbol{0}_{p\times 1}, \boldsymbol{I}_{p\times p}\right)^T E\left[{\mathbb{X}}_n^T{\mathbb{D}}_n{\mathbb{X}}_n\right]^{-1} \boldsymbol{R}^T
\left(\boldsymbol{R} E\left[{\mathbb{X}}_n^T{\mathbb{D}}_n{\mathbb{X}}_n\right]^{-1} \boldsymbol{R}^T \right)^{-1}
\boldsymbol{R} E\left[{\mathbb{X}}_n^T{\mathbb{D}}_n{\mathbb{X}}_n\right]^{-1}\left(\boldsymbol{0}_{p\times 1}, \boldsymbol{I}_{p\times p}\right)
\boldsymbol\Delta(\boldsymbol{d})
\\
&\quad -
\boldsymbol{V}_R(\boldsymbol{d})^T\boldsymbol{V}_R(\boldsymbol{d})\boldsymbol\Delta(\boldsymbol{d}, \boldsymbol{d})
E\left[\frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) ^T {\mathbb{B}}_{n,R} \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) \right] \\
&\quad + \boldsymbol V(\boldsymbol{d})^{-1}(1, \boldsymbol{0}_p^T)^T (1, \boldsymbol{0}_p^T)\boldsymbol V(\boldsymbol{d})^{-1}
\boldsymbol\Delta(\boldsymbol{d})
E\left[ \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) ^T {\mathbb{B}} _n \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) \right],
\end{aligned}
\end{equation*}
see \citet{johnson1995}. Determining critical values requires the estimation of the last expression. We can use the estimators stated in \eqref{variance_est} and for all other components, we simply replace the unknown expressions by their sample version, e.g. estimate ${\mathbb{B}}_{n,R}$ by
$$
\widehat{\mathbb{B}}_n + {\mathbb{D}}_n \widehat{\mathbb{X}}_n\left(\widehat{\mathbb{X}}_n^T{\mathbb{D}}_n\widehat{\mathbb{X}}_n\right)^{-1} \boldsymbol{R}^T
\left(\boldsymbol{R} \left(\widehat{\mathbb{X}}_n^T{\mathbb{D}}_n\widehat{\mathbb{X}}_n\right)^{-1} \boldsymbol{R}^T \right)^{-1}
\boldsymbol{R} \left(\widehat{\mathbb{X}}_n^T{\mathbb{D}}_n\widehat{\mathbb{X}}_n\right)^{-1}\widehat{\mathbb{X}}_n^T{\mathbb{D}}_n.
$$
\subsection{Testing the transformation parameter and the slope coefficients}
Finally, we consider the combined restrictions for $\boldsymbol\beta$ and $\lambda$. Suppose we want to test
\begin{align*}
H_0: \boldsymbol{R}\boldsymbol\beta_0 = \boldsymbol c \quad and \quad \lambda_0 = \lambda_R.
\end{align*}
In contrast to the hypothesis stated in \eqref{test_beta} we do not need to estimate $\widehat{\lambda}_R$. Therefore, the distance metric statistic is for this case given by
\begin{align*}
DM_{\boldsymbol\beta,\lambda} = \frac{1}{n}
\left( \widehat{\mathbb{Y}}_n(\lambda_R) - \widehat{\mathbb{X}}_n \widehat{\boldsymbol{\beta}}_R(\lambda_R)\right)^T
{\mathbb{D}}_n \,
\left( \widehat{\mathbb{Y}}_n(\lambda_R) - \widehat{\mathbb{X}}_n \widehat{\boldsymbol{\beta}}_R(\lambda_R)\right)
-
\frac{1}{n}\widehat{\mathbb{Y}}_n(\widehat\lambda)^T \; \widehat {\mathbb{B}} _n \widehat{\mathbb{Y}}_n(\widehat\lambda).
\end{align*}
In addition, $DM_{\boldsymbol\beta,\lambda}$ does not converge to the same expression as $DM_{\boldsymbol\beta}$ as $\lambda_R$ is fixed. Therefore, we state the following proposition.
\begin{prop} \label{prop_test_beta_lambda}
Assume that the conditions of Proposition \ref{AN_prop} hold true. Then, uniformly with respect to $h\in\mathcal{H}_{sc,n}$, $\boldsymbol d\in \mathcal{D}$ and $s\in S_n$,
\begin{align*}
DM_{\boldsymbol\beta, \lambda}&\\
-
&n^{-3/2}\boldsymbol{A}_n^T\Bigg(
\left(\boldsymbol{0}_{p\times 1}, \boldsymbol{I}_{p\times p}\right)^T E\left[{\mathbb{X}}_n^T{\mathbb{D}}_n{\mathbb{X}}_n\right]^{-1} \boldsymbol{R}^T
\left(\boldsymbol{R} E\left[{\mathbb{X}}_n^T{\mathbb{D}}_n{\mathbb{X}}_n\right]^{-1} \boldsymbol{R}^T \right)^{-1}
\boldsymbol{R} E\left[{\mathbb{X}}_n^T{\mathbb{D}}_n{\mathbb{X}}_n\right]^{-1}\left(\boldsymbol{0}_{p\times 1}, \boldsymbol{I}_{p\times p}\right)\\
&\quad + \boldsymbol V(\boldsymbol{d})^{-1}(1, \boldsymbol{0}_p^T)^T (1, \boldsymbol{0}_p^T)\boldsymbol V(\boldsymbol{d})^{-1}
E\left[ \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) ^T {\mathbb{B}} _n \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) \right]
\Bigg)\boldsymbol{A}_nn^{-3/2} = o_{\mathbb{P}}(1),
\end{align*}
under $H_0$ and $\mathbb{P}(n^{-1}DM_{\boldsymbol\beta, \lambda} > c) \rightarrow 1$ for any $c>0$ if $H_0$ does not hold.
\end{prop}
The process in Proposition \ref{prop_test_beta_lambda} is asymptotically tight and for each $\boldsymbol{d}$ behaves asymptotically as a weighted sum of $p - r$ independent chi-squares, where the weights are the positive eigenvalues of
\begin{equation*}\label{dist_beta_lambda}
\begin{aligned}
&\left(\boldsymbol{0}_{p\times 1}, \boldsymbol{I}_{p\times p}\right)^T E\left[{\mathbb{X}}_n^T{\mathbb{D}}_n{\mathbb{X}}_n\right]^{-1} \boldsymbol{R}^T
\left(\boldsymbol{R} E\left[{\mathbb{X}}_n^T{\mathbb{D}}_n{\mathbb{X}}_n\right]^{-1} \boldsymbol{R}^T \right)^{-1}
\boldsymbol{R} E\left[{\mathbb{X}}_n^T{\mathbb{D}}_n{\mathbb{X}}_n\right]^{-1}\left(\boldsymbol{0}_{p\times 1}, \boldsymbol{I}_{p\times p}\right)
\boldsymbol\Delta(\boldsymbol{d})
\\
&\quad + \boldsymbol V(\boldsymbol{d})^{-1}(1, \boldsymbol{0}_p^T)^T (1, \boldsymbol{0}_p^T)\boldsymbol V(\boldsymbol{d})^{-1}
\boldsymbol\Delta(\boldsymbol{d})
E\left[ \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) ^T {\mathbb{B}} _n \frac{\partial }{\partial \lambda} {\mathbb{Y}}_n(\lambda_0) \right],
\end{aligned}
\end{equation*}
see \citet{johnson1995}. Determining critical values requires the estimation of the last expression. We can use the estimators stated in \eqref{variance_est} and for all other components, we simply replace the unknown expressions by their sample version.
\begin{rmk}
The Propositions of Section \ref{Test} are also valid if the unknown parameters $\lambda$ and $\boldsymbol{\beta}$ are estimated without the intercept nuisance parameter $\gamma$. In that case ${\mathbb{D}}_{n}$ is replaced by $\boldsymbol \Omega_{n}$ in the statements. Moreover, when estimating the unknown variance ${\mathbb{D}}_{n,inf}$ has to be replaced by $\boldsymbol{I}_{n\times n}$.
\end{rmk}
\section{Small sample study and real data application} \label{small_sample_study}
In this section we consider the small sample behavior of our estimator. We conduct several simulation experiments to consider bias and standard deviation for the estimated parameters. In addition, we conduct hypothesis tests as discussed in Section \ref{Test}.
In the proof of Theorem \ref{AN} it was established that
\begin{align*}
\left((\widehat \lambda, \widehat {\boldsymbol \beta}(\widehat \lambda )^T)^T - (\lambda_0, \boldsymbol {\beta}_0 ^T)^T\right) &= - \boldsymbol{V}(\boldsymbol{d})^{-1}\Bigg(\frac{1}{n}\sum_{j=1}^{n}\varepsilon_j{f}_z(\boldsymbol Z_j) E\left[\boldsymbol{\tau}_i(\boldsymbol{d})\,\boldsymbol\Omega_{n,ij}(\boldsymbol{d})\mid \boldsymbol X_j, \boldsymbol Z_j\right]\\
&\quad -
\frac{1}{n} \sum_{k=1}^{n}\varepsilon_k f_z(\boldsymbol Z_k) E\left[\boldsymbol{\tau}_i(\boldsymbol{d}) \boldsymbol\Omega_{n,ik}^Z \boldsymbol\Omega_{n,ij}^X\mid \boldsymbol Z_k\right]
\Bigg) +
o_{\mathbb{P} }\left(n^{-1/2}\right).
\end{align*}
The second sum in this asymptotic representation is due to the estimation of $\boldsymbol \eta_\lambda$. In order to propose a simpler procedure, in our simulation experiments we also investigated what happens when one neglects the second part in the asymptotic representation, that is abusively consider
$
E\left[\boldsymbol{\tau}_i(\boldsymbol{d}) \boldsymbol\Omega_{n,ik}^Z \boldsymbol\Omega_{n,ij}^X\mid \boldsymbol Z_k\right]= \boldsymbol 0_{p+1} $. The estimator is labeled SmoothMD* in this section. The reason for this investigation is that $E\left[\boldsymbol{\tau}_i(\boldsymbol{d}) \boldsymbol\Omega_{n,ij}^Z \boldsymbol\Omega_{n,ij}^X\mid \boldsymbol Z_k\right]$ is indeed null. It appears that considering an intercept $\gamma$ compensates for this \emph{ad-hoc} simplification and allows reasonably accurate results to be obtained.\footnote{Note that when estimating the model without constant one has to replace $\boldsymbol{\tau}_i(\boldsymbol{d})$ by
$
\widetilde{\boldsymbol{\tau}}_i(\boldsymbol{d})
:=
\left(\frac{\partial}{\partial \lambda}\mathbb{Y}_{n,i},-\mathbb{X}_{n,i}^T \right)^{T},
$
but $E\left[\widetilde{\boldsymbol{\tau}}_i(\boldsymbol{d}) \boldsymbol\Omega_{n,ij}^Z \boldsymbol\Omega_{n,ij}^X\mid \boldsymbol Z_k\right] \neq \boldsymbol 0_{p+1} $.
}
\subsection{Simulation setup} \label{sim}
During the simulation, we consider four different models. The models are given by
\begin{enumerate}
\item[Model 1:] $T(Y,\lambda_0) = X\beta_0 + m(Z) + \varepsilon$,
$m(Z) = \frac{\exp\{Z\}}{1 + \exp\{Z\}} + \frac{1}{3}$ with $Z \sim N(1,1)$, $\lambda_0 = 0$ and $\beta_0 = 1$,
$X = -\frac{2}{3} Z + u$ with $u \sim N(0,1)$ and $\varepsilon = \sqrt{\frac{1 + X^2}{2}} \; \widetilde u$ with $\widetilde u \sim N\left(0,\frac{1}{13}\right)$.
\item[Model 2:] $T(Y,\lambda_0) = X\beta_0 + m(Z) + \varepsilon$, $m(Z) = \frac{\exp\{Z\}}{1 + \exp\{Z\}} + 3$ with $Z \sim N(1,1)$, $\lambda_0 = 0.5$
and $\beta_0 = 1$, $X = -\frac{2}{3} Z + u$ with $u \sim N(0,1)$ and $\varepsilon \sim N\left(0,\frac{1}{9}\right)$.
\item[Model 3:] $T(Y,\lambda_0) = X\beta_0 + m(Z) + \varepsilon$, $m(Z) = \frac{\exp\{Z\}}{1 + \exp\{Z\}} - 1$ with $Z \sim U(-3,-1)$, $\lambda_0 = -1$
and $\beta_0 = 1$, $X = \frac{2}{3} Z + u$ with $u \sim U(-1,1)$ and $\varepsilon \sim U\left(-\sqrt{1/9},\sqrt{1/9}\right)$.
\item[Model 4:] $T(Y,\lambda_0) = X_1\beta_{10} + \boldsymbol X_2\boldsymbol \beta_{20} + m(Z_1, Z_2) + \varepsilon$,
$m(Z_1, Z_2) = \frac{1}{3} + Z_1 + Z_2 + Z_1Z_2$ with $Z_1, Z_2 \sim N(0,1)$, $\lambda_0 = 0$, $\beta_{10} = 1$,
$X_1 = -\frac{1}{3} \left( Z_1 + Z_2\right) + u$ with $u \sim N(0,1)$, $X_{2,l} \overset{i.i.d.}{\sim} Ber(0.2)$ and $\beta_{2,l} \overset{i.i.d.}{\sim} U(-1,1)$ for $l=1,\ldots,30$ , $\varepsilon \sim N\left(0,\frac{1}{9}\right)$.
\end{enumerate}
The main difference in the models is the transformation parameter $\lambda$. Model 1 and Model 4 have $\lambda_0 = 0$, whereas Model 2 has $\lambda_0 = 0.5$ and Model 3 $\lambda_0 = -1$. To ensure that $Y > 0$ in Model 3 we draw the random variables from uniform distributions. In all other models positivity of $Y$ is ensured as well. Model 1 has heteroskedastic error terms which are captured by the developed theory. Model 4 contains 30 dummy variables, $\boldsymbol X_2$, which take the value 1 with probability $20\%$.
In addition, Model 4 has the same structure as the model of the application we consider in section \ref{real_data}.
The estimators are computed by employing a normal kernel for $K(\cdot)$. $\boldsymbol Z$ is standardized componentwise by the corresponding standard deviations and $h \propto n^{-1/3.5}$. This bandwidth choice satisfies the assumptions of Theorem \ref{AN}.
The components of $\boldsymbol d$ defining the diagonal matrix $\boldsymbol{D}$ in $\boldsymbol{\Omega}_n$, are set equal to the componentwise standard deviations of $X$ and $\boldsymbol Z$ when $X$ is continuous. In the case of the dummy variables $\boldsymbol{X}_2$, an indicator of the event that the observations have the same value, is employed. For Model 4, we ensure in the simulations that for every observation there exists at least 4 observations with the same dummy variable combination.
In the estimation, we define a grid for values of $\lambda$ that are considered during the optimization.
This optimization grid for $\lambda$ is given in our simulation by the grid $[\lambda_0 -0.8, \lambda_0 + 0.8]$ with step size $0.001$. We minimize $G_n^{-\lambda} \widehat{\mathbb{Y}}_n(\lambda)^T \; \widehat {\mathbb{B}} _n \; G_n^{-\lambda} \widehat{\mathbb{Y}}_n(\lambda) $ over the defined grid to get $\widehat \lambda$ and $\boldsymbol{\widehat\beta}(\widehat \lambda)$, where $G_n = \prod\limits_{i=1}^{n} Y_i^{1/n}$ is the sample geometric mean.
In the simulation, we compare the proposed estimator where $\gamma$ is employed with the estimator that does not use $\gamma$. Both estimators converge asymptotically to a normal distribution. The only difference is that we have to replace $\boldsymbol{\tau}_i(\boldsymbol{d})$ by
$
\widetilde{\boldsymbol{\tau}}_i(\boldsymbol{d})
:=
\left(\frac{\partial}{\partial \lambda}\mathbb{Y}_{n,i},-\mathbb{X}_{n,i}^T \right)^{T}
$
in the case of the estimator without $\gamma$. It is therefore interesting to compare both estimators.
We consider the bias and standard deviation of the estimators as well as the power and size of the distance metric statistics proposed in Section \ref{Test}. In addition, we test by a simple Z-Test, if the estimated parameters are significantly different from the true value. Therefore, we employ the variance estimator stated in equation \eqref{variance_est}, and the necessary adjustments for the estimator without $\gamma$ are replacing $\mathbb{D}_{n,inf}$ by $\boldsymbol I_{n\times n}$ and $\mathbb{D}_n$ by $\boldsymbol \Omega_n$.
To estimate the error variance, we employ the Eiker-White variance estimator. In order to see the influence of the estimated $\boldsymbol \eta_\lambda$ on the variance, we consider all tests also without taking the estimation error of $\boldsymbol \eta_\lambda$ into account. Therefore, we replace $\widehat{\boldsymbol \Phi}_{n}$ by $\boldsymbol{\Omega}_n$ in the variance estimator. As mentioned before, we label this estimator SmoothMD*.
In addition, the Nonlinear two-stage Least Squares (NL2SLS) estimator for the Box-Cox model introduced by \citet{amemiya1981comparison} is considered as a competitor. In order to be able to employ this estimator, it is assumed that the function $m(\cdot)$ is known and, thus,
$m(\boldsymbol Z)$ can be added as additional regressor. The instruments are given therefore by $\boldsymbol V_i = (1,\boldsymbol X_i,\boldsymbol X_i^2,m(\boldsymbol Z_i),m(\boldsymbol Z_i)^2)$. We consider the Z-Test for the NL2SLS estimator as well where we employ the Eiker-White variance estimator again.
\subsection{Simulation results}
Table \ref{Bias_Std_Model1} states the results for the bias and standard deviation for $\lambda$ and $\beta$ in Model 1. All three estimators have comparable results for the bias and the bias decreases with sample size for $\beta$, whereas it is the lowest for $n=500$ for the SmoothMD estimators in the case of $\lambda$. Surprisingly, the standard deviation is also comparable for all three estimators even though $m(\cdot)$ is given for the NL2SLS estimator.
\begin{table}[H]
\caption{\textit{Bias and Standard Deviation of the estimators for $\lambda$ and $\beta$ in Model 1. } }
\begin{tabular}{@{}lcd{3.5}d{3.5}d{3.5}d{3.5}d{3.5}d{3.5}}
\toprule\midrule
& $s$ & \multicolumn{3}{c}{Bias} &\multicolumn{3}{c}{St. dev.}\\\cmidrule(lr){3-5} \cmidrule(lr){6-8} $n$ & & \multicolumn{1}{c}{250} & \multicolumn{1}{c}{500} & \multicolumn{1}{c}{1000} & \multicolumn{1}{c}{250} & \multicolumn{1}{c}{500} & \multicolumn{1}{c}{1000} \\
\midrule
$\lambda$ estimator & & & & & & & \\
\midrule
SmoothMD with $\gamma$ & $G_n$ & 0.003 & 0.0001 & 0.001 & 0.042 & 0.03 & 0.021 \\
SmoothMD without $\gamma$ & $G_n$ & 0.002 & 0.0001 & 0.001 & 0.041 & 0.029 & 0.02 \\
NL2SLS & $G_n$ & -0.003 & -0.001 & 0.0001 & 0.042 & 0.029 & 0.02 \\
\midrule
$\beta$ estimator & & & & & & & \\
\midrule
SmoothMD with $\gamma$ & $G_n$ & -0.001 & -0.001 & 0.0004 & 0.036 & 0.025 & 0.017\\
SmoothMD without $\gamma$ & $G_n$ & -0.001 & -0.001 & 0.0004 & 0.035 & 0.024 & 0.017 \\
NL2SLS & $G_n$ & -0.002 & -0.001 & 0.0002 & 0.035 & 0.024 & 0.016 \\
\midrule\bottomrule
\end{tabular}
\vskip 0.1cm
\textit{Notes: For the SmoothMD estimators, $h \propto n^{-1/3.5}$. The components of $\boldsymbol d$ are set equal to the componentwise standard deviations for all variables. The grid for $\lambda$ is $[\lambda_0-0.8,\lambda_0+0.8]$. 2000 Monte Carlo samples were used for all simulations. }
\label{Bias_Std_Model1}
\end{table}
Table \ref{Bias_Std_Model4} states the results for the bias and standard deviation for $\lambda$, $\beta_1$ and $\beta_2$, one representative parameter out of the 30 parameters in $\boldsymbol{\beta}_2$, in Model 4. All three estimators have comparable results for the bias and the bias decreases with sample size for all three parameters. In contrast to the results for Model 1, the standard deviation is smaller in the case of the NL2SLS estimator for $\beta_1$ and $\beta_2$. This result should be expected as $m(\cdot)$ is given for the NL2SLS estimator. The standard deviations for the SmoothMD estimators with and without $\gamma$ are as in Model 1, nearly the same.
\begin{table}[H]
\caption{\textit{Bias and Standard Deviation of the estimators for $\lambda$, $\beta_1$ and $\beta_2$ in Model 4. } }
\begin{tabular}{@{}lcd{3.5}d{3.5}d{3.5}d{3.5}d{3.5}d{3.5}}
\toprule\midrule
& $s$ & \multicolumn{3}{c}{Bias} &\multicolumn{3}{c}{St. dev.}\\\cmidrule(lr){3-5} \cmidrule(lr){6-8} $n$ & & \multicolumn{1}{c}{250} & \multicolumn{1}{c}{500} & \multicolumn{1}{c}{1000} & \multicolumn{1}{c}{250} & \multicolumn{1}{c}{500} & \multicolumn{1}{c}{1000} \\
\midrule
$\lambda$ estimator & & & & & & & \\
\midrule
SmoothMD with $\gamma$ & $G_n$ & 0.0001 & -0.0002 & -0.0001 & 0.015 & 0.01 & 0.008 \\
SmoothMD without $\gamma$ & $G_n$ & 0.0001 & -0.0002 & -0.0001 & 0.015 & 0.01 & 0.008 \\
NL2SLS & $G_n$ & -0.0004 & -0.0002 & -0.0002 & 0.014 & 0.01 & 0.005 \\
\midrule
$\beta_1$ estimator & & & & & & & \\
\midrule
SmoothMD with $\gamma$ & $G_n$ & -0.002 & -0.002 & -0.001 & 0.036 & 0.023 & 0.017\\
SmoothMD without $\gamma$ & $G_n$ & -0.002 & -0.002 & -0.001 & 0.036 & 0.023 & 0.017 \\
NL2SLS & $G_n$ & 0.0004 & -0.0001 & -0.0001 & 0.025 & 0.015 & 0.011 \\
\midrule
$\beta_2$ estimator & & & & & & & \\
\midrule
SmoothMD with $\gamma$ & $G_n$ & 0.004 & -0.001 & -0.0002 & 0.133 & 0.065 & 0.042\\
SmoothMD without $\gamma$ & $G_n$ & 0.004 & -0.001 & -0.0002 & 0.133 & 0.065 & 0.042 \\
NL2SLS & $G_n$ & -0.001 & 0.002 & -0.0006 & 0.091 & 0.046 & 0.028 \\
\midrule\bottomrule
\end{tabular}
\vskip 0.1cm
\textit{Notes: For the SmoothMD estimators, $h \propto n^{-1/3.5}$. The components of $\boldsymbol d$ are set equal to the componentwise standard deviations for all continuous variables and for the dummy variables an indicator of the event that the observations have the same value is employed. The grid for $\lambda$ is $[\lambda_0-0.8,\lambda_0+0.8]$. 2000 Monte Carlo samples were used for all simulations. }
\label{Bias_Std_Model4}
\end{table}
\begin{table}[H]
\caption{\textit{Empirical Level for distance metric statistics of the estimators for $\lambda$ and $\beta$ in Model 2. } }
\begin{tabular}{@{}lcd{3.5}d{3.5}d{3.5}d{3.5}d{3.5}d{3.5}}
\toprule\midrule& $s$ & \multicolumn{3}{c}{5\% level} &\multicolumn{3}{c}{10\% level}\\\cmidrule(lr){3-5} \cmidrule(lr){6-8} $n$ & & \multicolumn{1}{c}{250} & \multicolumn{1}{c}{500} & \multicolumn{1}{c}{1000} & \multicolumn{1}{c}{250} & \multicolumn{1}{c}{500} & \multicolumn{1}{c}{1000} \\
\midrule
Test for $\lambda$ & & & & & & & \\
\midrule
SmoothMD with $\gamma$ & $G_n$ & 10.3 & 8.15 & 7.8 & 12.75 & 10.65 & 11.45 \\
SmoothMD* with $\gamma$ & $G_n$ & 10.15 & 7.95 & 7.7 & 12.75 & 10.6 & 11.45 \\
SmoothMD without $\gamma$ & $G_n$ & 9.55 & 7.0 & 6.0 & 12.45 & 10.85 & 10.35 \\
\midrule
Test for $\beta$ & & & & & & & \\
\midrule
SmoothMD with $\gamma$ & $G_n$ & 10.4 & 8.25 & 7.95 & 12.55 & 11.1 & 11.9 \\
SmoothMD* with $\gamma$ & $G_n$ & 10.2 & 8.3 & 7.7 & 12.55 & 11.3 & 11.1 \\
SmoothMD without $\gamma$ & $G_n$ & 9.7 & 6.85 & 6.7 & 12.55 & 10.2 & 11.2 \\
\midrule \bottomrule
\end{tabular}
\vskip 0.1cm
\textit{Notes: For the SmoothMD estimators, $h \propto n^{-1/3.5}$. The components of $\boldsymbol d$ are set equal to the componentwise standard deviations for all variables. The variances are estimated by the Eiker-White variance estimator. For SmoothMD* the additional variance part due to the estimation of $\boldsymbol{\eta_\lambda}$ is not taken into account. For SmoothMD the additional variance part is taken into account. 2000 Monte Carlo samples were used for all simulations.}
\label{Emp_DistM_Model2}
\end{table}
Table \ref{Emp_DistM_Model2} states the empirical level for distance metric statistics of the estimators for $\lambda$ and $\beta$ in Model 2. Here we state the results for the SmoothMD estimator with correctly estimated variance as well as with the variance estimate that does not account for the estimation of $\boldsymbol \eta_\lambda$. For all three estimators, the empirical levels converge to the nominal levels if the sample size increases and $\beta$ seems to need a larger sample size than $\lambda$ to get close to the nominal level. However, in this setup, the results of the SmoothMD estimators with and without $\gamma$ differ. In addition, it can be seen that the estimator SmoothMD* leads to almost the same results as SmoothMD.
Table \ref{Emp_LevZ_Model4} states the empirical level for the Z-Tests for $\lambda$, $\beta_1$ and $\beta_2$ in Model 4. The fact that we do not consider the estimation error has almost no influence on the results. In addition, both SmoothMD versions lead to similar results. However, in order to get close to the nominal level, the sample size needs to be large as only for $n = 1000$ do the SmoothMD estimators get close to the nominal level. The NL2SLS estimator gives more convincing results for smaller sample sizes. Note that the dummy variable coefficient $\beta_2$ seems to require a larger sample size than the other two parameters to get close to the nominal level when employing the SmoothMD estimators.
\begin{table}[H]
\caption{\textit{Empirical Level for Z-Tests of the estimators for $\lambda$, $\beta_1$ and $\beta_2$ in Model 4. } }
\begin{tabular}{@{}lcd{3.5}d{3.5}d{3.5}d{3.5}d{3.5}d{3.5}}
\toprule\midrule& $s$ & \multicolumn{3}{c}{5\% level} &\multicolumn{3}{c}{10\% level}\\\cmidrule(lr){3-5} \cmidrule(lr){6-8} $n$ & & \multicolumn{1}{c}{250} & \multicolumn{1}{c}{500} & \multicolumn{1}{c}{1000} & \multicolumn{1}{c}{250} & \multicolumn{1}{c}{500} & \multicolumn{1}{c}{1000} \\
\midrule
Test for $\lambda$ & & & & & & & \\
\midrule
SmoothMD with $\gamma$ & $G_n$ & 8.4 & 7.1 & 5.5 & 15.45 & 13.5 & 10.55 \\
SmoothMD* with $\gamma$ & $G_n$ & 8.5 & 7.1 & 5.45 & 15.3 & 13.45 & 10.55 \\
SmoothMD without $\gamma$ & $G_n$ & 9.6 & 8.15 & 5.2 & 16.45 & 12.6 & 11.65 \\
NL2SLS & $G_n$ & 9.6 & 7.6 & 6.0 & 16.8 & 12.5 & 11.1 \\
\midrule
Test for $\beta_1$ & & & & & & & \\
\midrule
SmoothMD with $\gamma$ & $G_n$ & 11.8 & 8.55 & 6.4 & 18.25 & 14.05 & 11.75 \\
SmoothMD* with $\gamma$ & $G_n$ & 11.8 & 8.55 & 6.45 & 18.25 & 14.05 & 11.8 \\
SmoothMD without $\gamma$ & $G_n$ & 11.75 & 8.9 & 6.4 & 18.35 & 14.8 & 12.1 \\
NL2SLS & $G_n$ & 8.25 & 4.95 & 5.1 & 13.95 & 10.4 & 10.35 \\
\midrule
Test for $\beta_2$ & & & & & & & \\
\midrule
SmoothMD with $\gamma$ & $G_n$ & 13.6 & 8.25 & 7.05 & 20.7 & 15.65 & 12.7 \\
SmoothMD* with $\gamma$ & $G_n$ & 13.7 & 8.25 & 7.05 & 20.75 & 15.55 & 12.75 \\
SmoothMD without $\gamma$ & $G_n$ & 13.75 & 8.35 & 6.65 & 20.6 & 14.55 & 12.25 \\
NL2SLS & $G_n$ & 7.65 & 6.55 & 4.7 & 13.45 & 13.05 & 9.55 \\
\midrule \bottomrule
\end{tabular}
\vskip 0.1cm
\textit{Notes: For the SmoothMD estimators, $h \propto n^{-1/3.5}$. The components of $\boldsymbol d$ are set equal to the componentwise standard deviations for all continuous variables and for the dummy variables, an indicator of the event that the observations have the same value is employed. The variances are estimated by the Eiker-White variance estimator. For SmoothMD* the additional variance part due to the estimation of $\boldsymbol{\eta_\lambda}$, is not taken into account. For SmoothMD, the additional variance part is taken into account. $\beta_2$ is one representative parameter out of the 30 parameters in $\boldsymbol{\beta}_2$. 2000 Monte Carlo samples were used for all simulations.}
\label{Emp_LevZ_Model4}
\end{table}
\begin{figure}[h]
\begin{minipage}{0.4\textwidth}
\caption{\textit{Power function of the distance metric \\ statistic
for $\lambda$ of Model 3 with $n=250$.}}
\includegraphics[scale=0.4]{{Power_Plot_lambda_10_alpha01_Model4_N=250H=3,5}}
\label{fig_power_4.1}
\end{minipage}
\hfill
\begin{minipage}{0.4\textwidth}
\caption{\textit{Power function of the distance metric \\ statistic
for $\lambda$ of Model 3 with $n=1000$.}}
\includegraphics[scale=0.4]{{Power_Plot_lambda_10_alpha01_Model4_N=1000H=3,5}}
\label{fig_power_4.2}
\end{minipage}
\begin{spacing}{0.8}
\flushleft{\textit{\footnotesize Notes: For the SmoothMD estimators, $h \propto n^{-1/3.5}$. The components of $\boldsymbol d$ are set equal to the componentwise standard deviations for all variables. The variances are estimated by the Eiker-White variance estimator. Only the SmoothMD estimators that take the additional variance part due to the estimation of $\boldsymbol{\eta_\lambda}$ into account are considered. 2000 Monte Carlo samples were used for all simulations. The nominal level is $ 10\%$.}}
\end{spacing}
\end{figure}
Figures \ref{fig_power_4.1} and \ref{fig_power_4.2} state the power functions of the distance metric statistic for $\lambda$ in Model 3 with $n=250$ and $n = 1000$. In the case of $n=250$, the power function is skewed and the power for values larger than $-1$ is small. In addition, the power function is smaller than the nominal value at $-0.85$ and $-0.7$.
For the SmoothMD estimator without $\gamma$, the power function is larger than for the SmoothMD estimator with $\gamma$ at values larger than $-1$. These issues disappear for the larger sample size $n=1000$.
Figures \ref{fig_power_2.1} and \ref{fig_power_2.2} state the power functions of the distance metric statistic for $\beta$ in Model 2 with $n=250$ and $n = 500$.
As in Figure \ref{fig_power_4.1} the power function for $n=250$ is skewed but the effect is less distinct. However, the power function is smaller than the nominal value at $0.8$. For the SmoothMD estimator without $\gamma$, the power function is larger than for the SmoothMD estimator with $\gamma$ at values smaller than $1$ for both sample sizes. For $n=500$, the skewness is less pronounced and the power function has no values lower than the nominal value. The main conclusion from both power functions is that the samples size should not be too small so that the tests have a reasonable power.
\begin{figure}[h!]
\begin{minipage}{0.4\textwidth}
\caption{\textit{Power function of the distance metric \\ statistic
for $\beta$ of Model 2 with $n=250$.}}
\includegraphics[scale=0.4]{{Power_Plot_beta_10_alpha01_Model2_N=250H=3,5}}
\label{fig_power_2.1}
\end{minipage}
\hfill
\begin{minipage}{0.4\textwidth}
\caption{\textit{Power function of the distance metric \\ statistic
for $\beta$ of Model 2 with $n=500$.}}
\includegraphics[scale=0.4]{{Power_Plot_beta_10_alpha01_Model2_N=500H=3,5}}
\label{fig_power_2.2}
\end{minipage}
\begin{spacing}{0.8}
\flushleft{\textit{\footnotesize Notes: For the SmoothMD estimators, $h \propto n^{-1/3.5}$. The components of $\boldsymbol d$ are set equal to the componentwise standard deviations for all variables. The variances are estimated by the Eiker-White variance estimator. Only the SmoothMD estimators that take the additional variance part due to the estimation of $\boldsymbol{\eta_\lambda}$ into account are considered. 2000 Monte Carlo samples were used for all simulations. The nominal level is $ 10\%$.}}
\end{spacing}
\end{figure}
Before we consider a real data application, we close the discussion with Figures \ref{fig_m_1.1} and \ref{fig_m_3.1}. The figures state the estimated $m(Z)$ for Model 1 with $n=250$ and for Model 3 with $n=500$. For the estimation, the NW estimator with same kernel and bandwidth as for the SmoothMD estimator was used. No matter if the SmoothMD estimator with or without $\gamma$ is employed, the results are very accurate. In practice one can of course use cross validation to choose the bandwidth or employ the local linear estimator instead of the NW estimator.
\begin{figure}[h!]
\begin{minipage}{0.4\textwidth}
\caption{\textit{Estimated $m(Z)$ for Model 1 with $n=250$.}}
\includegraphics[scale=0.4]{{Estimated_mZ_alpha01_Model6_N=250H=3,5}}
\label{fig_m_1.1}
\end{minipage}
\hfill
\begin{minipage}{0.4\textwidth}
\caption{\textit{Estimated $m(Z)$ for Model 3\\ with $n=500$.}}
\includegraphics[scale=0.4]{{Estimated_mZ_alpha01_Model4_N=500H=3,5}}
\label{fig_m_3.1}
\end{minipage}
\begin{spacing}{0.8}
\flushleft{\textit{\footnotesize Notes: For the estimation the NW estimator with normal kernel and $h \propto n^{-1/3.5}$ is employed. The $25\%$ and $75\%$ quantiles as well as the mean are reported. 2000 Monte Carlo samples were used for all simulations.}}
\end{spacing}
\end{figure}
\subsection{Real data application} \label{real_data}
We consider in this section an application of our estimator to investigate the returns of social and cognitive skills in the labor market. For this purpose, we apply the proposed transformation partially linear estimator to a dataset studied in \citet{deming2017growing}. In particular, we consider the regressions (4) and (5) in TABLE I of \citet{deming2017growing} that are based on the National Longitudinal Survey of Youth 1979 (NLSY79). NLSY79 is a nationally representative sample taken in the US, of young people aged from 14 to 22. The survey was conducted yearly from 1979 to 1993 and biannually from 1994 through 2012.
\citet{deming2017growing} estimates the model
\begin{equation}
\begin{aligned}
\log(wage_{ijt}) &= \alpha + \beta_1 \cdot COG_i + \beta_2 \cdot SS_i + \beta_3 \cdot COG_i \times SS_i + \beta_4 \cdot NCOG_i \\
& \quad + \boldsymbol{C}_{ijt}^T \boldsymbol{\rho} + \delta_j + \zeta_t + \varepsilon_{ijt},
\end{aligned}
\label{reg_dem}
\end{equation}
where $COG$, $SS$ and $NCOG$ denote measures of cognitive, social and noncognitive skills. The model includes controls $\boldsymbol{C}$ for race-by-gender indicators, indicators for region and urbanicity as well as age (indexed by $j$) and year (indexed by $t$) fixed effects.
In his paper, \citet{deming2017growing} develops a theoretical model that is written in levels instead of logs as in equation \eqref{reg_dem}. Nevertheless, he estimates the log-linearized model in his paper to follow standard practice in the literature, as he argues. Results for the model in levels are stated in an online appendix. Therefore, it makes sense to use the Box-Cox transformation for $wage$ and estimate the transformation parameter $\lambda$ together with the remaining model parameters to decide whether the model in logs or in levels is more appropriate.
Furthermore, we consider an unknown functional form for cognitive and social skills to see if the linear form, $\beta_1 \cdot COG + \beta_2 \cdot SS + \beta_3 \cdot COG \times SS$ used by \citet{deming2017growing}, is reasonable. The transformation partially linear model is, thus, given by
\begin{align}
T(wage_{ijt},\lambda) = m(COG_i, SS_i) + \beta \cdot NCOG_i + \boldsymbol{C}_{ijt}^T \boldsymbol{\rho} + \delta_j + \zeta_t + \varepsilon_{ijt},
\label{reg_smooth}
\end{align}
where $m(\cdot)$ is an unknown function. The model stated in \eqref{reg_smooth} is closely related to Model 4 of the simulations in section \ref{sim}, where $Y = wage$, $\boldsymbol Z = (COG, SS)^T$, $X_1 = NCOG$ and $X_2 = (\boldsymbol{C}^T,1,1)^T$.
As proxy for cognitive skills, the Armed Forces Qualifying Test (AFQT) was taken. \citet{deming2017growing} uses raw scores from \citet{altonji2012changes}
and normalizes them to have mean 0 and standard deviation 1. The social skill measure is constructed from the following four variables of the NLSY79:
\begin{enumerate}
\item Self-reported sociability in 1981 (extremely shy, somewhat shy, somewhat outgoing, extremely outgoing)
\item Self-reported sociability in 1981 at age 6 (retrospective)
\item The number of clubs in which the respondent participated in high school
\item Participation in high school sports (yes/no).
\end{enumerate}
Each variable is normalized to have mean 0 and standard deviation 1. The social skill measure is the average of these four normalized variables (also normalized to standard deviation 1). In addition to social and cognitive skill measures \citet{deming2017growing} includes a noncognitive skill measure in his regression. He uses the Rotter Locus of Control and the Rosenberg Self-Esteem Scale as also used by \citet{heckman2006effects}. In the following discussion we use these variables to estimate the models stated in \eqref{reg_dem} and \eqref{reg_smooth}.
\citet{deming2017growing} used a weighted log-linearized OLS estimator to estimate the returns of cognitive and social skills on wage and excluded respondents under the age of $23$ or who were enrolled in school. The weighting was necessary as in each survey year of the NLSY79 a set of sampling weights was constructed. These weights provided the researcher with an estimate of how many individuals in the United States each respondent's answers represented. We also employ these weights in our analysis.
Table \ref{base_reg} shows the regression results. The first column, (4), provides the results of \citet{deming2017growing} estimating equation \eqref{reg_dem}. The second and third column state the transformation partially linear estimator of equation \eqref{reg_smooth} with and without employing $\gamma$. The fourth and fifth column state the transformation partially linear estimator of equation \eqref{reg_smooth} with and without employing $\gamma$ imposing $\lambda=0$. This is a standard partially linear model as studied by \citet{robinson1988root} and \citet{li1996root}.
For the inner smoothing of the estimations in column 2-5, we use a normal kernel with $h \propto n^{-1/3.5}$, which is in line with the developed theory and usual bandwidth choices for bivariate smoothing with the Nadaraya-Watson estimator. The components of $\boldsymbol d$ defining the diagonal matrix $\boldsymbol{D}$ in $\boldsymbol{\Omega}_n$ are set equal to the componentwise standard deviations for all continuous variables. In the case of the controls and fixed effects, an indicator of the event that the observations have the same value is employed.
The optimization grid for $\lambda$ is given by $[-0.1, 0.1]$ with step size $0.001$.\footnote{We evaluated subsamples of the dataset before we made the final estimation. The estimated $\lambda$'s in the subsamples are contained in the used grid.} We minimize
$G_n^{-\lambda} \widehat{\mathbb{Y}}_n(\lambda)^T \; \widehat {\mathbb{B}} _n \; G_n^{-\lambda} \widehat{\mathbb{Y}}_n(\lambda) $ over the defined grid to get $\widehat \lambda$ and the estimates of the remaining coefficients, where $G_n$
is the sample geometric mean.
The results in the first column of Table \ref{base_reg} show that all \citet{deming2017growing} estimated coefficients are significantly different from $0$. In the remaining four columns, we cannot state parameter estimates for cognitive and social skills and the interaction of both as these variables are contained in $m(\cdot)$. However, the parameter estimates for noncognitive skills are comparable to the estimate from the first column. In addition, the estimates for $\lambda$ with and without $\gamma$ are equal and close to zero which would imply that a log-transformation of the dependent variable is appropriate. The estimated coefficient for noncognitive skills is significantly different from $0$ in all SmoothMD estimations whereas both estimates for $\lambda$ are not significantly different from $0$.
In order to check if the linear specification for cognitive and social skills employed by \citet{deming2017growing} is reasonable, we proceed as follows. We estimate the parameters of the transformation partially linear model as stated in \eqref{reg_smooth} to obtain the residuals
$$
\widehat\varepsilon_{ijt} = T(wage_{ijt},\widehat\lambda) - \widehat\beta \cdot NCOG_i - \boldsymbol{C}_{ijt}^T \widehat{\boldsymbol{\rho}}- \widehat \delta_j - \widehat \zeta_t.
$$
We now estimate the unknown function $m(\cdot)$ by smoothing $\widehat\varepsilon$ with the NW estimator. In addition, we also regress $\widehat\varepsilon$ on $COG$, $SS$ and $COG \times SS$. To see if the linear specification is appropriate, we compare the MSE of the linear and nonlinear estimates. We employ a normal density kernel for the NW estimator and let $h \propto n^{-1/6}$.
Table \ref{base_OLS} states the results where OLS indicates that we used the linear model to fit the residuals. The MSE of the linear and nonlinear estimates are identical no matter if we use the SmoothMD estimator with or without $\gamma$ to estimate the unknown model parameters. The same holds true for the SmoothMD estimator with or without $\gamma$ where $\lambda = 0$ is imposed. All results show that the linear representation of \citet{deming2017growing} seems to be reasonable.
\begin{table}[H]
\begin{adjustwidth}{+0.0cm}{+0.0cm}
\caption{\textit{Labor Market Returns to Cognitive and Social Skills in the NLSY79}}
\begin{tabular}{@{}lccccc}
\toprule
\midrule
Outcome: (log) hourly wage &(4) & SmoothMD & SmoothMD & SmoothMD & SmoothMD \\
\hskip 1.6cm(in 2012 dollars) & & with $\gamma$ & without $\gamma$ & with $\gamma$, $\lambda = 0$ & without $\gamma$, $\lambda = 0$ \\
\midrule
$\lambda$ & - & \multicolumn{1}{l}{\;-0.007} & \multicolumn{1}{l}{\;-0.007} & - & - \\
& & [0.005] & [0.005] & & \\
Cognitive skills & 0.189*** & - & - & - & - \\
& [0.007] & & & & \\
Social skills & 0.043*** & - & - & - & - \\
& [0.006] & & & & \\
Cognitive $\times$ Social & 0.019*** & - & - & - & - \\
& [0.006] & & & & \\
Noncognitive skills & 0.048*** & 0.047*** & 0.047*** & 0.048*** & 0.048*** \\
& [0.006] & [0.004] & [0.004] & [0.004] & [0.004]\\
Demographics and age/ & X & X & X & X & X \\
year fixed effects & & & & & \\
Number of Observations & 126191 & 126191 & 126191 & 126191 & 126191 \\
\midrule
\bottomrule
\end{tabular}
\vskip 0.1cm
\textit{Notes: The data source is the National Longitudinal Survey of Youth 1979 cohort (NLSY79). (4) denotes the OLS regression proposed by \citet{deming2017growing}. In all SmoothMD estimations, $h \propto n^{-1/3.5}$. The components of $\boldsymbol d$ are set equal to the componentwise standard deviations for all continuous variables and for controls and fixed effects, an indicator of the event that the observations have the same value is employed. The grid for $\lambda$ is $[-0.1,0.1]$ and $s = G_n$. Cognitive skills are measured by each NLSY79 respondent's score on the Armed Forces Qualifiying Test (AFQT) and are normalized to have mean 0 and standard deviation 1. The AFQT score crosswalk of \citet{altonji2012changes} is used. Social skill is a standardized composite of four variables, (i) sociability in childhood, (ii) sociability in adulthood, (iii) participation in high school clubs and (iv) participation in team sports; see the text and \citet{deming2017growing} for details on the construction of the social skills measure. The noncognitive skills measure is the normalized average of the Rotter and Rosenberg scores in the NLSY. The regressions also control for race-by-gender indicator variables, age, year, census region and urbanicity. Standard errors are in brackets and are clustered at the individual level for (4). The remaining standard errors are estimated by the Eiker-White variance estimator. ***$p < .01$, **$p < .05$, *$p < .1$}
\label{base_reg}
\end{adjustwidth}
\end{table}
\begin{table}[H]
\caption{\textit{MSE of estimated nonlinear part in the transformation partially linear model}}
\begin{tabular}{@{}lcccccccc}
\toprule \midrule
& \multicolumn{2}{c}{SmoothMD} &\multicolumn{2}{c}{SmoothMD} & \multicolumn{2}{c}{SmoothMD} &\multicolumn{2}{c}{SmoothMD}\\
& \multicolumn{2}{c}{with $\gamma$} &\multicolumn{2}{c}{without $\gamma$} & \multicolumn{2}{c}{with $\gamma$, $\lambda = 0$} &\multicolumn{2}{c}{without $\gamma$, $\lambda = 0$}\\
\cmidrule(lr){2-3} \cmidrule(lr){4-5} \cmidrule(lr){6-7} \cmidrule(lr){8-9}
& OLS & NW & OLS& LL & OLS & NW & OLS& LL\\
\midrule
MSE & 0.282 & 0.277 & 0.282 & 0.277 & 0.293 & 0.288 & 0.293 & 0.288 \\
Number of Observations & 126191 & 126191 & 126191 & 126191 & 126191 & 126191 & 126191 & 126191 \\
\midrule
\bottomrule
\end{tabular}
\vskip 0.1cm
\textit{Notes: For the NW estimator a normal kernel with $h \propto n^{-1/6}$ is employed. OLS indicates that the linear model is used to fit the residuals.}
\label{base_OLS}
\end{table}
In a second step, we include \textit{years of completed education} as an additional explanatory variable in the regression models. In one of his estimations \citet{deming2017growing} controls for years of education as well. Table \ref{edu_reg} states the regression results for all considered models.
The first column, (5), provides the results of \citet{deming2017growing} estimating equation \eqref{reg_dem} with \textit{years of completed education} as control. The results show that all \citet{deming2017growing} estimated coefficients are significantly different from $0$. However, the coefficients become smaller compared to the first specification. In addition, the coefficient of the interactive effect is only significant at the $10\%$ level. In the remaining four columns, the parameter estimates for noncognitive skills are comparable to the estimate from the first column. In addition, the estimates for $\lambda$ with and without $\gamma$ are equal and close to zero which would imply that a log-transformation of the dependent variable is appropriate. The estimated coefficient for noncognitive skills is significantly different from $0$ in all SmoothMD estimations whereas both estimates for $\lambda$ are not significantly different from $0$.
Table \ref{edu_OLS} states the MSE of the estimated nonlinear part in the transformation partially linear models. The MSE of the linear and nonlinear estimates are identical no matter if we use the SmoothMD estimator with or without $\gamma$ to estimate the unknown model parameters. The same holds true for the SmoothMD estimator with or without $\gamma$ where $\lambda = 0$ is imposed. All results show that the linear representation of \citet{deming2017growing} seems to be reasonable.
\begin{table}[H]
\begin{adjustwidth}{+0.0cm}{+0.0cm}
\caption{\textit{Labor Market Returns to Cognitive and Social Skills in the NLSY79 controlling for education}}
\begin{tabular}{@{}lccccc}
\toprule
\midrule
Outcome: (log) hourly wage &(5) & SmoothMD & SmoothMD & SmoothMD & SmoothMD \\
\hskip 1.6cm(in 2012 dollars) & & with $\gamma$ & without $\gamma$ & with $\gamma$, $\lambda = 0$ & without $\gamma$, $\lambda = 0$ \\
\midrule
$\lambda$ & - & 0.002~~~~~ & 0.002~~~~~ & - & - \\
& & [0.005] & [0.005] & & \\
Cognitive skills & 0.126*** & - & - & - & - \\
& [0.008] & & & & \\
Social skills & 0.029*** & - & - & - & - \\
& [0.006] & & & & \\
Cognitive $\times$ Social & 0.011*~~~& - & - & - & - \\
& [0.006] & & & & \\
Noncognitive skills & 0.040*** & 0.037*** & 0.037*** & 0.037*** & 0.037*** \\
& [0.006] & [0.004] & [0.004] & [0.004] & [0.004]\\
Demographics and age/ & X & X & X & X & X \\
year fixed effects & & & & & \\
Years of completed education & X & X & X & X & X \\
Number of Observations & 126191 & 126191 & 126191 & 126191 & 126191 \\
\midrule
\bottomrule
\end{tabular}
\vskip 0.1cm
\textit{Notes: The data source is the National Longitudinal Survey of Youth 1979 cohort (NLSY79). (5) denotes the OLS regression proposed by \citet{deming2017growing}. In all SmoothMD estimations, $h \propto n^{-1/3.5}$. The components of $\boldsymbol d$ are set equal to the componentwise standard deviations for all continuous variables and for controls and fixed effects, an indicator of the event that the observations have the same value is employed. The grid for $\lambda$ is $[-0.1,0.1]$ and $s = G_n$. Cognitive skills are measured by each NLSY79 respondent's score on the Armed Forces Qualifiying Test (AFQT) and are normalized to have mean 0 and standard deviation 1. The AFQT score crosswalk of \citet{altonji2012changes} is used. Social skill is a standardized composite of four variables (i) sociability in childhood, (ii) sociability in adulthood, (iii) participation in high school clubs and (iv) participation in team sports; see the text and \citet{deming2017growing} for details on the construction of the social skills measure. The noncognitive skills measure is the normalized average of the Rotter and Rosenberg scores in the NLSY. The regressions also control for race-by-gender indicator variables, age, year, census region, urbanicity and years of completed education. Standard errors are in brackets and are clustered at the individual level for (5). The remaining standard errors are estimated by the Eiker-White variance estimator. ***$p < .01$, **$p < .05$, *$p < .1$}
\label{edu_reg}
\end{adjustwidth}
\end{table}
Before we close the section, we plot the estimated labor market returns to cognitive and social skills of model \eqref{reg_smooth} with and without controlling for years of completed education. The returns are estimated with the NW estimator employing a normal kernel with $h \propto n^{-1/6}$. Figures \ref{fig_app_1.1} and \ref{fig_app_2.2} present the results. Note that the mean was subtracted. The return increases no matter if the social or cognitive indicator is increased. However, the cognitive effect seems to be stronger. In addition, if we control for education, it seems that being too social might sometimes lower the wage somewhat. Nevertheless, both plots confirm that the linear model used by \citet{deming2017growing} is reasonable.
\begin{table}[H]
\caption{\textit{MSE of estimated nonlinear part in the transformation partially linear model controlling for education}}
\begin{tabular}{@{}lcccccccc}
\toprule \midrule
& \multicolumn{2}{c}{SmoothMD} &\multicolumn{2}{c}{SmoothMD} & \multicolumn{2}{c}{SmoothMD} &\multicolumn{2}{c}{SmoothMD}\\
& \multicolumn{2}{c}{with $\gamma$} &\multicolumn{2}{c}{without $\gamma$} & \multicolumn{2}{c}{with $\gamma$, $\lambda = 0$} &\multicolumn{2}{c}{without $\gamma$, $\lambda = 0$}\\
\cmidrule(lr){2-3} \cmidrule(lr){4-5} \cmidrule(lr){6-7} \cmidrule(lr){8-9}
& OLS & NW & OLS& LL & OLS & NW & OLS& LL\\
\midrule
MSE & 0.288 & 0.283 & 0.288 & 0.283 & 0.284 & 0.280 & 0.284 & 0.280 \\
Number of Observations & 126191 & 126191 & 126191 & 126191 & 126191 & 126191 & 126191 & 126191 \\
\midrule
\bottomrule
\end{tabular}
\vskip 0.1cm
\textit{Notes: For the NW estimator a normal kernel with $h \propto n^{-1/6}$ is employed. OLS indicates that the linear model is used to fit the residuals.}
\label{edu_OLS}
\end{table}
\begin{figure}[h]
\begin{minipage}{0.4\textwidth}
\caption{\textit{Estimated Labor Market Returns to\\ Cognitive
and Social Skills in the NLSY79.}}
\includegraphics[scale=0.8]{{SmoothMD_with_gamma}}
\label{fig_app_1.1}
\end{minipage}
\hfill
\begin{minipage}{0.4\textwidth}
\caption{\textit{Estimated Labor Market Returns to \\ Cognitive
and Social Skills in the NLSY79 \\ controlling for education.}}
\includegraphics[scale=0.8]{{SmoothMD_with_gamma_edu}}
\label{fig_app_2.2}
\end{minipage}
\begin{spacing}{0.8}
\flushleft{\textit{\footnotesize Notes: The coefficients are estimated by the SmoothMD estimator with $\gamma$. For the NW estimator a normal kernel with $h \propto n^{-1/6}$ is employed.}}
\end{spacing}
\end{figure}
\section{Discussion} \label{discussion}
In this paper, we study the semiparametric partially linear model with Box-Cox transformed dependent variable. We follow the SmoothMD approach introduced by \citet{lavergne2013smooth}, which is based on conditional moment restrictions, and we extend it to the case with an infinite-dimensional nuisance parameter. Our results are new both for the transformation regression models and the semiparametric partially linear model. We establish model identification, consistency as well as $\sqrt{n}$-asymptotic normality. In addition, we proposed a distance metric statistic to test the model parameters. A Monte Carlo experiment showed the usefulness of the proposed estimator in finite samples. An application to a large real data sample already studied in Labor Economics is also reported.
The SmoothMD approach is a convenient approach for nonlinear regression models. It could be interpreted as a generalized least-squares method where the weights are given by a suitable positive-definite matrix with each entry representing a measure of discrepancy between a pair of covariate vectors. It is unnecessary to localize the measure of discrepancy between the observed covariate vectors, and this represents a significant advantage for the practitioner who thus avoids the choice of an additional tuning parameter. One price to pay is on the semiparametric efficiency for the estimators of the finite dimensional parameters. However, \citet{lavergne2013smooth} showed that a two-step procedure, where the first step allows nonparametrically suitable weights to be estimated for building an asymptotically optimal measure of discrepancy, allows to achieve semiparametric efficiency. An efficient estimator would also induce a distance metric test statistic with a usual asymptotic chi-square distribution under the null hypothesis. We expect that the same results extend to the present framework. However, with our semiparametric model, the two-step procedure would likely result in a numerically unstable, complex to calibrate, inference procedure. We believe this theoretical, and quite technical, refinement to be of little use for the applications, and we therefore do not consider it.
Another slight drawback for using SmoothMD with a fixed weighting matrix, is the condition that $\boldsymbol Z$ should be of dimension $q$ less than or equal to 3. This restriction would not be binding in most applications. However, if necessary, one could use higher-order kernels to diminish the bias induced by the nonparametric estimation of the nuisance parameter, and thus allow for larger $q$. It is noticeable that SmoothMD does not involve any denominator and therefore the higher-order kernels would not induce the numerical problems, due to division by zero, usually encountered in semiparametric methods.
\addcontentsline{toc}{section}{Acknowledgments}
\section*{Acknowledgments}
This research was supported by the DFG through KN 567/5-1 and the Hausdorff Center for Mathematics. We furthermore thank the Regional Computing Center of the
University of Cologne (RRZK) for providing computing time on the DFG-funded High Performance Computing (HPC) system CHEOPS as well as support. This work was started while the first author was visiting CREST-ENSAI. V. Patilea acknowledges support from ‘Models and mathematical processing of very large data’, a Joint Research Initiative under the aegis of Risk Foundation, with partnership of MEDIAMETRIE and GENES, and from the Romanian Minister of Education and Research, CNCS – UEFISCDI, project number PN-III-P4-ID-PCE-2020-1112, within PNCDI III.
\clearpage
\newpage
\begin{appendix}
\numberwithin{equation}{section}