Regularized Nonlinear Regression with Dependent Errors

and its Application to a Biomechanical Model

Hojun You1, Kyubaek Yoon3, Wei-Ying Wu4, Jongeun Choi3 and Chae Young Lim2

1University of Houston
2Seoul National University
3Yonsei University
4National Dong Hwa University

> Abstract: A biomechanical model often requires parameter estimation and
> selection in a known but complicated nonlinear function. Motivated by
> observing that data from a head-neck position tracking system, one of
> biomechanical models, show multiplicative time dependent errors, we develop
> a modified penalized weighted least squares estimator. The proposed method
> can be also applied to a model with non-zero mean time dependent additive
> errors. Asymptotic properties of the proposed estimator are investigated
> under mild conditions on a weight matrix and the error process. A simulation
> study demonstrates that the proposed estimation works well in both parameter
> estimation and selection with time dependent error. The analysis and
> comparison with an existing method for head-neck position tracking data show
> better performance of the proposed method in terms of the variance accounted
> for (VAF).
>
> Key words and phrases: nonlinear regression; temporal dependence;
> multiplicative error; local consistency and oracle property

##  1 Introduction

A nonlinear regression model has been widely used to describe complicated
relationships between variables (Wood, 2010; Baker and Foley, 2011; Paula et
al., 2015; Lim et al., 2014). In particular, various nonlinear problems are
considered in the field of machinery and biomechanical engineering (Moon et
al., 2012; Santos and Barreto, 2017). Among such nonlinear problems, a head-
neck position tracking model with neurophysiological parameters in
biomechanics motivated us to develop an estimation and selection method for a
nonlinear regression model in this work.

The head-neck position tracking application aims to figure out how
characteristics of the vestibulocollic and cervicocollic reflexes (VCR and
CCR) contribute to the head-neck system. The VCR activates neck muscles to
stabilize the head-in-space and the CCR acts to hold the head on the trunk. A
subject of the experiment follows a reference signal on a computer screen with
his or her head and a head rotation angle is measured during the experiment. A
reference signal is the input of the system and the measured head rotation
angle is the output. The parameters related to VCR and CCR in this nonlinear
system are of interest to understand the head-neck position tracking system.

The head-neck position problem has been widely studied in the literature of
biomechanics (Peng et al., 1996; Chen et al., 2002; Forbes et al., 2013;
Ramadan et al., 2018; Yoon et al., 2022). One of the prevalent issues in
biomechanics is that a model suffers from a relatively large number of
parameters and limited availability of data because the subjects in the
experiment cannot tolerate sufficient time without being fatigued. This leads
to overfitting as well as nonidentifiability of the parameters. To resolve
this issue, selection approaches via a penalized regression method have been
implemented to fix a subset of the parameters to the pre-specified values
while the remaining parameters are estimated (Ramadan et al., 2018; Yoon et
al., 2022).

Figure 1: The black curve represents the measured responses (the observations)
from the subject No. 8 in the head-neck position tracking experiment. The red
dashed curve represents the estimated responses (the fitted values) from the
nonlinear regression model with additive errors introduced in Yoon et al.
(2022).

Figure 2: Sample autocorrelation (left) and sample partial autocorrelation
(right) of the residuals $(\mbox{measured response}-\mbox{estimated
response})$ for the subject No. 8 from the head-neck position tracking
experiment. The estimated response is obtained from the method in Yoon et al.
(2022).

The existing approaches, however, have some limitations. The fitted values
from the penalized nonlinear regression with additive errors in Yoon et al.
(2022) show larger discrepancy from the observed values when the head is
turning in its direction. For example, Figure 2 shows the fitted values (the
estimated responses) from the method in Yoon et al. (2022) and the
observations (the measured responses) of the subject No. 8 in a head-neck
position tracking experiment. The detail description of the data from the
experiment is given in Section 4. We can see the tendency that the larger the
measured response in its magnitude is, the bigger the gap exists between the
measured response and the estimated response. Note that the largest measured
responses occur when the head is turning in its direction to follow the
reference signal’s direction changes. It is more difficult for some subjects
to track the end positions correctly, which can create more errors at each
end. Hence, the multiplicative errors rather than the additive errors would be
more suitable. Indeed, Figure 2 shows the model with additive errors did not
successfully accommodate this characteristic.

Previous studies did not also take into account temporal dependence while the
experimental data exhibit temporal correlation. For example, Figure 2 shows
sample autocorrelation function and sample partial autocorrelation function of
the residuals from the fitted model for the subject No. 8 by the method in
Yoon et al. (2022). It clearly shows temporal dependence in the residuals.
However, Yoon et al. (2022) studied a penalized nonlinear least squares
estimator and its asymptotic properties under the independent error
assumption. Lastly, Ramadan et al. (2018) and Yoon et al. (2022) restrict the
number of sensitive parameters to five, where the sensitive parameters refer
to the parameters whose estimates are not shrunk to the pre-specified values.
Not only may this restriction increase computational instability but also
reduce estimation and prediction performances. Provided the head-neck position
tracking task already suffers from computational challenges, additional
computational issues should be avoided.

To resolve the above-mentioned issues, we consider a nonlinear regression
model with multiplicative errors for the head-neck position tracking system,
which can be written as

$z_{t}=g(\bm{x}_{t};\bm{\theta})\times\varsigma_{t},$ (1.1)

where $z_{t}$ is a measured response, $g(\bm{x};\bm{\theta})$ is a known
nonlinear function with an input $\bm{x}$. $\bm{\theta}$ is a set of
parameters in $g$ and $\varsigma_{t}$ is a multiplicative error. The details
for $g$ and $\bm{\theta}$ for the head-neck position tracking system are
described in Ramadan et al. (2018) and Yoon et al. (2022). A typical approach
is to take the logarithm in both sides of (1.1) so that the resulting model
becomes a nonlinear model with additive errors:

$y_{t}=f(\bm{x}_{t};\bm{\theta})+\epsilon_{t},$ (1.2)

with $y_{t}=\log(z_{t})$,
$f(\bm{x}_{t};\bm{\theta})=\log(g(\bm{x}_{t};\bm{\theta}))$ and
$\epsilon_{t}=\log(\varsigma_{t})$ by assuming all components are positive.
Note that the additive error $\epsilon_{t}$ in (1.2) may not have zero-mean
anymore due to the log transformation.

Motivated by the head-neck position tracking application, we propose a
parameter estimation method for a nonlinear model with multiplicative error in
(1.1) by applying a modified weighted least squares estimation method to (1.2)
which accommodate temporal dependence as well as non-zero mean errors. Given
by the structure, the proposed estimation can also handle the nonlinear
regression model with non-zero mean additive errors in (1.2). The asymptotic
properties of the estimator obtained from the proposed method are studied
under the assumption of temporally correlated errors as we observed temporal
dependence in the head-neck position tracking system data.

A penalized estimator and its asymptotic properties are also investigated. In
the application of the head-neck position tracking system, a set of parameters
needs to be shrunk to the pre-specified values instead of the zero-values.
Thus, we use the penalized estimation approach to shrink estimates to the pre-
specified values. By doing so, we not only resolve the non-identifiability
issue but also keep all estimates meaningful in the head-neck position
tracking system. While the application of the head-neck position tracking
system motivated us to develop the proposed method, it is applicable to
general nonlinear regression models as well.

Least squares estimators and their theoretical properties for the nonlinear
regression models have been studied in the literature. Jennrich (1969) first
proved the strong consistency and asymptotic normality of the nonlinear least
squares estimator with independent errors. Wu (1981) provided a necessary
condition for the existence of any kind of weakly consistent estimator in a
nonlinear regression model with an additive error. In the same article, the
asymptotic properties of the least squares estimator were proved under weaker
conditions than those in Jennrich (1969). The results by Wu (1981) were
generalized in several other papers. Van de Geer (1990) studied three
estimation methods for a general regression model: least squares, least
absolute deviations, and penalized least squares. A stochastic nonlinear
function and martingale difference errors were considered in Lai (1994).
Pollard and Radchenko (2006) proved asymptotic properties of a least squares
estimator for a nonlinear regression model under the second moment assumptions
on the errors.

There are a few studies on the nonlinear regression with multiplicative
errors. Xu and Shimada (2000) studied least squares estimation for nonlinear
multiplicative noise models with independent errors. A weighted least squares
estimation method was proposed in Xu and Shimada (2000), but the estimator
induces a bias and needs correction. Lim et al. (2014) also investigated the
nonlinear multiplicative noise models with independent errors by the log
transformation and proposed the modified least square estimation by including
a sample mean component in the objective function. We also take a similar
approach but the least squares are modified with weights and temporally
dependent errors are assumed in addition to the penalization. Bhattacharyya et
al. (1992) showed that an ordinary least squares estimator for a nonlinear
regression model with additive errors may not possess strong consistency when
the true underlying model is a nonlinear regression model with multiplicative
errors.

Parameter selection via a penalty function has attracted a great deal of
attention since Tibshirani (1996) introduced the least absolute shrinkage and
selection operator (LASSO) and Fan and Li (2001) developed non-concave
penalized likelihood approach with the smoothly clipped absolute deviation
(SCAD). In particular, Fan and Li (2001) established asymptotic properties,
especially the oracle property for the penalized likelihood estimator under
mild regularity conditions. For a nonlinear regression model with independent
errors, asymptotic properties of the estimator from the penalized estimation
method were investigated in Jiang et al. (2012), Wu et al. (2014), and Lv et
al. (2014). Jiang et al. (2012) proposed a penalized weighted composite
quantile regression estimator and developed its asymptotic properties. They
highlighted the proposed method works as efficient as the oracle maximum
likelihood estimator with various error distributions and even works for
heavy-tailed error distributions. Wu et al. (2014) studied nonlinear
independence screening and nonnegative garrote for a high-dimensional
nonlinear additive model. Lv et al. (2014) investigated a nonlinear modal
regression model with the increasing number of variables. Asymptotic
properties such as oracle property and finite sample properties were also
studied in the paper.

To deal with a temporal dependence in the errors, mixing conditions are
considered: strong mixing ($\alpha$-mixing), $\phi$-mixing, and $\rho$-mixing.
These mixing conditions have been frequently studied for handling dependence
of temporal data (Machkouri et al., 2017; Geller and Neumann, 2018). Various
regression models with strong-mixing errors have been previously investigated
(Zhang and Liang, 2012; Yang et al., 2017). Zhang and Liang (2012) studied a
semi-parametric regression model with strong mixing errors and established the
asymptotic normality of a least squares estimator and a weighted least squares
estimator. Yang et al. (2017) developed probability inequalities about a least
squares estimator of nonlinear regression with strong mixing errors.

In this work, we propose a modified weighted least squares estimator for a
nonlinear regression model with non-zero mean errors, which is temporally
correlated. As we mentioned before, this approach can handle a nonlinear
regression model with multiplicative dependence errors after log
transformation. We investigate the asymptotic properties of the proposed
estimator and its penalized version. Specifically, we establish the local
consistency and asymptotic normality for both estimators and the oracle
property of the penalized estimator.

For a penalty function, we consider LASSO and SCAD in the numerical study. The
performance in the numerical study indicates the proposed penalized estimator
works well with finite samples. At last, we exhibit the results of adopting
the proposed method to the head-neck position tracking data with comparison to
those from the existing method.

