/*** Setup the global parameters ***/
/*The parameters below should be changed accordingly*/
*%LET filename="c:\zzy\mnarmediation\active1.txt"; * data file directory and name;
%LET varname=x m y;  *specify variable names in the data file. Please 
use x for the input variable, m for the mediation variable, and y for 
the output variable. a1 and a2 are two auxiliary variables in the data 
file. You can use any names except for x, m, and y for naming the 
auxiliary variables;
%LET missing=99999; *specify the missing data value;
%LET nboot = 1000;  *define the number of bootstraps B;
%LET alpha = 0.95;  *define the confidence level;
%LET seed = 2010;   *random number seed;
/*** End of setup of global parameters ***/
data all;
input col1-col4;
cards;
0 0 0 0
;
run;

%MACRO runsim(i);
data temp;
number = &i;
datastr=put(number, 4.);
datastrcomp=compress(datastr,'');
commd ='cp /pscratch/zzhang4/Private/mediation/data/simucond/data-noav-'||datastrcomp||'.txt data.txt';
run;

data _null_;
set temp;
file 'cp.sh';
put commd;
run;

data _null_;
X 'chmod 755 cp.sh';
X './cp.sh';
X 'rm -f cp.sh';
run;
quit;

/*In general, there is no need to change the codes below*/
/*Read data into sas*/
DATA dset;
  INFILE "data.txt";
  INPUT &varname;
  ARRAY nvarlist &varname;
  DO OVER nvarlist;
    IF nvarlist = &missing THEN nvarlist = .;  
  END;
RUN;

/*Calculate covariance matrix based on pairwise deletion*/
PROC CORR DATA=dset	cov OUTP=ocov nosimple nocorr noprint;
RUN; QUIT;

DATA ocov;
SET ocov;
KEEP x m y;
RUN;

/*** Bootstraping data to obtain standard errors and confidence 
intervals ***/
DATA bootsamp;
  DO sampnum = 1 to &nboot;
     DO i = 1 TO nobs;
        ran = ROUND(RANUNI(&seed) * nobs);
        SET dset
        nobs = nobs
        point = ran;
        OUTPUT;
     END;
  END;
  STOP;
  DROP i;
RUN; QUIT;

PROC CORR DATA=bootsamp cov OUTP=bcov nosimple nocorr noprint;
by sampnum;
RUN; QUIT;

DATA bcov;
SET bcov;
KEEP x m y;
RUN;

/*** Calculate the BC intervals based on the point estimates from 
different bootstrap samples and produce a table containing the points 
estimates, standard errors, confidence intervals in the output 
window.***/
PROC IML;
  START main;
  USE ocov;
  READ ALL INTO Y;
  USE bcov;
  READ ALL INTO X;

  n=NROW(X);
  m=n/6;

  oest=	Y[1,2]/Y[1,1]*(Y[2,3]*Y[1,1]-Y[1,2]*Y[1,3])/(Y[1,1]*Y[2,2]-Y[1,2]*Y[1,2]);

  best=J(m,1,0);
  a=J(1,m,0);
  b=J(1,m,0);
  temp=0;
  DO j=1 TO m;
     temp=1+(j-1)*6;
	 a[j]= X[temp,2]/X[temp,1];
	 b[j]= (X[temp+1,3]*X[temp,1]-X[temp,2]*X[temp,3])/(X[temp,1]*X[temp+1,2]-X[temp,2]*X[temp,2]);
     best[j]=a[j]*b[j];
  END;

  result=J(1,4,0);
  result[1] = oest;

  alphas=1-(1-&alpha)/2;
  zcrit = PROBIT(alphas);

	result[2]=SQRT((SSQ(best) -(SUM(best))**2/n)/(n-1)); 
	print result;

    number=0;
	DO i=1 TO m;
		IF oest[1]<best[i] THEN number=number+1;
	END;

	p=number/m;

	z0hat=PROBIT(p);

	q1=z0hat+(z0hat-zcrit);
	q2=z0hat+(z0hat+zcrit);
	alpha1=PROBNORM(q1);
	alpha2=PROBNORM(q2);

	vec=best;
	CALL SORT(best,{1});

	low=int(alpha1*(m+1));
	up=int(alpha2*(m+1));

	IF low<1 THEN low=1;
	IF up>m THEN up=m;
	result[3]=best[low];
	result[4]=best[up];

create bcci from result;
append from result;
  FINISH;
  RUN main;
QUIT;
/*Combine results*/
%MEND runsim;

%MACRO batch;
%LET r=1001;
%DO %WHILE (&r <=1100);
  %runsim(&r);
  %LET r = %eval(&r+1);
  data all;
   set all bcci;
run;
%END;
%MEND batch;

%batch;


data _NULL_;
set all;
file "all.txt";
put col1-col4;
run;
quit;