/*This file estimates several triple index models using Recursive Differencing
and Residual Controls */

storebRd   = {};
storebRdRc = {};

library pgraph;
Library Maxlik;

N = 2000;

ii =1;
do while ii .le 100;
Print ii;

seed1 =1234+ii;
seed2 =76+ii;
seed3 = 2043 + ii;
seed4 =5731 + ii;
seed5 = 8954 + ii;
seed6 = 58993 + ii;
seed7 = 14121412 + ii;

x1 = rndns(2*N,1,seed1);
X2 = rndns(2*N,1,seed2);
X3 = rndns(2*N,1,seed3);
x4  = rndns(2*N,1,seed5);

s1 = abs(x1) .le 3;
s2 = abs(x2) .le 3;
s3 = abs(x3) .le 3;
s4 = abs(x4) .le 3;
x1 = selif(x1,s1);
x2 = selif(x2,s2);
x3 = selif(x3,s3);
x4 = selif(x4,s4);
x1 = x1[1:N];
x2 = x2[1:N];
x3 = x3[1:N];	
x4 = x4[1:N];
u =   rndns(N,1,seed4);
v1 = x1 + x4;
v2 = X2 - X4;
v3 = X3 + X4;
True = 1|-1|1;
v1 = v1/sqrt(2);
v2 = v2/sqrt(2);
v3 = v3/sqrt(2);
X = X1~x2~x3~X4;
v = v1~v2~v3;


/*Quadratic*/
B = v1;
M = v1.*v3;
L = v2^2;
L = L/1.3447592;
E = B + M + L;
Y = E + u;


/*Cubic
B = v1;
M = v1.*v3;
L = v2^3;
L = L/3.6820149;
E = B + M + L;
E = 2*E/1.4;
Y = E + u;*/


/*Exponential
B = v1;
M = v1.*v3;
L = exp(v2)/2;
E = B + M + L;
E = 2*E/1.6300842 ;
Y = E + u;*/

/*Sin
B = v1;
M = v1.*v3;
L = sin(v2);
L = L/ 0.65702512; 
E = B + M + L;
E = 2*E/1.3929493;
Y = E + u;*/

deltasm=0.001; 
deltadjust = 1/30;
@t1  =  smtrm(x[.,1], deltasm,.98).*smtrm(x[.,2],deltasm,.98)
.*smtrm(x[.,3],deltasm,.98).*smtrm(x[.,4],deltasm,.98);@

t1  =  trimdexpan(X1~X2~X3~x4,deltasm,.98);
t2  =  trimd(X1~X2~X3~x4,.97);

D=Y~X; 

start1 = 0;
start2 = 0;
start3 = 0;
start=start1|start2|start3;
@gstart=gridmaxStart(3,.5,1,d,start);@
gstart = 1|-1|1;
maxset;
/*_max_linesearch = 3;
{bRd1,f,g,cRd,ret}  = maxlik(D,0,&slsRd1,start1);
ev1 = X[.,1] + X[.,4]*brd1;
EY1 = ce1(Y,ev1);
D2 = (Y-EY1)~X;
{bRd2,f,g,cRd,ret}  = maxlik(D2,0,&slsRd2,start2);
ev2 = X[.,2] + X[.,4]*brd2;
EY2 = ce1(Y-eY1,ev2);
D3 = (Y-EY2)~X;
{bRd3,f,g,cRd,ret}  = maxlik(D3,0,&slsRd3,start3);
Maxset;
maxset;
_max_maxiters = 100;
start = brd1|brd2|brd3;*/


{bRd,f,g,cRd,ret} = maxlik(D,0,&slsRd,gstart);

/*SLS with Recusive Differencing and  Residuial Control*/
vv1  = X[.,1] + X[.,4]*bRd[1];
vv2  = X[.,2] + X[.,4]*bRd[2];
vv3  = x[.,3] + x[.,4]*bRd[3];
evv  = vv1~vv2~vv3;
tv   = trimd(evv,.97);
@tvv1 = trimd(evv,.98);
tv1  = smtrm(evv[.,1], deltasm,.98).*smtrm(evv[.,2],deltasm,.98).*smtrm(evv[.,3],deltasm,.98);@
tv1 = trimdexpan(evv,deltasm,.98);
Gscal1   = Gscale1(evv);
Gscal2   = Gscale2(evv);


