function [score_mat,info_op,info_h] = score_info_GEV_xi_sigma_mu(Y,xsi,sigma,mu)
    % Compute Score Matrix with respect to xsi, sigma, mu for the GEV likelihood
    % Compute the informtion (covariance matrix of the score) using:
    % 1) the outer product of the scores: cov_mat_op
    % 2) the Hessian of the log-likelihood function: cov_mat_h
    % Compute centered numerical derivatives and compute covariance matrix

    % Y is the data matrix .. it is T x k
   
    T = size(Y,1);

    % Numerical derivatives
    % score_mat = score_mat_numerical(Y,[xsi sigma mu]);
    
    % Analytic derivatives
    score_mat = score_mat_analytic(Y,xsi,sigma,mu);
    
    % Outer product
    info_op = score_mat'*score_mat/T;

    % Hessian
    H = NaN(T,3,3);
    deriv_delta = 0.0001;
    small = 0.00001;
    g_p = NaN(T,3);
    g_m = NaN(T,3);
    delta_vec = NaN(1,3);
    x = [xsi sigma mu];
    for i = 1:3
      delta = max([abs(deriv_delta*x(i)) small]);
      x_p = x;
      x_m = x;
      x_p(i) = x(i) + delta;
      x_m(i) = x(i) - delta;
      score_mat_p = score_mat_analytic(Y,x_p(1),x_p(2),x_p(3));
      score_mat_m = score_mat_analytic(Y,x_m(1),x_m(2),x_m(3));
      H(:,:,i) = (score_mat_p - score_mat_m)/(2*delta);
   end
   info_h = -squeeze(mean(H,1));
   info_h = (info_h + info_h')/2;
  end


function score_mat = score_mat_analytic(Y,xsi, sigma, mu)
   T = size(Y,1);
   score_mat = NaN(T,3);
   for t = 1:T
      score_mat(t,:) = score_gev_xsi_sigma_mu(Y(t,:),xsi,sigma,mu)';
   end
  end

function [score] = score_gev_xsi_sigma_mu(Y, xsi, sigma, mu)
% score_gev_xsi_sigma_mu: This computes the score for a single vector observation Y from a GEV observation, where Y is k x 1 and are the order statistics.
% Recoding of UM's fortran code.

  k = length(Y);
  score = zeros(3,1);
  val = zeros(k,3);
  
  Z = (Y - mu)/sigma;
  iZ = (1+xsi*Z)<0;
  if sum(iZ) > 0
      score = zeros(3,1);
      return;
  end
  g = -log1p(xsi*Z)/xsi;
  dl1p = dlog1p(xsi*Z);
  dgxsi = -g/xsi-(Z/xsi).*dl1p;
  dldz=-(1+xsi)*dl1p;
  dldz(end) = dldz(end) + exp(g(end))*dl1p(end);
  val(:,1) = g + (1+xsi)*dgxsi;
  val(:,2) = -Z.*dldz/sigma;
  val(:,3) = -dldz/sigma;
  score = sum(val,1)';
  score(1) = score(1) - exp(g(end))*dgxsi(end);
  score(2) = score(2) - k/sigma;
end

function score_mat = score_mat_numerical(Y,x)
  % Parameters for derivatives 
  T = length(Y);
  deriv_delta = 0.0001;
  small = 0.00001;
  g_p = NaN(T,3);
  g_m = NaN(T,3);
  delta_vec = NaN(1,3);
  for i = 1:3
    delta_vec(i) = max([abs(deriv_delta*x(i)) small]);
    x_p = x;
    x_m = x;
    x_p(i) = x(i) + delta_vec(i);
    x_m(i) = x(i) - delta_vec(i);
    xsi_p = x_p(1);
    sigma_p = x_p(2);
    mu_p = x_p(3);
    xsi_m = x_m(1);
    sigma_m = x_m(2);
    mu_m = x_m(3);
    g_p(:,i) = log(gevpdf(Y,xsi_p,sigma_p,mu_p));
    g_m(:,i) = log(gevpdf(Y,xsi_m,sigma_m,mu_m));
   end
   score_mat = (g_p-g_m)./repmat((2*delta_vec),T,1);
  end