In Section 2, we demonstrate the proposed estimation method for a nonlinear
regression model and establish the asymptotic properties of the proposed
estimators. In Section 3, several simulation studies are conducted with
various settings. The analysis on head-neck position tracking data with the
proposed method is introduced in Section 4. At last, we provide a discussion
in Section 5. The technical proofs for the theorems and additional results of
the simulation studies are presented in a supplementary material.

##  2 Methods

###  2.1 Modified Weighted Least Squares

We consider a following nonlinear regression model

$y_{t}=f(\bm{x}_{t};\bm{\theta})+\epsilon_{t},$

for $t=1,\cdots,n$, where $f({\bm{x}};\bm{\theta})$ is a known nonlinear
function on ${\bm{x}}\in D\subset R^{d}$, which also depends on the parameter
vector $\bm{\theta}:=(\theta_{1},\theta_{2},...,\theta_{p})^{T}\in\Theta$.
${A}^{T}$ is the transpose of a matrix $A$. $\epsilon_{t}$ is temporally
correlated and its mean, $E(\epsilon_{t})=\mu$, may not be zero. The
assumption on a non-zero mean of the error comes from a nonlinear regression
model with multiplicative errors introduced in (1.1). We further assume that
only a few entries of the true parameter are non-zero. Without loss of
generality, we let the first $s$ entries of the true $\bm{\theta}_{0}$ be non-
zero. That is,
$\bm{\theta}_{0}=(\theta_{01},\theta_{02},...,\theta_{0s},\theta_{0s+1},...,\theta_{0p})^{T}$
and $\theta_{0t}\neq 0$ for $1\leq t\leq s$ and $\theta_{0t}=0$ for $s+1\leq
t\leq p$.

We begin with a modified least squares method using

$\displaystyle\displaystyle\bm{S}_{n}^{(ind)}(\bm{\theta})$
$\displaystyle=\sum^{n}_{t=1}\left(y_{t}-f({\bm{x}_{t}};\bm{\theta})-\frac{1}{n}\sum_{t^{\prime}=1}^{n}(y_{t^{\prime}}-f(\bm{x}_{t^{\prime}};\bm{\theta}))\right)^{2},$
$\displaystyle=\left(\bm{y}-\bm{f}(\bm{x},\bm{\theta})\right)^{T}\bm{\Sigma}_{n}\left(\bm{y}-\bm{f}(\bm{x},\bm{\theta})\right),$

where $\bm{y}=(y_{1},\ldots,y_{n})^{T},\
\bm{f}(\bm{x},\bm{\theta})=(f(\bm{x}_{1};\bm{\theta}),\ldots,f(\bm{x}_{n};\bm{\theta}))^{T}$,
and $\bm{\Sigma}_{n}=\bm{I}_{n}-n^{-1}\bm{1}\bm{1}^{T}$. $\bm{I}_{n}$ is the
identity matrix and $\bm{1}$ is the column vector of one’s. This objective
function is different from a typical least squares expression in that the
sample mean of the errors is subtracted from the error at each data point.
This is motivated by taking into account possibly non-zero mean of the errors
(Lim et al., 2014).

Since we consider temporally dependent data, we introduce a temporal weight
matrix $\bm{W}$ to account for the temporal dependence so that the modified
objective function is

$\displaystyle\bm{S}_{n}(\bm{\theta})$
$\displaystyle=\left(\bm{y}-\bm{f}(\bm{x},\bm{\theta})\right)^{T}\bm{W}^{T}\bm{\Sigma}_{n}\bm{W}\left(\bm{y}-\bm{f}(\bm{x},\bm{\theta})\right),$
$\displaystyle=\left(\bm{y}-\bm{f}(\bm{x},\bm{\theta})\right)^{T}\bm{\Sigma}_{w}\left(\bm{y}-\bm{f}(\bm{x},\bm{\theta})\right),$
(2.3)

where $\bm{\Sigma}_{w}=\bm{W}^{T}\bm{\Sigma}_{n}\bm{W}$. If we know the true
temporal dependence model of the error process, a natural choice of the weight
matrix is from the true covariance matrix of the error process. However, we
allow the weight matrix more flexible since we do not want to assume a
specific temporal dependence model for the error process. Conditions for
$\bm{W}$ will be introduced in the next section.

We can add a penalty function $p_{\tau_{n}}(\cdot)$ when our interest is to
detect the relevant parameters. Then, the penalized estimator is obtained by
minimizing

$\bm{Q}_{n}(\bm{\theta})=\bm{S}_{n}(\bm{\theta})+n\sum^{p}_{i=1}p_{\tau_{n}}(|\theta_{i}|).$
(2.4)

To investigate theoretical properties of the proposed estimators obtained by
minimizing $\bm{S}_{n}(\bm{\theta})$ and $\bm{Q}_{n}(\bm{\theta})$, we
introduce notations and assumptions in the next section.

###  2.2 Notations and Assumptions

We start with three mixing conditions for temporal dependence:
$\alpha$-mixing, $\phi$-mixing, and $\rho$-mixing.

######  Definition 1.

  1. Consider a sequence of random variables, $\\{\xi_{i},i\geq 1\\}$ and let $\mathcal{F}_{i}^{j}$ denote the $\sigma$-field generated by $\\{\xi_{i},\ldots,\xi_{j}\\}$. Then, $\\{\xi_{i},i\geq 1\\}$ is said to be

  2. (a)

strong mixing or $\alpha$-mixing if $\alpha(m)\rightarrow 0$ as
$m\rightarrow\infty$, where

$\displaystyle\alpha(m)$
$\displaystyle=\displaystyle\underset{n}{\sup}~{}\alpha(\mathcal{F}_{1}^{n},\mathcal{F}_{n+m}^{\infty})$
$\displaystyle\hbox{with}~{}~{}\alpha(\mathcal{F},\mathcal{G})$
$\displaystyle=\displaystyle\underset{A\in\mathcal{F},B\in\mathcal{G}}{\sup}~{}|P(AB)-P(A)P(B)|,$

  3. (b)

$\phi$-mixing if $\phi(m)\rightarrow 0$ as $m\rightarrow\infty$, where

$\displaystyle\phi(m)$
$\displaystyle=\displaystyle\underset{n}{\sup}~{}\phi(\mathcal{F}_{1}^{n},\mathcal{F}_{n+m}^{\infty})$
$\displaystyle\hbox{with}~{}~{}\phi(\mathcal{F},\mathcal{G})$
$\displaystyle=\displaystyle\sup_{A\in\mathcal{F},B\in\mathcal{G},P(A)>0}|P(B|A)-P(B)|,~{}~{}\hbox{and}$

  4. (c)

$\rho$-mixing if $\rho(m)\rightarrow 0$ as $m\rightarrow\infty$, where

$\displaystyle\rho(m)$
$\displaystyle=\displaystyle\underset{n}{\sup}~{}\rho(\mathcal{F}_{1}^{n},\mathcal{F}_{n+m}^{\infty}),$
$\displaystyle\hbox{with}~{}~{}\rho(\mathcal{F},\mathcal{G})$
$\displaystyle=\displaystyle\sup_{f\in\mathcal{L}^{2}_{real}(\mathcal{F}),\
g\in\mathcal{L}^{2}_{real}(\mathcal{G})}|corr(f,g)|.$

Here, $L^{2}_{real}(\mathcal{F})$ denotes the space of square-integrable,
$\mathcal{F}$-measurable real-valued random variables (Bradley, 2005).

These mixing conditions have been widely adopted to explain dependence of a
random sequence in the literature (Machkouri et al., 2017; Geller and Neumann,
2018). It is well-known that $\phi$-mixing implies $\rho$-mixing,
$\rho$-mixing implies $\alpha$-mixing, and the strong mixing condition is one
of the weakest conditions among many mixing conditions (Peligrad and Utev,
1997; Bradley, 2005). We assume an appropriate mixing condition for our
temporal data and derive the asymptotic properties of the proposed estimators
under such condition. The details appear in Assumption 1.

Next, we introduce notations and assumptions for theoretical results. Define
$d_{t}(\bm{\theta},\bm{\theta}^{\prime})=f({\bm{x}}_{t};\bm{\theta})-f({\bm{x}}_{t};\bm{\theta}^{\prime})$
and
$\bm{d}=(d_{1}(\bm{\theta},\bm{\theta}^{\prime}),d_{2}(\bm{\theta},\bm{\theta}^{\prime}),\ldots,d_{n}(\bm{\theta},\bm{\theta}^{\prime}))^{T}$.
When $f$ is twice differentiable with respect to $\bm{\theta}$, let
$\bm{f}_{k}=\left(\frac{\partial{f(\bm{x}_{1},\bm{\theta})}}{\partial\theta_{k}},\cdots,\frac{\partial
f(\bm{x}_{n},\bm{\theta})}{\partial\theta_{k}}\right)^{T}$ and
$\bm{f}_{kl}=\left(\frac{\partial^{2}\bm{f}(\bm{x}_{1},\bm{\theta})}{\partial\theta_{k}\partial\theta_{l}},...,\frac{\partial^{2}\bm{f}(\bm{x}_{n},\bm{\theta})}{\partial\theta_{k}\partial\theta_{l}}\right)^{T}$.
Using $\bm{f}_{k}$ and $\bm{f}_{kl}$, we define
${\bm{\dot{F}}}(\bm{\theta})=(\bm{f}_{1},...,\bm{f}_{p})$ and
$\bm{{\ddot{F}}}(\bm{\theta})$=Block($\bm{f}_{kl}$) so that
$\bm{\dot{F}}(\bm{\theta})$ is $n\times p$ matrix whose $k$th column is
$\bm{f}_{k}$ and $\bm{{\ddot{F}}}$ is a $pn\times p$ block matrix whose
$(k,l)$th block is $\bm{f}_{kl}$. $\mbox{E}(\bm{\epsilon})=\mu\bm{1}$ and
$\operatorname{var}(\bm{\epsilon})=\bm{\Sigma}_{\epsilon}$. Let $\lambda_{w}$
and $\lambda_{\epsilon}$ denote the maximum eigenvalues of $\bm{\Sigma}_{w}$
and $\bm{\Sigma}_{\epsilon}$, respectively. We consider a temporal weight
matrix satisfying $\bm{W}\bm{1}=\bm{1}$, i.e. the row-sums are 1’s. This
condition is to handle the nonzero mean of the errors. Let $\|\cdot\|$ for a
vector denote a euclidean norm and $\|\cdot\|_{1}$ and $\|\cdot\|_{\infty}$
for a matrix denote 1-norm and infinity norm, respectively.

The assumptions on the nonlinear function $f$, the errors $\epsilon_{i}$, the
weight matrix $\bm{W}$ and the penalty function $p_{\tau_{n}}(\cdot)$ to
investigate asymptotic properties are now introduced.

######  Assumption 1.

  1. (1)

The nonlinear function $f\in C^{2}$ on the compact set
$\mathcal{D}\times\Theta$ where $C^{2}$ is the set of twice continuously
differentiable functions.

  2. (2)

