function [score_mat,info_op,info_h] = score_info_Exceeedance_GEV_tvtau(Y,nobs_vec,tau_vec,xsi,sigma,mu);
    % Compute Score Matrix with respect to xsi, sigma, mu in the Exceedance Model with potentially time varying threshold (tau)
        % Compute the informattion (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
       
        T = length(Y);
    
        % Analytic derivatives
        % score_mat = score_mat_numerical(Y,nobs_vec,tau_vec,[xsi sigma mu]);
        score_mat = score_mat_analytic(Y,nobs_vec,tau_vec,[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_numerical(Y,nobs_vec,tau_vec,x_p);
          % score_mat_m = score_mat_numerical(Y,nobs_vec,tau_vec,x_m);
          score_mat_p = score_mat_analytic(Y,nobs_vec,tau_vec,x_p);
          score_mat_m = score_mat_analytic(Y,nobs_vec,tau_vec,x_m);
          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,nobs_vec,tau_vec,x)
        T = size(Y,1);
        score_mat = NaN(T,3);
        xsi = x(1);
        sigma = x(2);
        mu = x(3);
        for t = 1:T
           k = nobs_vec(t);
           threshold = tau_vec(t);
           score_mat(t,:) = score_exceedence_gev_xsi_sigma_mu(Y(t,:),k,threshold,xsi,sigma,mu)';
        end
    end
    
    function [score] = score_exceedence_gev_xsi_sigma_mu(Y, k, threshold, xsi, sigma, mu)
        % score_exceedanc_gev_xsi_sigma_mu: This computes the score for a single vector observation Y from a exceedabce GEV observation.
        
        
          score = zeros(3,1);
          val = zeros(k,3);

          Z = (threshold - mu)/sigma;
          g = -log1p(xsi*Z)/xsi;
          dl1p = dlog1p(xsi*Z);
          dgxsi = -g/xsi-(Z/xsi)*dl1p;
          dldz = dl1p*exp(g);
          score(1)= -exp(g)*dgxsi;
          score(2)= -Z*dldz/sigma;
          score(3)= -dldz/sigma;

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

    function score_mat = score_mat_numerical(Y,nobs_vec,tau_vec,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);
          g_p(:,i) = llf_exceedance_gev(Y,nobs_vec,tau_vec,x_p(1),x_p(2),x_p(3));
          g_m(:,i) = llf_exceedance_gev(Y,nobs_vec,tau_vec,x_m(1),x_m(2),x_m(3));
         end
         score_mat = (g_p-g_m)./repmat((2*delta_vec),T,1);
        end

        function llf_vec = llf_exceedance_gev(Y_data,nobs_vec,tau_vec,xi,sigma,mu)
            % Minus log-likelihood from exceedance likelihood using GEV parameterization

            T = size(Y_data,1);
            llf_vec = zeros(T,1);
            X = (tau_vec-mu)/sigma;
            xi_X = xi*X;
            % Compute log-likelihood
            for t = 1:T
                k = nobs_vec(t);
                ln_1 = -(log(factorial(k))+k*log(sigma));
                x  = (tau_vec(t)-mu)/sigma;
                x1 = -((1+xi*x).^(-1/xi));
                x2 = 0;
                for ik = 1:k
                    tmp = xi*((Y_data(t,ik)-mu)/sigma);
                    tmp = log1p(tmp);
                    x2 = x2-(1+1/xi)*tmp;
                end
                llf_vec(t) = ln_1+x1+x2;
            end
            
          end