fitting a GARCH@CARR Model with R

Question

I try to fit a GARCH@CARR (proposed here: https://www.sciencedirect.com/science/article/abs/pii/S1544612318302782) function with R. It's not part of the rugarch package, I m using for other GARCH Models. I tried modeling it myself with both optim and mle, but because I am pretty new to working with R I didn't get anything useful. I use SP500 log returns and log range, obtained from yahoo.finance.

The Model is stated as following:

getSymbols("^GSPC",from="2000-01-01", to="2020-08-01")
Data_SP500<-subset(GSPC, select = 6)
##compute daily log return
return_SP500<-Data_SP500[2:length(Data_SP500),0]
for (n in 2:length(Data_SP500)) {
  return_SP500[n]<-log(as.numeric(Data_SP500[n])/as.numeric(Data_SP500[n-1]))
}  

Range_Data_SP500<-subset(GSPC, select = c(2,3))

##compute Range
Range_SP500<-Data_SP500[1:length(Data_SP500),0]
for (n in 1:nrow(Range_Data_SP500)) {
  Range_SP500[n]<-log(as.numeric(Range_Data_SP500[n,1]))-log(as.numeric(Range_Data_SP500[n,2]))
}
#delete NA
return_SP500<-na.omit(return_SP500)
Range_SP500<-na.omit(Range_SP500)

Data_SP500_1<-merge.xts(Rendite_SP500,Range_SP500)
Data_SP500_1<-na.omit(Data_SP500_1)

#GARCH@CARR Model
#p=q=1
#r_t = roh * lambda_t * z_t, z_t~N(0,1) 
#ln(lambda_t) = omega + alpha_1 * ln(lambda_t-1) + beta_1 * ln(x_t-1) + pi_1 * z_t-1
#x_t = lambda_t * u_t, u_t~LN(-sigma^2/2,sigma)

I tried using altering GARCH Models, available in the rugarch package in a way to fit the GARCH@CARR Model, but it didn't work either. I failed to build anything useful from scratch myself as well. The closest I got so far is with mle:

#GARCH CARR
N <- nrow(Data_SP500_1)#length Timeseries
r <- Data_SP500_1[,1]#returnseries
x <- Data_SP500_1[,2]#Rangeseries

LL <- function(alpha0,beta1,gamma1,roh,mu1,sigma1,mu2,sigma2) {
  # Find residuals
  #
  lambda=mean(Data_SP500_1[,2])
  y1 = r - roh * lambda
  y2 <-for (i in 2:N) {
    lambda[i] - exp(alpha0 + beta1*log(lambda[i-1]) + gamma1*log(x[i-1]))
  } 
  y3 = x - lambda
  #
  # Calculate the likelihood for the residuals (with mu and sigma as parameters)
  #
  R1 = suppressWarnings(dnorm(y1, mu1, sigma1,log = T))
  R3 = suppressWarnings(dlnorm(y3, mu2, sigma2,log = T))
  #
  # Sum the log likelihoods for all of the data points
  #
  -sum(R1+R3)
}
fit <- mle(LL, start = list(alpha0=0.2, beta1=0.4, gamma1=0.2, roh=1, mu1=0, sigma1=1, mu2=0, sigma2=lambda[1]))

But I still get a Error similar to MLE error in R: initial value in 'vmmin' is not finite . Using log() also didn't help.

Error in optim(start, f, method = method, hessian = TRUE, ...) : 
  initial value in 'vmmin' is not finite

Here there are 3 different possible reasons for this error listed, I can't see which one is relevant for this problem, or how to fix that. https://quant.stackexchange.com/questions/44098/mle-error-in-r-initial-value-in-vmmin-is-not-finite Going for a standard method from one of the packages is unfortunately not an option, because this model isn't part in any package, I've seen so far.