% PROGRAM stoc_log_lam
%	Calculates stochastic log growth rate
% 	by simulation and by Tuljapurkar's approximation;
% 	penv allows some matrices to be chosen more often
% 	than others.
%   REQUIRES THE FUNCTION eigenall.m

%************* 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];	
% number of time intervals to simulate - should be large
maxt=50000;			
%********************************************************

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

% Calculate mean matrix Abar and covariance matrix C, 
% taking into account the fact that all matrices may not 
% be equally likely; uses the fact that 
% Cov(x,y)=E(x*y)-E(x)*E(y), where E() denotes
% expectation
Abar=zeros(s^2,1);		
Exy=zeros(s^4,1);
for i=1:numenvts 
   A=matrices(:,i);
   Exy=Exy+penv(i)*kron(A,A); % kron is Kronecker tensor product function
   Abar=Abar+penv(i)*A;       % Abar is sum of A(i)'s weighted by penv(i)
end
C=(Exy-kron(Abar,Abar))*numenvts/(numenvts-1); % Cov(x,y)=E(x*y)-E(x)*E(y)
C=reshape(C,s^2,s^2);
Abar=reshape(Abar,s,s);

% Calculate dominant eigenvalue, lambdabar, and
% matrix of eigenvalue sensitivities, S, for mean matrix Abar;
% requires user-defined function "eigenall"
[lambdas,lambdabar,W,w,V,v]=eigenall(Abar);			
S=v*w'/(v'*w);								

% Calculate stochastic log lambda by simulation
n=w;				% Start at stable distribution of Abar, with 
					%      total population size = 1
for t=1:maxt                        % For each year,
   x=sum(rand>=cumdist)+1;			% choose a matrix at random,
   A=reshape(matrices(:,x),s,s);	% extract the chosen matrix,
   n=A*n;							% project 1 year ahead,
   N=sum(n);		% sum n to get new total pop. size,
   r(t)=log(N);	    % calculate log growth rate, and
   n=n/N;			% renormalize so sum(n)=1 to avoid pop. 
   					%   sizes too small or too large for computer to handle.
end
loglsim=mean(r);			   % simulated stochastic log growth rate
dse=1.96*sqrt(var(r)/maxt);	   % standard error of loglsim
CL1=[loglsim-dse loglsim+dse]; % approx. 95% confidence interval - 
							   %    see Caswell 2001, eq. 14.62
lamsim=exp(loglsim);  		   % simulated stochastic growth rate
CL2=exp(CL1);				   % confidence limits on lamsim

% Calculate stochastic log lambda using Tuljapurkar's approximation
Svec=S(:);		
tau2=Svec'*C*Svec;
loglams=log(lambdabar)-tau2/(2*lambdabar^2);  
lams=exp(loglams);

disp('Stochastic log growth rate computed by simulation');
disp('(followed by 95% confidence interval):');
disp(loglsim);
disp(CL1);
disp('Stochastic growth rate computed by simulation');
disp('(followed by 95% confidence interval):');
disp(lamsim);
disp(CL2);
disp('Stochastic log growth rate computed by Tuljapurkar''s approximation:');
disp(loglams);
disp('Stochastic growth rate computed by Tuljapurkar''s approximation:');
disp(lams);