0

I am using R package costrOptim.nl.

I need to minimize a function with the following constraints:

Alpha < sqrt(2*omega) and omega > 0

In my code expressed as:

theta[3] < sqrt(2*theta[1]) and theta[1] > 0

I write these conditions as:

Image

But when I call optimizer and run it.

I'm getting the following problem: 

1: In sqrt(2 * theta[1]) : NaNs produced

Why? Did I set the proper conditions?

Mikhail Kholodkov
  • 23,642
  • 17
  • 61
  • 78
Francesco Campigli
  • 1
  • 1
  • Welcome to Stack Overflow. Can you please modify your code so it's clear what exactly is being executed? – wibeasley Jun 25 '18 at 16:14
  • of course. i have a very long code. thank you for the suggestion – Francesco Campigli Jun 25 '18 at 16:22
  • Thanks, that's better, but I suggest including some more work, such a [reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example?rq=1) where we can run your code an get the same error message you're getting on your machine. Your first two lines of code still look like pseudo-code. (If not --that's your problem --R uses the `&` for the `and` operator). – wibeasley Jun 25 '18 at 16:27

1 Answers1

0
This is my whole code. 

data <- read.delim(file = file, header = FALSE)
ind <- seq(from = 1, to = NROW(data), by = 1)
data <- data.frame(ind = ind, Ret = data$V1, Ret2 = data$V1^2)

colnames(data)[1] <- "Ind"
colnames(data)[2] <- "Ret"
colnames(data)[3] <- "Ret2"
T <- length(data$Ret)

m <- arima(x = data$Ret2, order = c(3,0,0), include.mean = TRUE, method = c("ML"))

b_not <- m$coef

omega <- 0.1
alpha <- 0.005
beta <-  0.9

theta <- c(omega,beta,alpha) # "some" value of theta

s0 <- theta[1]/(1-theta[2])


theta[3] < sqrt(2*theta[1]) # check whether the Feller condition is verified 

N <- 30000 
reps <- 1
rho <- -0.8
n <- 100

heston.II <- function(theta){
set.seed(5)
u <- rnorm(n = N*reps,mean = 0, sd = 1)
u1 <- rnorm(n = N*reps,mean = 0, sd = 1)
u2 <- rho*u + sqrt((1-rho^2))*u1

sigma <- matrix(0, nrow = N*reps, ncol = 1)
ret.int <- matrix(0, nrow = N*reps, ncol = 1)
sigma[1,1] <- s0

for (i in 2:(N*reps)) {
sigma[i,1] <- theta[1] + theta[2]*sigma[i-1,1] + theta[3]*sqrt(sigma[i-1,1])*u1[i]
 # if(sigma[i,1] < 0.00000001){ sigma[i,1] =  s0}
}

for (i in 1:(N*reps)) {
ret.int[i,1] <- sqrt(sigma[i,1])*u2[i]
}


ret <- matrix(0, nrow = N*reps/n, ncol = 1)
ret[1,1] <- sum(ret.int[1:n],1)

for (i in 2:((N*reps)/n)) {
ret[i,] <- sum(ret.int[(n*i):(n*(i+1))])
ret[((N*reps)/n),] <- sum(ret.int[(n*(i-1)):(n*i)])
}

ret2 <- ret^2 

model <- arima(x = ret2, order = c(3,0,0), include.mean = TRUE)

beta_hat <- model$coef
m1 <- beta_hat[1] - b_not[1]
m2 <- beta_hat[2] - b_not[2]
m3 <- beta_hat[3] - b_not[3]
m4 <- beta_hat[4] - b_not[4]

D <- cbind(m1,m2,m3,m4)

DD <- (D)%*%t(D)/1000
DD <- as.numeric(DD)
return(DD)

}

heston.sim <- heston.II(theta)

hin <- function(theta){
h <- rep(NA, 2)
h[1] <- theta[1] 
h[2] <- sqrt(2*theta[1]) - theta[3]
return(h)
}

hin(theta = theta)

.opt <- constrOptim.nl(par = theta, fn = heston.II, hin = hin)
.opt
Francesco Campigli
  • 1
  • 1
  • Two things: (a) move this code from an 'answer' to the 'question' (you can edit questions and (b) this example is probably more complicated than you need to demonstrate the problem you're facing. Read and follow the advice of [How to make a great R reproducible example?](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example?rq=1). – wibeasley Jun 25 '18 at 18:22