R cannot replicate results from Stata for log-likelihood maximization

Question

I have been working on this algorithm all week, and I do not seem to find the problem.

I employ Stata's NLSUR command for a simple QUAIDS maximization. Stata requires me to write a program evaluator in which I parametrize my system of equations as well as imposing parametric restrictions. Those with experience in such implementations will find the code below familiar. Then, Stata'NLSUR uses MLE to find the parameters. This is not problem, so I use these results (which are correct) to check my log-likelihood optimization problem in R.

R is a bit trickier, because it requires me to write a similar script to that in Stata, but I need to additionally specified the log-likelihood and so the variance-covariance matrix in a function. This function is shown below. I have tried the BBoptim() solver and the results for the parameters do not match, only those in the variance-covariance matrix. I also wrote my own gradient descent algorithm and check it with other toy examples. Then, I check it with this function. My gradient descent algorithm works with the toy examples, but not with my function. I believe the problem is my function, which is similar to the one employ in Stata's NLSUR.

I have to also point out that I use the Stata's estimated parameters in my own function with the same data and I find the same log-likelihood value. I also predict values with Stata and my R's function and I find the same values. I have also tested if the parameter restrictions have been imposed, and it seems this is the case.

Could you help me figure it out what I am doing wrong? The reason that I need to write this specification is that I will be incorporating truncations in the log-likelihood, but I started step-by-step by first trying to solve the simpler problem without the truncations, but only using the log-likelihood.

quaids.loglike <- function(param, s1, s2, s3, lnp1, lnp2, lnp3, lnp4, lnw){
  # alphas
  a1 <- param[1]
  a2 <- param[2]
  a3 <- param[3]
  a4 <- (1 - a1 - a2 - a3)
  
  # betas
  b1 <- param[4]
  b2 <- param[5]
  b3 <- param[6]
  b4 <- (-b1 - b2 - b3)
  
  # gammas
  g11 <- param[7]  
  g12 <- param[8]
  g13 <- param[9]
  g14 <- (-g11 - g12 - g13)
  
  g21 <- g12
  g22 <- param[10]
  g23 <- param[11]
  g24 <- (-g21 - g22 - g23)
  
  g31 <- g13
  g32 <- g23
  g33 <- param[12]
  g34 <- (-g31 - g32 - g33)
  
  g41 <- g14
  g42 <- g24
  g43 <- g34
  g44 <- (-g41 - g42 - g43)
  
  # lambdas
  l1 <- param[13]
  l2 <- param[14]
  l3 <- param[15]
  
  # Sigmas
  sig11 <- param[16]
  sig12 <- param[17]
  sig13 <- param[19]
  
  sig21 <- sig12
  sig22 <- param[18]
  sig23 <- param[20]
  
  sig31 <- sig13
  sig32 <- sig23
  sig33 <- param[21]
  
  # Lnpindex (Ln a(p) where a0 = 5)
  lnpindex <- 5 + a1*lnp1 + a2*lnp2 + a3*lnp3 + a4*lnp4
  lnpindex <- lnpindex + 0.5*(g11*lnp1*lnp1 + g12*lnp1*lnp2 + g13*lnp1*lnp3 + g14*lnp1*lnp4)
  lnpindex <- lnpindex + 0.5*(g21*lnp2*lnp1 + g22*lnp2*lnp2 + g23*lnp2*lnp3 + g24*lnp2*lnp4)
  lnpindex <- lnpindex + 0.5*(g31*lnp3*lnp1 + g32*lnp3*lnp2 + g33*lnp3*lnp3 + g34*lnp3*lnp4)
  lnpindex <- lnpindex + 0.5*(g41*lnp4*lnp1 + g42*lnp4*lnp2 + g43*lnp4*lnp3 + g44*lnp4*lnp4)
  #for(i in as.character(1:4)) {
  #  for(j in as.character(1:4)) {
  
  #    gij <- get(paste0("g", i, j))
  #    lnpi  <- get(paste0("lnp", i))
  #    lnpj  <- get(paste0("lnp", j))
  #    lnpindex <- lnpindex + 0.5*gij*lnpi*lnpj
  
  #  }
  #}
  
  # b(p) price index
  bofp <- lnp1*b1 + lnp2*b2 + lnp3*b3 + lnp4*b4 
  #for(i in as.character(1:4)) {
  
  #  lnpi <- get(paste0("lnp", i))
  #  bi  <- get(paste0("b", i))
  #  bofp <- bofp + lnpi*bi
  
  #}
  bofp <- exp(bofp)
  
  # The parametric shares
  u1 <- a1 + g11*lnp1 + g12*lnp2 + g13*lnp3 + g14*lnp4 + b1*(lnw - lnpindex) + (l1/bofp)*(lnw - lnpindex)^2
  u2 <- a2 + g21*lnp1 + g22*lnp2 + g23*lnp3 + g24*lnp4 + b2*(lnw - lnpindex) + (l2/bofp)*(lnw - lnpindex)^2
  u3 <- a3 + g31*lnp1 + g32*lnp2 + g33*lnp3 + g34*lnp4 + b3*(lnw - lnpindex) + (l3/bofp)*(lnw - lnpindex)^2
  U <- c(mean(u1, na.rm = TRUE), mean(u2, na.rm = TRUE), mean(u3, na.rm = TRUE))
  
  # The vcov matrix:
  sigma <- c(sig11, sig12, sig13, sig21, sig22, sig23, sig31, sig32, sig33)
  sigma <- matrix(sigma, 3, 3)
  
  # The shares
  S <- cbind(s1, s2, s3)
  
  # the individual log-likelihood
  ll <- dmvnorm(S, U, sigma = sigma, log = TRUE)
  return(sum(ll))
}