*Wouter Gelade, Vincenzo Verardi and Catherine Vermandele
*Version 1.1
cap program drop sqn
program define sqn, eclass


* Qn

version 10.0


if replay()& "`e(cmd)'"=="sqn" {
	ereturn  display
exit
}

syntax varname [if] [in]
marksample touse

tempname  m b V b2 V2 ones

qui count if `touse'
local N=r(N)

ereturn clear
mata: _qn("`varlist'","`touse'")
di ""
di ""
di "The Qn robust scale estimator of Rousseeuw and Croux (1993) is: " `m'
di ""

ereturn post, esample(`touse')
ereturn scalar qn=`m'
ereturn scalar N=`N'

end

mata:
void _qn(string scalar varlist,string scalar touse)
{
real vector X, L, R, P, Q, M, W, Mp, D
real scalar n, diff, C, kw, i, j, m, Z, PT, QT, M0, d

st_view(X=.,.,tokens(varlist),touse)
X=sort(X,-1)

n=rows(X)

L=J(n,1,1)
R=J(n,1,n)

P=J(n,1,.)
Q=J(n,1,.)

M=L
W=J(n,1,1)

diff=1
C=n^2

kw=comb((round(n/2)+1),2)+comb(n,2)+n

while ((C>n)&(diff>0)) {

	for (i=1; i<=n; i++) {
		j=floor(L[i]+(R[i]-L[i])/2)	
		
		if (L[i]>R[i]) {
		M[i]=0
		W[i]=0
		}
		
		else {
		W[i]= R[i] - L[i] + 1
		if (W[i]!=0){
		if (W[i]/2!=round(W[i]/2)) {
		M[i]= At(X,i,j)
		}
		else {
		M[i]= (At(X,i,j)+At(X,i,j+1))/2
		}
		}		
	}
}

m=mm_median(M,W)
	for (i=1;i<=n;i++) {
		Z=PQ(L[i],R[i],m,i,X)
		Q[i]=Z[1,1]	
		P[i]=Z[1,2]
	}
PT=sum(P)
QT=sum(Q:-1)

	if (kw<=PT) {
		if (mreldif(R,P)==0) {
		diff=0
		}
	R=P 
	}

	else if (kw>QT) {
		if (mreldif(L,Q)==0) {
		diff=0
		}
	L=Q
	}

	else {
	diff=0
	}

C=sum(R-L:+1)

} 

//Create a vector of remaining candidates for the median
Mp=0

for (i=1;i<=n;i++) {
	if (W[i,1]:>0) {
		for (j=L[i];j<=R[i]; j++) {
		M0=At(X,i,j)
		Mp=(Mp\M0)
		}
	}
}

Mp=Mp[2..rows(Mp),1]
m=max(Mp)

	for (i=1;i<=n;i++) {
		Z=PQ(L[i],R[i],m,i,X)
		Q[i]=Z[1,1]	
		P[i]=Z[1,2]
	}
Mp=sort(Mp,-1)
maxindex(Mp,1,r1,r2)
PT=sum(Q:-1)-rows(r1)	
	
//PT=sum(P)
//Mp=sort(Mp,-1)
m=Mp[kw-PT]


if (n<=9) {
D=(0.399,0.994,0.512,0.844,0.611,0.857,0.669,0.872)
d=D[1,n-1]
}

else {
	if (round(n/2)==n/2) {
	d=n/(n+3.8)
	}

	else {
	d=n/(n+1/4)
	}
}

m=-m*2.2219
m=m*d

st_local("m",strofreal(m))
}
end

mata:
//Do the binary search to find the bounds for the median
real matrix PQ(real scalar L, real scalar R, real scalar m, real scalar i, real matrix X) {
real scalar P, Q, A1, An, k, Ak

P=0
Q=0

A1=At(X,i,1)
An=At(X,i,rows(X))


if (A1<=m) {
P=0
}
else {
	L=1
	R=rows(X)
	while (R>L+1) {
		k=floor(L+(R-L)/2)
		Ak=At(X,i,k)

		if (Ak<=m) {
		R=k
		}
		
		else if (Ak>m) {
		L=k
		}
	}	
	P=L
}
	
if (An>=m) {
	Q=rows(X)+1
		if(An>m){
		P=Q-1
		}
	}
else if (A1 < m) {
	Q = 1
}
else {
	L=1
	R=rows(X)
		while (R>L+1) {
		k=floor(L+(R-L)/2)
		Ak=At(X,i,k)

		if (Ak<m) {
		R=k
		}
		
		else if (Ak>=m) {
		L=k
		}

		}	
	Q=R
}
	
return((Q,P))
}


mata:
// Implicitly create the MC matrix
real scalar At(matrix X, scalar i, scalar j) {
real scalar A

A=(X[i]-X[rows(X)-j+1])
return(A)
}

end