[R] Problem with numerical integration and optimization with BFGS

From: Deepankar Basu <basu.15_at_osu.edu>
Date: Thu, 24 May 2007 18:13:53 -0400


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

#############################
# 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
} 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



R-help_at_stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Thu 24 May 2007 - 22:16:57 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 25 May 2007 - 04:31:19 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.