As $\|\bm{\theta}-\bm{\theta}_{0}\|\rightarrow 0$,
$\left({\bm{\dot{F}}}(\bm{\theta}_{0})^{T}\bm{\Sigma}_{w}{\bm{\dot{F}}}(\bm{\theta}_{0})\right)^{-1}{\bm{\dot{F}}}(\bm{\theta})^{T}\bm{\Sigma}_{w}{\bm{\dot{F}}}(\bm{\theta})\rightarrow
I_{p}$, elementwisely and uniformly in $\bm{\theta}$.

  3. (3)

There exist symmetric positive definite matrices $\bm{\Gamma}$ and
$\bm{\Gamma}_{\epsilon}$ such that

$\displaystyle\frac{1}{n\lambda_{w}}\bm{\dot{F}}(\bm{\theta}_{0})^{T}\bm{\Sigma}_{w}\bm{\dot{F}}(\bm{\theta}_{0})$
$\displaystyle\rightarrow\bm{\Gamma}$
$\displaystyle\frac{1}{n\lambda_{\epsilon}\lambda_{w}^{2}}{\bm{\dot{F}}}(\bm{\theta}_{0})^{T}\bm{\Sigma}_{w}\bm{\Sigma}_{\epsilon}\bm{\Sigma}_{w}{\bm{\dot{F}}}(\bm{\theta}_{0})$
$\displaystyle\rightarrow\bm{\Gamma}_{\epsilon}.$

  4. (4)

$\frac{\|\bm{W}\|_{1}\cdot\|\bm{W}\|_{\infty}}{\|\bm{W}^{T}\bm{\Sigma}_{n}\bm{W}\|_{2}}=o(n^{1/2}\lambda_{\epsilon}^{1/2})$.

  5. (5)

$O(1)\leq\lambda_{\epsilon}\leq o(n)$ and $\lambda_{w}\geq O(1)$.

  6. (6)

$\\{\epsilon_{i}^{2}\\}$ is uniformly integrable.

  7. (7)

One of the following conditions is satisfied for $\epsilon_{i}$.

     * $(a)$

$\\{\epsilon_{i}\\}$ is a $\phi$-mixing.

     * $(b)$

$\\{\epsilon_{i}\\}$ is a $\rho$-mixing and
$\sum_{j\in\mathcal{N}}\rho(2^{j})<\infty$.

     * $(c)$

For $\delta>0$, $\\{\epsilon_{i}\\}$ is a $\alpha$-mixing,
$\\{|\epsilon_{i}|^{2+\delta}\\}$ is uniformly integrable, and
$\sum_{j\in\mathcal{N}}n^{2/\delta}\alpha(n)<\infty$.

######  Assumption 2.

The first derivative of a penalty function $p_{\tau_{n}}(\cdot)$ denoted by
$q_{\tau_{n}}(\cdot)$, has the following properties:

  1. (1)

$c_{n}=\max_{i\in\\{1,\ldots,s\\}}\left\\{|q_{\tau_{n}}(|\theta_{0i}|)|\right\\}=O\left(\left(\lambda_{\epsilon}/n\right)^{1/2}\right)$

  2. (2)

$q_{\tau_{n}}(\cdot)$ is Lipschitz continuous given $\tau_{n}$

  3. (3)

$n^{1/2}\lambda_{\epsilon}^{-1/2}\lambda_{w}^{-1}\tau_{n}\rightarrow\infty$

  4. (4)

For any $C>0$,
$\displaystyle\liminf_{n\rightarrow\infty}\inf_{\theta\in\left(0,C(\lambda_{\epsilon}/n)^{1/2}\right)}\tau_{n}^{-1}q_{\tau_{n}}(\theta)>0$

Assumption 1 imposes mild conditions on the nonlinear function, its domain,
the weight matrix, and the error process. Assumption 1-(1), (2), and (3)
introduce reasonably weak conditions for the nonlinear function and the domain
of data and parameters. The first condition in Assumption 1-(3) is a modified
version of Grenander condition for our objective function (Grenander, 1954;
Wang and Zhu, 2009; Lim et al., 2014). The second condition in Assumption
1-(3) is required to derive the variance of the asymptotic distribution.
Remark 1 explains the plausibility of these conditions by addressing that
slightly weaker conditions can be easily satisfied. We impose a weak condition
on the temporal weight matrix in Assumption 1-(4) so that flexible weight
matrices are allowed. Remark 2 further discusses on the condition for the
temporal weight matrix. In Assumption 1-(5), a lower bound for
$\lambda_{\epsilon}$ can be attained if the error process is stationary with a
bounded spectral density. Assumption 1 contains additional conditions for the
asymptotic normality of the unpenalized estimator from
$\bm{S}_{n}(\bm{\theta})$. Assumptions 1-(6) and (7) refer to Peligrad and
Utev (1997), which studied central limit theorems for linear processes.
Assumption 1-(6) implies uniform boundedness of the second moment for the
errors. Assumption 1-(7) provides weak conditions for temporal dependence of
the errors. The detailed discussion on Assumption 1-(7) is given in Remark 3.

Assumption 2 demonstrates typical conditions for a penalty function. The first
two conditions guarantee that the penalized least squares estimator possess
consistency with the same order as the modified weighted least squares
estimator. The other two conditions contribute to the sparsity of the
penalized estimator. The conditions in Assumption 2 are similar to those in
Fan and Li (2001) and Wang and Zhu (2009). Typically, LASSO and SCAD penalty
functions are considered. The former satisfies only the first two conditions
in Assumption 2 while the latter satisfies all conditions in Assumption 2 with
properly chosen $\tau_{n}$. This means LASSO fails to correctly identify
significant parameters while an estimator using the SCAD penalty function
possesses selection consistency as well as estimation consistency.

######  Remark 1.

We discuss the positive definiteness of $\bm{\Gamma}$ and
$\bm{\Gamma}_{\epsilon}$ and the boundedness of the sequences of the matrices.
First, $\bm{\Sigma}_{w}$ is a symmetric and semi-positive definite matrix
since $\bm{\Sigma}_{w}=\bm{W}^{T}\bm{\Sigma}_{n}\bm{W}$ and $\bm{\Sigma}_{n}$
has rank of $n-1$. Despite the rank deficiency, the sequences of the matrices
are $p\times p$ matrices with $p<n$, so we believe that the limits of the
sequences are likely to acquire positive definiteness. Next, the first
sequence of the matrices in Assumption 1-(3) are clearly bounded above. Since
$\bm{\Sigma}_{w}$ is a semi-positive definite matrix,
$(n\lambda_{w})^{-1}\bm{\dot{F}}(\bm{\theta}_{0})^{T}\bm{\Sigma}_{w}\bm{\dot{F}}(\bm{\theta}_{0})\leq
n^{-1}\bm{\dot{F}}(\bm{\theta}_{0})^{T}\bm{\dot{F}}(\bm{\theta}_{0})=O(1)$ by
the compactness of the domain (Assumption 1-(1)). With
$\lambda_{max}(\bm{\Sigma}_{w}\bm{\Sigma}_{\epsilon}\bm{\Sigma}_{w})\leq\lambda_{\epsilon}\lambda_{w}^{2}$,
we obtain the same result for the second sequence.

######  Remark 2.

We give detailed justification for assumptions on the temporal weight matrix.
By H$\ddot{\mbox{o}}$lder’s inequality,
$\|\bm{W}\|_{2}^{2}\leq\|\bm{W}\|_{1}\|\bm{W}\|_{\infty}$. Hence, with
$\lambda_{w}=\|\bm{W}^{T}\bm{\Sigma}_{n}\bm{W}\|_{2}$ Assumption 1-(4) leads
to $\|\bm{W}\|_{2}\leq o(n^{1/4}\lambda_{\epsilon}^{1/4}\lambda_{w}^{1/2})$.
Recall $O(1)\leq\lambda_{\epsilon}\leq o(n)$ and $\lambda_{w}\geq O(1)$ from
Assumption 1-(5). Thus, the upper bound is sufficiently large for
$\|\bm{W}\|_{2}$ to allow flexible $\bm{W}$. In addition, since product
matrices $\bm{A}\bm{B}$ and $\bm{B}\bm{A}$ share their eigenvalues,
$\lambda_{w}=\lambda_{max}(\bm{W}^{T}\bm{\Sigma}_{n}\bm{W})=\lambda_{max}(\bm{\Sigma}_{n}\bm{W}\bm{W}^{T})\leq\|\bm{W}\|_{2}^{2}$.
In summary, we obtain $O(1)\leq\|\bm{W}\|_{2}^{2}\leq
o(n^{1/2}\lambda_{\epsilon}^{1/2}\lambda_{w})$, so Assumptions 1-(4) and (5)
together provide a flexible upperbound and lowerbound for $\bm{W}$.

######  Remark 3.

There exist many familiar time series processes that satisfy Assumption 1-(7).
Autoregressive (AR) processes and moving average (MA) processes are strongly
mixing under mild conditions (Athreya and Pantula, 1986). Athreya and Pantula
(1986) also mentions that finite order autoregressive moving average (ARMA)
processes are $\phi$-mixing under mild conditions. Furthermore, ARMA processes
and bilinear processes are strong mixing with $\alpha(n)=O(e^{-n\rho})$ with
some $\rho>0$ (Roussas et al., 1992).

###  2.3 Theoretical results

First, we construct the existence and the consistency of the modified weighted
least squares estimator and the penalized least squares estimator.

######  Theorem 1.

For any $\varepsilon>0$ and $a_{n}=(\lambda_{\epsilon}/n)^{1/2}$, under
Assumption 1-(1), (2), (3), and (5), there exists a positive constant $C$ such
that

$P\left(\inf_{\|\bm{v}\|=C}\bm{S}_{n}(\bm{\theta}_{0}+a_{n}\bm{v})-\bm{S}_{n}(\bm{\theta}_{0})>0\right)>1-\varepsilon$

for large enough $n$. Therefore, with probability tending to 1, there exists a
local minimizer of $\bm{S}_{n}(\bm{\theta})$, say $\hat{\bm{\theta}}^{(s)}$,
in the ball centered at $\bm{\theta}_{0}$ with the radius $a_{n}\bm{v}$. Since
$a_{n}=o(1)$ by Assumption 1-(5), we have the consistency of
$\hat{\bm{\theta}}^{(s)}$.

######  Theorem 2.

For any $\varepsilon>0$ and $b_{n}=(\lambda_{\epsilon}/n)^{1/2}+c_{n}$, under
Assumptions in theorem 1 and 2-(1),(2), there exists a positive constant $C$
such that

$P\left(\inf_{\|\bm{v}\|=C}\bm{Q}_{n}(\bm{\theta}_{0}+b_{n}\bm{v})-\bm{Q}_{n}(\bm{\theta}_{0})>0\right)>1-\varepsilon$

for large enough $n$. Therefore, with probability tending to 1, there exists a
local minimizer of $\bm{Q}_{n}(\bm{\theta})$, say $\hat{\bm{\theta}}$, in the
ball centered at $\bm{\theta}_{0}$ with the radius $b_{n}\bm{v}$. By
Assumptions 1-(5) and 2-(1), $b_{n}=o(1)$, which leads to the consistency of
$\hat{\bm{\theta}}$.

The next two theorems establish the asymptotic normality of the modified
weighted least squares estimator from Theorem 1 and the oracle property of the
penalized least squares estimator from Theorem 2.

######  Theorem 3 (Asymptotic normality).

Under Assumption 1,

