function [mu,sigma2]=dennisholmes(RawCounts,L,taumax)
% dennisholmes(RawCounts,L,taumax) computes mu and sigma^2 
% using the Dennis-Holmes method;
% RawCounts=vector of census counts
% L=length of running sum
% taumax=maximum time lag in regressions;
% See Holmes, Proceeding of the National Academy of
% Sciences USA 98: 5072-5077 (2001).

NumCounts=length(RawCounts);

% First, compute running sums

RunSum=[];		
for i=1:NumCounts-L+1
   RunSum=[RunSum sum(RawCounts(i:i+L-1))]; 
end;

% Second, compute means and variances of log growth rate at
% each time lag, tau

taus=[];		% array to store time lags
means=[];		% array to store means
vars=[];	    % array to store variances   
for tau=1:taumax                          % For each time lag,
   R=[];
   for i=1:NumCounts-L+1-tau              % compute log growth
      R=[R log(RunSum(i+tau)/RunSum(i))]; % rates using running 
   end;                                   % sums, then compute
   taus=[taus; tau];
   means=[means; mean(R)];                % their means and
   vars=[vars; var(R)];                   % variances
end;

% Third, compute regression slopes of means and vars. vs. tau

% Step 1: center tau, means, and vars by subtracting 
% their means
taus=taus-mean(taus);
means=means-mean(means);
vars=vars-mean(vars);

% Step 2: with centered dependent and independent variables,
% the right-hand sides of the following two lines
% give the slopes of linear regressions of the means and
% variances of the log growth rates on tau, which are 
% our estimates of mu and sigma^2
mu=(taus'*means)/(taus'*taus);      
sigma2=(taus'*vars)/(taus'*taus);