Maxset;
_max_maxiters = 50;
{bRdRc,f,g,cRdRc,ret} = maxlik(D,0,&SLSRdRc,bRd);

storebRd    = storebRd|bRd';
storebRdRc  = storebRdRc|bRdRc';
save storebRd,storebRdRc ;

ii = ii + 1;
endo;
/*Results*/
erd   = storebrd   - True';
erdrc = storebrdrc - True';



/*Programs for RD as Sole Bias Control*/

proc gridmaxStart(n,fine,step,d,bstart);
local grid,b,stepi,i,l,bval,Lval;
grid=gengamma(-n,n,fine,3);
stepi=1;
do while stepi<=step;
b=grid+bstart;
bval = b';
save bval;
i=1;
l=ones(cols(b),1); 
do while i<=cols(b);
l[i]=sumc(SLSstart(b[.,i],d)); 
i=i+1;
endo;
Lval = L;
save Lval;
bstart=b[.,maxindc(l)];
stepi=stepi+1;
endo;
retp(bstart);
endp;

proc(1)=gengamma(gamstart,gamend,inc,gamdim);  /*create gamma matrix [gamstart,gamend]^gamdim with increment*/
local seq,gamseqn,gam,i;
gamseqn=1+(gamend-gamstart)/inc;
seq=seqa(gamstart,inc,gamseqn);
gam=ones(gamdim,gamseqn^gamdim);
i=1;
do while i<=gamdim;
gam[i,.]=reshape(seq*ones(1,gamseqn^(gamdim-i)),1,gamseqn^gamdim);
i=i+1;
endo;
retp(gam);
endp;

Proc slsRd1(b,D);
Local Y,X,v1,v2,v3,v,E,Ehat,e1,e2,e3,f1,f2,f3,g3,g1,g2,r2,s0,s1;
Y = D[.,1];
X = D[.,2:cols(D)];
v1 = X[.,1] + X[.,4]*b;
v2 = X[.,2]; 
v3 = x[.,3];
v = v1~v2~v3;
Ehat = ce1(Y,v);
r2 = (Y-Ehat)^2;
retp(-t2.*r2);
endp;

Proc slsRd2(b,D);
Local Y,X,v1,v2,v3,v,E,Ehat,e1,e2,e3,f1,f2,f3,g3,g1,g2,r2,s0,s1;
Y = D[.,1];
X = D[.,2:cols(D)];
v1 = X[.,1] ;
v2 = X[.,2]*b; 
v3 = x[.,3];
v = v1~v2~v3;
Ehat = ce1(Y,v);
r2 = (Y-Ehat)^2;
retp(-t2.*r2);
endp;

Proc slsRd3(b,D);
Local Y,X,v1,v2,v3,v,E,Ehat,e1,e2,e3,f1,f2,f3,g3,g1,g2,r2,s0,s1;
Y = D[.,1];
X = D[.,2:cols(D)];
v1 = X[.,1] ;
v2 = X[.,2]; 
v3 = x[.,3] + X[.,4]*b;
v = v1~v2~v3;
Ehat = ce1(Y,v);
r2 = (Y-Ehat)^2;
retp(-t2.*r2);
endp;

