clear all
load cook

n = 11944;
bs = 10000;

replica = zeros(n,bs);

disp('generate random numbers')
for i=1:bs,
    vrand = random('unid',n,n,1);
    replica(:,i) = cook(vrand);
end

%rolling standard deviation, skewness and first-order autocorrelation
disp('generate sample moments')
vstdev = zeros(n,bs);
vskew = zeros(n,bs);
vautocor = zeros(n,bs);

w = 250; %window half-width

for i=1:bs,
    disp(i)
    for j= w+1:n-w,
        vstdev(j,i) = std(replica(j-w:j+w,i));
        vskew(j,i) = skewness(replica(j-w:j+w,i));
    end
    for j= w+2:n-w,
        vcc = corrcoef(replica(j-w:j+w,i),replica(j-1-w:j-1+w,i));
        vautocor(j,i) = vcc(2,1);
    end
end

cookstd = zeros(n,1);
cookskw = zeros(n,1);
cookcor = zeros(n,1);

stdavg = zeros(n,1);
std025 = zeros(n,1);
std975 = zeros(n,1);
skwavg = zeros(n,1);
skw025 = zeros(n,1);
skw975 = zeros(n,1);
coravg = zeros(n,1);
cor025 = zeros(n,1);
cor975 = zeros(n,1);

%average and 95% confidence interval
disp('generate mean and bounds')
for i= w+1:n-w,
    cookstd(i) = std(cook(i-w:i+w));
    cookskw(i) = skewness(cook(i-w:i+w));
end
for i= w+2:n-w,
    cookcor(i) = corr(cook(i-w:i+w),cook(i-1-w:i-1+w));
end

for i = w:n-w,
    stdavg(i) = mean(vstdev(i,:));
    std025(i) = prctile(vstdev(i,:),2.5);
    std975(i) = prctile(vstdev(i,:),97.5);
    skwavg(i) = mean(vskew(i,:));
    skw025(i) = prctile(vskew(i,:),2.5);
    skw975(i) = prctile(vskew(i,:),97.5);
end

for i = w+1:n-w,
    coravg(i) = mean(vautocor(i,:));
    cor025(i) = prctile(vautocor(i,:),2.5);
    cor975(i) = prctile(vautocor(i,:),97.5);
end

%autocorrelation function
nac = 20;
cookacf = zeros(nac,1);
vact = zeros(nac,bs);

for j=1:nac,
    cookacf(j) = corr(cook(1:n-j),cook(1+j:n));
end

for i=1:bs,
    for j=1:nac,
        vacf(j,i) = corr(replica(1:n-j,i),replica(1+j:n,i));
    end
end

for i = 1:nac,
    acfavg(i) = mean(vacf(i,:));
    acf025(i) = prctile(vacf(i,:),2.5);
    acf975(i) = prctile(vacf(i,:),97.5);
end

clear v*
save cook