mata:
 real matrix nwest (real vector x,real scalar lag)
{
real scalar t, xom,xomm,temp11
   t = rows(x)
   xom= (x'*x)/t
   xomm = 0;
  for (i=1; i<=lag; i++) { 
      temp11 = sum(x[1+i::t-1,.]:*x[1::t-1-i,.])
      xomm = (1-(i/(lag+1)))*(temp11/t) + xomm
                            }  
     xom = xom +2*xomm

return(xom)
}
mata mosave nwest(), replace
end

mata:
real matrix logtm (real matrix y, real scalar mk)
{
real scalar t, trat, n
real matrix xx, be, re, sig, my, y1, yy, ly1, rht, trd, trd1, re1
t = rows(y);
n=cols(y);
lag=3

yy=y;

my = mean(yy');
my = my'
my=J(1,n,my)
y1=yy:/my;

y11 = (y1 - J(t,n,1)):*(y1 - J(t,n,1));
ly1 = mean(y11')
ly1 = ly1[.,1]:/ly1;
ly1 = ln(ly1');

rht = ly1[mk+1::t];
trd = 1::t;
trd1 = ln(trd[mk+1::t,.]);

rht = rht - 2*ln(trd1);
xx= (J(t-mk,1,1), trd1);

be = invsym(xx'xx)*xx'*rht;
re = rht - xx*be;
re1=vec(re)
sig = nwest(re1,lag);
sig = sig*invsym(xx'xx);
sig = sig[2,2];
trat = be[2,1]/sqrt(sig);
return(trat)
}
mata mosave logtm(), replace

end 

mata:

real matrix fixedcore(real matrix x, real scalar mk)
{
x2=x

T=rows(x2)
lasty = x2[T,.]'
idx=order(-1*lasty,1)
idx=idx'
x1= x2[2::T,idx]
T=rows(x1)

ij=1
ii=1
while(ij<999){
    ww=(x1[.,ii], x1[.,ii+1])
	ti1=logtm(ww,mk)
	if(ti1<-1.65){
	    ij=ij+1
		ii=ij+1
	}
	else {
	    ij=1000
	}
}

tt1=J(10-ii,2,0)

for(ii2=ii+1; ii2<=10; ii2++){
    ww=x1[.,ii::ii2]
	tt1[ii2-ii,1]=ii2
	tt1[ii2-ii,2]=logtm(ww,mk)
}

v = J(0,1,.) ; w = J(0,2,.)
maxindex(tt1[,2],1,v,w)
maxt=tt1[v,1]
coreidx=idx[.,ii::maxt]

return(coreidx)
}
mata mosave fixedcore(), replace
end

mata:
real matrix setdiff2 (real rowvector a, real rowvector b)
{
    sa=cols(a)
	sb=cols(b)
    j=0
	p=1
	q=1
    output = J(1, sa,0)
    for(mm=1; mm<=sa; mm++){
	    for(nn=1; nn<=sb; nn++)
		{
		    if(a[1,mm]!=b[1,nn]){
			    j=j+1
			}
			else{
			    j=0
				nn=sb
				q=q+1
			}
			if(j==sb){
			    output[1,p]=a[1,mm]
				p=p+1
				q=q+1
				j=0
			}
		}
	}
	output=select(output,(output:>0))
return(output)
}
mata mosave setdiff2(), replace
end

mata

real matrix sortclub_(real matrix x2, real scalar mk, real scalar ccc, real scalar rd, real rowvector fc)
{

T=rows(x2)
n=cols(x2)
idx=x2[1,.]
x1 = x2[2::T,.]


if(rd==1){
	ii=1
	coreidx=fc
	core=x1[,coreidx]
	dd=x2
	T=rows(dd)
	lasty=dd[T,.]'
	idx=order(-1*lasty,1)
	idx=idx'
	icore=setdiff2(idx,coreidx)
	ccore=x1[,icore]
	T=rows(x1)
	xx=x1[,idx]
}
else{
T=rows(x1)

ij=1
ii=1
while(ij<999){
    ww=(x1[.,ii], x1[.,ii+1])
	ti1=logtm(ww,mk)
	if(ti1<-1.65){
	    ij=ij+1
		ii=ij+1
	}
	else {
	    ij=1000
	}
}

tt1=J(10-ii,2,0)

for(ii2=ii+1; ii2<=10; ii2++){
    ww=x1[.,ii::ii2]
	tt1[ii2-ii,1]=ii2
	tt1[ii2-ii,2]=logtm(ww,mk)
	if(tt1[ii2-ii,2]<-1.65){
	    tt1=tt1[1::ii2-ii,]
	    break
	}
}

v = J(0,1,.) ; w = J(0,2,.)
maxindex(tt1[,2],1,v,w)
maxt=tt1[v,1]
coreidx=idx[.,ii::maxt]
core=x1[.,ii::maxt]

ncore = cols(coreidx)
icore = J(1, n-ncore,0)
ccn=cols(icore)
ccore = J(T,ccn,0)

if(ii>1) {
	for(kk=1; kk<=ii-1; kk++){
		icore[1,kk]=idx[1,kk]
		ccore[.,kk]=x1[.,kk]
	}
	for(kk=ii; kk<=ccn; kk++){
		icore[1,kk]=idx[1,kk+ncore]
		ccore[.,kk]=x1[.,kk+ncore]
	}
}
if(ii==1) {
	icore = idx[1,ncore+1::n]
	ccore=x1[.,ncore+1::n]
}
xx=x1
}

ccn=cols(icore)
idxclub=coreidx
club = core
for(jj=1; jj<=ccn; jj++){
	yy = (core, ccore[.,jj])
	t5 = logtm(yy,mk)
	if(t5 > ccc){
		idxclub = (idxclub, icore[1,jj])
		club = (club, ccore[.,jj])
	}
}

sizeclub = 	cols(idxclub)
sizecomp= n-sizeclub
idxcomp = J(1,sizecomp,0)
ccomp = J(T,sizecomp,0)

j=0
p=1
q=1
for(mm=1; mm<=n; mm++){
	for(nn=1;nn<=sizeclub;nn++){
		if(idxclub[1,nn]!=idx[1,mm]){
			j=j+1
		}
		else{
			j=0
			nn=sizeclub
			q=q+1
		}
		if(j==sizeclub){
			idxcomp[1,p]=idx[1,mm]
			ccomp[,p]=xx[,q]
			p=p+1
			q=q+1
			j=0
			}
	}
}

tr=logtm(club,mk)
tr1=logtm(ccomp,mk)
xcomp=(idxcomp\ccomp)

all=J(5+T,n,.)
all[1,1::sizeclub]=idxclub
all[2,.]=J(1,n,tr)
all[3,.]=J(1,n,sizeclub)
all[4,.]=J(1,n,tr1)
if(sizecomp>=1){
all[5::T+5,1::sizecomp]=xcomp
}
return(all)
}
mata mosave sortclub_(), replace
end

mata

real matrix sortall_(real matrix y, real scalar mk, real scalar ccc, real rowvector fc)
{

T=rows(y)
n=cols(y)
icm=J(30,1,0)
x1=y

idx1=J(1,1,1::n)'
x2=(idx1\x1)

for(jj=1;jj<=30;jj++){
	alldata=sortclub_(x2,mk,ccc,jj, fc)
	nclub= alldata[3,1]
	ncomp=cols(alldata)-nclub
	ico = alldata[1,1::nclub]
	tc = alldata[2,1]
	to = alldata[4,1]
	iout1 = alldata[5,1::ncomp]
	xcomp = alldata[5::T+5,1::ncomp]

	if(tc<-1.65){
		printf("Error: Club %g does not converge. \n", jj)
		printf("Increase critical value \n")
		break
	}
		
	if (nclub==0){
		break
	}
	
	icm[jj,.] = nclub

	if(jj==1){
		ic=ico
	}
	else{
	ic = (ic, ico)
	}

	if(to==.){
		to=-999
	}
	
	if(to>-1.65){
		printf("There are no divergent members at T=%g \n \n", T)
		icm[jj+1,1]=ncomp
		ic = (ic, iout1)
		break
	}
	
	if(ncomp<10){
		for(ll=1; ll<=ncomp; ll++){
			printf("Divergent member at T=%g : \n", T) 
			printf("%g \n", iout1[1,ll])
		}
		break
	}
	x2=xcomp
}

n_icm = sum(icm:!=0) 
n_ic = cols(ic)
comb=J(2,n_ic,0)
comb[1,1::n_icm]=icm[1::n_icm]'
comb[2,1::n_ic]=ic	

return(comb)
}
mata mosave sortall_(), replace
end

mata

void autosort_(real matrix xx1, real scalar mk, real scalar ccc, real rowvector fc)
{

x2=xx1
T=rows(x2)

trat=logtm(x2,mk)
if(trat>-1.65){
    printf("The entire dataset converges, t=1:%g \n \n", T)
	printf("t-ratio: %g \n", trat)
}
else{

comb = sortall_(x2,mk,ccc,fc)
icm=comb[1,.]'
ic=comb[2,.]

maxc=max(icm)
n_icm = sum(icm:!=0) 
w=J(n_icm,maxc,0)

j=1
for(i=1;i<=n_icm;i++){
    bb=icm[i,1]
    w[i,1::bb]=ic[1,j::bb+j-1]
	j=icm[i,1]+j
}

if(n_icm==1){
	rcon=x2[,w]
    tra = logtm(rcon,mk)
    printf("The entire dataset, minus the divergent members, converges, t=1:%g \n",T)
	printf("t-ratio: %g \n \n", tra)	
}
else{
	
for(jj=1; jj<=6; jj++){
	printf("Merging clubs round %g, t=1:%g \n \n", jj, T)
	printf("--------------------------------------------------------------------")
	printf("--------------------------------------------------------------------\n")
	printf("--------------------------------------------------------------------")
	printf("--------------------------------------------------------------------\n")
	nic=ic


	for(il=1;il<=n_icm;il++){
		in=w[il,]
		in=select(in, (in :> 0))
		cc=x2[,in]
		t=logtm(cc,mk)
		
		printf("convergent club %g members: \n \n", il)
		sort(in',1)'
		printf("t-ratio: %g \n \n", t)
	}
	
	v=2
	w1=w[1,]
	w1=select(w1,(w1:>0))
	w2=w[2,]
	w2=select(w2, (w2:>0))
	z=(w1, w2)
	n=cols(x2)
	mm=J(n_icm,n,0)
	p=1
	mm[,1::maxc]=w[,]
	
	for(i=1; i<=n_icm-1;i++){
	    zz=x2[,z]
		tz=logtm(zz,mk)
		nic=setdiff2(nic,z)
		
		if(tz>-1.65){
		    mm[i+1,1::cols(z)]=z[,]
			p=p+1
			
			if(v>2){
			    if(v==3){
				    printf("\n club %g converges with club %g and club %g, \n", i+1, i-1, i)
					printf(" \n merged club %g, %g, and %g members: \n", i-1, i, i+1)
					sort(z',1)'
					printf("\n t-ratio: %g \n \n", tz)
				}
				if(v==4){
				    printf("\n club %g converges with club %g, club %g, and club %g, \n", 
					i+1, i-2, i-1, i)
					printf(" \n merged club %g, %g, %g, and %g members: \n", i-2, i-1, i, i+1)
					sort(z',1)'
					printf("\n t-ratio: %g \n \n", tz)
				}
				if(v==5){
				printf("\n club %g converges with club %g, club %g, club %g and club %g, \n", 
					i+1, i-3, i-2, i-1, i)
					printf(" \n merged club %g, %g, %g, %g and %g members: \n", 
					i-3, i-2, i-1, i, i+1)
					sort(z',1)'
					printf("\n t-ratio: %g \n \n", tz)
				}
				if(v==6){
				printf("\n club %g converges with club %g, club %g, club %g club %g and	club %g, \n", 	i+1, i-4, i-3, i-2, i-1, i)
					printf(" \n merged club %g, %g, %g, %g and %g members: \n", 
					i-4, i-3, i-2, i-1, i, i+1)
					sort(z',1)'
					printf("\n t-ratio: %g \n \n", tz)
				}
			}
			else{
			    printf("\n club %g converges with club %g \n", i, i+1)
				printf("\n merged club %g and %g members: \n", i, i+1)
				sort(z',1)'
				printf("\n t-ratio: %g \n \n", tz)
			}
			
			if(cols(nic)>0){
			    zc=x2[,nic]
				tc=logtm(zc,mk)
				
				if(tc>-1.65){
				    printf("the remaining clubs can be merged \n \n")
					printf("\n merged club members:\n")
					sort(nic',1)'
					printf("\n t-ratio: %g \n \n", tc)
					break
				}
				else{
				    printf("the remaining clubs do not converge \n \n")
				}
			}
			
			if(i+1<=n_icm-1){
			    wi=w[i+2,]
				wi=select(wi, (wi :> 0))
			    z=(z, wi)
			}
			v=v+1
		}
		else{
		    if(v>2){
			    if(v==3){
				    printf("\n club %g does not converge with clubs %g and %g \n", i+1, i-1, i)
				}
				if(v==4){
				    printf("\n club %g does not converge with clubs %g, %g, and %g \n", i+1, i-2, i-1, i)
				}
				if(v==5){
				    printf("\n club %g does not converge with clubs %g, %g, %g, and %g \n", i+1, i-3, i-2, i-1, i)
				}
				if(v==6){
				    printf("\n club %g does not converge with clubs %g, %g, %g, %g, and %g \n", i+1, i-4, i-3, i-2, i-1, i)
				}
			}
			else{
			    printf("\n club %g does not converge with club %g \n", i, i+1)
			}
			if(i+1<=n_icm-1){
			    wi1=w[i+1,]
				wi1=select(wi1,(wi1:>0))
				wi2=w[i+2,]
				wi2=select(wi2,(wi2:>0))
			    z=(wi1, wi2)
			}
			v=2
		}
}
if(p==1){
	printf("\n no other clubs can be merged \n \n")
	break
}

qq=mm
for(i=1; i<=n_icm-1; i++){
    if(mm[i,1]==mm[i+1,1]){
	    qq[i,]=J(1,n,0)
	}
	}

qq=select(qq, rowsum(abs(qq)))
icm=rowsum(qq:!=0)
w=qq[,1::max(icm)]
n_icm=rows(icm)
maxc=max(icm)

}
}
}
}
mata mosave autosort_(), replace

end