close all

start = tic;

% loading data

WTI = xlsread('../../Data/WTI/WTI_Best.xlsx','WTI_Best','B2:B1002');
Brent = xlsread('../../Data/Brent/Brent_Best.xlsx','Brent_Best','B2:B1002');

% initialisation

n = length(WTI);
norm = 2;
steps = n;
ns = floor(linspace(0,1,steps).*n);
sopt = floor(linspace(100,250,150));
pos = linspace(0,0,length(sopt));
pvalue = linspace(0,0,length(sopt));
teststat = linspace(0,0,length(sopt));
result = linspace(0,0,length(sopt));
LowTDC = linspace(0,0,length(sopt));

% calculating pseudo-observations

rcl = zeros(n,1);
[hv,ar] = sort(WTI);
for l = 1:n
    rcl(ar(l)) = l;
end
rcl = rcl/(n+1);

rhh = zeros(n,1);
[hv,ar] = sort(Brent);
for l = 1:n
    rhh(ar(l)) = l;
end
rhh = rhh/(n+1);

% test statistic and p-values

for r = 1:length(sopt)

    EmpLtdc = linspace(0,0,length(ns));
    for j = 1:length(ns)                
        EmpLtdc(j) = 1/sopt(r)*sum((rcl(1:ns(j)) <= sopt(r)/(n+1)).*(rhh(1:ns(j)) <= sopt(r)/(n+1)));           
    end
    
    stat = sqrt(sopt(r)).*(EmpLtdc - EmpLtdc(length(ns))*ns/n);

    stat = stat.*(EmpLtdc(length(ns)).^(-1/2));    

    if norm == 1
        [normstat pos(r)] = max(abs(stat));
        result(r) = normstat > 1.358; % 0.05 critical value
    else
        normstat = mean(stat.^2);
        result(r) = normstat > 0.4633; % 0.05 critical value
    end
    
    teststat(r) = normstat;
    pvalue(r) = pvaluebrownianbridge(teststat(r),0,1000,1000,norm);
    LowTDC(r) = EmpLtdc(length(ns));
       
end

% additional functions to be copied in a separate m-file

function [pbb] = pvaluebrownianbridge(x, s, numsim, steps, norm)

% finds the p-value of the supremum of the Brownian Bridge
% s ratio, numsim number of simulations, steps discretisation for the
% interval [0,1], norm applies either the supremumsnorm or euclidean norm

maxbb = zeros(steps,1);
for i = 1:numsim
    [bb] = brownianbridge(steps);
    if norm ==1
        maxbb(i,1) = max(abs(bb((floor(s*steps)+1):steps,1)));
    else 
        maxbb(i,1) = mean(bb((floor(s*steps)+1):steps,1).^2);
    end
end
pbb = mean(maxbb(:,1) >= x);
end

function [bb] =  brownianbridge(steps)

%simulates a Brownian Bridge between 0 and 1

bm = 1/sqrt(steps)*cumsum(randn(steps,1)); %column vector
bb = bm - (1:steps)'./steps.*bm(steps,1);
%bb = [0; bb]; %normal Brownian Bridge which starts in zero

end

