# [R] Problem with numerical integration and optimization with BFGS

Deepankar Basu basu.15 at osu.edu
Fri May 25 00:13:53 CEST 2007

```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

> 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:

> 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

```