**[R] nl** -- Nonlinear least-squares estimation

__Syntax__

Interactive version

**nl** **(***depvar***=**<*sexp*>**)** [*if*] [*in*] [*weight*] [**,** *options*]

Programmed substitutable expression version

**nl ***sexp_prog* **:** *depvar* [*varlist*] [*if*] [*in*] [*weight*] [**,** *options*]

Function evaluator program version

**nl** *func_prog* **@** *depvar* [*varlist*] [*if*] [*in*] [*weight*]**,**
{__param__**eters(***namelist***)**|__nparam__**eters(***#***)**} [*options*]

where

*depvar* is the dependent variable;
*<sexp>* is a substitutable expression;
*sexp_prog* is a substitutable expression program; and
*func_prog* is a function evaluator program.

*options* Description
-------------------------------------------------------------------------
Model
__va__**riables(***varlist***)** variables in model
__in__**itial(***initial_values***)** initial values for parameters
* __param__**eters(***namelist***)** parameters in model (function evaluator
program version only)
* __nparam__**eters(***#***)** number of parameters in model (function
evaluator program version only)
*sexp_options* options for substitutable expression program
*func_options* options for function evaluator program

Model 2
__ln__**lsq(***#***)** use log least-squares where ln(*depvar - #*)
is assumed to be normally distributed
__nocons__**tant** the model has no constant term; seldom used
__h__**asconstant(***name***)** use *name* as constant term; seldom used

SE/Robust
**vce(***vcetype***)** *vcetype* may be **gnr**, __r__**obust**, __cl__**uster**
*clustvar*, __boot__**strap**, __jack__**knife**, **hac**
*kernel*, **hc2**, or **hc3**

Reporting
__l__**evel(***#***)** set confidence level; default is **level(95)**
__lea__**ve** create variables containing derivative of
*E*(y)
**title(***string***)** display *string* as title above the table of
parameter estimates
**title2(***string***)** display *string* as subtitle
*display_options* control column formats and line width

Optimization
*optimization_options* control the optimization process; seldom
used
**eps(***#***)** specify *#* for convergence criterion; default
is **eps(1e-5)**
__del__**ta(***#***)** specify *#* for computing derivatives; default
is **delta(4e-7)**

__coefl__**egend** display legend instead of statistics
-------------------------------------------------------------------------
* For function evaluator program version, you must specify
**parameters(***namelist***)** or **nparameters(***#***)**, or both.
**bootstrap**, **by**, **jackknife**, **rolling**, **statsby**, and **svy** are allowed; see
prefix.
Weights are not allowed with the **bootstrap** prefix.
**aweight**s are not allowed with the **jackknife** prefix.
**vce()**, **leave**, and weights are not allowed with the **svy** prefix.
**aweight**s, **fweight**s, and **iweight**s are allowed; see weight.
**coeflegend** does not appear in the dialog box.
See **[R] nl postestimation** for features available after estimation.

__Menu__

**Statistics > Linear models and related > Nonlinear least-squares**
**estimation**

__Description__

**nl** fits an arbitrary nonlinear regression function by least squares.
With the interactive version of the command, you enter the function
directly on the command line or in the dialog box by using a
substitutable expression. If you have a function that you use regularly,
you can write a substitutable expression program and use the second
syntax to avoid having to reenter the function every time. The function
evaluator program version gives you the most flexibility in exchange for
increased complexity; with this version, your program is given a vector
of parameters and a variable list, and your program computes the
regression function.

When you write a substitutable expression program or function evaluator
program, the first two letters of the name must be **nl**. *sexp_prog* and
*func_prog* refer to the name of the program without the first two letters.
For example, if you wrote a function evaluator program named **nlregss**, you
would type **nl regss @ ...** to estimate the parameters.

__Options__

+-------+
----+ Model +------------------------------------------------------------

**variables(***varlist***)** specifies the variables in the model. **nl** ignores
observations for which any of these variables have missing values.
If you do not specify **variables()**, then **nl** issues an error message
with return code 480 if the estimation sample contains any missing
values.

**initial(***initial_values***)** specifies the initial values to begin the
estimation. You can specify a 1 x k matrix, where k is the number of
parameters in the model, or you can specify a parameter name, its
initial value, another parameter name, its initial value, and so on.
For example, to initialize **alpha** to 1.23 and **delta** to 4.57, you would
type

**nl ... , initial(alpha 1.23 delta 4.57)...**

Initial values declared using this option override any that are
declared within substitutable expressions. If you specify a
parameter that does not appear in your model, **nl** exits with error
code 480. If you specify a matrix, the values must be in the same
order that the parameters are declared in your model. **nl** ignores the
row and column names of the matrix.

**parameters(***namelist***)** specifies the names of the parameters in the model.
The names of the parameters must adhere to the naming conventions of
Stata's variables; see **[U] 11.3 Naming conventions**. If you specify
both **parameters()** and **nparameters()**, the number of names in the
former must match the number specified in the latter; if not, **nl**
issues an error message with return code 198.

**nparameters(***#***)** specifies the number of parameters in the model. If you
do not specify names with the **parameters()** option, **nl** names them **b1**,
**b2**, ..., **b***#*. If you specify both **parameters()** and **nparameters()**, the
number of names in the former must match the number specified in the
latter; if not, **nl** issues an error message with return code 198.

*sexp_options* refer to any options allowed by your *sexp_prog*.

*func_options* refer to any options allowed by your *func_prog*.

+---------+
----+ Model 2 +----------------------------------------------------------

**lnlsq(***#***)** fits the model by using log least-squares, which we define as
least squares with shifted lognormal errors. In other words, ln(
*depvar*-*#*) is assumed to be normally distributed. Sums of squares and
deviance are adjusted to the same scale as *depvar*.

**noconstant** indicates that the function does not include a constant term.
This option is generally not needed, even if there is no constant
term in the model, unless the coefficient of variation (over
observations) of the partial derivative of the function with respect
to a parameter is less than **eps()** and that parameter is not a
constant term.

**hasconstant(***name***)** indicates that parameter *name* be treated as the
constant term in the model and that **nl** should not use its algorithm
to find a constant term. As with **noconstant**, this option is seldom
used.

+-----------+
----+ SE/Robust +--------------------------------------------------------

**vce(***vcetype***)** specifies the type of standard error reported, which
includes types that are derived from asymptotic theory (**gnr**), that
are robust to some kinds of misspecification (**robust**), that allow for
intragroup correlation (**cluster** *clustvar*), and that use bootstrap or
jackknife methods (**bootstrap**, **jackknife**); see **[R] ***vce_option*.

**vce(gnr)**, the default, uses the conventionally derived variance
estimator for nonlinear models fit using Gauss-Newton regression.

**nl** also allows the following:

**vce(hac** *kernel* [*#*]**)** specifies that a heteroskedasticity- and
autocorrelation-consistent (HAC) variance estimate be used. HAC
refers to the general form for combining weighted matrices to
form the variance estimate. There are three kernels available
for **nl**:

__nw__**est** | __ga__**llant** | __an__**derson**

*#* specifies the number of lags. If *#* is not specified, N - 2 is
assumed.

**vce(hac** *kernel* [*#*]**)** is not allowed if weights are specified.

**vce(hc2)** and **vce(hc3)** specify alternative bias corrections for the
robust variance calculation. **vce(hc2)** and **vce(hc3)** may not be
specified with the **svy** prefix. By default, **vce(robust)** uses
sigma_j^2 = {n/(n-k)} u_j^2 as an estimate of the variance of the
jth observation, where u_j is the calculated residual and n/(n-k)
is included to improve the overall estimate's small-sample
properties.

**vce(hc2)** instead uses u_j^2/(1-h_jj) as the observation's
variance estimate, where h_jj is the jth diagonal element of the
hat (projection) matrix. This produces an unbiased estimate of
the covariance matrix if the model is homoskedastic. **vce(hc2)**
tends to produce slightly more conservative confidence intervals
than **vce(robust)**.

**vce(hc3)** uses u_j^2/(1-h_jj)^2 as suggested by Davidson and
MacKinnon (1993 and 2004), who report that this often produces
better results when the model is heteroskedastic. **vce(hc3)**
produces confidence intervals that tend to be even more
conservative.

See, in particular, Davidson and MacKinnon (2004, 239), who
advocate the use of **vce(hc2)** or **vce(hc3)** instead of the plain
robust estimator for nonlinear least squares.

+-----------+
----+ Reporting +--------------------------------------------------------

**level(***#***)**; see **[R] estimation options**.

**leave** leaves behind after estimation a set of new variables with the same
names as the estimated parameters containing the derivatives of *E*(y)
with respect to the parameters. If the dataset contains an existing
variable with the same name as a parameter, then using **leave** causes
**nl** to issue an error message with return code 110.

**leave** may not be specified with **vce(cluster** *clustvar***)** or the **svy**
prefix.

**title(***string***)** specifies an optional title that will be displayed just
above the table of parameter estimates.

**title2(***string***)** specifies an optional subtitle that will be displayed
between the title specified in **title()** and the table of parameter
estimates. If **title2()** is specified but **title()** is not, **title2()** has
the same effect as **title()**.

*display_options*: **cformat(***%fmt***)**, **pformat(%***fmt***)**, **sformat(%***fmt***)**, and
**nolstretch**; see **[R] estimation options**.

+--------------+
----+ Optimization +-----------------------------------------------------

*optimization_options*: __iter__**ate(***#***)**, [__no__]__lo__**g**, __tr__**ace**. **iterate(***#***)** specifies
the maximum number of iterations, **log**/**nolog** specifies whether to show
the iteration log, and **trace** specifies that the iteration log should
include the current parameter vector. These options are seldom used.

**eps(***#***)** specifies the convergence criterion for successive parameter
estimates and for the residual sum of squares. The default is
**eps(1e-5)**.

**delta(***#***)** specifies the relative change in a parameter to be used in
computing the numeric derivatives. The derivative for parameter b_i
is computed as {f(X,b_1,b_2,...,b_i + d, b_[i+1],...) - f(X,
b_1,b_2,...,b_i,b_[i+1],...)}/d, where d is delta(b_i + delta). The
default is **delta(4e-7)**.

The following option is available with **nl** but is not shown in the dialog
box:

**coeflegend**; see **[R] estimation options**.

__Remarks__

Remarks are presented under the following headings:

Substitutable expressions
Examples
Some commonly used models
Substitutable expression programs
Example
Function evaluator programs
Example

__Substitutable expressions__

Using a substitutable expression is the easiest way to define your
nonlinear function. Substitutable expressions are just like any other
mathematical expression in Stata, except that the parameters of your
model are bound in braces. There are three rules to follow:

1. Parameters of the model are bound in braces: **{b0}**, **{param}**, etc.

2. Initial values for parameters are given by including an equal
sign and the initial value inside the braces: **{b1=1.267}**,
**{gamma=3}**, etc. If you do not specify an initial value, that
parameter is initialized to zero. The **initial()** option overrides
initial values provided in substitutable expressions.

3. Linear combinations can be included using the notation
**{***eqname***:***varlist***}**:

**{xb:mpg price weight}** is equivalent to
**{xb_mpg}*mpg +** **{xb_price}*price +** **{xb_weight}*weight**

__Examples__

1. To fit the model

y = alpha + beta*x^gamma

where alpha, beta, and gamma are parameters and a starting value
for gamma is one, you would type

**. nl (y = {alpha} + {beta}*x^{gamma=1})**

2. To regress y on a constant and the reciprocal of x you could do

**. nl (y = {b0} + {b1} / x), initial(b0 2 b1 3)**

which obviates the need to generate a new variable equal to 1/x
before calling **regress**. Here b0 is initialized to two and b1 is
initialized to three.

__Some commonly used models__

The following models are used so often that they are built into **nl**.

Exponential regression with one asymptote:

**exp3** y = b0 + b1*b2^x
**exp2** y = b1*b2^x
**exp2a** y = b1*(1-b2^x)

Logistic function (symmetric sigmoid shape)(*):

**log4** y = b0 + b1/(1 + exp(-b2*(x-b3)))
**log3** y = b1/(1 + exp(-b2*(x-b3)))

Gompertz function (asymmetric sigmoid shape):

**gom4** y = b0 + b1*exp(-exp(-b2*(x-b3)))
**gom3** y = b1*exp(-exp(-b2*(x-b3)))

(*) not to be confused with logistic regression

To use any of these, you type

**. nl** *model* **:** *depvar* *indepvar*

For example,

**. nl exp3: y x**
**. nl gom3: response dosage**

Initial values are chosen automatically, though you can override the
defaults by using the **initial()** option.

__Substitutable expression programs -- ____sexp_prog____s__

If you use the same nonlinear function repeatedly, then you can write a
substitutable expression program so that you do not have to retype the
expression every time. The first two letters of the program name must by
**nl**. The **nl***sexp_prog* is to accept a *varlist*, an **if** *exp*, and, optionally,
weights. The program will then return a substitutable expression in the
r-class macro **r(eq)** and, optionally, a title in **r(title)**.

The outline of an **nl***sexp_prog* program is

**program nl***sexp_prog***, rclass**
**version 15.1**
**syntax ***varlist* **[aw fw iw]** **if**
*(obtain initial parameters if desired)*
**return local eq "***<sexp>***"**
**return local title "***title***"**
**end**

__Example__

Returning to the model

y = alpha + beta*x^gamma

one way to obtain initial values is to let gamma = 1 and then run a
regression of x on y to obtain alpha and beta. The substitutable
expression program is

**program nlmyreg, rclass**
**version 15.1**
**syntax varlist(min=2 max=2) [aw fw iw] if**
**local lhs: word 1 of `varlist'**
**local rhs: word 2 of `varlist'**
**regress `lhs' `rhs' [`weight'`exp'] `if'**
**tempname a b**
**scalar `a' = _b[_cons]**
**scalar `b' = _b[`rhs']**
**return local eq "`lhs' =**
**{alpha=`=`a''}+{beta=`=`b''}*`rhs'^{gamma=1}"**
**return local title "`lhs' = alpha+beta*`rhs'^gamma"**
**end**

To fit your model, you type

**. nl myreg: y x**

(There is a space between **nl** and **myreg**, even though the program is
named **nlmyreg**.)

The substitutable expression does not need to account for weights or the
**if** *exp*, though you do need to use them in obtaining initial values.
Also, the substitutable expression is not bound in parentheses, unlike
when typing it in interactively.

__Function evaluator programs -- ____func_prog____s__

If your function is particularly complex, then you may find that writing
one substitutable expression is impractical. In those cases, you can
write a function evaluator program instead. Whenever **nl** needs to
evaluate your function, it calls your program with a vector of
parameters. Your program then fills in the dependent variable with
function values.

Function evaluator programs must accept a *varlist*, an **if** *exp*, and an
option named **at()** that accepts the name of a matrix. It may optionally
accept weights as well. Unlike substitutable expression programs, a
function evaluator program is not declared to be r-class. The outline of
a **nl***func_prog* program is

**program nl***func_prog*
**version 15.1**
**syntax** *varlist* **[aw fw iw] if, at(name)**
**local lhs: word 1 of `varlist'**
**local rhs: subinstr local varlist "`lhs'" "", word**
*(evaluate the function at matrix)*
**replace `lhs' =*** <the function values>* **`if'**
**end**

When evaluating your function, remember to restrict the estimation sample
by using **`if'**. Also, remember to include the weights if using commands
such as **summarize** or **regress** if you intend to do weighted estimation.

__Example__

The CES production function can be written

ln Q = b0 - 1/rho*ln{delta*K^-rho + (1-delta)*L^-rho}

where Q denotes output and b0, rho, and delta are parameters to be
estimated. The function evaluator program is

**program nlces**
**version 15.1**
**syntax varlist(min=3 max=3) [aw fw iw] if, at(name)**
**local logout: word 1 of `varlist'**
**local capital: word 2 of `varlist'**
**local labor: word 3 of `varlist'**
**// Retrieve parameters out of at matrix**
**tempname b0 rho delta**
**scalar `b0' = `at'[1,1]**
**scalar `rho' = `at'[1,2]**
**scalar `delta' = `at'[1,3]**
**// Some temporary variables**
**tempvar kterm lterm**
**generate double `kterm' = `delta'*`capital'^(-1*`rho') `if'**
**generate double `lterm' = (1-`delta')*`labor'^(-1*`rho') `if'**
**// Now fill in dependent variable**
**replace `logout' = `b0' - 1/`rho'*ln(`kterm'+`lterm') `if'**
**end**

If your variables are **logout**, **capital**, and **labor**, then any of the
following methods can be used to estimate the parameters:

1. This method uses b0 = 0 as an initial value by default:

**. nl ces @ logout capital labor, parameters(b0 rho delta)**
**initial(rho 1 delta 0.5)**

2. This method initializes b0 to 2, rho to 1, and delta to 0.5.
Because we do not give parameter names, nl names them b1, b2,
and b3:

**. nl ces @ logout capital labor, nparameters(3) initial(b1 2**
**b2 1 b3 0.5)**

3. This method sets up a vector holding the initial values:

**. matrix ivals = (2, 1, 0.5)**
**. nl ces @ logout capital labor, parameters(b0 rho delta)**
**initial(ivals)**

or

**. nl ces @ logout capital labor, nparameters(3)**
**initial(ivals)**

__Example__

Setup
**. webuse production**

Fit CES production function with initial values **rho**=1 and **delta**=.5
**. nl (lnoutput = {b0} - 1/{rho=1}*ln({delta=0.5}*capital^(-1*{rho}) +**
**(1-{delta})*labor^(-1*{rho})))**

__Stored results__

**nl** stores the following in **e()**:

Scalars
**e(N)** number of observations
**e(k)** number of parameters
**e(k_eq_model)** number of equations in overall model test; always **0**
**e(df_m)** model degrees of freedom
**e(df_r)** residual degrees of freedom
**e(df_t)** total degrees of freedom
**e(mss)** model sum of squares
**e(rss)** residual sum of squares
**e(tss)** total sum of squares
**e(mms)** model mean square
**e(msr)** residual mean square
**e(ll)** log likelihood assuming i.i.d. normal errors
**e(r2)** R-squared
**e(r2_a)** adjusted R-squared
**e(rmse)** root mean squared error
**e(dev)** residual deviance
**e(N_clust)** number of clusters
**e(lnlsq)** value of **lnlsq** if specified
**e(log_t)** **1** if **lnlsq** specified, **0** otherwise
**e(gm_2)** square of geometric mean of (y-k) if **lnlsq**, **1**
otherwise
**e(cj)** position of constant in **e(b)** or **0** if no constant
**e(delta)** relative change used to compute derivatives
**e(rank)** rank of **e(V)**
**e(ic)** number of iterations
**e(converged)** **1** if converged, **0** otherwise

Macros
**e(cmd)** **nl**
**e(cmdline)** command as typed
**e(depvar)** name of dependent variable
**e(wtype)** weight type
**e(wexp)** weight expression
**e(title)** title in estimation output
**e(title_2)** secondary title in estimation output
**e(clustvar)** name of cluster variable
**e(hac_kernel)** HAC kernel
**e(hac_lag)** HAC lag
**e(vce)** *vcetype* specified in **vce()**
**e(vcetype)** title used to label Std. Err.
**e(type)** **1** = interactively entered expression
**2** = substitutable expression program
**3** = function evaluator program
**e(sexp)** substitutable expression
**e(params)** names of parameters
**e(funcprog)** function evaluator program
**e(rhs)** contents of **variables()**
**e(properties)** **b V**
**e(predict)** program used to implement **predict**
**e(marginsok)** predictions allowed by **margins**
**e(marginsnotok)** predictions disallowed by **margins**

Matrices
**e(b)** coefficient vector
**e(init)** initial values vector
**e(V)** variance-covariance matrix of the estimators
**e(V_modelbased)** model-based variance

Functions
**e(sample)** marks estimation sample

__References__

Davidson, R., and J. G. MacKinnon. 1993. *Estimation and Inference in*
*Econometrics*. New York: Oxford University Press.

------. 2004. *Econometric Theory and Methods*. New York: Oxford
University Press.