Identification from implicit function Suppose my observed data $y$ and $x$ is generated by the following relationship for each observation $i$:
$$ y_i = h(y_i,\theta) + x_i + \varepsilon_i$$
where $x_i$ is a strictly exogenous variable , $\varepsilon_i$ is an i.i.d. error term and $\theta$ the vector or parameters of interest. Assume that the solution to this equation is unique.
When $h$ is sufficiently simple (e.g. linear), I can manually solve for $y_i$ and consistently estimate $\theta$ with OLS.
However, when $h$ is more involved, e.g., $h(y_i, \theta) = (1+\exp(\theta_1 - \theta_2 y_i)) / \theta_2$, then there is no closed-form solution. It feels tempting to estimate $\theta$ using non-linear least squares:
$$ \arg \min_\theta \sum_i(y_i - h(y_i, \theta) - x_i)^2$$
but there are obvious endogeneity concerns with respect to the right-hand side of this expression. 
Is there any instrument that allows dealing with these endogeneity concerns? Or how else is such a non-linear relationship without closed-form solution approached?
 A: 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 - ...
A: 

Or how else is such a non-linear relationship without closed-form
  solution approached?

For maximum likelihood we want the density of $y$ conditional on $x$ (because this is our sample), $f_{y|x}(y|x)$.
We start with making an assumption on the density of $\varepsilon$ conditional on $x$, $f_{\varepsilon|x}(\varepsilon|x)$.
When the relation is not implicit, say $y = ax + \varepsilon$ the change-of-variables method has a Jacobian determinant equal to unity so we simply have
$$f_{y|x}(y|x) = f_{\varepsilon|x}(y-ax|x)$$
and we can proceed as usual.
But when the relation is implicit we have
$$\varepsilon = y-h(y,\theta)-x \implies \frac{\partial \varepsilon}{\partial y} = 1- \frac{\partial h(y,\theta)}{\partial y}$$
so here, the observation density will be
$$f_{y|x}(y|x) = \left|1- \frac{\partial h(y,\theta)}{\partial y}\right| \cdot f_{\varepsilon|x}(y-h(y,\theta)-x|x)$$
