% Simulation of PVA using an unstructured random growth model
% and parameter values for species in Table 2 of Brook et al.
% Supplementary Data
%
% After drawing a set of simulated data, the model parameters
% are estimated, and the model with estimated parameters is 
% simulated to estimate the risk of the population decreasing
% to 10% or less of its initial abundance. This is compared to
% the true risk, estimated by simulating the model with the
% true parameters.

% Select a species by setting the value of jspecies
% Possible values are 1,2,3,...21
jspecies=4; 

% Read in parameter values for the species
load Brook2.txt; 

% Pull out the values for the selected species
nyears=Brook2(jspecies,1); 
mu=Brook2(jspecies,2); 
sigma=Brook2(jspecies,4);

% simulate nyears of data
rdata=mu+sigma*randn(nyears,1);

% calculate sample estimates of mu and sigma
muhat=mean(rdata);
sigmahat=std(rdata); 

disp([mu muhat sigma sigmahat]);

% Simulate 5000 replicates of the population for 50 years
% using the estimated and true values of mu and sigma
%
% Note how this is vectorized to speed things up: each column
% of logxhat is a replicate, so there is no loop over replicates 

logxhat=ones(50,5000);
logx=ones(50,5000);
for j=2:50
	rvals=muhat+sigmahat*randn(1,5000);
	logxhat(j,:)=rvals+logxhat((j-1),:);
	rvals=mu+sigma*randn(1,5000);
	logx(j,:)=rvals+logx((j-1),:);
end

% find the smallest population size in each replicate
minxhat=min(logxhat);
minx=min(logx);

% find fraction below 10% of initial size
estrisk=sum(minxhat<log(0.1))/5000;
truerisk=sum(minx<log(0.1))/5000;
disp([truerisk estrisk]);