I cannot replicate Stata's results in R log-likelihood maximization

[Noé J Nava]

Dear Community, 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))
    }

[Nick Cox]

Cross-posted at https://stackoverflow.com/questions/...d-maximization