Maximum Likelihood for multivariable non-linear normal distrubution

Question

I am trying to write a code to determine g, sigma and lambda using maximum likelihood for the following Function Equation. I have needed to restrict the intervals to ensure that the equation is not taking the log of a negative number, however, it is coming up with the error "L-BFGS-B needs finite values". What can I do to fix this?

lnQs<- c(0.4452211,  3.2828926,  3.0752400,  2.6613305,  5.8122312,  0.5629881,  0.8445112,  4.4336806,  3.8253957,  0.9336889,  1.5188934,  4.3915304, 2.5368227,  0.3729370,  2.3683679,  2.3555777,  0.5985054,  0.6240360,  0.9462143,  0.5440311,  2.6390102,  0.4921728,  0.2971820,  0.1939826, 0.5621709,  0.7839881,  2.2367834, 10.4101839,  5.8010886,  6.0974008)
Rf<- c(0.0484, 0.0537, 0.0598, 0.0590, 0.0590, 0.0571, 0.0497, 0.0497, 0.0492, 0.0588, 0.0586, 0.0492, 0.0468, 0.0480, 0.0405, 0.0405, 0.0488, 0.0452, 0.0675, 0.0395, 0.0395, 0.0385, 0.0400, 0.0432, 0.0394, 0.0397, 0.0397, 0.0407, 0.0436, 0.0436)
S<- c(0.8,  2.3,  2.2,  3.3,  4.9,  8.7,  0.8,  0.9,  1.8,  2.3,  1.3,  1.9,  5.7,  9.3,  4.9, 18.7,  3.6,  2.4, 15.1, 10.3,  0.8,  2.9, 12.0,  0.8,  9.9,  1.3,  8.9, 12.3,  4.2,  4.2)
T<- c(1.0,  5.0,  5.0,  7.0,  5.1, 21.0, 14.1,  5.0,  5.0, 12.0,  7.0,  3.0,  7.0, 21.0,  7.0, 21.0, 21.0, 21.0, 21.0, 20.0,  5.0, 21.0, 21.0, 21.1, 21.0, 14.0,  12.0, 14.0,  5.0,  5.0)

LL<- function(g,sigma,lambda){
  R=(dnorm(lnQs,(log(1+Rf+lambda)-sigma^2/2)*S-log(((1+Rf+lambda)/(1+g))^T-1),sigma*S^0.5))
#
    -sum(log(R), log=TRUE)
}
fit<-mle(minuslogl=LL, start=list(g=.05, sigma=.2, lambda=.1), method = "L-BFGS-B",lower=c(0,0,0.0915),upper=c(0.13,Inf,Inf))
summary(fit)
#Criteria required 
#   lambda>-(1+Rf)  - easily done with restriction lambda>0
#   lambda>(g-Rf)   - NOT SURE HOW TO DEAL WITH lowest Rf=0.0385, tried putting upper limit on g and lower limit on lambda for now
#   sigma>0         - easily done with restriction sigma>0
#   Problem that L-BFGS-B needs finite values offit<-mle(minuslogl=LL, start=list(g=.077, sigma=.256, lambda=.110), method = "BFGS")

Answer 1

Do you need the Inf limits for your sigma and lambda parameters for MLE? for example the following works without error (although the lambda estimate is pretty bad with high std error):

set.seed(123)
T<-0.5
S<-1
Rf<-2
g <-.1 
sigma <- .5
lambda <- 1
lnQs <- rnorm(100,(log(1+Rf+lambda)-sigma^2/2)*S-log(((1+Rf+lambda)/(1+g))^T-1),sigma*S^0.5)

# negative ll fn  
LL<- function(g,sigma,lambda){
  R <- dnorm(lnQs,(log(1+Rf+lambda)-sigma^2/2)*S-log(((1+Rf+lambda)/(1+g))^T-1),sigma*S^0.5)
  -sum(log(R))
}
fit<-mle(minuslogl=LL, start=list(g=.05, sigma=.2, lambda=.1), method = "L-BFGS-B",lower=c(0,0,0.0915),upper=c(0.13,5,5))
summary(fit)
Maximum likelihood estimation

Call:
mle(minuslogl = LL, start = list(g = 0.05, sigma = 0.2, lambda = 0.1), 
    method = "L-BFGS-B", lower = c(0, 0, 0.0915), upper = c(0.13, 
        5, 5))

Coefficients:
       Estimate Std. Error
g       0.10583  16.869598
sigma   0.45412   0.032111
lambda  0.39076 363.867321

With the data you provided, the LL evaluates to infinity at many points with many different values of your parameters, for example using a simple grid search on [0,1]x[0,1]x[0,1] you can find the following points where LL evaluates to infinity (and sometimes NAN too), that's why L-BFGS-G fails:

g <- seq(0,1,length=10)
sigma <- seq(0,1,length=10)
lambda <- seq(0,1,length=10)
# grid search
for (i in 1:10) {
  for (j in 1:10) {
    for (k in 1:10) {
      if (LL(g[i],sigma[j],lambda[k]) == Inf) {
        print(paste(g[i],sigma[j],lambda[k]))
      }
    }
  }
}

Some points where LL evaluates to Inf

[1] "0 0 0"
[1] "0 0 0.111111111111111"
[1] "0 0 0.222222222222222"
[1] "0 0 0.333333333333333"
[1] "0 0 0.444444444444444"
[1] "0 0 0.555555555555556"
[1] "0 0 0.666666666666667"
[1] "0 0 0.777777777777778"
[1] "0 0 0.888888888888889"
[1] "0 0 1"
[1] "0 0.111111111111111 0.111111111111111"
[1] "0 0.111111111111111 0.222222222222222"
[1] "0 0.111111111111111 0.333333333333333"
[1] "0 0.111111111111111 0.444444444444444"

There is a simple hack that may fix the problem:

LL<- function(g,sigma,lambda){
      R=(dnorm(lnQs,(log(1+Rf+lambda)-sigma^2/2)*S-log(((1+Rf+lambda)/(1+g))^T-1),sigma*S^0.5))
      #
      ll <- -sum(log(R), log=TRUE)
      ifelse(ll == Inf, 9999999, ll) # return a large enough number if Inf
}

fit<-mle(minuslogl=LL, start=list(g=.05, sigma=.2, lambda=.1), method = "L-BFGS-B",lower=c(0,0,0.0915),upper=c(0.13,Inf,Inf))
summary(fit)


Maximum likelihood estimation

Call:
mle(minuslogl = LL, start = list(g = 0.05, sigma = 0.2, lambda = 0.1), 
    method = "L-BFGS-B", lower = c(0, 0, 0.0915), upper = c(0.13, 
        Inf, Inf))

Coefficients:
        Estimate Std. Error
g      0.1249300 0.07295480
sigma  0.8795551 0.07349924
lambda 0.0915000 0.07459289

-2 log L: 160.1491