% PROGRAM ricker_corr
%    Uses the stochastic Ricker model N(t+1)=N(t)exp{r(1-N(t)/K)+e(t)}
%    with correlation rho between environmental deviations e(t) in 
%    adjacent years to calculate (by simulation) the probability that
%    a population starting at Nc will fall below the quasi-extinction
%    threshold Nx at or before time tmax

%****************** SIMULATION PARAMETERS ***********************
tmax=50;			% time horizon 
NumReps=10000;		% number of replicate populations to simulate
Nc=10;				% current population size
Nx=5;				% quasi-extinction threshold
r=0.8;				% intrinsic rate of increase
K=10;				% equilibrium population size
sigma2=0.05;		% total environmental variance
rho=0.1;			% first-order environmental autocorrelation
%****************************************************************

sigma=sqrt(sigma2);	
randn('state',sum(100*clock));	% Seed the random number generator.
beta=sqrt(1-rho^2);	% Beta scales new random numbers so that
                    % total environmental variance=sigma2.	
NumExt=0;           % No populations are extinct initially.
for i=1:NumReps,    % For each replicate population, 
    N=Nc;           % start at current population size,
    eold=randn;     % generate an initial environmental deviation,
    for j=1:tmax,                         % then looping over time,
        enew=rho*eold + sigma*beta*randn; % make new (auto-correlated)
										  % environmental deviation and 
        N=N*exp(r*(1-N/K)+enew);	      % use the Ricker model to 
                                          % project one year ahead.
		eold=enew;              % Save old environmental deviation.          
        if N<Nx break; end;     % Stop projecting if threshold is hit.          
	end;
    if N<Nx 
        NumExt=NumExt+1;        % Store number of extinct populations. 
    end;
end;

disp('Probability of hitting quasi-extinction threshold Nx at or before tmax:');
disp(NumExt/NumReps);