/*   This is file GaiaQuai.GAU: Given p and q matrices and QUAIDS parameters (for a particular
           value of NSystem), it calculates Geary and GAIA real incomes, and writes results to disk.
	This uses NLSYS to solve the budget share equation directly for the GAIA world prices 

(c) J.P. Neary 1999

*/
new;  /*************************************************************************/
library gearylib, pgraph, nlsys;      #include nlsys.ext; nlset;  /* Set up NLSYS files */
     /* Note: pgraph is not used but myst still be set up, since otherwise gearylib conflicts */
/*************************************************************************/
/*          USER-SPECIFIED VARIABLES: ECONOMIC DATA AND PARAMETERS  */
n=11;                       /* n = # goods        N.B. This and next two overridden if TESTIF=1 */
nyear=1980;
/*************************************************************************/
/*    USER-SPECIFIED VARIABLES WHICH CONTROL OPERATION OF PROGRAM  */
NSystem=2;   /* NSystem = 1 for HAIDS; =2 for AIDS; =3 for QUAIDS */
TESTIF=0;     /* TESTIF=1 on a test run (using small data set only): overrides n, nyear, NSystem  */
kstart =0; /* Starting value for k: must be >=0 and <=klast */
klast=10;   /*   Highest value of k for which GAIA estimates are calculated (0 <= klast <= n-1)  */
printsor=0;      /* printsor= 0 for no printing of data in PROC SORTDATA; 1 for printing */
printif=1;          /* printif= 1 for printing of results in each round */
printif9=1;       /* printif9=1 for printing of summary table */
/*************************************************************************/
/* Read Demand Parameter Matrices (for k=0,n-1), created by Qk_It.GAU from disk:  
       alphaALL, betaALL, and lamALL of dimension n-by-(n-1); gammaALL is (n*(n-1))-by(n-1).
    First, set location of files: depends on which system is being considered */

IF NSYSTEM==1;
   load path = c:\me\research\geary2\empirics\gaussout\HAIDS ;
   save path = c:\me\research\geary2\empirics\gaussout\HAIDS ;
   screen on; print "ML HAIDS parameters, n-1 goods";
   ELSEIF NSYSTEM==2;
      load path = c:\me\research\geary2\empirics\gaussout\AIDS ;
      save path = c:\me\research\geary2\empirics\gaussout\AIDS ;
      screen on; print "ML AIDS parameters, n-1 goods";
   ELSEIF NSYSTEM==3;
      load path = c:\me\research\geary2\empirics\gaussout\QUAIDS ;
      save path = c:\me\research\geary2\empirics\gaussout\QUAIDS ;
      screen on; print "ML QUAIDS parameters, n-1 goods";
ENDIF;

load  alphaALL betaALL lamALL gammaALL;
screen OFF; format /rd 7,4; print  "k:          alpha's:"; print seqa(0,1,n)~alphaALL;
                                                     print  "k:          beta's:"; print seqa(0,1,n)~betaALL;
                                                     print  "k:          lambda's:"; print seqa(0,1,n)~lamALL;
/*************************************************************************/
{ m,p,q } = readdata(n,nyear,testif);  /* Call PROC READDATA to read raw data and print
   heading:  m=#countries; p and q are n-by-m matrices of prices and quantities */
countrynum=seqa(1,1,m);
/**************************************************************************/
            /*  INITIALISE MATRICES  */
mbase=1;        /* mbase = index of base country:  =1 for richest;  =m for poorest country  
     WARNING:  Geary section only works for mbase=1 ;  24.2.98 */
mrel=m;    /* Index of country relative to which the results are printed */
omegastr = zeros (n,m); /* n-by-m matrix of virtual budget shares */
qstar = zeros (n,m); /* n-by-m matrix of virtual quantities */
ZGAIAALL=ones(m,n);  /*  column kk gives the calculated GAIA index for k=kk-1 */ 
piall=zeros(n,n);  /*  column kk gives the calculated world prices for k=kk-1 */ 
/***********************************************************************/
/* Call PROC SORTDATA to sort p and q by z, scale p and z by sample mean 
     and calculate: pprimeq, z, w, w_trim, poverz and poz_trim */
{p, q, pprimeq, z, zraw, w, w_trim, poverz, poz_trim}=SORTDATA(p,q,printsor);
/***********************************************************************/
/*  CALCULATE GEARY INCOMES and WORLD PRICES  */
{ ZGeary, ppiGeary } = Gearyinc(w,q);   /* Call PROC Gearyinc to calculate Geary incomes */
      screen on; format /rd 6,3;  print "zG:          " ZGeary';    print "pi:          "  ppiGeary';
/*************************************************************************/
/*  START UTILITY-BASED DERIVATIONS, LOOPING OVER VALUES OF k */
k=kstart-1;
DO WHILE k<=klast-1;   /*  START BIG LOOP: over k=kstart, klast  */
   k=k+1;
