clc
losses=load('sp500.txt'); %Data for S&P 500, 1950 - 2011, losses only (with positive sign)

q = quantile(losses,.99); % 99% quantile
y = losses(losses>q) - q; % exceedances

[parameters,confidence_interval]=gpfit(y);
xi=parameters(1)    % shape parameter
alpha=1/xi          % tail index
beta=parameters(2)  % scale parameter


[F,yi] = ecdf(y);   % empirical CDF
plot(yi,gpcdf(yi,xi,beta),'-');
hold on; stairs(yi,F,'r'); hold off;
legend('Fitted Generalized Pareto CDF','Empirical CDF','location','southeast');
pause

bins = 0:.005:0.1;
h = bar(bins,histc(y,bins)/(length(y)*.005),'histc');
set(h,'FaceColor',[.9 .9 .9]);
ygrid = linspace(0,max(y),100);
line(ygrid,gppdf(ygrid,xi,beta));
xlim([0,0.1]); xlabel('Exceedance'); ylabel('Probability Density');
pause

p=0.0002;
N=15519; % length of original data (both profits and losses)
VaR10EVTml=(10^xi)*(q+beta/xi*((p*N/length(y))^(-xi)-1))

% Mean excess function, Hill and Tail estimator

N=15519; % length of original data (both profits and losses)

for i=1:length(losses)
    MEF(i)=sum(losses(1:i))/i-losses(i);
    Hill(i)=1/(sum(log(losses(1:i)))/i-log(losses(i)));
    Tail(i)=losses(i)*(i/(N*p))^(1/Hill(i));
    VaR10EVT(i)=10^(1/Hill(i))*Tail(i);
end    

subplot(1,3,1)
scatter(losses(20:end),MEF(20:end))

subplot(1,3,2)
plot(Hill(1:1000))
axis([0 400 0 10])

subplot(1,3,3)
plot(VaR10EVT(1:1000))
axis([0 400 0 2])

mean(VaR10EVT(100:400))