% Function to calculate the Hessian of R0 to the s by 1 parameter vector
% theta where there is only one kind of offspring (stage 1).
%
% Input:
%    U = the transition matrix (n by n matrix)
%    F = the fertility matrix (n by n matrix)
%    s = length of parameter vector theta
%    DUt = D[vecU, theta] (n^2 by s matrix)
%    DFt = D[vecF, theta] (n^2 by s matrix)
%    HUt = H[vecU, theta] (n^2*s by s matrix)
%    HFt = H[vecF, theta] (n^2*s by s matrix)
%
% Function outputs result of Equation (46) (s by s matrix).
%

function H = HR0_theta_1(U, F, DUt, DFt, HUt, HFt)
A = U + F;
[n2,s] = size(DUt);
[n, n] = size(A);
In = eye(n); 
In2 = eye(n^2);
Is = eye(s);
Ins = eye(n*s);

e1 = zeros(n,1);
e1(1) = 1;

F = A-U;
N = inv(In - U);
R = F*N;
Rt = R';

DRt = kron(N', In)*DFt + kron(N', R)*DUt; % Eq (C-27)
DNt = kron(N', N)*DUt; % Eq (C-33)
Dvec1 = kron(kron(In,Kmn(n,n)),In) * kron(In2, In(:))*DNt; % Eq (C-31)
Dvec2 = kron(kron(In,Kmn(n,n)),In) * ( kron(In2, Rt(:))*DNt + kron(N(:), In2)*Kmn(n,n)*DRt ); % Eq (C-32)

HRt = kron(N', Ins)*HFt + kron(In2, DFt')*Dvec1 + kron(kron(N',R),Is)*HUt + kron(In2, DUt')*Dvec2; % Eq (C-30)

H = kron(kron(e1', e1'),Is) * HRt; % Eq (46)

end