/*   This is file Basics.GAU: Calculates basic indexes.
	This program reads p and q matrices from disk; then it calculates Laspeyres, EKS, 
        CCD and Geary real incomes and writes results to disk.  

(c) J.P. Neary 1999

*/
/*	It also reads GAIA real incomes for all 3 systems (created by GaiaQuai.gau) from disk.
	Then it calculates all correlations etc. and prints Figures 4-6.  Started 8.11.99   */
new;  /*************************************************************************/
/*          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  */
library gearylib, pgraph;
TESTIF=0;            /* TESTIF=1 on a test run (using small data set only)  */
printsor=0;              /* printsor= 0 for no printing of data in PROC SORTDATA; 1 for printing */
printifS=1;       /* printifS=1 for printing of summary table */
GRAPHIF=1; /* graphif=0: no graph; >1: print graphs */
/*************************************************************************/
 
save path = c:\me\research\geary2\empirics\gaussout\ALL ;  /* Output path here */
/*************************************************************************/
/* These variables are not needed; but for the present they should be defined, since
      otherwise there are conflicts with GEARY.SRC (where they are assumed as globals).
      Good idea to fix this sometime.   */
alpha=zeros(n-1,1);    /* (n-1)-by-1 vector of ML alpha's:  n'th to be estimated residually */
beta=zeros(n-1,1);      /* (n-1)-by-1 vector of ML beta's:  n'th to be estimated residually */
lambda=zeros(n-1,1);  /* (n-1)-by-1 vector of ML lambdas: n'th to be estimated residually */
ggamma=zeros(n-1,n-1);   /* (n-1)-by-(n-1) matrix of ML gamma's:  n'th to be estimated residually */
/*************************************************************************/
{ 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 */
mrel=m;    /* Index of country relative to which the results are printed */
/***********************************************************************/
/* 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 LASPEYRES AND FISHER INDEXES */
L=pprimeq./z;         /* calculate Laspeyres Q-Matrix ;  N.B. NOT transpose of pprimeq */
Star = L./L[.,mbase];   /* Relative Laspeyres Matrix */
F = sqrt(L.*(1./L'));	/* calculate Fisher q-Matrix: {j,k} gives income of k relative to j */
/***********************************************************************/
/*  Calculate Maximum and Minimum of Laspeyres Star Matrix */  
star=star./star[.,mrel];  /* N.B. Only calculate Starmax AFTER deflation by mrel */
StarMax=Maxc(star);  /* Maximum of Star Matrix */
StarMin=Minc(star);  /* Minimum of Star Matrix */
StarMaxI=Maxindc(star);  /* Index of Maximum of Star Matrix */
StarMinI=Minindc(star);    /* Index of Minimum of Star Matrix */
/*************************************************************************/
/*                              CALCULATE EKS RELATIVE INCOME INDEXES   */
FRel= F ./ F[.,mrel] ;   /* m-by-m Matrix of log of Fisher indexes relative to country mrel */
EKS=exp(sumc( ln(FRel) )/m);     /* m-by-1 vector of EKS indexes */
QRel = ( ln(Q) - ln(Q[.,mbase]) ) ;
/*************************************************************************/
/*           CALCULATE CCD RELATIVE INCOME INDEXES   */
/*  First, calculate m-by-m Tornqvist matrix: {j,k} gives country j's income relative to k's */
TT = zeros(m,m);
jj=0; DO UNTIL jj==m;  jj=jj+1;
             kk=0; DO UNTIL kk==m; kk=kk+1;
                              ii=0; DO UNTIL ii==n; ii=ii+1;
			             TT[jj,kk] = TT[jj,kk] + 0.5* (W[ii,jj]+W[ii,kk]) * ln (Q[ii,jj]/Q[ii,kk]);
	                            ENDO;
                        ENDO;
ENDO;
TT = exp(TT);
TTRel = TT./TT[mrel,.];
CCD = exp ( sumc (ln(TTRel')/m) );
/*************************************************************************/
/*  CALCULATE GEARY INCOMES  */
{ ZGeary, ppi } = Gearyinc(w,q);   /* Call PROC Gearyinc to calculate Geary incomes */
/*************************************************************************/
/*  SAVE RESULTS TO DISK */
save z;
save StarMax;
save StarMin;
save EKS;
save CCD;
save ZGeary;
save ppi;
/*************************************************************************/
/*                      PRINT FINAL SUMMARY INFORMATION               */
/*   Assemble the key matrix "results" to be printed and charted  */
/* Decomment whichever of the next lines is desired:
results = z~EKS~ZGeary;
results = CCD;
results = Star' ~EKS~CCD~ZGeary;   /* Print all Laspeyres Star matrix (N.B. Transpose) */
*/
results = z~StarMax~StarMin~EKS~CCD~ZGeary;   /* Figure 1 */
/*************************************************************************/
results = results./results[mrel,.];   /*  Reexpress all results relative to country indexed "mrel"  */
call printout (results, EKS, ZGeary, mrel, printifS,graphif);
/*************************************************************************/
finish:
end;