Proc  ce1(Y,v);
local h,E1,Temp,n,r,j,tob,tj,Z,K;
n=rows(v);
h=(n^(-1/7))*stdc(v)';
E1 = {};
j=1;
tob = seqa(1,1,N);
do while j .le n;
tj = (tob .ne j);
Z=(v-v[j,.]);
k=pdfn(Z./h)./h;
k=prodc(k').*tj; 
Temp = meanc(Y.*K)/meanc(K);
E1 = E1|Temp;
j=j+1;
endo;
retp(E1);
endp;

Proc slsStart(b,D);
Local Y,X,v1,v2,v3,v,E,Ehat,e1,e2,e3,f1,f2,f3,g3,g1,g2,r2,s0,s1;
Y = D[.,1];
X = D[.,2:cols(D)];
v1 = X[.,1] + X[.,4]*b[1];
v2 = X[.,2] + X[.,4]*b[2];
v3 = x[.,3] + x[.,4]*b[3];
v = v1~v2~v3;
Ehat = MhatStart(Y~v);
r2 = (Y-Ehat)^2;
retp(-t2.*r2);
endp;

Proc MhatStart(D);
Local Y,X,v1,v2,v,E,e1,e2,e3,f1,f2,f3,g3,g1,g2,r2;
Y = D[.,1];
X = D[.,2:cols(D)];
e1  =  CE1Rd(Y,x,x,1/7);
retp(e1); 
endp;


Proc slsRd(b,D);
Local Y,X,v1,v2,v3,v,E,Ehat,e1,e2,e3,f1,f2,f3,g3,g1,g2,r2,s0,s1;
Y = D[.,1];
X = D[.,2:cols(D)];
v1 = X[.,1] + X[.,4]*b[1];
v2 = X[.,2] + X[.,4]*b[2];
v3 = x[.,3] + x[.,4]*b[3];
v = v1~v2~v3;

Ehat = MhatRd(Y~v);
r2 = (Y-Ehat)^2;
retp(-t2.*r2);
endp;

Proc MhatRd(D);
Local Y,X,v1,v2,v,E,e1,e2,e3,f1,f2,f3,g3,g1,g2,r2;
Y = D[.,1];
X = D[.,2:cols(D)];
e1  =  CE1Rd(Y,x,x,1/11.99  );
e2  =  CEsRd(y,e1,x,x,1/11.99);
e3  =  CEsRD(y,e2,x,x,1/11.99);
retp(e3); 
endp;

proc ce1Rd(y, x, arg, r);
local e1,h,narg,n,j,tob,tj,Xj,Yj,k,gscale,tI,delta,hh,Z,Z1,mm,m,m1,m2,d,Zbar;
narg=rows(arg);
n=rows(x);
h=(n^(-r))*stdc(x)';
m = {};
j=1;
tob = seqa(1,1,N);
do while j .le narg;
tj = (tob .ne j);
Z=(x-arg[j,.]);
k=pdfn(Z./h)./h;
k=prodc(k').*tj; 
m1 = meanc(Y.*K)/meanc(K);
Zbar =meanc(z.*K)./meanc(K) ;
K = K;
d  = inv( (Z.*K)'Z )*(Z.*K)'(Y-m1);
m2 = m1 - Zbar'd;
m  = m|M2;
j=j+1;
endo;
retp(m);
endp;

proc CEsRd(y, e3,x, arg, r);
local YY,e1,f3,g3,e,h,narg,n,j,tob,tj,Xj,Yj,k,A;
narg=rows(arg);
n=rows(x);
h=(n^(-r))*stdc(x)';
e={};
f3 = {};
g3 = {};
j=1;
tob = seqa(1,1,N);
do while j .le narg;
tj = (tob .ne j).*t1;
YY = Y - (e3 - e3[j]);
k=pdfn((arg[j,.]-x)./h)./h;
k=prodc(k');
g3 = g3|meanc(tj.*k);
f3 = f3|meanc(tj.*k.*yy);
j=j+1;
endo;
e = f3./g3;
retp(e);
endp;

/*Programs for RD and Newey's Residual Control*/

Proc SlsRdRc(b,D);
Local Y,X,v1,v2,v3,v,E,Ehat,e1,e2,e3,f1,f2,f3,g3,g1,g2,r2,s0,s1;
Y = D[.,1];
X = D[.,2:cols(D)];
v1 = X[.,1] + X[.,4]*b[1];
v2 = X[.,2] + X[.,4]*b[2];
v3 = x[.,3] + x[.,4]*b[3];
v = v1~v2~v3;
Ehat = MhatRdRc(Y~v);
r2 = (Y-Ehat)^2;
retp(-tv.*r2);
endp;

Proc MhatRdRc(D);
Local Y,X,v1,v2,v,E,e1,e2,e3,f1,f2,f3,g3,g1,g2,r2;
Y = D[.,1];
X = D[.,2:cols(D)];
e1  =  CE1RDRc(Y,x,x,1/15); 
e2  =  CEsRdRc(y, e1,x, x, 1/15);
e3  =  CEsRdRc(y, e2,x, x, 1/15);
retp(e3);
endp;

proc ce1RdRc(y, x, arg, r);
local mz, g1,h,narg,n,cx,dent,j,tob,tj,Xj,Yj,k,delta,hh,Z,Z1,mm,m,m1,m2,d,Zbar,tI;
narg=rows(arg);
n=rows(x);
    cx = cols(x);
    dent = eye(cx);
h=(n^(-r))*stdc(x)';
hh =  N^-deltadjust;
m = {};
g1 = {};
j=1;
tob = seqa(1,1,N);
do while j .le narg;
tj = (tob .ne j);
Z=(x-arg[j,.]);
k=pdfn(Z./h)./h;
k=prodc(k').*tj; 
 g1 = g1|meanc(K);
j=j+1;
endo;
tI = smdentrm(tden1(x),hh);
delta=(1-tI).*hh.*gscal1;
g1 = g1 + delta;
tob = seqa(1,1,N);
j = 1;
do while j .le Narg;
tj = (tob .ne j);
Z=(x-arg[j,.]);
k=pdfn(Z./h)./h;
k=prodc(k').*tj; 
m1   = meanc(Y.*K)/g1[j];
Zbar = meanc(Z.*K)/g1[j];
K = K;  
mz = (Z.*K)'Z + delta[j]*dent; 
d  = inv( mz )*(Z.*K)'(Y-m1);
m2 = m1 - Zbar'd;
m  = m|M2;
j=j+1;
endo;

retp(m);
endp;

proc CEsRdRc(y, e3,x, arg, r);
local YY,e1,f3,g3,e,h,narg,n,j,tob,tj,Xj,Yj,k,A,delta,hh,tI;
narg=rows(arg);
n=rows(x);
h=(n^(-r))*stdc(x)';
@hh=(n^(-r/2)); @
 hh = N^-deltadjust;
e={};
f3 = {};
g3 = {};
j=1;
tob = seqa(1,1,N);
do while j .le narg;
tj = (tob .ne j).*tv1;
YY = Y - (e3 - e3[j]);
k=pdfn((arg[j,.]-x)./h)./h;
k=prodc(k');
g3 = g3|meanc(tj.*k);
f3 = f3|meanc(tj.*k.*yy);
j=j+1;
endo;
@tI  = insidetrimv(x,.01,.99); @
tI = smdentrm(tden1(x),hh);
delta=(1-tI).*(hh^.5).*gscal2; 
e=( f3./(g3+delta) ) ;
retp(e);
endp;

/*Trimming and Adjustment Programs*/

 
Proc  tden1(v);
local h,g1,n,r,j,tob,tj,Z,K;
n=rows(v);
h=(n^(-1/7))*stdc(v)';
g1 = {};
j=1;
tob = seqa(1,1,N);
do while j .le n;
tj = (tob .ne j);
Z=(v-v[j,.]);
k=pdfn(Z./h)./h;
k=prodc(k').*tj; 
g1 = g1|meanc(K);
j=j+1;
endo;
retp(g1);
endp;


Proc  Gscale1(v);
local h,g1,n,r,j,tob,tj,Z,K;
n=rows(v);
h=(n^(-1/7))*stdc(v)';
g1 = {};
j=1;
tob = seqa(1,1,N);
do while j .le n;
tj = (tob .ne j);
Z=(v-v[j,.]);
k=pdfn(Z./h)./h;
k=prodc(k').*tj; 
g1 = g1|meanc(K);
j=j+1;
endo;
retp(quantile(g1,.01));
endp;


Proc tden2(v);
local h,g1,n,r,j,tob,tj,Z,K;
n=rows(v);
h=(n^(-1/7))*stdc(v)';
g1 = {};
j=1;
tob = seqa(1,1,N);
do while j .le n;
tj = (tob .ne j).*tv1;
Z=(v-v[j,.]);
k=pdfn(Z./h)./h;
k=prodc(k').*tj; 
g1 = g1|meanc(K);
j=j+1;
endo;
retp(g1);
endp;


Proc Gscale2(v);
local h,g1,n,r,j,tob,tj,Z,K;
n=rows(v);
h=(n^(-1/7))*stdc(v)';
g1 = {};
j=1;
tob = seqa(1,1,N);
do while j .le n;
tj = (tob .ne j).*tv1;
Z=(v-v[j,.]);
k=pdfn(Z./h)./h;
k=prodc(k').*tj; 
g1 = g1|meanc(K);
j=j+1;
endo;
retp(quantile(g1,.01));
endp;


proc trimd(X,f);
local N,vx,svx,tx,Low,High,t,cx,i;
N = rows(x);
cx = cols(x);
tx = zeros(cx,N);
i = 1;
do while i .le cx;
	vx = X[.,i];
	svx = sortc(vx,1);
	Low  =  quantile(svx,1-f);
	High =  quantile(svx,f);
	tx[i,.] = ((vx .ge low) .and (vx .le high))';
i = i + 1;
endo;
t = prodc(tx);
retp(t);
endp;

/* indicator trimming with expansion */
proc trimdexpanc(X,delta,fi,fo);
local N,vx,svx,tx,Low,High,t,cx,i,mi,ma,tt;
N = rows(x);
cx = cols(x);
tt = {};
i = 1;
do while i .le cx;
	vx = X[.,i];
	svx = sortc(vx,1);
    mi= quantile(svx,1-fi);
    ma= quantile(svx,fi);
	Low  =  (quantile(svx,1-fo)+mi*delta*ln(N))./(1+delta*ln(N));
	High =  (quantile(svx,fo)+ma*delta*ln(N))./(1+delta*ln(N));  
    
	tt = tt~(pilotrm(vx,Low).*pilotrm(High,vx));
i = i + 1;
endo;
t = prodc(tt');
retp(t);
endp;


proc trimdexpan(X,delta,f);
local N,vx,svx,tx,Low,High,t,cx,i,mi,ma;
N = rows(x);
cx = cols(x);
tx = {};
i = 1;
do while i .le cx;
	vx = X[.,i];
	svx = sortc(vx,1);
    mi=minc(svx);
    ma=maxc(svx);
	Low  =  (quantile(svx,1-f)+mi*delta*ln(N))./(1+delta*ln(N));
	High =  (quantile(svx,f)+ma*delta*ln(N))./(1+delta*ln(N));    
	tx = tx~((vx .ge low) .and (vx .le high));
i = i + 1;
endo;
t = prodc(tx');
retp(t);
endp;

proc  smdentrm(g,hh);
local mminv,mmaxv,trmL,trmU,tt,i;
trmL = pilotrm(g,hh); 
retp(trmL);
endp;

@proc  insidetrim(g,hh);
local mminv,mmaxv,trmL,trmU,tt,i;
trmL = pilotrm(g,hh); 
retp(trmL);
endp;@

proc  insidetrimg(g,hh);
local mminv,mmaxv,trmL,trmU,tt,i;
trmL = pilotrm(g,hh); 
retp(trmL);
endp;

proc  insidetrimv(ev,delta,f);
local mminv,mmaxv,trmL,trmU,tt,i,mi,ma;
tt = {};
i = 1;
do while i .le cols(v);
mi=minc(ev[.,i]);
ma=maxc(ev[.,i]);
@mminv = quantile(v[.,i],0.01);@
mminv = (quantile(ev,1-f)+mi*delta*ln(N))./(1+delta*ln(N));
mmaxv = (quantile(ev,f)+ma*delta*ln(N))./(1+delta*ln(N));
trmL = pilotrm(v[.,i],  mminv); 
trmU = pilotrm(mmaxv,  v[.,i]); 
tt = tt~(trmL.*trmU);
i = i + 1;
endo;
retp(prodc(tt'));
endp;


Proc Smtrm(ev,delta,f);
local q1,q2,mi,ma,t;
mi=minc(ev);
ma=maxc(ev);
q1 = (quantile(ev,1-f)+mi*delta*ln(N))./(1+delta*ln(N));
q2 = (quantile(ev,f)+ma*delta*ln(N))./(1+delta*ln(N));
t = pilotrm(ev,q1).*pilotrm(q2,ev);
retp(t);
endp;
proc  Pilotrm(g,delt);
local N, tp,  temp;
N = maxc(rows(g)|rows(delt));
temp=  1 + exp( -(ln(N)*ln(N))*(g-delt) );
tp =  1./ temp;
retp(tp);
endp;