function [ST, realvarT, S, v]=Fang(k, r, S0, v0, kappa, theta, eta, lambda0, kappalambda, thetalambda, etalambda, mu, delta, rho, T) 

    % Simulates a price process for the Fang model. The stochastic jump rate follows a square root process. 
    % Number of jumps conditional the jump rate is Poisson. Jump sizes are normally distributed. Between jumps
    % a continuous stochastic volatility path of heston type is simulated. The algorithm returns the asset price 
    % and the variance realized at maturity, and the whole asset and stochastic volatility paths in [0,T].
    
    % k= number of mesh points; 
    % r= interest rate;
    % S0= initial stock price;
    % v0= initial volatility;  
    % kappa theta, eta= stochastic variance diffusion parameters; 
    % kappalambda, thetalambda, etalambda= stochastic jump rate diffusion parameters;   
    % lambda0= initial jump rate;
    % delta= variance of the jumps;
    % mu= jump percentage/compensator
    % rho= correlation between stock and volatility,
    % T= maturity
    
     step=T/k;
     S=zeros(k,1);
     v=zeros(k,1);
     y=zeros(k,1);
     tau=zeros(k,1);
     lambda=zeros(k,1);
     realvar=zeros(k,1);                        
     v(1)=v0;
     S(1)=S0;
     y(1)=lambda0;
     realvar(1)=v0*step;
     lambda(1)=y(1)*step;
    
      for j=1:k                                                                                           %Simulates a path for the jump rate
        y(j+1)=abs(y(j))+kappalambda*(thetalambda-abs(y(j)))*step+etalambda*sqrt(abs(y(j))*step)*randn;
         lambda(j+1)=lambda(j)+(step/2)*(y(j)+y(j+1));
         abs(lambda(j+1));
       end

     R=poissrnd(lambda(k));                                         %Simulates a Poisson distributed number of jumps of parameter lambda(k)
  
     times=[0,sort( rand(1,R))*T , T];                              %Jump times given the number of jumbs are uniformly distributed
    while(T-times(R+1)<step)
           times=[0,sort( rand(1,R))*T , T];
    end
         
  
   tau(1)=1;
     for j=1:R+1
     for t=tau(j):k
        if( times(j+1)-t*step<step)
         tau(j+1)=t+1;
        break;
        end
     end
      tau(R+2)=k;
      
        for i=tau(j):tau(j+1)-1                                  %Between jumps the stock price follows an Heston type diffusion
         x=randn;
         v(i+1)=abs(v(i))+kappa*(theta-abs(v(i)))*step+eta*sqrt(abs(v(i))*step)*x;
         S(i+1)=abs(S(i))+(r-y(i)*mu)*abs(S(i))*step+abs(S(i))*sqrt(abs(v(i))*step)*(x*rho+randn*(sqrt(1-rho^2)));
         v(i+1)=abs(v(i+1));
         realvar(i+1)=realvar(i)+(step/2)*(v(i)+v(i+1));
        end

      if(j ~=R+1)                                                   %Increments of the path when a jump occurs
         x=randn;
         J=log(mu+1)-delta^2/2+delta*randn;
         jump=(exp(J)-1);
      S(tau(j+1))=abs(S(tau(j+1)-1))+(r-y(j)*mu)*abs(S(tau(j+1)-1))*step+abs(S(tau(j+1)-1))*sqrt(abs(v(tau(j+1)-1))*step)*(x*rho+randn*(sqrt(1-rho^2)))+S(tau(j+1)-1)*jump;
      v(tau(j+1))=abs(v(tau(j+1)-1))+kappa*(theta-abs(v(tau(j+1)-1)))*step+eta*sqrt(abs(v(tau(j+1)-1))*step)*randn;
       realvar(tau(j+1))=realvar(tau(j+1))+J^2;
      end
     end;   
  
      realvarT=realvar(k);
      ST=S(k);
     
end
