How to estimate a nonlinear equation system in R?

Question

I have trouble finding the right packages and methods to estimate a system of three nonlinear equations with cross-equation restrictions using R.

I want to estimate the parameters of a CES production functions using the methodology developed by Miguel A. León-Ledesma, Peter McAdam, and Alpo Willman (2010):

The estimator used for the system is a nonlinear feasible generalized least squares (FGLS) method which accounts for possible cross equation error correlation (much like an SUR model in linear contexts). The estimator performs NLLS on each individual equation and uses the estimated errors to build a variance-covariance (VCV) matrix and then estimates the system by GLS, completing one iteration. The estimated VCV matrix will be updated with each iteration until the system converges to a predetermined criterion.

I have completed the first step of a nonlinear estimation of the individual equations using the nlsLM function and build the vcv matrix of the estimated residuals but I'm struggling as to how I should implement the iterated FGLS procedure.

An example of my code is below:

psi = 1
sigma = 0.6
alpha.e = 0.01
alpha.l = 0.01

These are my starting values

es.f = les ~ (sigma-1)/sigma * (log(psi) + alpha.e*t)+ 
(1- sigma)/sigma * ly + (1+sigma)*log(ebar) + log(1-normals[6])

nl.e.lm = nlsLM(es.f,start = list(psi = psi, sigma = sigma, alpha.e = alpha.e))

les, ly, ebar, tbar are from the data, the elements of normalsare known constants. I only have a limited number of observations (23) which is why I am opting for the systems approach, which offers better identification especially with small sample sizes compared to other approaches such as the Kmenta approximation, which is implemented in the package micEconCES.

I get reasonable results in all three equations for at least 1 or 2 of the parameters, which is encouraging, but of course so far the cross-equation restrictions are not implemented because the equations are just estimated individually.

I would be glad for any help, just pointing me in the right direction would be helpful. I was thinking of trying the BBpackage but I am confused as to how I should implement what I want to do there.

P.S.: Side-question I just remembered: How come the dat$year reference to my data isn't working inside the call to nlsLM? I had to remove all references with $ and predefine those variables separately before my estimations.

edit: I have tried implementing the system using nlsystemfitbut it did not work.

Answer 1

After scouring the web for more information as to how a system of nonlinear equations can efficiently be solved in R, I stumbled upon the nlsur method available in Stata. I know not everyone is able to use Stata, but for those that have the opportunity, it's helpful to know that the implementation in Stata is very easy. I have reused all my datapreparation work and exported only the transformed variables as a .csv

This is the code which solves the following equation system:

#delimit ; /*sets line end indicator to semicolon*/
nlsur (les = ({sigma}-1)/{sigma} * (ln({psi}) + {alphae}*t)+ 
    (1- {sigma})/{sigma} * ly + (1+{sigma})*ln(ebar) + ln(1-gv))
    (lks = ({sigma}-1)/{sigma} * (ln({psi}) + {alphal}*beta*t + beta *  ln(lbar) + (1-beta) * ln(kbar)) + 
    (1- {sigma})/{sigma} * ly +  ln(gv* (1-beta)))
    (ly = ln({psi}) + {sigma}/({sigma} - 1) * ln(gv*(exp({alphal} * beta*t)*(lbar)^beta* 
    (kbar)^(1- beta))^(({sigma}-1)/{sigma}) + (1-gv)* (exp({alphae}*t)*ebar) ^(({sigma}-1)/{sigma}))), 
    initial(sigma 0.6 psi 1 alphae 0.05 alphal 0.05);
#delimit cr /*sets line end indicator back to carriage return*/

This uses a 2stage FGLS to solve the following three equations(if an iterated FGLS is wanted the flag ifgnls can be set after the initials): $$ \log S_E =\frac{\sigma -1}{\sigma}\left(\log{\psi} + \alpha_E \tilde{t}\right) + \frac{1-\sigma}{\sigma}\log{\tilde{Y}}+ (\sigma + 1)\log{\tilde{E}} + \log{\gamma_E} $$ $$ \log{S_L} = \frac{\sigma -1}{\sigma}\left(\log{\psi} + \alpha_L \beta \tilde{t}+\beta \log{\tilde{L}} + (1-\beta)\log{\tilde{K}}\right) + \frac{1-\sigma}{\sigma}\log{{\tilde{Y}}} + \log{\gamma_V \beta} $$ $$ \log\tilde{Y} = \log\psi + \frac{\sigma}{\sigma - 1} \log ( \gamma_{V}\left(e^{\alpha_{L}\tilde{t}\beta}\left(\tilde{L}\right)^{\beta} \left(\tilde{K}\right)^{1-\beta}\right)^{\frac{\sigma-1}{\sigma}} +\gamma_{E}\left(e^{\alpha_{E}\tilde{t}} \tilde{E}\right)^{\frac{\sigma-1}{\sigma}}) $$

The cross equation restrictions on my parameters are automatically recognized and the convergence behavior seems good. You have to use curly braces to explicitly declare the parameters. I personally am not a good enough programmer, but I hope that a package analogous to nlsurin Stata will become available in R.

Edit: One big draw-back is that this method does NOT support Newey-West standard errors which is a big issue for my estimation. When writing Stata support they have answered that development for this is currently not planned. This is a shame because hac-robust estimation is available for the normal nonlinear estimation, guess I will have to find something more suitable after all.