% PROGRAM demstoch;
%    Simulates the theta-logistic model with environmental and 
%    demographic stochasticity.  Assumes the population growth  
%    rate at time t is lognormally distributed with mean 
%    exp(r(1-(Nt/K)^theta)) and variance Vd/Nt + sigma^2, 
%    where Nt is population size at time t; r, K, and theta are
%    model parameters; and both the demographic variance Vd and 
%    the environmental variance sigma^2 may be density-
%    dependent, as Vd=ad*Nt+bd and sigma^2=ae*Nt+be

%***************** SIMULATION PARAMETERS ***************
NumPops=10; % number of trajectories to simulate
Nc=10;      % initial population density
r=0.1;      % finite rate of increase
K=15;       % carrying capacity 
theta=1;    % pattern of density dependence
ad=0;       % slope of demographic variance vs. Nt
bd=1;       % intercept of demographic variance vs. Nt
ae=0;       % slope of environmental variance vs. Nt
be=.1;      % intercept of environmental variance vs. Nt
tmax=50;    % duration of simulation
%********************************************************

randn('state',sum(100*clock));	% seed the random number generator
N=Nc*ones(1,NumPops);           % start all trajectories at Nc

Ns=[N];                         % store initial densities
for t=1:tmax                         % for each future time,    
    Mlam=exp(r*(1-(N/K).^theta));    % compute means of lambdas,
    Vlam=(ad*N+bd)./N + ae*N+be;	 % variances of lambdas,
    CV2=Vlam./(Mlam.^2);		     % squared coefficients of 
                                     %    variation of lambdas,
    Mx=log(Mlam)-0.5*log(CV2+1);     % means of Xt's, and 
    SDx=sqrt(log(CV2+1));		     % standard deviations of Xt's, 
    N=N.*exp(Mx+SDx.*randn(1,NumPops));  % then project populations 
                                         % one year ahead using
                                         % lambda=exp(Xt),
    Ns=[Ns; N];                          % and store the results.
end;

% For all trajectories, plot Nt vs t, with a log scale for the y axis
semilogy(Ns)        
xlabel('Time (years)');
ylabel('Population density');
axis([1 tmax .01 100]);