% Function to calculate the Hessian of R0 to the s by 1 parameter vector
% theta.
%
% 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 (44) (s by s matrix).
%

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

F = A-U;
N = inv(In - U);
R = F*N;
Rt = R';
[R0, wr, vr] = domeig(R);

DRt = kron(N', In)*DFt + kron(N', R)*DUt; % Eq (C-27)
HR0R = Hlambda_A(R); % H[R0,vecR], using Eq (24)

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)

B = kron(kron(wr', vr'), Is)*HRt + DRt'*HR0R*DRt; % Eq (57)
H = 1/2*(B + B');

end