/*************************************************************************/
/* Retrieve n-1 parameter values for current value of k */
alpha = alphaALL[k+1,.]';      /* Note: ROW # k+1 of alphaALL becomes COL. vector alpha */
beta = betaALL[k+1,.]';     
lambda = lamALL[k+1,.]';
nfirstga = k*(n-1) + 1;
nlastga = (k+1)*(n-1);
ggamma = gammaALL[nfirstga:nlastga,.] ;
/*************************************************************************/
/* (Given current value of k) calculate parameters for good n by subtraction */
alphan=alpha | 1-sumc(alpha);   /* calculate alpha for good #n by subtraction */
betan= beta | -sumc(beta);   /* calculate beta for good #n by subtraction */
lambdan= lambda | -sumc(lambda);   /* calculate lambda for good #n by subtraction */
temp = ggamma | (-sumc(ggamma))';
ggamman= temp ~ (-sumc(temp'));   /* calculate gamma's for good #n by subtraction */
/*************************************************************************/
/* Print all parameter values  */        screen on; 
format /rd 2,0; print "Params. for Semiflexible QUAIDS [nsystem =" nsystem "] of Rank " k ;
format /rd 6,4; print "alpha, beta, lambda and gamma: " alphan~betan~lambdan~ggamman; 
/*************************************************************************/
/*  CALCULATE m-by-1 VECTOR OF ACTUAL COUNTRY UTILITY LEVELS   */
logofy = log(z) - logalfap(p); 
u=1./( lambofp(p)+ bbetaofp(p)./logofy); 
/*************************************************************************/

/*   SKIP THIS.

/*  If NSYSTEM=1:  Calculate GAIA index for HAIDS case directly */
IF NSYSTEM==1;  /* Special HAIDS case: No longer needed, but a useful check on GaiaQuai */
  ZGAIA = 10^(u-u[mrel,.]);  /* Under HAIDS, relative GAIA index is just antilog of utility differences */
  ZGAIAALL[.,k+1] = ZGAIA;   /* Read ZGAIA into appropriate col. (N.B.) of summary ZGAIAALL matrix */
ENDIF;  

*/

/*************************************************************************/
/* CALCULATE GEARY-KONUS PRICES and GAIA REAL INCOMES */
/***********************************************************************/
TIMETAKEN=HSEC;   /*  Note current time */
/**************************************************************************/
/* Call NLSYS to solve non-linearly: Find the square root of PI (to avoid negative prices) */
ppi = ppiGeary;  /* Starting value for ppi is Geary vector in each round? */
/* Using ppi from previous loop [k-1] as starting value doesn't work for Q[0] */
x0 = sqrt(ppi[1:n-1,1]);   /* Starting vector is sqrt of (n-1)-by-1 vector of Geary world prices */

IF PRINTIF==0;
   screen off; { xxx, fvc, jc, tcode } = NLSYS(&func1, x0) ;   /* Call NLSYS with no printing */ 
ELSEIF PRINTIF==1;
   screen On; format /rd 2,0;
   PRINT "\LGAIA Prices and Incomes Calculated using NLSYS: nsystem= " nsystem ", k= " k;
   { xxx, fvc, jc, tcode } = NLPRT (NLSYS(&func1, x0)) ;   /* Call NLSYS, print all iterations */ 
ENDIF;
/*************************************************************************/
proc func1(xxx);   /* PROC to calculate deviations of budget share equations  */
local ppi, logofy, omegastr, zstar, errorvec;   /* globals: u, alphan, betan, ggamman, lambdan, m, W */
   /*  CALCULATE n-by-m Matrices of QUAIDS Virtual Budget Shares and Quantities */
ppi = xxx^2 | eye(1);  /* Input vector is the square root of the price vector for n-1 goods */
logofy = logystar(ppi,u);
omegastr = (alphan+ggamman*log(ppi)) *(ones(m,1))' + betan*(logofy)' + 
     lambdan*(((logofy)^2)/bbetaofp(ppi))' ; /* n-by-m matrix of virtual budget shares */ 
zstar = 10^(ones(m,1)*logalfap(ppi) + logofy);
errorvec  =  trimr ( ( omegastr - W )* zstar, 0, 1 );  /* Remove bottom row */
retp (errorvec);
endp;
/*************************************************************************/

ppi = xxx^2 | eye(1);
logofy = logystar(ppi,u);
ZGAIA = 10^(ones(m,1)*logalfap(ppi) + logofy);  /*  Formula is the same as for zstar */
ZGAIA = ZGAIA./ZGAIA[mrel,1];

/***********************************************************************/
/*      FINISH GEARY LOOP FOR THIS VALUE OF k */
IF PRINTIF==1;        /* Print results */
      screen on; format /rd 6,3;  print "z*:          " ZGAIA';    print "pi:          "  ppi';
ENDIF;
TIMETAKEN=HSEC-TIMETAKEN;   /* Calculate time taken on GAIA calculations */
screen off; print "\LTime taken for GAIA calculations:  " TIMETAKEN " 100ths of a second \L";
/*************************************************************************/
ZGAIAALL[.,k+1] = ZGAIA;   /* Read ZGAIA into appropriate column (N.B.) of summary matrix */
piall[.,k+1] = ppi;   /* Read ppi into appropriate column (N.B.) of summary matrix */
save ZGAIAALL;  /* Summary matrix saved each round, in case of crashes at higher values of k */
ENDO;     /*  END OF BIG DO LOOP  */
/***********************************************************************/
/*                      PRINT FINAL SUMMARY INFORMATION               */
IF PRINTIF9==1;
   screen on; format /rd 3,0;
   format /rd 6,3; print "Matrix of ppi:\L    k:    " seqa(0,1,n)' seqa(1,1,n)~piall;
   print "               SUMMARY TABLE (Relative to country number " mrel ")";
   print "               Expend.    Geary   GAIA[k=0,10], NSystem =" NSystem; 
   results=z~ZGeary~ZGAIAALL;  
   results = results./results[mrel,.];   /*  Reexpress all results relative to country # "mrel"  */
   format /rd 6,3;  print countrynum~results; 
   format /rd 3,0; print "            Expend.    Geary   GAIA[k=0,10], NSystem =" NSystem;
ENDIF;
/************************************************************************/
finish:
end;