$\left(\frac{n}{\lambda_{\epsilon}}\right)^{1/2}\left(\hat{\bf{\bm{\theta}}}^{(s)}-\bm{\theta}_{0}\right)~{}\overset{d}{\longrightarrow}~{}N\left(0,\bm{\Gamma}^{-1}\bm{\Gamma}_{\epsilon}\bm{\Gamma}^{-1}\right),$

where $\hat{\bm{\theta}}^{(s)}$ is a consistent estimator introduced in
Theorem 1 using $\bm{S}_{n}(\bm{\theta})$.

######  Theorem 4 (Oracle property).

With $\hat{\bm{\theta}}$, a consistent estimator introduced in Theorem 2 using
$\bm{Q}_{n}(\bm{\theta})$, if Assumptions 1 and 2 are satisfied,

  1. (i)

$P\left(\hat{\theta}_{i}=0\right)\rightarrow 1,$ for $i\in\\{s+1,\ldots,p\\}$.

  2. (ii)

Also,

$\left(\frac{n}{\lambda_{\epsilon}}\right)^{1/2}\left(\hat{\bm{\theta}}_{1}-\bm{\theta}_{01}+\left((2\lambda_{w}\bm{\Gamma})^{-1}\right)_{11}\bm{\beta}_{n,s}\right)~{}\overset{d}{\longrightarrow}~{}N\left(0,\left(\bm{\Gamma}^{-1}\bm{\Gamma}_{\epsilon}\bm{\Gamma}^{-1}\right)_{11}\right),$

where $\hat{\bm{\theta}}_{1}=(\hat{\theta}_{1},\ldots,\hat{\theta}_{s})^{T},\
\bm{\theta}_{01}=(\theta_{01},\ldots,\theta_{0s})^{T},\
\bm{\beta}_{n,s}=(q_{\tau_{n}}({|\theta_{01}|)sgn(\theta}_{01}),\ldots,q_{\tau_{n}}(|\theta_{0s}|)sgn(\theta_{0s}))^{T}$
and $\bm{A}_{11}$ is the $s\times s$ upper-left matrix of $\bm{A}$.

Note that the estimators from theorems 1 and 2 have convergence orders of
$a_{n}$ and $b_{n}$, respectively. Also, $a_{n}$ and $b_{n}$ eventually have
the same order of $(\lambda_{\epsilon}/n)^{1/2}$ by Assumption 2-(1). One may
think $\lambda_{w}$, information of $\bm{W}$, makes no contribution to both
theorems, even though we have $\bm{W}$ in the objective functions. Recall that
$\lambda_{w}$ contributes to $\tau_{n}$ via Assumption 2-(2) and (4), and
$c_{n}$, which appears in $b_{n}$, is related to $\tau_{n}$ by Assumption
2-(1). This is where $\lambda_{w}$ implicitly comes into the theorems. We
could impose different conditions to make $\lambda_{w}$ explicitly appear in
the theorems. However, such conditions restrict the flexibility of
$\bm{\Sigma}_{\epsilon}$ so that the applicability of the proposed methods
becomes limited. Instead, we decide to keep the current assumptions to allow
flexible $\bm{\Sigma}_{\epsilon}$ and attain the implicit involvement of
$\lambda_{w}$ in the theorems. The proofs for Theorems 1-4 are given in the
supplementary material.

##  3 Simulation Study

We investigate performance of the proposed estimator, in particular, the
penalized version with two different penalty functions, LASSO and SCAD using
simulated data sets. First, we consider the data generated from a nonlinear
additive error model:

$\displaystyle
y_{t}=\frac{1}{1+\exp(-\bm{x}_{t}^{T}\bm{\theta}_{0})}+\epsilon_{t},$ (3.5)

where $\bm{\theta}_{0}=(\theta_{01},\theta_{02},\ldots,\theta_{0,20})^{T}$
with $\theta_{01}=1,\theta_{02}=1.2,\theta_{03}=0.6$, and the others being
zero. The first component of the covariate $\bm{x}$ comes from $U[-1,1]$, a
uniform distribution on $[-1,1]$, and the other components of $\bm{x}$ are
simulated from a joint normal distribution with the zero mean, the variance
being 0.6 and pairwise covariance being 0.1. This is a slight modification of
the covariates setting used in Jiang et al. (2012). For $\epsilon_{t}$, the
AR(1) and ARMA(1,1) with the non-zero mean are considered since these
processes not only represent typical time series processes but also possess
the strong mixing property. The choices of the AR(1) coefficient ($\rho$) are
0.5 and 0.9. For the ARMA process, the parameters for the AR and MA parts are
fixed as 0.8 ($\rho$) and 0.4 ($\phi$), respectively. For the non-zero mean,
$\mu$, the choices are 0.1 and 0.5 and the same for the standard deviation,
$\sigma$. We only show the results when $\sigma=0.5$ to highlight findings as
it is more difficult settings due to a larger variance. We consider three
sample sizes; $n=50,100$ and $200$ to investigate improvements according to
the increasing sample sizes.

A coordinate descent (CD) algorithm, in particular, a cyclic CD algorithm
(Breheny and Huang, 2011), was implemented to calculate the minimizer of the
objective function given in (2.4). Although Breheny and Huang (2011)
considered the convergence of the CD algorithm in a linear model, the cyclic
CD algorithm worked well for our penalized nonlinear regression problem as
well.

To select the tuning parameter, $\tau_{n}$, of the penalty function, a BIC-
type criterion (Wang et al., 2007) was used. The tuning parameter $a$ in the
SCAD penalty was fixed at 3.7 as recommended in Fan and Li (2001). The BIC-
type criterion we consider is
$\mbox{BIC}=\log(\hat{\sigma}^{2})+\log(n)\cdot\widehat{df}/n,$ where
$\widehat{df}$ is the number of significant estimates. For $\hat{\sigma}^{2}$,
we used $\hat{\sigma}^{2}=\overline{r^{2}}-\bar{r}^{2}$ where
$\bar{r}=n^{-1}\sum_{t=1}^{n}r_{t}$ and
$\overline{r^{2}}=n^{-1}\sum_{t=1}^{n}r_{t}^{2}$ with
$r_{t}=y_{t}-f(\bm{x}_{t};\bm{\hat{\theta}})$.

Our proposed method is denoted as penalized modified weighted least squares
(PMWLS). For the simulation study, the square roots of the inverse of the
covariance matrices from the AR(1) process with the AR coefficient
$\rho=0.5,0.9$ and the ARMA(1,1) process with AR and MA coefficients
$(\rho,\phi)=(0.8,0.4)$ are considered for the weight matrices after scaling
to have the row-sums 1. We compare the results of the proposed PMWLS method
with the results from a penalized least squares with a weight matrix (PWLS).
For fair comparison with the proposed method, a temporal weight matrix is also
considered for the PWLS method. In the PWLS method, we introduce an intercept
term to account for a possible non-zero mean of the error in the equation
(3.5), which is a straightforward way to handle non-zero mean of the error.
That is, PWLS minimizes

$\left(\bm{y}-\beta_{0}-\bm{f}(\bm{x};\bm{\theta})\right)^{T}\widetilde{\bm{W}}(\bm{y}-\beta_{0}-\bm{f}(\bm{x};\bm{\theta}))+n\sum_{i=1}^{p}p_{\tau_{n}}(|\theta_{i}|),$

where $\beta_{0}$ is the intercept term and $\widetilde{\bm{W}}$ is a weight
matrix. For weight matrices, we used the inverse of the covariance matrices
from the same models considered for the PMWLS method such as AR(1) and
ARMA(1,1) processes.

Tables 1-3 report the values of mean squared error (MSE) with standard
deviation of squared error (SD) in parenthesis for the estimates with a SCAD
penalty from 100 repetitions of data. MSE and SD are calculated as

$\displaystyle MSE$
$\displaystyle=\frac{1}{R\,p}\sum_{j=1}^{R}\left\|\hat{\bm{\theta}}^{(j)}-\bm{\theta}_{0}\right\|^{2},$
$\displaystyle SD$
$\displaystyle=\displaystyle\sqrt{\frac{1}{R-1}\sum_{j=1}^{R}\left\|\hat{\bm{\theta}}^{(j)}-{\bm{\bar{\theta}}}\right\|^{2}},$

where $\hat{\bm{\theta}}^{(j)}$ stands for the estimate from the $j$-th
repetition and ${\bm{\bar{\theta}}}$ is the sample mean of
$\hat{\bm{\theta}}^{(j)}$ for $j=1,\ldots,100$. The results for the estimates
with a LASSO penalty are provided in the supplementary material.

The MSE and SD results in Tables 1-3 show good performances in estimating the
true parameters. Overall, the estimation results are robust over various
weight matrices. While the results with the LASSO penalty for both PMWLS and
PWLS are comparable (Tables S1-S3 in the supplementary material), the proposed
method (PMWLS) outperforms the PWLS method for most cases with the SCAD
penalty. In particular, the improvement by the PMWLS method with the SCAD
penalty is more apparent when the error process has stronger dependence
(Tables 2 and 3). This performance improvement would be from estimation
efficiency in finite sample since the PMWLS method estimates one less number
of parameters, an intercept, compared to the PWLS method. Weight matrices for
both PMWLS and PWLS methods help to reduce MSE and the improvement is more
clear when the dependence is strong. On the other hand, the choice of the
weight matrices do not make much difference except that the performance is
better when the weight matrix is introduced.

