% Script file for Math 305 - lab assignment 8

% Sample paths and expected duration of a bounded simple random walk
 
set(0,'DefaultAxesFontSize',18) % Increases axes labels
sim=1000; % Number of simulations, 
t=zeros(1,sim);
q=0.55;  % Probability of losing
for j=1:sim
  j %  Simulation number; prints on the computer screen, 
  clear r 
  r(1)=50; 
  i=1;
  while r(i)>0 & r(i)<100 
      y=rand; 
      if y<=q
          r(i+1)=r(i)-1;
      else
          r(i+1)=r(i)+1; 
      end % End of if/else 
      i=i+1; 
  end % End of while loop   
  t(j)=i; % Time until absorption  
  if j<=3 % Plots three sample paths.
      l1=stairs([0:1:i-1],r) % Draws stairstep graph, 
      set(l1,'LineWidth',2) % Thickens the line width, 
      hold on % Holds the current plot  
  end % End of sample path loop  
end % end of j loop
mindur=min(t)  % Minimum duration; printed on screen 
maxdur=max(t)  % Maximum duration; printed on screen 
meandur=mean(t) % Mean duration; printed on screen  
stdevdur=std(t) % Standard deviation; printed on screen, 
xlabel('Games') 
ylabel('Capital')
hold off % Erases previous plots before drawing new ones. 

%  Looking at the figure, which shows the first 3 simulations, what was the
%  approximate shortest time and longest time to bankruptcy for these 3?
%  For all the 1000 simulations, what was the shortest, longest and mean?
 
% The next program generates three realizations 
% for the simple birth process 
clear 
set(0,'DefaultAxesFontSize',18); % Increases axes labels. 
b=.7;                    % birthrate
x=linspace(0,100,101); % type help linspace to see what this command does
                       % in this case it is equivalent to x=0:100
x0=1                 % initial condition
y=x0*exp(b*x);             % exponential growth, y(0)=1,
n=linspace(1,100,100); % Defines the population vector, 
for j=1:3; t(1)=0; 
    for i=1:99;
        t(i+1)=t(i)-log(rand)/(abs(b)*n(i)); 
    end % End of i loop.
    s=stairs(t,n); % Draws stairstep graph. 
    set(s,'LineWidth',2); % Thickens the line width, 
    hold on % Holds current plot, 
end % end of j loop
plot(x,y,'k--','LineWidth',2); % Plots the exponential.
axis([0,8,0,100]);
xlabel('Time');
ylabel('Population Size');
hold off % Erases previous plots before drawing new ones.

%  Exercise 1:  Change the number of realizations t0 10.
%  Exercise 2: Change this to a death process; 
%  replace the birthrate .7 with a death rate -.5;
%  change the initial condition to 100;
%  change the population vector n appropriately