new;

@ computes boostrapped critical value for AR(1), 0.8 @
@ based on the paper version with one new series one and multiple blocks  @
@ REWRITTEN @


nobs=200;
sobs=220;

/*
nobs = 400;
sobs=600;
*/

d1=zeros(100,1);
d2=zeros(100,1);
d3 = zeros(100,1);

g1=zeros(102,1);
g2=zeros(102,1);
g3 = zeros(102,1);



k2 = 1/(2*sqrt(pi));

k1  = 1/(sqrt(2*pi));

term = 1/sqrt(2*pi);

rho = 0.5;



mu = 0; @ location @
c = 1; @ scale @

@ simulate the process 100 times @


rep=1;

do until rep>100;


uu = rndu(sobs,1);

cou = c* tan(pi*(uu - 0.5)) + mu;


y = zeros(sobs,1);
@ initial value 0a @

i = 0;

do until i > (sobs-2);

k = sobs-i-1;
kak = sobs - i;

y[k,1] = rho*y[kak,1] + cou[k,1];

i = i+1;

endo;

@ cut off 100 first and 100 last @

y = trimr(y,100,100);
/*
output file = check.out reset;

print rows(y);

print ;
print nobs;

*/


sy = sortc(y,1);


@ bandwith @

  hh1=1;

num = 0;
den= 0;

met= zeros(5,1);
eta = zeros(5,1);
fu =  zeros(5,1);


gri = zeros(5,1);

gri[1] = sy[60,1];
gri[2] = sy[80,1];
gri[3] = sy[100,1];
gri[4] = sy[120,1];
gri[5] = sy[140,1];

seu=sy[100,1];


@ mean : @

i = 1;
do until i>5;

a = gri[i];

num = 0;
den= 0;

j = 1;
do until j > (nobs-1);
x1 = ((a-y[j]))^2/(2*hh1);
num = num + (exp(-x1)) * y[j+1];
x2 = (a-y[j])^2/(2*hh1);
den = den + exp(-x2);

j = j+1;
endo;
  

met[i] = num/den;

i = i+1;
endo;


@ variance: @

i = 1;
do until i>5;

a = gri[i];

num = 0;
den= 0;

j = 1;
do until j > (nobs-1);
x1 = ((a-y[j]))^2/(2*hh1);
num = num + (exp(-x1)) * ((y[j+1]-met[i])^2);
x2 = (a-y[j])^2/(2*hh1);
den = den + exp(-x2);

j = j+1;
endo;

eta[i] = num/den;

i = i+1;
endo;

@ density: @



i = 1;
do until i> 5;

a = gri[i];

num = 0;

j = 1;
do until j > nobs;

xx1 = (y[j]-a)^2/(2*hh1);
num = num + k1*(exp(-xx1));

j = j+1;
endo;

fu[i] = (1/(nobs*hh1)) * num;

i = i+1;
endo;


res = met-gri;

u1 = k2*eta;
u2 = nobs*hh1* fu;

sed = u1./u2;
rs = sqrt(sed);

st1 = res./rs;


g1[1,1] = st1[2]^2+st1[3]^2+st1[4]^2;

g2[1,1] = sumc(st1^2);

g3[1,1] = st1[3]^2;




ny = zeros(nobs,1);


@ under H0 process is ny @

ny[1]=y[1];

i = 2;

do until i>nobs;



a = y[i-1];

num = 0;
den= 0;

j = 1;
do until j > (nobs-1);
x1 = ((a-y[j]))^2/(2*hh1);

num = num + (exp(-x1)) * y[j+1];
x2 = (a-y[j])^2/(2*hh1);
den = den + exp(-x2);

j = j+1;
endo;

vmet = num/den;



ny[i] = y[i] + y[i-1] - vmet;

i = i+1;
endo;



@ statistic for the H0 process @ 

sy = sortc(ny,1);


@ bandwith @

  hh1=0.5;

num = 0;
den= 0;

met= zeros(5,1);
eta = zeros(5,1);
fu =  zeros(5,1);

gri = zeros(5,1);

gri[1] = sy[60,1];
gri[2] = sy[80,1];
gri[3] = sy[100,1];
gri[4] = sy[120,1];
gri[5] = sy[140,1];




@ mean : @

i = 1;
do until i>5;

a = gri[i];

num = 0;
den= 0;

j = 1;
do until j > (nobs-1);
x1 = ((a-ny[j]))^2/(2*hh1);
num = num + (exp(-x1)) * ny[j+1];
x2 = (a-ny[j])^2/(2*hh1);
den = den + exp(-x2);

j = j+1;
endo;


met[i] = num/den;

i = i+1;
endo;


@ variance: @

i = 1;
do until i>5;

a = gri[i];

num = 0;
den= 0;

j = 1;
do until j > (nobs-1);
x1 = ((a-ny[j]))^2/(2*hh1);
num = num + (exp(-x1)) * ((ny[j+1]-met[i])^2);
x2 = (a-ny[j])^2/(2*hh1);
den = den + exp(-x2);

j = j+1;
endo;

eta[i] = num/den;

i = i+1;
endo;

@ density: @



i = 1;
do until i> 5;

a = gri[i];

num = 0;

j = 1;
do until j > nobs;

xx1 = (ny[j]-a)^2/(2*hh1);
num = num + k1*(exp(-xx1));

j = j+1;
endo;

fu[i] = (1/(nobs*hh1)) * num;

