simhawkes1(t,par)

function[h] = simhawkes1(t,par)
% SIMHAWKES1  Simulate a univariate Hawkes process 
%             (Constant-unconditional-intensity first-order exponential-
%             decay univariate Hawkes process)
% Inputs      t   - end time (Zero starting time assumed)
%             par - process-parameters structure containing fields
%                   'mu'    - M*1 vector (unconditional intensities)
%                   'alpha' - M*M matrix (multiply the exponent terms)
%                   'beta'  - M*M matrix (degrees of the exponents)
% Outputs     h   - process history, 1*1 cell array; vector H{1} stores
%                   event-occurrence times between 0 and t
% Notes       The code is a horribly inefficient literal translation of 
%             Algorithm 1 of Bravaccino (2004, p. 80)
% Example     See HAWKESDEMO
% Author      Dimitri Shvorob, [e-mail], 12/12/07
h = -log(rand)/par.mu;
if h < t
   do5 = 1;
   i = 0;
   j = 0;
   k = 0;
   n = 1;
   while true
      if do5;
         k = k + 1;
         Lstar(k) = inthawkesm(1,h(n),{h},par);            %#ok
      end
      j = j + 1;
      U(j) = rand;                                         %#ok
      i = i + 1;
      u(i) = -log(U(j))/Lstar(k);                          %#ok
      if i == 1
         Tstar(i) = u(i);                                  %#ok
      else
         Tstar(i) = Tstar(i-1) + u(i);
      end   
      if Tstar(i) > t
         h = {h};
         break
      end   
      j = j + 1;
      U(j) = rand;
      if U(j) <= inthawkesm(1,Tstar(i),{h},par)/Lstar(k)
         n = n + 1;
         h = [h Tstar(i)];
         do5 = true;
      else
         k = k + 1;
         Lstar(k) = inthawkesm(1,Tstar(i),{h},par);
         do5 = false;
      end
   end   
else
   h = {[]};
end