*This file ammends /homes/trefler/hov/test/feb94/gqopt.sas;
*Trefler, Nov 12, 1994;

options nocenter ls=172;
options compress=no;

libname data '/public/trefler/pub-dat/test';
* label at='factor content of trade, F=A*T'
         v='endowment'
	 y='gdp, world bank, y=p*Q'
	 b='trade balance, world bank b=p*T'
	 ypc='gdp per capita, PWT'
	 country='alphabetic country name'
	 factor='alphabetic factor name'
	 indexc='country identifier'
	 indexf='factor identifier'
	 pop='population, wb';
* keep indexc country indexf factor at v y b ypc;

*data subscripted by both factors and countries;
data one;
   set data.hov_pub;
   keep country indexc factor indexf at v;

*data subscripted by countries only;
data two;
   set data.hov_pub;
   *output all the data for a single factor;
   if indexf=1 then output;
   keep country y b ypc;
   	 
proc iml;
   use one;
   read all var {indexc} into indexc;
   read all var {indexf} into indexf;
   read all var {at}     into at1;
   read all var {v}      into v1;
   use two;
   read all var {y}     into y;
   read all var {b}     into b;
   read all var {ypc}  into ypc;
   read all var {country}  into country;

   v  =j(9,33,.);
   at=j(9,33,.);
   do i=1 to 297;
      indc=indexc[i];
      indf=indexf[i];
      v[indf,indc]   =   v1[9*(indc-1)+indf];
      at[indf,indc] = at1[9*(indc-1)+indf];
      end;
   ct=country`;
*  print v[colname=ct];
*  print at[colname=ct];
   
   yw=0;
   bw=0;
   dvw=j(9,1,0);
   vw =j(9,1,0);   
   delta=j(33,1,1);
   do ic=1 to 33;
      ypc[ic]=ypc[ic]/ypc[33];
*     delta[ic]=ypc[ic];
      yw=yw+y[ic];
      bw=bw+b[ic];
      dvw=dvw+delta[ic]*v[,ic];
      vw=vw+v[,ic];
      end;
      
   s  = (y-b)/(yw-bw);
   rho= y/yw;
*  print country y b s rho yw bw      ;
*  print vw dvw;

   *compute scaling;
   e=j(9,33,.);
   p=j(9,33,.);
   sigma=j(9,1,.);
   do if=1 to 9;
      sigma0=0;
      emean=0;
      do ic=1 to 33;
         p[if,ic]=v[if,ic]-s[ic]*vw[if];
         e[if,ic]=at[if,ic] - p[if,ic];
	 emean=emean+e[if,ic]/33;
	 end;
      do ic=1 to 33;
	 sigma0=sigma0+(e[if,ic]-emean)**2;
	 end;
      sigma[if]=sqrt( sigma0 / 32);
      end;
*  print sigma;
   w=j(9,33,.);
   do if=1 to 9;
      do ic=1 to 33;
         w[if,ic]=sigma[if] * sqrt(s[ic]);
	 end;
      end;
*  print indexc indexf w;


****************************
 Repeat JPE as data check
****************************;
   cs=(y-b)/yw;
   cx=j(32,32,.);
   cy=j(32,1,.);
   cbf=j(32,1,0);
   do if=1 to 9;
      do ic=1 to 32;
         cy[ic]=at[if,ic]+cs[ic]*v[if,33];
         do icc=1 to 32;
            if ic=icc then cx[ic,icc]=(1-cs[ic])*v[if,ic];
            else           cx[ic,icc]=-cs[ic]*v[if,icc];
            end;
         end;
      cb=inv( (cx`)*cx ) * (cx`)*cy;
      cbf=cbf||cb;
      end;
      cbff=cbf[,2:10];
      print country cbff[format=5.2];


****************************
 Do SIGN tests
****************************;
   sign =j(9,33,0);
   pred =j(9,33,0);
   signa=j(9,33,0);
   preda=j(9,33,0);
   ata=j(9,33,0);
   msign=j(33,1,.);
   msigna=j(33,1,.);
   do c=1 to 33;
      do f=1 to 9;
         pred[f,c] = v[f,c] - s[c]*vw[f];
	 if at[f,c]*pred[f,c]>=0 then sign[f,c]=1;
	 ata[f,c]  = at[f,c] - (b[c]/(yw-bw)) * vw[f];
	 preda[f,c] = v[f,c]   - (y[c]/(yw-bw)) * vw[f];
         if ata[f,c]*preda[f,c]>=0 then signa[f,c]=1;
	 end;
      end;
   msign=t(sign[+,]/9);
   msigna=t(signa[+,]/9);
   print country msign[format=5.2] msigna[format=5.2];

****************************
 Do RANK tests
****************************;
NC=33;
NF=9;
count = j(nc,1,0);
prop  = j(nc,1,.);
counta = j(nc,1,0);
propa  = j(nc,1,.);
*scale because cross-factor comparisons will be made;
do c=1 to nc;
   do f=1 to nf;
      at[f,c]=at[f,c]/w[f,c];
      ata[f,c]=ata[f,c]/w[f,c];
      pred[f,c]=pred[f,c]/w[f,c];
      preda[f,c]=preda[f,c]/w[f,c];
      end;
   end;
do c=1 to nc;
   do f1=2 to nf;
      f1m = f1-1;
      do f2=1 to f1m;
         if (pred[f1,c]-pred[f2,c]) * (at[f1,c]-at[f2,c]) >= 0
	    then count[c] = count[c]+1;
         if (preda[f1,c]-preda[f2,c]) * (ata[f1,c]-ata[f2,c]) >= 0
	    then counta[c] = counta[c]+1;
         end;
      end;
   prop[c] = count[c]/(nf*(nf-1)/2);
   propa[c] = counta[c]/(nf*(nf-1)/2);
   end;
   print country ypc[format=5.2]
       msign[format=5.2] msigna[format=5.2] prop[format=5.2] propa[format=5.2];

*Compute statstics for all countries together;
c=count[+];ca=counta[+];
d=nc*nf*(nf-1)/2;
p=c/d; pa=ca/d;
print c ca p pa;


* prepare graphs;
e = at - pred;
esign=j(nf,nc,0);
predsign=j(nf,nc,0);
do c=1 to nc;
   do f=1 to nf;
      if e[f,c]<0 then esign[f,c]=1;
      if pred[f,c]>=0 then predsign[f,c]=1;
      end;
   end;
mes=t(esign[+,]);
mps=t(predsign[+,]);
print country mes mps;