% PROGRAM simext 
%   Estimates by simulation the quasi-extinction time cumulative
%   distribution function (CDF) for a structured population in an iid
%   stochastic environment

%************* USER-SPECIFIED PARAMETERS ***************
% Enter name of m file containing matrices; e.g., the file  
% hudmats.m contains 4 matrices (A85, A86, A87, A88),
% one for each year of Frost's demographic study
hudmats;						
% Change names of matrices below to correspond with those 
% in the file named in the preceding command:
matrices=[A85(:) A86(:) A87(:) A88(:)];
% Enter probabilities of choosing each matrix below;
% NOTE: sum of penv must equal 1
penv=[.25 .25 .25 .25];	
% Enter farthest future time horizon
tmax=50;						
% Enter number of times to simulate CDF
maxruns=10;
% Enter number of realization of population growth 
% to simulate in each run
numreps=5000;				
% Enter initial population vector
n0=[4264; 3; 30; 16; 25; 5];			
% Enter quasi-extinction threshold, expressed as a density
Nx=500;
% Enter 0 to omit a class and 1 to include it when computing the 
% summed density to compare to the quasi-extinction threshold
sumweight=[1 1 1 1 1 1];
%********************************************************

rand('state',sum(100*clock));	% seed random number generator
s=sqrt(size(matrices,1));		% s=number of classes
cumdist=cumsum(penv);			% CDF for environmental states

Results=[];					% initialize array to store 
                            %     extinction time CDFs
for i=1:maxruns				% calculate ext. CDF "maxruns" times,
	PrExt=zeros(tmax,1);    % initialize counter to record extinctions
	for j=1:numreps			% with "numreps" realizations per run,
   	    n=n0;               % starting at the initial population vector,
   	    for t=1:tmax						% for each future time, 
        	x=sum(rand>=cumdist)+1;			% draw a matrix at random,
        	A=reshape(matrices(:,x),s,s);	% extract the chosen matrix 
                                            %    from "matrices", & 
        	n=A*n;							% project the population 
                                            %    1 year ahead
         	N=sumweight*n;					% compute weighted sum of 
                                            %    current densities
        	if N<Nx							% if quasi-extinct, 
           	    PrExt(t)=PrExt(t)+1;		% update counter
                break;					    % and start new realization
      	    end;
   	    end;  % for t
	end;  % for j
	PrExt=cumsum(PrExt/numreps);	% sum extinctions at each 
                                    %     time to get CDF
	Results=[Results PrExt];		% store the result
end;  % for i

plot(Results)					% plot all "maxruns" CDFs
xlabel('Years into the future')
ylabel('Cumulative probability of quasi-extinction')
meancdf=mean(Results');
figure                          % in a new figure,
plot(meancdf)					% plot mean CDF
xlabel('Years into the future')
ylabel('Cumulative probability of quasi-extinction')