% Function to calculate the Hessian of lambda to the vector of stage
% specific survival probabilties, sigma.
%
% Input:
%    U = the transition matrix (n by n matrix)
%    F = the fertility matrix (n by n matrix)
%
% Function outputs result of Equation (68) (n by n matrix).
%

function H = Hlambda_sigma(U, F)
A = U + F;
[n,n] = size(A);
In = eye(n);

sigma = sum(U)'; % stage specific survival vector
G = U*diag(sigma)^-1; %transitions conditional on survival, Eq (66)
DAs = kron(In, G) * diag(In(:)) * kron(ones(n, 1), In); % D[vecA, sigma], Eq (72)

B = DAs'*Hlambda_A(A)*DAs; % Eq (73)
H = 1/2*(B+B'); 

end