new;
cls;
clear all;

library pgraph;

/*load x[101,1] = c:\Users\jbhill\Desktop\QML_TTVT\sp500_2005_2009_FZ.txt;*/
/*load x[1041,1] = c:\Users\jbhill\Desktop\QML_TTVT\sp500_2005_2009_FZ_rts.txt;*/
/*load x[1511,1] = c:\Users\jbhill\Desktop\QML_TTVT\sp500_2005_2010_rts.txt;*/
/*load x[1511,1] = c:\Users\jbhill\Desktop\QML_TTVT\nasdaq_2005_2010_rts.txt;*/
/*load x[1469,1] = c:\Users\jbhill\Desktop\QML_TTVT\nikkei_2005_2010_rts.txt;*/
/*load x[1531,1] = c:\Users\jbhill\Desktop\QML_TTVT\dax_2005_2010_rts.txt;*/
load x[1510,1] = c:\Users\jbhill\Desktop\QML_TTVT\dji_2005_2010_rts.txt;

/*load x[1511,1] = c:\Users\jbhill\Desktop\QML_TTVT\sp500_2005_2010_ln_abs.txt;*/
/*load x[1469,1] = c:\Users\jbhill\Desktop\QML_TTVT\nikkei_2005_2010_ln_abs.txt;*/


band_exp = .225;                    /* kernel var. bandwidth = [m^band_exp]; band_exp <=.25  	*/
n = rows(x);
m_fract = .2;

ci = .10;                           /* confidence interval quantile */

m = seqa(5,1,round(m_fract*n));		/* the set of feasible sample tail fractiles {m(n)} 		*/

	/* ESTIMATION */

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

z_quant = cdfni(1-ci/2);


	/* OUTPUT */

"       Hill Estimator Output: k(m) = sig(m)/sqrt(m) where sig(m)^2 = kernel estimator of asymptotic variance";"";
"           m            a(m)-k(m)           a(m)          a(m)+k(m)      1/a(m)-k(m)         1/a(m)         1/a(m)+k(m)";
m~a_hat-z_quant*sqrt(var)~a_hat~a_hat+z_quant*sqrt(var)~a_inv-z_quant*sqrt(var_a_inv)~a_inv~a_inv+z_quant*sqrt(var_a_inv);
"           m            a(m)-k(m)           a(m)          a(m)+k(m)      1/a(m)-k(m)         1/a(m)         1/a(m)+k(m)";

graphset;
begwind;
window (2,1,1);
setwind(1);
title("Inverse Hill Estimator with 95% Kernel Confidence Bands");
xtics(minc(m),maxc(m),1,1);
xy(m,a_inv-z_quant*sqrt(var_a_inv)~a_inv~a_inv+z_quant*sqrt(var_a_inv));
nextwind;
title("Hill Estimator with 95% Kernel Confidence Bands");
xy(m,a_hat-z_quant*sqrt(var)~a_hat~a_hat+z_quant*sqrt(var));
endwind;
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;

/*____________________________________________________________________________________*/