% Stochsens
% This program was modified on 10/28/02 to correct an error in the 
% estimation of the variance and covariance sensitivity and elasticity
% values.

% This is a set of MATLAB commands to estimate the
% sensitivities and elasticities of stochastic growth rate,
% stochastic lambda. 
% You need to have computed and have available the mean matrix
% and also the covariance matrix for the matrix elements:
% the covariance matrix should be arranged in order of, 
% for example, a11, a21, a12, a22. It will be an n by n matrix, 
% where n is the total number of elements in the population
% matrix. 
% This program uses the function secder.m that calculates 
% second derivatives of the deterministic matrix lambda. 

load tortmxs % load the stored population and covariance matrix
mx = meanmx; % set mx equal to the name of stored pop'n matrix 
% if the stored cov matrix is not called covmx, rename it here
n = length(mx);
nn=length(covmx); 
[lambdas,lam,W,w,V,v]= eigenall(mx); % use eigenall.m (BOX 7.1)
sensmx = v*w'/(v'*w); % to get the sensitivities of mx elements
elastmx = (sensmx.*mx)/lam; % find mx element elasticities
vecSS = reshape(sensmx,nn,1); % make the sensitivities a vector

% use Tuljapurkar's approximation for stochastic-lambda
stochlam=exp(log(lam)-((vecSS')*covmx*vecSS)/(2*(lam^2)));

% following loops calculate the sensitivities of
% log(stochlambda) to the mean matrix elements:
  for kk=1:n
     for ll =1:n
        DD=secder(mx,kk,ll);
        vecDD = reshape(DD,nn,1);
        p1 = ...
     (sensmx(kk,ll)/lam)*(1+ ((vecSS')*covmx*vecSS)/(lam^2));
        p2 = ...
(((vecDD')*covmx*vecSS) + ((vecSS')*covmx*vecDD))/(2*(lam^2));
        mnsens(kk,ll) = stochlam*(p1 - p2);
		end; % ll
	end; % kk

% converting to sens. and elast. of stochastic lambda
mnsens = mnsens*stochlam; 
mnelast = mnsens.*mx/stochlam;

disp('sensitivities of mean matrix elements are:')
disp(mnsens)
disp('elasticities of mean matrix elements are:')
disp(mnelast)

% results for variances and covariances of matrix elements
covsens = -stochlam*repmat(vecSS,1,nn).*repmat(vecSS',nn,1)/(lam*lam);
covsens = covsens./(eye(nn)+ones(nn)); %adjust variances (Equ 9.15)
covelast = covsens.*covmx/stochlam;
% the following displays can be used, but the covariance
% matrices are 64 by 64 for desert tortoise and will overrun the screen:
%disp('Sensitivities of element covariances are:')
%disp(covsens)
%disp('Elasticities of element covariances are:')
%disp(covelast)

% the following just shows elast./sens. results for variances:
varsens = reshape( diag(covsens),8,8);
varelast = reshape(diag(covelast),8,8);
disp('sensitivities to the matrix element variances are:')
disp(varsens)
disp('elasticities for matrix element variances are:')
disp(varelast)

covs = reshape( diag(covmx),8,8)