AR(1) with $\rho=0.5$
---
$(\mu,\sigma)$ | Methods | $n=50$ | $n=100$ | $n=200$
$(0.1,0.5)$ | PMWLS | 15.34 (1.67) | 4.98 (0.97) | 2.38 (0.66)
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 11.36 (1.38) | 4.17 (0.85) | 2.36 (0.64)
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 10.69 (1.30) | 4.22 (0.84) | 2.22 (0.61)
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 10.94 (1.29) | 4.42 (0.86) | 2.49 (0.64)
PWLS | 15.00 (1.64) | 4.89 (0.96) | 2.51 (0.67)
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 11.88 (1.40) | 4.22 (0.85) | 2.29 (0.63)
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 11.15 (1.34) | 4.09 (0.83) | 2.42 (0.63)
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 11.07 (1.33) | 4.53 (0.86) | 2.38 (0.63)
$(0.5,0.5)$ | PMWLS | 9.44 (1.32) | 5.15 (0.98) | 2.60 (0.70)
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 9.84 (1.31) | 4.19 (0.84) | 2.02 (0.59)
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 9.58 (1.27) | 4.07 (0.82) | 1.97 (0.58)
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 11.71 (1.39) | 4.73 (0.87) | 2.24 (0.62)
PWLS | 10.17 (1.34) | 6.78 (1.10) | 3.81 (0.83)
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 9.17 (1.24) | 4.20 (0.84) | 2.17 (0.60)
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 10.66 (1.34) | 4.37 (0.83) | 2.16 (0.59)
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 12.17 (1.43) | 5.10 (0.90) | 2.67 (0.66)
$\ast$ The actual MSE values are $0.01\times$ the reported values.
Table 1: Estimation results with SCAD for the equation (3.5) when the error
process is AR(1) with the AR coefficient $\rho=0.5$. Mean squared error (MSE)
with standard deviation (SD) values in the parenthesis are presented. PMWLS
and PWLS stand for penalized modified weighted least squares (our proposed
method) and penalized least squares with a weight matrix, respectively. In the
Methods column, the value of $\rho$ indicates that the weight matrix from the
AR(1) process with $\rho$ as the AR coefficient is considered for estimation.
The value of $(\rho,\phi)$ indicates that the weight matrix from the ARMA(1,1)
process with $\rho$ and $\phi$ as AR and MA coefficients are considered for
estimation. The rows with no $\rho$ or $(\rho,\phi)$ indicate that no weight
matrix is used. Additionally, $\mu$ and $\sigma$ are the mean and standard
deviation values of the error process. Three samples sizes $(n=50,100,\mbox{
and }200)$ are considered to verify sampling properties of the estimators.
AR(1) with $\rho=0.9$
---
$(\mu,\sigma)$ | Methods | $n=50$ | $n=100$ | $n=200$
$(0.1,0.5)$ | PMWLS | 9.42 (1.31) | 4.14 (0.85) | 2.49 (0.68)
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 4.90 (0.83) | 3.65 (0.69) | 1.37 (0.47)
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 4.96 (0.79) | 3.76 (0.67) | 1.49 (0.47)
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 5.29 (0.82) | 3.76 (0.69) | 1.42 (0.46)
PWLS | 14.59 (1.64) | 4.58 (0.89) | 2.47 (0.68)
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 5.14 (0.86) | 2.93 (0.66) | 1.41 (0.48)
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 4.95 (0.77) | 3.56 (0.65) | 1.41 (0.46)
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 5.08 (0.80) | 3.48 (0.65) | 1.44 (0.47)
$(0.5,0.5)$ | PMWLS | 9.38 (1.35) | 4.51 (0.91) | 2.33 (0.67)
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 5.13 (0.84) | 2.95 (0.65) | 1.42 (0.47)
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 5.19 (0.82) | 2.84 (0.63) | 1.41 (0.47)
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 5.76 (0.84) | 3.22 (0.67) | 1.32 (0.45)
PWLS | 12.02 (1.53) | 7.19 (1.17) | 2.93 (0.72)
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 5.53 (0.90) | 3.31 (0.68) | 1.65 (0.50)
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 5.25 (0.82) | 3.36 (0.67) | 1.69 (0.50)
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 5.75 (0.83) | 3.38 (0.67) | 1.60 (0.48)
$\ast$ The actual MSE values are $0.01\times$ the reported values.
Table 2: Estimation results with SCAD for the equation (3.5) when the error
process is AR(1) with the AR coefficient $\rho=0.9$. The other configurations
are identical to Table 1. ARMA(1,1) with $\rho=0.8,\phi=0.4$
---
$(\mu,\sigma)$ | Methods | $n=50$ | $n=100$ | $n=200$
$(0.1,0.5)$ | PMWLS | 5.83 (1.07) | 2.62 (0.71) | 0.83 (0.41)
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 3.53 (0.78) | 1.94 (0.57) | 0.50 (0.31)
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 3.38 (0.75) | 1.93 (0.55) | 0.50 (0.30)
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 3.40 (0.74) | 1.78 (0.54) | 0.46 (0.29)
PWLS | 5.73 (1.05) | 2.66 (0.71) | 0.90 (0.43)
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 3.52 (0.78) | 1.95 (0.57) | 0.43 (0.29)
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 3.22 (0.73) | 1.93 (0.56) | 0.59 (0.33)
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 3.60 (0.75) | 1.89 (0.55) | 0.51 (0.31)
$(0.5,0.5)$ | PMWLS | 4.86 (0.91) | 2.35 (0.67) | 0.85 (0.40)
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 4.09 (0.76) | 1.49 (0.50) | 0.58 (0.32)
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 4.16 (0.74) | 1.53 (0.50) | 0.58 (0.31)
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 4.46 (0.77) | 1.58 (0.51) | 0.55 (0.30)
PWLS | 9.08 (1.30) | 3.07 (0.77) | 1.01 (0.44)
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 4.68 (0.80) | 1.86 (0.55) | 0.88 (0.37)
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 4.29 (0.73) | 1.81 (0.51) | 0.79 (0.35)
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 4.56 (0.75) | 1.94 (0.55) | 0.79 (0.35)
$\ast$ The actual MSE values are $0.01\times$ the reported values.
Table 3: Estimation results with SCAD for the equation (3.5) when the error
process is ARMA(1,1) with AR and MA coefficients $\rho=0.8$ and $\phi=0.4$.
The other configurations are identical to Table 1.

Tables 4-6 demonstrate selection results of PMWLS and PWLS methods with the
SCAD penalty. The results for the LASSO penalty are provided in Tables S4-S6
in the supplementary material. True positive (TP) counts the number of
significant estimates among the significant true parameters and true negative
(TN) counts the number of insignificant estimates among the insignificant true
parameters. As the sample size increases, the values of TP approaches the true
value. The performance in terms of TN is better for SCAD compared to LASSO.
These results correspond to Theorem 4 since LASSO does not fulfill all
conditions in Assumption 2 as discussed. Lastly, selection performance between
PMWLS and PWLS methods are comparable. Different from the estimation
performance, the results are comparable over the choice of weight matrices
including no weight matrix.

AR(1) with $\rho=0.5$
---
$(\mu,\sigma)$ | Methods | TP |  | TN
50 | 100 | 200 |  | 50 | 100 | 200
$(0.1,0.5)$ | PMWLS | 1.35 | 2.04 | 2.40 |  | 16.99 | 17.00 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 1.31 | 2.02 | 2.33 |  | 17.00 | 17.00 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 1.26 | 1.97 | 2.35 |  | 16.97 | 17.00 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 1.11 | 1.93 | 2.33 |  | 17.00 | 17.00 | 17.00
PWLS | 1.27 | 2.07 | 2.38 |  | 17.00 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 1.28 | 1.98 | 2.35 |  | 16.98 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 1.26 | 2.00 | 2.31 |  | 17.00 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 1.17 | 1.92 | 2.36 |  | 16.98 | 17.00 | 17.00
$(0.5,0.5)$ | PMWLS | 1.51 | 2.03 | 2.39 |  | 16.99 | 17.00 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 1.38 | 1.99 | 2.39 |  | 17.00 | 16.99 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 1.36 | 2.00 | 2.41 |  | 16.96 | 17.00 | 16.99
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 1.19 | 1.90 | 2.39 |  | 17.00 | 17.00 | 16.99
PWLS | 1.35 | 1.81 | 2.24 |  | 17.00 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 1.35 | 1.98 | 2.32 |  | 17.00 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 1.25 | 1.90 | 2.32 |  | 17.00 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 1.20 | 1.80 | 2.25 |  | 17.00 | 17.00 | 17.00
|  |  |  |  |  |  |  |
Table 4: Selection results with SCAD for the equation (3.5) when the error
process is AR(1) with the AR coefficient $\rho=0.5$. TP counts the number of
significant estimates among the significant true parameters and TN counts the
number of insignificant estimates among the insignificant true parameters. The
other configurations are identical to Table 1.  AR(1) with $\rho=0.9$
---
$(\mu,\sigma)$ | Methods | TP |  | TN
50 | 100 | 200 |  | 50 | 100 | 200
$(0.1,0.5)$ | PMWLS | 1.69 | 2.04 | 2.50 |  | 16.97 | 17.00 | 16.99
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 1.69 | 1.91 | 2.50 |  | 17.00 | 17.00 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 1.64 | 1.88 | 2.46 |  | 17.00 | 17.00 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 1.58 | 1.89 | 2.47 |  | 17.00 | 17.00 | 17.00
PWLS | 1.62 | 1.99 | 2.49 |  | 16.98 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 1.64 | 1.93 | 2.49 |  | 17.00 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 1.62 | 1.92 | 2.47 |  | 17.00 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 1.62 | 1.93 | 2.49 |  | 17.00 | 17.00 | 17.00
$(0.5,0.5)$ | PMWLS | 1.83 | 2.05 | 2.49 |  | 16.98 | 17.00 | 16.99
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 1.71 | 2.09 | 2.50 |  | 17.00 | 17.00 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 1.65 | 2.09 | 2.49 |  | 16.98 | 17.00 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 1.54 | 2.01 | 2.51 |  | 16.98 | 17.00 | 17.00
PWLS | 1.68 | 1.98 | 2.35 |  | 16.95 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 1.66 | 1.98 | 2.41 |  | 17.00 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 1.62 | 1.96 | 2.39 |  | 17.00 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 1.52 | 1.95 | 2.40 |  | 17.00 | 17.00 | 17.00
|  |  |  |  |  |  |  |
Table 5: Selection results with SCAD for the equation (3.5) when the error
process is AR(1) with the AR coefficient $\rho=0.9$. The other configurations
are identical to Table 4. ARMA(1,1) with $\rho=0.8$, $\phi=0.4$
---
$(\mu,\sigma)$ | Methods | TP |  | TN
50 | 100 | 200 |  | 50 | 100 | 200
$(0.1,0.5)$ | PMWLS | 2.15 | 2.43 | 2.85 |  | 17.00 | 16.98 | 16.98
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 2.15 | 2.40 | 2.83 |  | 16.99 | 17.00 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 2.17 | 2.39 | 2.82 |  | 16.99 | 17.00 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 2.15 | 2.44 | 2.84 |  | 17.00 | 17.00 | 17.00
PWLS | 2.14 | 2.44 | 2.85 |  | 16.98 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 2.16 | 2.40 | 2.85 |  | 17.00 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 2.16 | 2.40 | 2.79 |  | 17.00 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 2.09 | 2.42 | 2.82 |  | 17.00 | 17.00 | 17.00
$(0.5,0.5)$ | PMWLS | 1.93 | 2.53 | 2.79 |  | 16.98 | 16.98 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 1.88 | 2.50 | 2.76 |  | 16.99 | 17.00 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 1.85 | 2.48 | 2.75 |  | 17.00 | 17.00 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 1.82 | 2.49 | 2.77 |  | 17.00 | 17.00 | 17.00
PWLS | 1.80 | 2.47 | 2.73 |  | 17.00 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 1.74 | 2.39 | 2.62 |  | 17.00 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 1.76 | 2.34 | 2.65 |  | 17.00 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 1.72 | 2.36 | 2.66 |  | 17.00 | 17.00 | 17.00
|  |  |  |  |  |  |  |
Table 6: Selection results with SCAD for the equation (3.5) when the error
process is ARMA(1,1) with AR and MA coefficients $\rho=0.8$ and $\phi=0.4$.
The other configurations are identical to Table 4.

Next, we consider a following nonlinear multiplicative model:

$\displaystyle
y_{t}=\frac{1}{1+\exp(-\bm{x}_{t}^{T}\bm{\theta}_{0})}\times\epsilon_{t}.$
(3.6)

For $\epsilon_{t}$, the exponentiated AR processes or an ARMA process are
considered since the $\epsilon_{t}$’s in the equation (3.6) are allowed to
have only positive values. The AR and ARMA coefficients and the parameter
setting of $\bm{\theta}$ are the same as the one in the model (3.5). We
transformed the model in the log scale and apply our approach. Then, we
compare the results with the PWLS method and an ‘additive’ method, where the
estimator from the additive method is calculated as if the data are from a
nonlinear additive model without log transformation. For this simulation, we
provide the results using the SCAD penalty.

