/************************************************************************************************/
/* B. Hill (1975) tail index estimator with nonparametric kernel asymptotic variance estimation */
/* 												*/
/* By Jonathan B. Hill (2010), based on my paper "On Tail Index Estimation for Dependent,	*/
/*					          Heterogeneous Data", Econometric Theory 26    */								
/************************************************************************************************/

/* The variable of interest is assumed to be loaded and called "x", an n-by-1 vector. 		*/

/* In general financial and macroeconomics time series will have to be differenced.  		*/

/* See "On Tail Index Estimation for Dependent, Heterogeneous Data" (2010) by Jonathan B. Hill
   for details on underlying assumptions. In general extremes of {x} must be Near Epoch Dependent,
   covering strong and uniform mixing and NED data each with short or long memor, including ARFIMA,
   FIGARCH, regime switching, nonlinear GARCH, nonlinear ARMA-nonlinear GARCH, Explosive GARCH, 
   IGARCH, and so on. 										*/

/* The code is written only for a Bartlett kernel since extensive simulation study reveals this
   kernel leads to a variance estimator that nearly matches simulation dispersion for a large
   array of time series. 									*/

/* There does not exist a published theory of tail fractile selection m(n) for tail index 
   estimation for dependent data. In the following a somewhat arbitrary window of fractiles
   is used.											*/

new;
cls;
clear all;


band_exp = .225;                 /* kernel var. bandwidth = [m^band_exp]; band_exp <=.25  	*/
n = rows(x);
m = seqa(5,1,n/2);		 /* the set of feasible sample tail fractiles {m(n)} 		*/


	/* ESTIMATION */

{a_hat,var,a_inv,var_a_inv} = a_hat_kern_var(x,m);


	/* OUTPUT */

"       Hill Estimator Output: k(m) = sig(m)/sqrt(m) where sig(m)^2 = kernel estimator of asymptotic variance";"";
"           m           1/a(m)-k(m)        1/a(m)         1/a(m)+k(m)";
m~a_inv-1.96*sqrt(var_a_inv)~a_inv~a_inv+1.96*sqrt(var_a_inv);
graphset;
title("Hill Estimator with 95% Kernel Confidence Bands");
xtics(minc(m),maxc(m),1,1);
xy(m,a_inv-1.96*sqrt(var_a_inv)~a_inv~a_inv+1.96*sqrt(var_a_inv));
end;

 
end;

/* @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ */


/*-------------------------------------------------------------------------
**  a_hat_kern_var.prc
**
**  Purpose:  Derive B. Hill's (1975) tail index estimator a, and Newey-West variance estimator
**                  
**  Usage:     {a_hat,var,a_inv,var_a_inv} = a_hat_kern_var_mse(x1,m);
**
**		a_hat = tail index esimtate over fractile window m
**		var = kernel estimator of asymptotic variance of a_hat
**		a_inv = 1/a_hat
**		var_a_inv = kernel estimator of asymptotic variance of a_inv
**
**  Input:      x1           nx1 time series
**              m           matrix of order indices used for tail index derivation
**
**  Output      a_hat       tail index estimator
**              var         asymptotic kernel variance estimator
**              var_b       bias correct small sample kernel variance estimator
**
**  Jonathan B. Hill (2010)
-----------------------------------------------------------------------*/

proc (4) = a_hat_kern_var(x1,m);
    local   s,t,n,i,j,k,w,a_hat,z,x_m_a,x_ord,n_pos;
    local   var_kern,var,a_inv,ln_xx,ln_xx_lag;
    local   rep,x_m,d1,d2,sq_m,bn;
    local   c_hat1,c_hat2,d_hat,kern_bart,t_seq;

    x1=abs(x1);           /* a_hat for all x: |x| */
    x1=packr(miss(x1,0));
    n = rows(x1);
    x_ord = sortc(x1,1);

    a_hat = zeros(rows(m),cols(m));
    a_inv = zeros(rows(m),cols(m));
    var = zeros(rows(m),cols(m));
    var_a_inv = zeros(rows(m),cols(m));

        /* compute a_hat */

    i=1;
    do until i > rows(m);
        j=1;
        do until j > cols(m);

            if m[i,j] .eq 0;
                goto skip_a;         /* reset m = 1 if ranked m is simply "0" */
            endif;

            a_inv[i,j] = meanc(ln(abs(x_ord[n-m[i,j]+1:n,1])))-ln(abs(x_ord[n-m[i,j],1]));

                /* compute kernal variance estimator for a_hat^(-1) for het., dep. data */
            a_hat[i,j] = 1./a_inv[i,j];
            bn = round(m[i,j].^band_exp);       /* bandwidth */
            ln_xx = (ln(x1[.,1])-ln(x_ord[n-m[i,j],1]));
            ln_xx = ln_xx.*(ln_xx .> 0) - (m[i,j]/n)*a_inv[i,j];
            x_m_a = x_ord[n-m[i,j],1].^a_hat[i,j];

            var_a_inv[i,j] = sumc(ln_xx[.,1].^2)./m[i,j];       /* asymptotic variance: square components */

            k=1;
            do until k > bn;
                ln_xx_lag = packr(lagn(ln_xx,(0|k)));
                w = 1-k/bn;
                var_a_inv[i,j] = var_a_inv[i,j] + 2*w*(ln_xx_lag[.,1]'ln_xx_lag[.,2])./m[i,j];
                k=k+1;
            endo;

            var[i,j] = abs((var_a_inv[i,j]).*(a_hat[i,j]^4));
            var_a_inv[i,j] = var_a_inv[i,j]./m[i,j];        /* scaled for interval estimation*/
            var[i,j] = var[i,j]./m[i,j];               

            skip_a:
            j=j+1;
        endo;
        i=i+1;
    endo;


    retp(a_hat,var,a_inv,var_a_inv);

endp;

/*____________________________________________________________________________________*/