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):
the integral is probably divergent
Below, I am providing the code that I use. Any help will be greatly appreciated.
############################# # COMPUTING THE LOGLIKELIHOOD
#############################
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 than0.01
########################################## # THE THREE FUNCTIONS TO INTEGRATE
########################################## 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}
> 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
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.