Tables 7-9 show estimation performances of PMWLS and PWLS methods for the
model given in (3.6) and Tables 10-12 describe selection performances. For
most data generation settings, our proposed method (PMWLS) shows better
results. In particular, when $(\mu,\sigma)=(0.5,0.5)$ and the sample size is
small, the difference in performance between PMWLS and PWLS methods becomes
more evident. Hence, we argue that PMWLS is preferred over PWLS in practice
since PMWLS shows better finite sample performance. In terms of choice of
weight matrices, the results are similar to those in the first simulation
study. That is, estimation performance is better when we use a weight matrix
for both approaches while the selection performances are comparable with and
without a weight matrix.

AR(1) with $\rho=0.5$
---
$(\mu,\sigma)$ | Methods | 50 | 100 | 200
$(0.1,0.5)$ | PMWLS | 0.88 (0.42) | 0.32 (0.26) | 0.15 (0.17)
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 0.79 (0.39) | 0.20 (0.20) | 0.08 (0.12)
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 0.85 (0.41) | 0.20 (0.20) | 0.08 (0.13)
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 0.87 (0.41) | 0.24 (0.22) | 0.08 (0.13)
PWLS | 0.91 (0.42) | 0.32 (0.25) | 0.13 (0.15)
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 0.73 (0.38) | 0.18 (0.19) | 0.07 (0.12)
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 0.77 (0.39) | 0.21 (0.20) | 0.07 (0.12)
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 0.83 (0.40) | 0.23 (0.21) | 0.08 (0.13)
$(0.5,0.5)$ | PMWLS | 0.68 (0.37) | 0.30 (0.25) | 0.14 (0.16)
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 0.56 (0.33) | 0.19 (0.19) | 0.08 (0.12)
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 0.64 (0.36) | 0.18 (0.19) | 0.09 (0.14)
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 0.75 (0.38) | 0.22 (0.21) | 0.11 (0.15)
PWLS | 0.78 (0.39) | 0.30 (0.24) | 0.14 (0.17)
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 0.71 (0.36) | 0.25 (0.21) | 0.10 (0.13)
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 0.74 (0.37) | 0.25 (0.22) | 0.08 (0.13)
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 0.93 (0.41) | 0.26 (0.23) | 0.12 (0.15)
|  |  |  |
Table 7: Estimation results with SCAD for the equation (3.6) when the error
process is the exponentiated AR(1) with $\rho=0.5$. Mean squared error (MSE)
with standard deviation (SD) values in the parenthesis are presented. The
other configurations are identical to Table 1. AR(1) with $\rho=0.9$
---
$(\mu,\sigma)$ | Methods | 50 | 100 | 200
$(0.1,0.5)$ | PMWLS | 0.61 (0.35) | 0.24 (0.22) | 0.13 (0.16)
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 0.25 (0.22) | 0.03 (0.08) | 0.02 (0.06)
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 0.32 (0.25) | 0.09 (0.14) | 0.01 (0.05)
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 0.40 (0.28) | 0.03 (0.08) | 0.01 (0.05)
PWLS | 0.76 (0.39) | 0.22 (0.21) | 0.12 (0.16)
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 0.27 (0.23) | 0.11 (0.15) | 0.02 (0.06)
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 0.39 (0.27) | 0.09 (0.14) | 0.01 (0.05)
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 0.40 (0.28) | 0.10 (0.14) | 0.01 (0.05)
$(0.5,0.5)$ | PMWLS | 0.55 (0.33) | 0.28 (0.24) | 0.13 (0.16)
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 0.16 (0.18) | 0.08 (0.13) | 0.02 (0.06)
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 0.17 (0.18) | 0.07 (0.12) | 0.01 (0.05)
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 0.19 (0.19) | 0.05 (0.10) | 0.02 (0.06)
PWLS | 0.47 (0.31) | 0.27 (0.23) | 0.14 (0.17)
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 0.22 (0.20) | 0.10 (0.13) | 0.03 (0.06)
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 0.23 (0.21) | 0.11 (0.14) | 0.01 (0.05)
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 0.22 (0.21) | 0.07 (0.12) | 0.02 (0.06)
|  |  |  |
Table 8: Estimation results with SCAD for the equation (3.6) when the error
process is the exponentiated AR(1) with $\rho=0.9$. The other configurations
are identical to Table 7. ARMA(1,1) with $\rho=0.8,\phi=0.4$
---
$(\mu,\sigma)$ | Methods | 50 | 100 | 200
$(0.1,0.5)$ | PMWLS | 0.37 (0.27) | 0.14 (0.17) | 0.09 (0.14)
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 0.18 (0.19) | 0.04 (0.08) | 0.02 (0.06)
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 0.17 (0.19) | 0.03 (0.08) | 0.01 (0.05)
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 0.12 (0.15) | 0.04 (0.09) | 0.01 (0.05)
PWLS | 0.29 (0.24) | 0.16 (0.17) | 0.09 (0.13)
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 0.18 (0.19) | 0.04 (0.09) | 0.02 (0.06)
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 0.17 (0.19) | 0.03 (0.08) | 0.01 (0.05)
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 0.12 (0.16) | 0.04 (0.09) | 0.01 (0.05)
$(0.5,0.5)$ | PMWLS | 0.44 (0.30) | 0.15 (0.18) | 0.07 (0.12)
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 0.19 (0.19) | 0.04 (0.09) | 0.02 (0.06)
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 0.18 (0.19) | 0.03 (0.08) | 0.01 (0.05)
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 0.21 (0.20) | 0.04 (0.08) | 0.02 (0.06)
PWLS | 0.42 (0.29) | 0.15 (0.17) | 0.06 (0.11)
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 0.25 (0.21) | 0.06 (0.10) | 0.04 (0.08)
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 0.23 (0.21) | 0.05 (0.10) | 0.02 (0.06)
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 0.32 (0.25) | 0.03 (0.08) | 0.02 (0.06)
|  |  |  |
Table 9: Estimation results with SCAD for the equation (3.6) when the error
process is the exponentiated ARMA(1,1) with AR and MA coefficients $\rho=0.8$
and $\phi=0.4$. The other configurations are identical to Table 7. AR(1) with
$\rho=0.5$
---
$(\mu,\sigma)$ | Methods | TP |  | TN
50 | 100 | 200 |  | 50 | 100 | 200
$(0.1,0.5)$ | PMWLS | 2.88 | 2.99 | 3.00 |  | 16.83 | 16.92 | 16.94
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 2.85 | 2.99 | 3.00 |  | 16.90 | 16.93 | 16.99
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 2.83 | 2.98 | 3.00 |  | 16.89 | 16.99 | 16.99
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 2.84 | 2.98 | 3.00 |  | 16.95 | 16.98 | 16.99
PWLS | 2.87 | 2.99 | 3.00 |  | 16.84 | 16.94 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 2.87 | 2.99 | 3.00 |  | 16.89 | 16.98 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 2.84 | 2.98 | 3.00 |  | 16.98 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 2.84 | 2.98 | 3.00 |  | 17.00 | 17.00 | 17.00
$(0.5,0.5)$ | PMWLS | 2.94 | 2.99 | 3.00 |  | 16.84 | 16.94 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 2.92 | 2.99 | 3.00 |  | 16.85 | 16.97 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 2.91 | 2.99 | 2.99 |  | 16.93 | 16.97 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 2.89 | 2.99 | 2.99 |  | 16.85 | 16.97 | 16.99
PWLS | 2.92 | 2.99 | 3.00 |  | 16.82 | 16.94 | 16.96
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 2.82 | 2.96 | 2.99 |  | 16.91 | 16.98 | 16.99
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 2.79 | 2.95 | 3.00 |  | 16.95 | 17.00 | 16.99
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 2.76 | 2.96 | 2.99 |  | 17.00 | 16.99 | 16.99
|  |  |  |  |  |  |  |
Table 10: Selection results with SCAD for the equation (3.6) when the error
process is the exponentiated AR(1) with $\rho=0.5$. TP counts the number of
significant estimates among the significant true parameters and TN counts the
number of insignificant estimates among the insignificant true parameters. The
other configurations are identical to Table 4. AR(1) with $\rho=0.9$
---
$(\mu,\sigma)$ | Methods | TP |  | TN
50 | 100 | 200 |  | 50 | 100 | 200
$(0.1,0.5)$ | PMWLS | 2.94 | 3.00 | 3.00 |  | 16.85 | 16.92 | 16.96
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 2.93 | 3.00 | 3.00 |  | 16.94 | 16.99 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 2.90 | 2.98 | 3.00 |  | 16.99 | 17.00 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 2.88 | 3.00 | 3.00 |  | 17.00 | 17.00 | 17.00
PWLS | 2.91 | 3.00 | 3.00 |  | 16.75 | 16.98 | 16.97
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 2.93 | 2.98 | 3.00 |  | 16.98 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 2.88 | 2.98 | 3.00 |  | 17.00 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 2.88 | 2.98 | 3.00 |  | 17.00 | 17.00 | 17.00
$(0.5,0.5)$ | PMWLS | 2.97 | 2.99 | 3.00 |  | 16.74 | 16.90 | 16.97
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 2.96 | 2.98 | 3.00 |  | 16.99 | 16.98 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 2.95 | 2.98 | 3.00 |  | 16.97 | 17.00 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 2.95 | 2.99 | 3.00 |  | 16.96 | 17.00 | 17.00
PWLS | 2.97 | 2.99 | 3.00 |  | 16.88 | 16.97 | 16.96
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 2.94 | 2.98 | 3.00 |  | 16.97 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 2.92 | 2.96 | 3.00 |  | 17.00 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 2.93 | 2.98 | 3.00 |  | 17.00 | 17.00 | 17.00
|  |  |  |  |  |  |  |
Table 11: Selection results with SCAD for the equation (3.6) when the error
process is the exponentiated AR(1) with $\rho=0.9$. The other configurations
are identical to Table 10. ARMA(1,1) with $\rho=0.8,\phi=0.4$
---
$(\mu,\sigma)$ | Methods | TP |  | TN
50 | 100 | 200 |  | 50 | 100 | 200
$(0.1,0.5)$ | PMWLS | 2.99 | 3.00 | 3.00 |  | 16.81 | 16.96 | 16.92
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 2.96 | 3.00 | 3.00 |  | 16.94 | 16.99 | 16.98
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 2.96 | 3.00 | 3.00 |  | 17.00 | 16.99 | 16.99
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 2.98 | 3.00 | 3.00 |  | 17.00 | 17.00 | 16.99
PWLS | 2.99 | 3.00 | 3.00 |  | 16.93 | 16.98 | 16.92
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 2.96 | 3.00 | 3.00 |  | 16.98 | 16.99 | 16.99
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 2.96 | 3.00 | 3.00 |  | 17.00 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 2.98 | 3.00 | 3.00 |  | 17.00 | 17.00 | 17.00
$(0.5,0.5)$ | PMWLS | 2.98 | 3.00 | 3.00 |  | 16.70 | 16.94 | 16.95
$\mbox{PMWLS }{\tiny(\rho=0.5)}$ | 2.95 | 3.00 | 3.00 |  | 16.98 | 16.98 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.9)}$ | 2.95 | 3.00 | 3.00 |  | 16.99 | 16.99 | 17.00
$\mbox{PMWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 2.94 | 3.00 | 3.00 |  | 17.00 | 16.99 | 16.99
PWLS | 2.97 | 3.00 | 3.00 |  | 16.76 | 16.95 | 16.99
$\mbox{PWLS }{\tiny(\rho=0.5)}$ | 2.93 | 3.00 | 3.00 |  | 16.91 | 16.98 | 16.99
$\mbox{PWLS }{\tiny(\rho=0.9)}$ | 2.92 | 2.99 | 3.00 |  | 16.97 | 17.00 | 17.00
$\mbox{PWLS }{\tiny(\rho=0.8,\phi=0.4)}$ | 2.90 | 3.00 | 3.00 |  | 16.97 | 17.00 | 17.00
|  |  |  |  |  |  |  |
Table 12: Selection results with SCAD for the equation (3.6) when the error
process is the exponentiated ARMA(1,1) with AR and MA coefficients $\rho=0.8$
and $\phi=0.4$. The other configurations are identical to Table 10.

For the comparison between PMWLS and the additive method, where the estimator
of the additive method is calculated as if the data are from a nonlinear
additive model without log transformation, the MSE values using the additive
method are large. The detailed results are provided in Tables S7 and S8 in the
supplementary material. This has been pointed out by Bhattacharyya et al.
(1992) that a mis-specified additive nonlinear model while the true model is a
multiplicative nonlinear model may lead to an inconsistent estimator. On the
other hand, both PMWLS and the additive method successfully discriminate
significant and insignificant parameters but PMWLS outperformed the additive
method in terms of TP.

##  4 Application to a head-neck position tracking system

In this section, we apply our PMWLS method to the parametric nonlinear model
of the head-neck position tracking task which is introduced in detail in
Ramadan et al. (2018). The penalization is done with the SCAD penalty. The
weight matrix is chosen by fitting the residuals from the PMWLS method without
the weight matrix to the ARMA(1,1) process. The weight matrix for the PWLS
method is also constructed in a similar way. The choice of the weight matrix
in real data analysis can be flexible, but fitting residuals from the
unweighted method to ARMA(1,1) worked well even for the data with slowly
decreasing autocorrelation.

There are ten subjects in total, and each subject participated in three trials
of an experiment. In each trial, the subject’s head-neck movement as an angle
(radian) for 30 seconds was collected with measurement frequency of 60Hz, i.e.
1800 observations per each trial. Reference signals as an input guided the
subjects to follow with their eyes. Since the neurophysiological parameters
are different by the subjects, the model is separately fitted to each subject.

Parameters | Max | Min | Description
---|---|---|---
$K_{vis}\left[\frac{Nm}{rad}\right]$ | $10^{3}$ | 50 | Visual feedback gain
$K_{vcr}\left[\frac{Nms^{2}}{rad}\right]$ | $10^{4}$ | 500 | Vestibular feedback gain
$K_{ccr}\left[\frac{Nm}{rad}\right]$ | 300 | 1 | Proprioceptive feedback gain
$\tau[s]$ | 0.4 | 0.1 | Visual feedback delay
$\tau_{1A}[s]$ | 0.2 | 0.01 | Lead time constant of the irregular vestibular afferent neurons
$\tau_{CNS1}[s]$ | 1 | 0.05 | Lead time constant of the central nervous system
$\tau_{C}[s]$ | 5 | 0.1 | Lag time constant of the irregular vestibular afferent neurons
$\tau_{CNS2}[s]$ | 60 | 5 | Lag time constant of the central nervous system
$\tau_{MS1}[s]$ | 1 | 0.01 | First lead time constant of the neck muscle spindle
$\tau_{MS2}[s]$ | 1 | 0.01 | Second lead time constant of the neck muscle spindle
$B\left[\frac{Nms}{rad}\right]$ | 5 | 0.1 | Intrinsic damping
$K\left[\frac{Nm}{rad}\right]$ | 5 | 0.1 | Intrinsic stiffness
Table 13: The neurophysiological parameters of the head-neck position tracking
model. The notation and description are adopted from Ramadan et al. (2018).
Max and Min are the range of the parameter values. The values in brackets are
units of the parameters.

The parametric model of the head-neck position tracking system involves a
highly nonlinear structure with 12 parameters to be estimated. Table 13 shows
a list of the parameters for the head-neck position tracking system. Each
parameter measures a neurophysiological function such as visual feedback gain,
vestibular feedback gain, and so on. The parametric model of the head-neck
position tracking task suffers from its complicated structure and limited data
availability, which could bring overfitting and non-identifiability. Thus, a
penalized method has been considered to stabilize parameter estimation and
build a sparse model for identifiability (Ramadan et al., 2018; Yoon et al.,
2022). The penalization method is applied to let the parameter values be
shrunk to the pre-specified values, called the typical values, instead of the
zero value since the parameters in the model have their physical meanings. For
the determination of the typical values, we pre-estimated the parameters 10
times and set the average as the typical values, $\bm{\tilde{\theta}}$. When
overfitting is highly concerned, the obtained estimates tend to possess non-
ignorable gaps with the true optimum. i.e., small changes in the initial
points could bring relatively large deviations among the fitted estimates.
Hence, averaging over the overfitted estimates may lead to a new estimate much
closer to the true optimum. Thus, we chose the average of pre-estimated values
as our typical values. We illustrate this in Figure 3. The overfitted
estimates
($\bm{\tilde{\theta}}_{1},\bm{\tilde{\theta}}_{2},\bm{\tilde{\theta}}_{3},\
\mbox{and}\ \bm{\tilde{\theta}}_{4}$) in Figure 3 are spread out widely around
the true optimum, $\bm{\theta}_{0}$ but the average, $\bm{\tilde{\theta}}$, is
closer to $\bm{\theta}_{0}$.

Figure 3: Illustration of the overfitted estimates. The background pattern
describes an exemplifying objective function to maximize. The lighter colors
mean higher values of the objective function. $\bm{\theta}_{0}$ indicates the
true maximum and
$\bm{\tilde{\theta}}_{1},\bm{\tilde{\theta}}_{2},\bm{\tilde{\theta}}_{3},\
\mbox{and}\ \bm{\tilde{\theta}}_{4}$ stand for the obtained estimates from 4
different pre-estimations and $\bm{\tilde{\theta}}$ stands for the average of
$\bm{\tilde{\theta}}_{1},\bm{\tilde{\theta}}_{2},\bm{\tilde{\theta}}_{3},\
\mbox{and}\ \bm{\tilde{\theta}}_{4}$.

As described in the introduction, the fitted values from the additive error
model with the penalized ordinary least square method used in Yoon et al.
(2022) indicate a possibility of the multiplicative errors and their residuals
still exhibit temporal dependence. Thus, we apply our approach in estimating
the neurophysiological parameters, $\bm{\theta}$, of the head-neck position
tracking model. Then we compare our PMWLS method to an additive approach
without log transformation (Yoon et al., 2022) and PWLS method. Note that both
PMWLS and PWLS methods are applied after log transformation by assuming
multiplicative errors. For the comparison, we consider variance accounted for
(VAF), which is defined as
${\text{VAF}}(\hat{\bm{\theta}})(\%)=\left[1-\frac{\sum^{n}_{t=1}(y_{t}-\hat{y}_{t})^{2}}{\sum^{n}_{t=1}y_{t}^{2}}\right]\times
100,$ where $\hat{y}_{t}$ is the $t$-th component of
$\bm{f}(\bm{x},\hat{\bm{\theta}})$. VAF has been frequently used to assess the
fit from the obtained estimates in biomechanics (Van Drunen et al., 2013;
Ramadan et al., 2018; Yoon et al., 2022). As observed in the expression of
VAF, the estimates with higher VAF values are translated into the estimates
with lower MSE values.

Recall that there are three sets of measurements (three trials) per subject.
We used the measurements from one trial (train set) to fit the model and the
measurements from two other trials (test set) were used to test the fitted
model. Thus, we have one VAF value using the train set and two VAF values
using the test set. We repeat this for each trial as a train set and report
the average VAF values.

No. | Train
---|---
Additive | PWLS | PMWLS
1 | 8453 | 8472 | 8471
2 | 6896 | 6927 | 7139
3 | 8229 | 8229 | 8229
4 | 8476 | 8464 | 8459
5 | 8424 | 8516 | 8647
6 | 8824 | 8795 | 8846
7 | 9308 | 9308 | 9308
8 | 7900 | 8624 | 8771
9 | 8857 | 8858 | 8842
10 | 7672 | 7672 | 7668
Average | 8304 | 8386 | 8438
No. | Test
Additive | PWLS | PMWLS
1 | 8815 | 8810 | 8823
2 | 5785 | 5827 | 6079
3 | 7680 | 7681 | 7680
4 | 8064 | 8058 | 8066
5 | 8350 | 8417 | 8502
6 | 8821 | 8816 | 8819
7 | 9114 | 9114 | 9114
8 | 8187 | 8975 | 9089
9 | 8276 | 8248 | 8327
10 | 7842 | 7842 | 7837
Average | 8093 | 8179 | 8234
$\ast$ The actual VAF values are $10^{-4}\times$ the reported values.
Table 14: VAFs for 10 subjects. ‘No.’ refers the subject number. ‘Additive’
refers to the additive method studied in Yoon et al. (2022), ‘PWLS’ refers to
the method that considers an intercept term, and ‘PMWLS’ refers to our
proposed method. In both PWLS and PMWLS methods, weight matrices in the
objective function are considered. The Train column is for the train set and
the Test column is for the test set.

Figure 4: Estimation and prediction results for the subject No.2 with measured
responses (black line) and fitted values from the additive approach studied in
Yoon et al. (2022) (blue dashed line, $--$), PWLS (yellow dot-dashed line,
$\cdot-\cdot$), and PMWLS (red dotted line, $\cdots$). The upper plot shows
the estimation result using the measurements from one exemplary trial and the
lower plot shows the prediction result for another trial using the fitted
model.

Figure 5: Estimation and prediction results for the subject No. 8 with
measured responses (black line) and fitted values from the additive approach
studied in Yoon et al. (2022) (blue dashed line, $--$), PWLS (yellow dot-
dashed line, $\cdot-\cdot$), and PMWLS (red dotted line, $\cdots$). The upper
plot shows the estimation result using the measurements from one exemplary
trial and the lower plot shows the prediction result for another trial using
the fitted model.

Table 14 shows VAF values of all ten subjects for the additive approach, PWLS
method and PMWLS method. The VAF values among different methods for some
subjects (No. 1, 3, 4, 7, 9 and 10) are similar and the differences are small.
On the other hand, the VAF values from our PMWLS method are higher than the
VAF values from the other methods for the subjects No. 2, 5, 6 and 8. The
difference among the three approaches is more clear when VAF values are
averaged over all subjects. The averages of VAF values over all subjects when
all the parameter values were set to the typical values, i.e. no penalized
estimation method, are 0.8149 for the train set and 0.7928 for the test set.
Hence, the improvement by our approach is larger compared to the improvements
by the other methods as well. We believe these results are originated from the
fact that multiplicative error assumption is valid and our PMWLS method has
successfully captured the error structure underneath the data.

The estimation and prediction results for the subjects No.2 and No.8 are
provided in Figures 4 and 5, respectively. The upper plots in both Figures 4
and 5 exhibit estimation results from one trial out of three trials per
subject as an example. The lower plots in both Figures 4 and 5 show prediction
results for another trial using the fitted model. In both plots, our PMWLS
method (the red dotted line) is better at capturing peaks of the measured
response than the additive method (the blue dashed line), which motivated this
study at the beginning. We believe that the reason our approach outperforms
the additive approach is that the data implicitly have multiplicative
structure. In addition, the PMWLS method also slightly excels the PWLS method
in both estimation and prediction. Therefore, one can benefit from adopting
our proposed method for data with multiplicative structure. The number of
selected parameters, i.e. sensitive parameters, on average were 2.33
(additive), 3.23 (PWLS) and 3.63 (PMWLS) per subject. This might imply that a
smaller number of selected parameters in the additive approach causes poor
performance in estimation and prediction.

##  5 Conclusion

We proposed an estimation and selection method for the parameters in a
nonlinear model when the errors have non-zero mean and temporal dependence.
Our approach can also handle the multiplicative error model as shown in the
simulation study and the real data example. One can consider simply adding an
intercept term to a nonlinear model and estimating the parameters using the
least squares to handle possible non-zero mean. Simulation results show that
both approaches are overall comparable but our approach performs better
compared to the approach with an intercept term when non-zero mean is larger
and a sample size is rather small. We provided asymptotic properties of the
proposed estimator and numerical studies supported our theoretical results.

Introducing a weight matrix to reflect the temporal dependence improves
estimation. One could assume a parametric model for temporal dependence of the
error process and construct a weight matrix from the covariance matrix of the
error process. However, the class of error processes for our theoretical
results is then confined to a class of some parametric models. Also,
simulation results show there is not much difference by a choice of weight
matrices. Hence, we keep freedom to choose a weight matrix without assuming
the exact temporal dependence model.

The proposed method was successfully applied to the head-neck position
tracking model and produced the state-of-the-art performance in both
estimation and prediction. Since the multiplicative structure with temporally
correlated errors is frequently observed in other fields such as finance,
signal processing and image processing, the use of proposed method is not
limited to the application we considered in this study.

The proposed method assumes a fixed number of parameters. Hence, it is natural
to think about an extension to an increasing number of parameters as the
sample size grows, which we leave as a future study. Since we shrink the
estimates toward the typical values, obtaining good typical values can play a
critical role in the final results. Therefore, a better way to locate the
typical values, though the task is very challenging due to the complicated
structure of the nonlinear function, would contribute greatly to solving the
problem in the head-neck position tracking task. In our result, we used the
same typical values as Yoon et al. (2022) and Ramadan et al. (2018) so that
the comparison was carried out fairly.

## Supplementary Materials

The supplementary materials contain proofs of theoretical results discussed in
the text and additional simulation studies.

## Acknowledgements

Wu was supported by Ministry of Science and Technology of Taiwan under grants
(MOST 111-2118-M-259 -002 -). Choi and Lim were supported by National Research
Foundation of Korea (NRF) grant funded by the Korea government (MSIT)
(NRF-2019R1A2C1002213, 2020R1A4A1018207, 2021R1A2B5B01002620).

## References

  * Athreya and Pantula (1986) Athreya, K. B. and S. G. Pantula (1986).  A note on strong mixing of arma processes.  Statistics $\&$ Probability Letters 4, 187–190.
  * Baker and Foley (2011) Baker, K. R. and K. M. Foley (2011).  A nonlinear regression model estimating single source concentrations of primary and secondarily formed pm2.5.  Atmospheric Environment 45, 3758–3767.
  * Bhattacharyya et al. (1992) Bhattacharyya, B., T. Khoshgoftaar, and G. Richardson (1992).  Inconsistent m-estimators: Nonlinear regression with multiplicative error.  Statistics $\&$ Probability Letters 14, 407–411.
  * Bradley (2005) Bradley, R. C. (2005).  Basic properties of strong mixing conditions. a survey and some open questions.  Probability Surveys 2, 107–144.
  * Breheny and Huang (2011) Breheny, P. and J. Huang (2011).  Coordinate descent algorithms for nonconvex penalized regression, with applications to biological feature selection.  The Annals of Applied Statistics 5(1), 232–253.
  * Chen et al. (2002) Chen, K. J., E. Keshner, B. Peterson, and T. Hain (2002).  Modeling head tracking of visual targets.  Journal of Vestibular Research 12(1), 25–33.
  * Fan and Li (2001) Fan, J. and R. Li (2001).  Variable selection via nonconcave penlized likelihood and its oracla properties.  Journal of the American Statistical Association 96(456), 1348–1360.
  * Forbes et al. (2013) Forbes, P. A., E. de Bruijn, A. C. Schouten, F. C. van der Helm, and R. Happee (2013).  Dependency of human neck reflex responses on the bandwidth of pseudorandom anterior-posterior torso perturbations.  Experimental Brain Research 226(1), 1–14.
  * Geller and Neumann (2018) Geller, J. and M. H. Neumann (2018).  Improved local polynomial estimation in time series regression.  Journal of Nonparametric Statistics 30(1), 1–27.
  * Grenander (1954) Grenander, U. (1954).  On the estimation of regression coefficients in the case of an autocorrelated disturbance.  The Annals of Mathematical Statistics 25(2), 252–272.
  * Jennrich (1969) Jennrich, R. I. (1969).  Asymptotic properties of nonlinear least squares estimators.  The Annals of Mathematical Statistics 40(2), 633–643.
  * Jiang et al. (2012) Jiang, X., J. Jiang, and X. Song (2012).  Oracle model selection for nonlinear models based on weighted composite quantile regression.  Statistica Sinica 22, 1479–1506.
  * Lai (1994) Lai, T. L. (1994).  Asymptotic properties of nonlinear least squares estimates in stochastic regression models.  The Annals of Statistics 22(4), 1917–1930.
  * Lim et al. (2014) Lim, C., M. Meerschaert, and H.-P. Scheffler (2014).  Parameter estimation for operator scaling random fields.  Journal of Multivariate Analysis 123, 172–183.
  * Lv et al. (2014) Lv, Z., H. Zhu, and K. Yu (2014).  Robust variable selection for nonlinear models with diverging number of parameters.  Statistics and Probability Letters 91, 90–97.
  * Machkouri et al. (2017) Machkouri, M. E., K. Es-Sebaiy, and I. Ouassou (2017).  On local linear regression for strongly mixing random fields.  Journal of Multivariate Analysis 156, 103–115.
  * Moon et al. (2012) Moon, K.-H., S. W. Han, T. S. Lee, and S. W. Seok (2012).  Approximate mpa-based method for performing incremental dynamic analysis.  Nonlinear Dynamics 67, 2865–2888.
  * Paula et al. (2015) Paula, D., B. Linard, D. A. Andow, E. R. Sujii, C. S. S. Pires, and A. P. Vogler (2015).  Detection and decay rates of prey and prey symbionts in the gut of a predator through metagenomics.  Molecular Ecology Resources 15, 880–892.
  * Peligrad and Utev (1997) Peligrad, M. and S. Utev (1997).  Central limit theorem for linear processes.  The Annals of Probability 25(1), 443–456.
  * Peng et al. (1996) Peng, G., T. Hain, and B. Peterson (1996).  A dynamical model for reflex activated head movements in the horizontal plane.  Biological Cybernetics 75(4), 309–319.
  * Pollard and Radchenko (2006) Pollard, D. and P. Radchenko (2006).  Nonlinear least-squares estimation.  Journal of Multivariate Analysis 97, 548–562.
  * Ramadan et al. (2018) Ramadan, A., C. Boss, J. Choi, N. P. Reeves, J. Cholewicki, J. M. Popovich, and C. J. Radcliffe (2018).  Selecting sensitive parameter subsets in dynamical models with application to biomechanical system identification.  Journal of Biomechanical Engineering 140(7), 074503.
  * Roussas et al. (1992) Roussas, G. G., L. T. Tran, and D. Ioannides (1992).  Fixed design regression for time series: Asymptotic normality.  Journal of Multivariate Analysis 40, 262–291.
  * Santos and Barreto (2017) Santos, J. D. A. and G. A. Barreto (2017).  An outlier-robust kernel rls algorithm for nonlinear system identification.  Nonlinear Dynamics 90, 1707–1726.
  * Tibshirani (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society. Series B (Methodological) 58(1), 267–288.
  * Van de Geer (1990) Van de Geer, S. (1990).  Estimating a regression function.  The Annals of Statistics 18(2), 907–924.
  * Van Drunen et al. (2013) Van Drunen, P., E. Maaswinkel, F. Van der Helm, J. Van Dieën, and R. Happee (2013).  Identifying intrinsic and reflexive contributions to low-back stabilization.  Journal of Biomechanics 46(8), 1440–1446.
  * Wang et al. (2007) Wang, H., R. Li, and C.-L. Tsai (2007).  Tuning parameter selectors for the smoothly clipped absolute deviation method.  Biometrika 94(3), 553–568.
  * Wang and Zhu (2009) Wang, H. and J. Zhu (2009).  Variable selection in spatial regression via penalized least squares.  The Canadian Jounal of Statistics 37(4), 607–624.
  * Wood (2010) Wood, S. N. (2010).  Statistical inference for noisy nonlinear ecological dynamic systems.  Nature 466, 1102–1107.
  * Wu (1981) Wu, C.-F. (1981).  Asymptotic theory of nonlinear least squares estimation.  The Annals of Statistics, 501–513.
  * Wu et al. (2014) Wu, S., H. Xue, Y. Wu, and H. Wu (2014).  Variable selection for sparse high-dimensional nonlinear regression models by combining nonnegative garrote and sure independence screening.  Statistica Sinica 24, 1365–1387.
  * Xu and Shimada (2000) Xu, P. and S. Shimada (2000).  Least squares parameter estimation in multiplicative noise models.  Communications in Statistics - Simulation and Computation 29, 83–96.
  * Yang et al. (2017) Yang, W., Y. Wang, and S. Hu (2017).  Some probability inequalities of least-squares estimator in non linear regression model with strong mixing errors.  Communications in Statistics – Theory and Methods 46, 165–175.
  * Yoon et al. (2022) Yoon, K., H. You, W.-Y. Wu, C. Y. Lim, J. Choi, C. Boss, A. Ramadan, J. M. Popovich Jr, J. Cholewicki, N. P. Reeves, et al. (2022).  Regularized nonlinear regression for simultaneously selecting and estimating key model parameters: Application to head-neck position tracking.  Engineering Applications of Artificial Intelligence 113, 104974.
  * Zhang and Liang (2012) Zhang, J.-J. and H.-Y. Liang (2012).  Asymptotic normality of estimators in heteroscedastic semi-parametric model with strong mixing errors.  Communications in Statistics – Theory and Methods 41, 2172–2201.

Hojun You, Department of Mathematics, University of Houston, TX, 77056, United
States of America E-mail<EMAIL_ADDRESS>
Kyubaek Yoon, School of Mechanical Engineering, Yonsei University E-mail:
<EMAIL_ADDRESS>Wei-Ying Wu, Department of Applied Mathematics,
National Dong Hwa University E-mail<EMAIL_ADDRESS>Jongeun
Choi, School of Mechanical Engineering, Yonsei University E-mail:
<EMAIL_ADDRESS>Chae Young Lim , Department of Statistics, Seoul
National University E-mail<EMAIL_ADDRESS>