% PROGRAM extremes
%    Simulates effect of catastrophes and/or bonanzas on extinction risk.  
%    Calculates the probability that a population starting at Nc will
%    hit the extinction threshold Nx at or before t=tmax, using a density-
%    independent model.  With probability "probout" each year, the 
%    population growth rate is drawn at random from the vector "outliers", 
%    which contains growth rates during catastrophes and/or bonanzas 

% ************** SIMULATION PARAMETERS *****************************
tmax=20;			% time horizon 
numreps=1E5;        % number of realizations to simulate
Nc=99;				% initial population size
Nx=20;				% extinction threshold
mu=0.01467;			% mean log population growth rate
sigma2=0.01167;		% environmental variance in log pop. growth rate
probout=0.0263;		% annual probability of an outlying growth rate

% vector of observed extreme growth rates
outliers=[0.7875, 1.3077];
% *******************************************************************

rand('state',sum(100*clock));  % Seed uniform random number generator 
randn('state',sum(100*clock)); % Seed normal random number generator 
sigma=sqrt(sigma2);	
numout=length(outliers);       % number of outlying growth rates

% Fill the matrix "lambdas" with "tmax" rows and "numreps" columns with 
% 'typical' population growth rates
lambdas=exp(mu+sigma*randn(tmax,numreps));

% In 'atypical' years, replace entry in "lambdas" with a growth rate
% chosen at random from "outliers"
for i=1:tmax           
   for j=1:numreps      
      if rand<probout					
         lambdas(i,j)=outliers(floor(numout*rand)+1);
      end;                                           
   end;
end;

N=Nc*ones(tmax,numreps);	% Start all realizations at Nc
N=cumprod(lambdas).*N;      % Project all realizations for all years
N=N.*cumprod(N>Nx);         % Zero realizations when they hit the
                            % quasi-extinction threshold and for all
                            % subsequent years
                            
% Compute probability of quasi-extinction at or before tmax                            
ProbExt=sum( (N(tmax,:)==0) )/numreps