Suppose the relationship between $y_{i,t}\equiv\log w_{i,t}$ and $z_{i,t}\equiv\left(s_{i,t},j_{i,t},e_{i,t},x_{i,t}\right)^{\top}$ is given by
$$y_{i,t}=g\left(z_{i,t}\right)+\varepsilon_{i,t},$$
where $g\left(\cdot\right)$ can be nonlinear. There can be two cases for this problem
$g\left(\cdot\right)$ is unknown. For example, if one use sieve estimation, which may lead to the high dimension problem you mentioned, to address this problem, he can apply some regulation methods like LASSO. However, this is not the case in this question, so I won't talk too much.
$g\left(\cdot\right)=m\left(z_{i,t};\beta\right)$ where the form of $m\left(\cdot\right)$ is known, like in this question, and $\beta=\left(D,\Pi,\lambda,\theta,\alpha\right)^{\top}$. This happens when there exists some economic theory or model that yields this kind of relationship, or the user needs to identify and analyze how the covariates affects the outcome. In this case, the NLS estimator is easier to analyze since it can be seen a special case of extreme estimator:
$$\hat{\beta}=\mathop{\arg\min}_{\beta\in\mathcal{B}}\hat{S}\left(\beta\right),$$
where $\hat{S}\left(\beta\right)\equiv\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\left(y_{i,t}-m\left(z_{i,t};\beta\right)\right)^{2}$. Following the similar discussion as in extreme estimators, if $\mathcal{B}$ is compact and $\hat{S}\left(\beta\right)$ converge uniformly in probability to some nonrandom function which has a unique minimizer, then $\hat{\beta}$ is consistent, and will be asymptotically normal under additional conditions.
To solve $\hat{\beta}$, you can first write the FOC of this problem
$$\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\left(y_{i,t}-m\left(z_{i,t};\hat{\beta}\right)\right)\frac{\partial m\left(z_{i,t};\hat{\beta}\right)}{\partial \beta}=0.$$
Generally there is no closed form solution for $\hat{\beta}$ hence it must be found by numerical methods. There are two iteration algorithms that can be used.
$$\beta_{r+1}=\beta_{r}-\frac{1}{2}\left(\sum_{i=1}^{N}\sum_{t=1}^{T}\frac{\partial m\left(z_{i,t};\beta_{r}\right)}{\partial\beta}\frac{\partial m\left(z_{i,t};\beta_{r}\right)}{\partial\beta^{\top}}\right)^{-1}\frac{\partial\hat{S}\left(\beta_{r}\right)}{\partial \beta}.$$
- Newton-Raphson Algorithm
$$\beta_{r+1}=\beta_{r}-\left(\frac{\partial^{2}\hat{S}\left(\beta_{r}\right)}{\partial\beta\partial\beta^{\top}}\right)^{-1}\frac{\partial \hat{S}\left(\beta_{r}\right)}{\partial\beta},$$
which requires the matrix in bracket being positive definite to converge.
As for programming, I'm not sure if any statistics software has NLS function. You can google it and check the help document if have one. If there's no existing package, you can use optimize to solve the origin minimization problem, or fsolve to find numerical solution to above FOC.
Update for std. Suppose $\sqrt{NT}\left(\hat{\beta}-\beta_{0}\right)\overset{d}{\to}\mathcal{N}\left(0,V\right)$, then $V$ can be estimated by
$$\begin{align*}\hat{V}&\equiv\left(\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\frac{\partial m\left(z_{i,t};\hat{\beta}\right)}{\partial\beta}\frac{\partial m\left(z_{i,t};\hat{\beta}\right)}{\partial\beta^{\top}}\right)^{-1}\\
&\left(\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\frac{\partial m\left(z_{i,t};\hat{\beta}\right)}{\partial\beta}\frac{\partial m\left(z_{i,t};\hat{\beta}\right)}{\partial\beta^{\top}}\hat{\varepsilon}_{i,t}^{2}\right)\left(\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\frac{\partial m\left(z_{i,t};\hat{\beta}\right)}{\partial\beta}\frac{\partial m\left(z_{i,t};\hat{\beta}\right)}{\partial\beta^{\top}}\right)^{-1}.\end{align*}$$
, where $D, \Pi$, all $\lambda$s $\theta$s, and $\alpha$ are coefficients, $j$, $e$ are age and experience, and $s$, $x$ means sex and education group.
