% Matlab Program  
% Sample paths of Continuous Time Markov Chain and 
% Stochastic Differential Equation
% SIS Epidemic Model with births and natural deaths
% Population size is constant
% ODE Solution and Three Sample Paths for CTMC and SDE
clear all
%Parameters and Initial Condition
beta=1; b=0.25;
gam=.25;N=100;
init=2;

%Number of simulations, time step for SDE, and end time
dt=0.01;
sim=3;
time=20;

% Euler's method ODE Model
y(1)=init;
for k=1:time/dt
  y(k+1)=y(k)+dt*(beta*(N-y(k))*y(k)/N-(b+gam)*y(k));
end
subplot(2,1,1)
plot([0:dt:time],y,'k-','LineWidth',2);
hold on

%CTMC Model
for k=1:sim
clear i t
t(1)=0;i(1)=init;
j=1;
while i(j)>0 & t(j)<time
  u1=rand;  u2=rand; % two uniform random numbers
  a=(beta/N)*i(j)*(N-i(j))+(b+gam)*i(j);
  probi=beta*(N-i(j))/N/(beta*(N-i(j))/N+b+gam);
  t(j+1)=t(j)-log(u1)/a;
  if u2 <= probi
    i(j+1)=i(j)+1;
  else
    i(j+1)=i(j)-1;
  end
  j=j+1;
end  
subplot(2,1,1)
plot(t,i,'r-','LineWidth',2)
hold on
end
axis([0,time,0,80]);
title('CTMC and ODE')
xlabel('Time t');
ylabel('I(t)');
hold off

%SDE Model
subplot(2,1,2)
plot([0:dt:time],y,'k-','LineWidth',2);
hold on
for k=1:sim
clear i t
i(1)=init; t(1)=0;
j=1;
while i(j)>0 & t(j)<time
   rn=randn; % standard normal random number
   mu=beta*i(j)*(N-i(j))/N-(b+gam)*i(j);
   sigma=sqrt(beta*i(j)*(N-i(j))/N+(b+gam)*i(j));
   i(j+1)=i(j)+mu*dt+sigma*sqrt(dt)*rn;
   t(j+1)=t(j)+dt;
   j=j+1;
   if i(j)<0.
       i(j)==0.;
   end
end
subplot(2,1,2)
plot(t,i,'r-','Linewidth',2);
hold on
end
axis([0,time,0,80]);
title('SDE and ODE')
xlabel('Time t');
ylabel('I(t)');
hold off