Matlab code for Calculating short-term elasticities for California sea lions

%This code calculates short-term stochastic elasticities for Ca. sea lions at the site 
%Granito: basic data is from Table 2 of Wielgus et al (2008).  

%First Get average demographic matrix aa from Wielgus et al. (2008)

aa = [0 0 0.530; 0.856 0.200 0; 0 0.389 0.645];

% Generate Beta variables for G_p S_j G_j and S_a and F with above means 

mean_Gp = aa(2,1); mean_Sj = aa(2,2); mean_Gj = aa(3,2); mean_Sa = aa(3,3); mean_FF = aa(1,3);

%We vary variances and get parameters of Beta 
var_Gp = [0.015:0.05:0.25]; var_Sj = [0.083:0.038:0.24]; var_Gj = [0.099:0.074:0.40]; var_Sa = [0.048:0.0312:0.2]; var_FF = [0.021:0.0125:0.071];

a_Gp = mean_Gp*(((mean_Gp*(1-mean_Gp))./var_Gp) - 1);
b_Gp = (1-mean_Gp)*(((mean_Gp*(1-mean_Gp))./var_Gp) - 1);

a_Sj = mean_Sj*(((mean_Sj*(1-mean_Sj))./var_Sj) - 1);
b_Sj = (1-mean_Sj)*(((mean_Sj*(1-mean_Sj))./var_Sj) - 1);

a_Gj = mean_Gj*(((mean_Gj*(1-mean_Gj))./var_Gj) - 1);
b_Gj = (1-mean_Gj)*(((mean_Gj*(1-mean_Gj))./var_Gj) - 1);

a_Sa = mean_Sa*(((mean_Sa*(1-mean_Sa))./var_Sa) - 1);
b_Sa = (1-mean_Sa)*(((mean_Sa*(1-mean_Sa))./var_Sa) - 1);

a_FF = mean_FF*(((mean_FF*(1-mean_FF))./var_FF) - 1);
b_FF = (1-mean_FF)*(((mean_FF*(1-mean_FF))./var_FF) - 1);

% Stochastic Simulations

N = 8000; %Number of iterations
N1= 2; %length(a_Gp); 
N2 = 2; %length(b_Gp);
for i = 1:N1
    for j = 1:N2          
Gp(i,j,:) = betarnd(a_Gp(i),b_Gp(j),[1,N+1]); 
Gj(i,j,:) = betarnd(a_Gj(i),b_Gj(j),[1,N+1]); 
Sj(i,j,:) = betarnd(a_Sj(i),b_Sj(j),[1,N+1]); 
Sa(i,j,:) = betarnd(a_Sa(i),b_Sa(j),[1,N+1]); 
FF(i,j,:) = betarnd(a_FF(i),b_FF(j),[1,N+1]); 
end
end

S = 3;    %Number of life stages
AVE = aa; 

%Initialize Age-vecs := UAGE2 

vec1 = ones(1,S);
vec1 = vec1/sum(vec1);
vec1 = vec1';
UAGE2(:,1) = vec1; 

%A = zeros(S,S);
C1=zeros(9,9); C2=C1;C3=C1;

%Initialize Age-vecs := UAGE11 
vec1 = ones(1,3);  stoch_2 = zeros(3,3,N1,N2,N+1); ave_2 = stoch_2; Sig_2 = stoch_2; 
stoch_3 = zeros(3,3,N1,N2); ave_3 = stoch_3; sig_3 = stoch_3;

vec1 = vec1/sum(vec1);
vec1 = vec1'; unitvec = ones(1,3)';
for jj=1:N1
    for kk=1:N2        
        UAGE11(:,jj,kk,1) = vec1;
        GROW(jj,kk,1) = 1;
    end
end

eeGR_1 = zeros(3,3,N1,N2);  eeaveGR_1 = eeGR_1; eesigGR_1 = eeGR_1;

for jj=1:N1
    for kk=1:N2
for i= 2:N+1        
 
 
A(:,:,jj,kk) = [0 0 mean_FF;Gp(jj,kk,i), Sj(jj,kk,i), 0;0 Gj(jj,kk,i),  mean_Sa];
vec1 = A(:,:,jj,kk) * vec1;    
GROW(jj,kk,i) = sum(vec1);
vec1 = vec1/GROW(jj,kk,i);       
UAGE11(:,jj,kk,i) = vec1;
UU1 = UAGE11(:,jj,kk,i-1);
UU2 = UAGE11(:,jj,kk,i);

for k=1:S
    for l=1:S
  C1((S+1)*k-S,(S+1)*l-S) = A(k,l,jj,kk);
C2((S+1)*k-S,(S+1)*l-S) = AVE(k,l);
C3((S+1)*k-S,(S+1)*l-S) = A(k,l,jj,kk)-AVE(k,l);

%Element perturbation (stoch elas)
stoch_1(:,:,k,l,i) = (unitvec*UU1').* C1(S*k-(S-1):S*k,S*l-(S-1):S*l)/GROW(jj,kk,i);
stoch_2(:,:,jj,kk,i) = stoch_2(:,:,jj,kk,i) + squeeze(stoch_1(:,:,k,l,i));

%Average perturbation (E_mu)
ave_1(:,:,k,l,i) = (unitvec*UU1').* C2(S*k-(S-1):S*k,S*l-(S-1):S*l)/GROW(jj,kk,i);
ave_2(:,:,jj,kk,i) = ave_2(:,:,jj,kk,i) + squeeze(ave_1(:,:,k,l,i));

%Variance perturbation (E_sigma)
Sig_1(:,:,k,l,i) = (unitvec*UU1').* C3(S*k-(S-1):S*k,S*l-(S-1):S*l)/GROW(jj,kk,i);
Sig_2(:,:,jj,kk,i) = Sig_2(:,:,jj,kk,i) + squeeze(Sig_1(:,:,k,l,i));
  
end
end

stoch_3(:,:,jj,kk) = stoch_3(:,:,jj,kk) + squeeze(stoch_2(:,:,jj,kk,i));
ave_3(:,:,jj,kk) = ave_3(:,:,jj,kk) + squeeze(ave_2(:,:,jj,kk,i));
sig_3(:,:,jj,kk) = sig_3(:,:,jj,kk) + squeeze(Sig_2(:,:,jj,kk,i));

eeGR_1(:,:,jj,kk) = stoch_3(:,:,jj,kk)/N; % this is the short-term E^S elasticity 
eeaveGR_1(:,:,jj,kk) = ave_3(:,:,jj,kk)/N; % this is the short-term E^Smu elasticity 
eesigGR_1(:,:,jj,kk) = sig_3(:,:,jj,kk)/N; % this is the short-term E^Ssigma elasticity 
end
 growth_GR(jj,kk) = mean(log(GROW(jj,kk,:))); %Stochastic growth rates
end
end