[R] Problem with numerical integration and optimization with BFGS

[Deepankar Basu]

Hi R users,

I have a couple of questions about some problems that I am facing with
regard to numerical integration and optimization of likelihood
functions. Let me provide a little background information: I am trying
to do maximum likelihood estimation of an econometric model that I have
developed recently. I estimate the parameters of the model using the
monthly US unemployment rate series obtained from the Federal Reserve
Bank of St. Louis. (The data is freely available from their web-based
database called FRED-II).

For my model, the likelihood function for each observation is the sum of
three integrals. The integrand in each of these integrals is of the
following form:

A*exp(B+C*x-D*x^2)

where A, B, C and D are constants, exp() is the exponential function and
x is the variable of integration. The constants A and D are always
positive; B is always negative, while there is no a priori knowledge
about the sign of C. All the constants are finite.

Of the three integrals, one has finite limits while the other two have
the following limits:

lower = -Inf
upper = some finite number (details can be found in the code below)

My problem is the following: when I try to maximize the log-likelihood
function using "optim" with method "BFGS", I get the following error
message (about the second integral):

> out <- optim(alpha.start, LLK, gr=NULL, method="BFGS", y=urate$y)
Error in integrate(f3, lower = -Inf, upper = upr2) :
        the integral is probably divergent

Since I know that all the three integrals are convergent, I do not
understand why I am getting this error message. My first question: can
someone explain what mistake I am making?

What is even more intriguing is that when I use the default method
(Nelder-Mead) in "optim" instead of BFGS, I do not get any such error
message. Since both methods (Nelder-Mead and BFGS) will need to evaluate
the integrals, my second question is: why this difference? 

Below, I am providing the code that I use. Any help will be greatly
appreciated.


Deepankar


************ CODE START *******************



#############################
# COMPUTING THE LOGLIKELIHOOD
# USING NUMERICAL INTEGRALS
#############################

LLK <- function(alpha, y) {

  n <- length(y)
  lglik <- numeric(n) # TO BE SUMMED LATER TO GET THE LOGLIKELIHOOD

           lambda <- numeric(n-1)    # GENERATING *lstar*
           for (i in 1:(n-1)) {      # TO USE IN THE 
           lambda[i] <- y[i+1]/y[i]  # RE-PARAMETRIZATION BELOW
           } 
           lstar <- (min(lambda)-0.01)


# NOTE RE-PARAMETRIZATION
# THESE RESTRICTIONS EMERGE FROM THE MODEL

  muep <- alpha[1]                                      # NO RESTRICTION
  sigep <-  0.01 + exp(alpha[2])                        # greater than
0.01
  sigeta <- 0.01 + exp(alpha[3])                        # greater than
0.01
  rho2 <- 0.8*sin(alpha[4])                             # between -0.8
and 0.8
  rho1 <- lstar*abs(alpha[5])/(1+abs(alpha[5]))         # between 0 and
lstar
  delta <- 0.01 + exp(alpha[6])                         # greater than
0.01


##########################################
# THE THREE FUNCTIONS TO INTEGRATE
# FOR COMPUTING THE LOGLIKELIHOOD
##########################################

  denom <- 2*pi*sigep*sigeta*(sqrt(1-rho2^2)) # A CONSTANT TO BE USED 
                                              # FOR DEFINING THE 
                                              # THREE FUNCTIONS


  f1 <- function(z1) {  # FIRST FUNCTION

       b11 <- ((z1-muep)^2)/((-2)*(1-rho2^2)*(sigep^2))  
       b12 <-
(rho2*(z1-muep)*(y[i]-y[i-1]+delta))/((1-rho2^2)*sigep*sigeta)
       b13 <- ((y[i]-y[i-1]+delta)^2)/((-2)*(1-rho2^2)*(sigeta^2))
     
     return((exp(b11+b12+b13))/denom)
  }

  f3 <- function(z3) { # SECOND FUNCTION

       b31 <- ((y[i]-rho1*y[i-1]-muep)^2)/((-2)*(1-rho2^2)*(sigep^2))  
       b32 <-
(rho2*(y[i]-rho1*y[i-1]-muep)*z3)/((1-rho2^2)*sigep*sigeta)
       b33 <- ((z3)^2)/((-2)*(1-rho2^2)*(sigeta^2))
     
     return((exp(b31+b32+b33))/denom)
  }

  f5 <- function(z5) { # THIRD FUNCTION

       b51 <- ((-y[i]+rho1*y[i-1]-muep)^2)/((-2)*(1-rho2^2)*sigep^2)  
       b52 <-
(rho2*(-y[i]+rho1*y[i-1]-muep)*(z5))/((1-rho2^2)*sigep*sigeta)
       b53 <- ((z5)^2)/((-2)*(1-rho2^2)*(sigeta^2))
      
     return((exp(b51+b52+b53))/denom)
  }

 
  for (i in 2:n) {   # START FOR LOOP

        upr1 <- (y[i]-rho1*y[i-1])      
        upr2 <- (y[i]-y[i-1]+delta)

     # INTEGRATING THE THREE FUNCTIONS
      out1 <- integrate(f1, lower = (-1)*upr1, upper = upr1) 
      out3 <- integrate(f3, lower = -Inf, upper = upr2)
      out5 <- integrate(f5, lower= -Inf, upper = upr2)
  
       pdf <- (out1$val + out3$val + out5$val)
 
     lglik[i] <- log(pdf) # LOGLIKELIHOOD FOR OBSERVATION i

     }               # END FOR LOOP

 return(-sum(lglik)) # RETURNING NEGATIVE OF THE LOGLIKELIHOOD
                     # BECAUSE optim DOES MINIMIZATION BY DEFAULT
}

***************** CODE ENDS *********************************

Then I use:

> urate <- read.table("~/Desktop/UNRATE1.txt", header=TRUE) # DATA
> alpha.start <- c(0.5, -1, -1, 0, 1, -1) # STARTING VALUES
> out <- optim(alpha.start, LLK, gr=NULL, y=urate$y) # THIS GIVES NO
ERROR

or

> out <- optim(alpha.start, LLK, gr=NULL, method="BFGS", y=urate$y)
Error in integrate(f3, lower = -Inf, upper = upr2) :
        the integral is probably divergent