Log likelihood of (marked) Hawkes Process

Question

I need to calculate the log likelihood of a Hawkes Process and also of a Marked Hawkes Process, where the jump in intensity is not the same for every event.
More specifically, I have some datasets consisting of around 100 events spread over 1 million years (each event happens in one of the years, no fractions) and I want to use likelihood maximization to estimate the parameters of the process.
I tried coding myself using Matlab's fmincon, but it doesn't seem to be working properly. It is too sensitive to the initial guess and to the dataset. Sometimes it throws the 'close to singular or badly scaled matrix' warning. If I change the initial guess for the first parameter I eventually find one that makes it work. I guess this happens because the values of the parameters are on the scale of 10^(-4)?
Here is the code I used for the Hawkes Process. 'times' stores the the year of each event and 'T' is the total length of the process.
Any thoughts on how to deal with this issue?
Thanks

func = @(x)hawkes_likelihood(times, T, x(1), x(2), x(3));
options = optimoptions(@fmincon,'Display','off');
[params, likelihood] = fmincon(func, [length(times)/T 0.001 0.001], ...
                              [], [], [], [0 0 0], [], [],options);

function likelihood = hawkes_likelihood(times, T, c, alpha, beta)
    % Compute the negative log-likelihood of the Hawkes Process
    
    % Integral of the intensity function    
    likelihood = c*T + sum(alpha.*(1-exp(-beta*(T-times)))/beta);
    
    for i = 1:length(times)
        % Intensity function at i-th event
        lambda_t = c + sum(alpha*exp(-beta*(times(i)-times(1:i-1))));
        likelihood = likelihood - log(lambda_t);
    end
end