I have to modify some R code to take jumps into account in my modelisation.

First I build a simulation of a Garch 1,1:

$r_t=\mu+\sqrt{h_t}\epsilon_t$

$h_t=w_0 +\alpha(r_{t-1}-\mu)^2+\beta h_{t-1} $

where $\epsilon$ follows a normal distribution.

I then try to estimate the parameters, so I build and optimise the log likelihood that way:

    garch_loglik<-function(para,x,mu){
# Parameters
omega0=para[1]
alpha=para[2]
beta=para[3]
# Volatility and loglik initialisation
loglik=0
h=var(x)
# Start of the loop
vol=c()
for (i in 2:length(x)){
h=omega0+alpha*(x[i-1]-mu)^2+beta*h
loglik=loglik+dnorm(x[i],mu,sqrt(h),log=TRUE)
}

Now I have a modified GARCH:

$r_t=\mu+\sqrt{h_t}(\epsilon_t+\sum_{i=0}^{N_{t}} x_{i,t})$

$h_t=w_0 +\alpha(r_{t-1}-\mu)^2+\beta h_{t-1}$

Where epsilon follows a normal distribution, N a poisson(0.1) distribution, $x_{i,t}$ a normal distribution(-0.1,0.3).

I have modified the simulation of my GARCH process consequently.

pois=rpois(n,0.1)
eps=rnorm(n,0,1)+ sum(rnorm(pois[i],-0.1,0.3))

But I strugle to modify the loglikelihood.

From what I understand I have to modify the line with dnorm in the loglikelihood, but I don't know how to find the distribution of a normal law + a sum of normal depending on a poisson process.

Is it possible to compute directly ? (a function to 'mix' distributions) Do I have to build it manually ? (simulating multiple variables to find the distribution.)