clear;
clc;

% This program computes the order of summability estimator and its
% corresponding subsampling confidence intervals. The partial demeaning 
% case that considers the presence of a constant term is considered. */


% /* Load Data */

fedfunds=importdata('C:\Users\Vanessa\Dropbox\Vanessa\MatlabV\MatlabVV\Data\TaylorRule\fedfunds.txt');

z0t=fedfunds;
[rowsz0t,colsz0t]=size(z0t);
T=rowsz0t;

% /* Deterministics */

% /* Constant Case: Partial Demeaning */

zt=zeros(T,1);
van=1;
while van <= T
zt(van)=z0t(van)-(sum(z0t(1:van))/van);
van=van+1;
end

zt=zt(2:T);
[rowszt,colszt]=size(zt);

T=rowszt;
b=floor(sqrt(T))+1; 

%/* Estimator of the Full Sample */

xdtn=log((1:1:T)');
ydtn=log(cumsum(zt).^2);

ydtn=ydtn-ydtn(1);

bmcon=(xdtn'*xdtn)\(xdtn'*ydtn);

%/******* Subsampling samples generator *******/

numsubsam=T-b+1;

subsamples=zeros(b,T-b+1);

i=1;
while i <= T-b+1;

subsamples(:,i)=zt(i:b+i-1);

i=i+1;
end

[rowsub,colsub]=size(subsamples);

%/*** Subsamples Estimation ***/

Nn=colsub;
bmco=zeros(Nn,1);

j=1;
while j <= Nn

xt=log((1:1:b)');
yt=log(cumsum(subsamples(:,j)).^2);

yt=yt-yt(1);
bmco(j)=log(b)*(abs(((xt'*xt)\(xt'*yt))-bmcon));

j=j+1;
end

%/*** Empirical CDF ***/

[c,y] = ecdf(bmco);
c=c(2:length(c));
y=y(2:length(y));

%/*** Critical Values ***/

alphasig=0.05;
confidence=1-alphasig;
confivector=ones(length(c),1).*confidence;
difconfic=abs(c-confivector);
[mincdifconfic,minindccdifconfic]=min(difconfic);
position=minindccdifconfic;
cnb1malphahalph=y(position,1);


%/*** Confidence Intervals ***/

Ilow=bmcon-((1/log(T))*cnb1malphahalph);
Iup=bmcon+((1/log(T))*cnb1malphahalph);

deltaIlow=(Ilow-1)/2;
deltamcon=(bmcon-1)/2;
deltaIup=(Iup-1)/2;

[deltaIlow deltamcon deltaIup]