function m_H = fn_hess_uncloglk_N2psi_NKbeta_Nsgm(v_theta, lambda, m_y, m_ysm, m_ysp, a_x, m_Wm, m_Wp)
[N T K] = size(a_x);

v_psim = v_theta(1:N, 1);
v_psip = v_theta((N + 1):(2 * N), 1);
v_beta = v_theta(((2 * N) + 1):((2 + K) * N), 1); % KN x 1: (v_beta_1', ..., v_beta_N')'
v_sgm  = v_theta((end - N + 1):end, 1);

m_beta_tr = reshape(v_beta, K, N); % K x N: (v_beta_1:...:v_beta_N)
m_beta = m_beta_tr';

m_Psim = diag(v_psim);
m_Psip = diag(v_psip);
v_sgmsq = v_sgm.^2;
v_sgm3h = v_sgm.^3;
v_sgm4h = v_sgm.^4;

m_A = eye(N) - (m_Psim * m_Wm) - (m_Psip * m_Wp);
invA = inv(m_A);
m_Qm = m_Wm * invA;
m_Qp = m_Wp * invA;

% compute residuals
a_beta = zeros(N, T, K);
for k = 1:K
     v_betak = m_beta(:, k);
     a_beta(:, :, k) = repmat(v_betak, 1, T);
end
m_psim = repmat(v_psim, 1, T);
m_psip = repmat(v_psip, 1, T);
m_eps = m_y - (m_psim .* m_ysm) - (m_psip .* m_ysp) - squeeze(sum(a_beta .* a_x, 3));

m_H11mm = (m_Qm .* m_Qm') + diag(sum(m_ysm .* m_ysm, 2) ./ T ./ v_sgmsq);
m_H11mp = (m_Qm .* m_Qp') + diag(sum(m_ysm .* m_ysp, 2) ./ T ./ v_sgmsq);
m_H11pp = (m_Qp .* m_Qp') + diag(sum(m_ysp .* m_ysp, 2) ./ T ./ v_sgmsq);
m_H13m = 2 .* diag(sum(m_ysm .* m_eps, 2) ./ T ./ v_sgm3h);
m_H13p = 2 .* diag(sum(m_ysp .* m_eps, 2) ./ T ./ v_sgm3h);
m_H33 = diag(-1 ./ v_sgmsq + (3 .* sum(m_eps .* m_eps, 2) ./ T  ./ v_sgm4h));

m_H12m = zeros(N, (N * K));
m_H12p = zeros(N, (N * K));
m_H22  = zeros(N * K);
m_H23 = zeros((N * K), N);

for i = 1:N
     ind = ((i - 1) * K + 1):(i * K);
     v_ysmi = m_ysm(i, :)';
     v_yspi = m_ysp(i, :)';
     m_Xi = permute(a_x(i, :, :), [2 3 1]); % T x K
     v_epsi = m_eps(i, :)';
     sgmsqi = v_sgmsq(i, 1);
     sgm3hi = v_sgm3h(i, 1);

     m_H12m(i, ind) = (v_ysmi' * m_Xi) ./ T ./ sgmsqi;
     m_H12p(i, ind) = (v_yspi' * m_Xi) ./ T ./ sgmsqi;
     m_H22 (ind, ind) = (m_Xi' * m_Xi) ./ T  ./ sgmsqi;
     m_H23(ind, i) = 2 .* (m_Xi' * v_epsi) ./ T  ./ sgm3hi;
end
m_H11 = [m_H11mm  m_H11mp; ...
         m_H11mp' m_H11pp];

m_H12 = [m_H12m; ...
         m_H12p];

m_H13 = [m_H13m; ...
         m_H13p];

m_H = -[[m_H11  m_H12  m_H13]; ...
        [m_H12' m_H22  m_H23]; ...
        [m_H13' m_H23' m_H33]];