i = i+1;
endo;


res = met-gri;

u1 = k2*eta;
u2 = nobs*hh1* fu;

sed = u1./u2;
rs = sqrt(sed);


st1 = res./rs;


g1[2,1] = st1[2]^2+st1[3]^2+st1[4]^2;

g2[2,1] = sumc(st1^2);

g3[2,1] = st1[3]^2;







@ BOOTSTRAP STARTS HERE @




coup=0.0;

grom=zeros(nobs,5);

i = 1;

do until i>nobs;



if ny[i] > seu;
grom[i,3] = 1;
elseif ny[i]<seu;
grom[i,3] = 2;
endif;

grom[i,2] = ny[i];

i = i+1;
endo;



i = 2;

do until i>nobs;

ni = i-1;

a=grom[i,3];
b=grom[ni,3];

if a /= b;

grom[i,4] = 1;
coup = coup+1;
grom[i,1] = coup;
endif;

i = i+1;
endo;




@ number of blocks is coup @
@ print "coup" coup; @


i=2;

do until i>nobs;

if grom[i,4] == 1;
gag=grom[i,1];
grom[i,5] = gag;
elseif grom[i,4] ==0;
grom[i,5] = grom[i-1,5];
endif;

i=i+1;
endo;


@ print grom[1:100,.]; @



ror = seqa(1,1,coup);


@ PERMUTATION STARTS HERE @

kek = 3;

do until kek > 102;




idx = ceil(coup*rndu(coup,1));

@ idx is the permuted order of blocks @
@ print idx; @

fky=zeros(1,5);

i = 1;

do until i> coup;

ab=idx[i];
aby = selif(grom, grom[.,5] .eq ab);
fky=fky|aby;
i = i+1;
endo;

fky = trimr(fky,1,0);

@ print fky[1:100,.]; @

@ print "rows" rows(fky);@



fay = fky[.,2];

rfk= rows(fky);

if rfk > nobs;
di= rfk - nobs;
fay = trimr(fay,0,di);
elseif rfk < nobs;
du = nobs-rfk;
fay= fay|y[1:du,1];

endif;





sy = sortc(fay,1);


@ bandwith @

  hh1=0.5; 


num = 0;
den= 0;

met= zeros(5,1);
eta = zeros(5,1);
fu =  zeros(5,1);


gri = zeros(5,1);

gri[1] = sy[60,1];
gri[2] = sy[80,1];
gri[3] = sy[100,1];
gri[4] = sy[120,1];
gri[5] = sy[140,1];


/*
 print gri; 
print rows(y);
*/


@ mean : @

i = 1;
do until i>5;

a = gri[i];

num = 0;
den= 0;

j = 1;
do until j > (nobs-1);
x1 = ((a-fay[j]))^2/(2*hh1);

num = num + (exp(-x1)) * fay[j+1];
x2 = (a-fay[j])^2/(2*hh1);
den = den + exp(-x2);

j = j+1;
endo;

met[i] = num/den;

i = i+1;
endo;


@ variance: @

i = 1;
do until i>5;

a = gri[i];

num = 0;
den= 0;

j = 1;
do until j > (nobs-1);
x1 = ((a-fay[j]))^2/(2*hh1);
num = num + (exp(-x1)) * ((fay[j+1]-met[i])^2);
x2 = (a-fay[j])^2/(2*hh1);
den = den + exp(-x2);

j = j+1;
endo;

eta[i] = num/den;

i = i+1;
endo;

@ density: @



i = 1;
do until i> 5;

a = gri[i];

num = 0;

j = 1;
do until j > nobs;

xx1 = (fay[j]-a)^2/(2*hh1);
num = num + k1*(exp(-xx1));

j = j+1;
endo;

fu[i] = (1/(nobs*hh1)) * num;

i = i+1;
endo;


res = met-gri;

u1 = k2*eta;
u2 = nobs*hh1* fu;

sed = u1./u2;
rs = sqrt(sed);

st1 = res./rs;


ake = kek+2;

g1[kek,1] = st1[2]^2+st1[3]^2+st1[4]^2;

g2[kek,1] = sumc(st1^2);

g3[kek,1] = st1[3]^2;






kek = kek +1;

endo;

output file = RBBPdatastat5.out on;
print g1[1,1]~g2[1,1]~g3[1,1];

output file = RBBPdatastat5.out off;


output file = RBBPH0stat5.out on;

print g1[2,1]~g2[2,1]~g3[2,1];

output file = RBBPH0stat5.out off;

gg1=g1[3:102,1];
gg2=g2[3:102,1]; 
gg3=g3[3:102,1];

sg1=sortc(gg1,1);

ci1= sg1[95,1];



sg2 = sortc(gg2,1);
ci2=sg2[95];


sg3=sortc(gg3,1);
ci3=sg3[95];

output file = RBBPcritnew5.out on;
print ci1~ci2~ci3;
output file = RBBPcritnew5.out off;




if g1[1,1] > ci1;
d1[rep,1] =1;
else;
d1[rep,1]=0;
endif;




 if g2[1,1] > ci2;
d2[rep,1] =1;
else;
d2[rep,1]=0;
endif;



if g3[1,1] > ci3;
d3[rep,1] =1;
else;
d3[rep,1]=0;
endif;



rep = rep+1;

endo;



output file = RBPPsize5.out reset;





print sumc(d1);
print sumc(d2);
print sumc(d3);







end;