options pageno = 1 ps=55;

/* This program calculates the probability that a county was chosen */
/* for a sample of size 2 given the other counties in the strata.   */
/* Note that this version does not require the sample size as one   */
/* of the macro arguments.                                          */

%macro strata(strata);

data a;
*Change the directory if necessary;
infile "/home/grego1/STAT 540/&strata..csv" dlm="," end=last;
input county $:11. aad;
/* StratumSize should be numeric, hence the put statement; I    */
/* initially forced it to have the right length--1.--to get the  */
/* ARRAY statements to work.  That seemed artificial, so I      */
/* switched to this syntax, which is a little clumsy.           */
/* SYMPUTN would work for StratumSize alone, but we would still */
/* need to handle ColIndex and PropIndex as shown below         */
if last=1 then do; call symputx("StratumSize",put(_n_,2.)); 
                   call symputx("ColIndex","col"||left(put(_n_,2.)));
		       call symputx("PropIndex","Prop"||left(put(_n_,2.)));
end; 
run;

/* Set up calculation of sampling probabilities */
data c;
   set a;
keep aad;
proc transpose data=c out=d;

/* The probability of selection for a sample of size 1 */
/* is proportional to traffic intensity  */ 
data e;
   set d;      **Set the probabilities into an array**;
total = sum(of col1-&ColIndex);

array b col1 - &ColIndex;
array c prop1 - &PropIndex;
do i = 1 to &StratumSize;
   c(i) = b(i)/total;
end;

/* Discard variables we don't need         */
data one;
   set e;
keep prop1 - &PropIndex;
data two;
   set one;

array p prop1-&PropIndex;

**CALCULATE ZJ TERMS**;      **Begin calculation PIE-i s**;
do i=1 to &StratumSize;
   sum1j=0;
   ind=i;
   ptempi=p(ind);

   do j=1 to &StratumSize;
      if j ne i then do;
      ind=j;
      ptempj=p(ind);
      zj=ptempj/(1-ptempj);
      sum1j=sum1j+zj;
      end;
    end;
/* Output the probability of selection for a */
/* sample of size 2                          */
final1=ptempi*(1+sum1j);
output;
end;

/* Clean-up again          */
data three; set two;
   keep final1;

/* This step doesn't make alot of sense anymore  */
/* We used to compute these values only for the  */
/* two counties we were actually going to sample */
/* from (it saved computational time)            */
data four;
   set a(firstobs=1 obs=&StratumSize); keep county;

data final (rename=(final1=Prob county=County strata=Stratum));
   merge three four;
   strata="&strata";
   keep county strata final1;

proc print noobs; var county prob; 
title "Sampling weights for &strata counties";

libname weights "/home/grego1/STAT 540/";

data weights.&strata;          **Write appropriate info to file**;
   set final;

%mend;

%strata(upurban)

run;