Here is some rough code in R implementing the non-linear two stage least squares estimator under the assumption that $\mathbb E[\epsilon\lvert x] = 0$. Just to see if it works ... it seems to be working up to the point of standard over-underflow issues there might be. But the non-linear two stage least squares estimator should be implemented in most statistical software.
N <- 10000
x <- rnorm(N)
e <- rnorm(N)
x <- c(0,x)
e <- c(0,e)
theta_1 <- 1
theta_2 <- 0.5
# function to solve for y for simulation of data given x and error e
makefunction <- function(x,e)
{
g <- function(y)
{
out <- y - (1 + exp(theta_1 - theta_2*y)/theta_2) - x - e
return(out)
}
return(g)
}
y <- rep(NA,N+1)
for (i in 1:(N+1))
{
g <- makefunction(x[i],e[i])
y[i] <- uniroot(g,lower=-100,upper=100)$root
}
# y is endogenous ... non-zero covariance
cov(y,e)
# x is exogenous
cov(x,e)
# The error-function
r <- function(theta_1,theta_2)
{
out <- y - (1 + exp(theta_1 - theta_2*y)/theta_2) - x
return(out)
}
# Instruments are function of x
X <- cbind(x,x^2,x^3)
objective <- function(theta)
{
a1 <- theta[1]
a2 <- theta[2]
out <- rbind(r(a1,a2))%*%X%*%solve(t(X)%*%X)%*%t(X)%*%cbind(r(a1,a2))
return(out)
}
optim(c(2,3),objective)
based on Amemiya (1983) in the Handbook of Econometrics formula (5.10).
I divide with $\theta_2$ within a prentheses instead of outside, but the code still illustrates how the estimator works so I'm not gonna alter it.
The residual function is
$$r(y,x,\theta) = y -h(y,\theta) - x$$
and
$$\mathbb E[r(y,x,\theta)\lvert x]=0$$ implying that
$$\mathbb E[t(x)r(y,x,\theta)]=0$$
for any function $t(x)$. Basically the two stage non linear least squares estimator is as such a general method of moments estimator.
A more recent text is Wooldridge, J.M. (1996) "Estimating systems of equations with different instruments for different equations" in Journal of Econometrics an article the gist of which is recapitulated in his (2010) Econometric Analysis of Cross Section and Panel Data in the chapter on GMM estimation page 530 - ...