#############################################################################################################
# MODEL PARAMETERS UNDER P AND Q

# Parameters under P
P  <- list()
# Log process
P$r           <- 0.0
P$q           <- 0.0
P$rho         <- 0
P$lambda.y0   <- 15
P$lambda.y1   <- 0.00001
P$mu.y        <- -0.02
P$sigma.y     <- 0.02
P$y0          <- log(1000)

# variance process
P$kappa       <- 4
P$theta       <- 0.02
P$sigma       <- 0.00005
P$lambda.v0   <- 0.00001
P$mu.v        <- 0.03
P$v0          <- 0.02

# Risk premiums
RP <- list()
RP$eta.y       <- 0.5
RP$gamma.y.min <- 0.005
RP$eta.v       <- -0.5
RP$gamma.v.maj <- 1

# Moment generating functions
MGF   <- function(x) exp(P$mu.y * x + 0.5 * P$sigma.y^2 * x^2)
fun   <- function(y) (MGF(1) - 1 ) + MGF(y) - MGF(1+y) - RP$gamma.y.min
RP$gamma.y.maj  <- uniroot(f = fun, interval = c(-50,50))$root

MGFV   <- function(x) 1/(1 - x * P$mu.v)

# Parameters under Q
Q <- list()
# Log process
Q$r         <- P$r
Q$q         <- P$q
Q$rho       <- P$rho
Q$lambda.y0 <- MGF(RP$gamma.y.maj) * P$lambda.y0 #
Q$lambda.y1 <- MGF(RP$gamma.y.maj) * P$lambda.y1 #
Q$mu.y      <- P$mu.y + RP$gamma.y.maj * P$sigma.y^2 #
Q$sigma.y   <- P$sigma.y
Q$y0        <- P$y0

# variance process
Q$kappa     <- P$kappa + P$sigma * RP$eta.v #
Q$theta     <- (P$kappa * P$theta)/(Q$kappa + P$sigma * RP$eta.v) #
Q$sigma     <- P$sigma
Q$lambda.v0 <- MGFV(RP$gamma.v.maj) * P$lambda.v0 #
Q$mu.v      <- MGFV(RP$gamma.v.maj) * P$mu.v #
Q$v0        <- P$v0

#############################################################################################################
# SIMULATION LOOP

for (j in 1:(S$n.days*S$n.step.obs)){


  # Variance jumps
  # ----------------------------------------------------------------------- #

  # is there a jump in the variance for this step?
  prob.jump.v        <- 1 - exp(-Q$lambda.v0 * S$h.obs)
  unif.jump.v        <- runif(S$n.traj.obs,0,1)

  # for experiment #2 and experiment #3
    mat               <- cbind(unif.jump.v, prob.jump.v)
    idx.T             <- which(mat[,1] < mat[,2])
    ind.jump.v        <- rep(0, S$n.traj.obs)
    ind.jump.v[idx.T] <- 1

  # which trajectory jumped for this step? which did not?
  idx.T              <- which(ind.jump.v == 1)
  idx.F              <- which(ind.jump.v == 0)

  # Jump size: for trajectories that jumped, generate a exponential
  Z.v[idx.T,j]       <- rexp(n=length(idx.T),rate=1/Q$mu.v)
  Z.v[idx.F,j]       <- 0

  # Variance
  # ----------------------------------------------------------------------- #

  # pre-jump variance
  V_[,j+1]        <- V[,j] + Q$kappa * (Q$theta - V[,j]) * S$h.obs +
    Q$sigma * sqrt(V[,j]) * W.v[,j] *sqrt(S$h.obs)

  # If pre-jump variance is negative, set to 0
  idx.T           <- which(V_[,j+1] < exp(-12))
  V_[idx.T,j+1]   <- exp(-12)

  # add jump
  V[,j+1]         <- V_[,j+1] + Z.v[,j]

  # Log price jumps
  # ----------------------------------------------------------------------- #

  if ( j == 1){
    lambda.y <- Q$lambda.y0
  }else {
    if ( j == 2){
      lambda.y <- Q$lambda.y1}
    else{
      lambda.y <- Q$lambda.y0 + Q$lambda.y1 * V_[,j+1]}
  }


  # is there a jump in the variance for this step?
  prob.jump.y     <- 1 - exp(-lambda.y * S$h.obs)
  unif.jump.y     <- runif(S$n.traj.obs,0,1)

  # If unif < jump prob, jump indicator = 1
    mat               <- cbind(unif.jump.y, prob.jump.y)
    idx.T             <- which(mat[,1] < mat[,2])
    ind.jump.y        <- rep(0, S$n.traj.obs)
    ind.jump.y[idx.T] <- 1

  # which trajectory jumped for this step? which did not?
  idx.T           <- which(ind.jump.y == 1)
  idx.F           <- which(ind.jump.y == 0)

  # Jump size: for trajectories that jumped, generate a exponential
  Z.y[idx.T,j]    <- rnorm(n=length(idx.T), mean=Q$mu.y, sd=Q$sigma.y)
  Z.y[idx.F,j]    <- 0

  # Log price
  # ----------------------------------------------------------------------- #


  alpha <- Q$r - Q$q +  - 0.5*V_[,j+1]  - (MGF(1) - 1)*lambda.y

  Y_[,j+1] <- Y[,j] +  alpha *S$h.obs + sqrt(V_[,j+1])*W.y[,j] *sqrt(S$h.obs)

  Y[,j+1]  <- Y_[,j+1] + Z.y[,j]


}