clc
load spatial

% ------------------------------------
% BCa-bootstrap interval
% ------------------------------------

alpha = 0.10;
B = 2000;
[blo, bhi, bvals, z0, ahat] = csbootbca(spatial' , 'mom', B, alpha);
disp('BCa-bootstrap intervals')
blo
bhi
range = bhi - blo

% ------------------------------------
% Bootstrap-t intervals
% ------------------------------------

% get sample size
n = length(spatial);

alpha = 0.10;
% number of bootstrap replications
B = 2000;
% get statistic of interest
thetahat = var(spatial);

[bootreps, bootsam] = bootstrp(B, 'var', spatial);
sehats = zeros(size(bootreps));

for i=1:B
    xstar = spatial(bootsam(:, i));
    bvals(i) = var(xstar);
    sehats(i) = std(bootstrp(50, 'var', xstar));
end
zvals = (bootreps - thetahat) ./ sehats;

SE = std(bootreps);

k = B*alpha/2;
szval = sort(zvals);
tlo = szval(k);
thi = szval(B - k);
disp('Bootstrap-t intervals')
blo_tci = thetahat - thi * SE
bhi_tci = thetahat - tlo * SE
range = bhi_tci - blo_tci

% ------------------------------------
% Bootstrap-percentile interval
% ------------------------------------

% get sample size
n = length(spatial);

alpha = 0.10;
% number of bootstrap replications
B =2000;
% get statistic of interest
thetahat = mean(spatial);

bvals = bootstrp(B, 'var', spatial);
k = B * alpha / 2;
sbval = sort(bvals);
disp('Bootstrap-percentile interval')
blo_percent = sbval(k)
bhi_percent = sbval(B - k)
range = bhi_percent - blo_percent