/* This is file GearyLib.SRC: It contains the definitions of procedures used in Geary project;
they are listed in my library C:\GAUSS\LIB\GEARYLIB.LCG.
BBETAOFP: For AIDS: Given an n-by-m price matrix p and an (n-1)-by-1 vector beta,
compute the m-by-1 vector of beta(p) functions [typical element: PI_i { p_i ^ beta_i }]
in each country
CHECKERR: Check two matrices for conformability and equality; print both if they differ.
CORRXY: Calculate h-by-1 vector of squared corr. coefficients between a given m-by-h
matrix x (which gives the values of h indexes for m countries) and a given m-by-1 vector y
GEARYINC: Solve M.theta = theta for Geary real incomes and calculate world prices.
GINICOEF: Given an m-by-h matrix x, which gives the values of h indexes for m countries,
calculate the h-by-1 vector of Gini coefficients for the columns of x.
LAMOFP: For QUAIDS: Given an n-by-m price matrix p, compute the m-by-1 vector of
lambda(p) functions [typical element: sigma_i { lambda_i * logp_ij }] in each country.
LOGALFAP: Calculates QUAIDS "subsistence" price index log{ alpha(p) }.
LOGPREL: Given an n-by-m price matrix p, compute the (n-1)-by-m matrix of log
prices, all relative to the price of good n.
LOGYSTAR: For QUAIDS: Given an n-by-m2 price matrix p, and an m1-by-1 vector of
utility levels u, compute the m1-by-m2 vector of u.beta/(1-u.lambda) functions.
PRINTOUT: Given a h-by-m matrix of h indexes "results", calculate summary statistics,
and (optionally) print summary table and Figures 4-6.
READDATA : Read n-by-m p (price) and q (quantity) matrices from files.
SORTDATA: Calculate expenditure vector z and budget share matrix w; sort p and q
matrices by expenditure; and scale p and z by sample mean.
TRANSFER: Calculate the h-by-1 vector of transfers implied by the UN benchmark
for foreign aid (0.7% of GNP) for the columns of a given m-by-h matrix x, which gives
the values of h indexes for m countries.
WINDOW1: Procedure (with zero returns) which sets up 1-window graph
WINDOW4: Procedure (with zero returns) which sets up 4-window graph
(c) J.P. Neary 1999
*/
/* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX */
/* PROC BBETAOFP: For AIDS: Given an n-by-m price matrix p and an (n-1)-by-1 vector beta,
compute the m-by-1 vector of beta(p) functions [typical element:
PI_i { p_i ^ beta_i }] in each country
N.B. Double-b "bbetaofp" to distinguish it from betaofp vector in StQua_It program */
PROC BBETAOFP(p);
local betan;
betan = beta | - sumc(beta); /* Calculate beta for good #n by subtraction: Repeat here to avoid
errors in QUAIDS program */
retp( prodc(p.^betan) );
endp;
/* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX */
/* PROC CHECKERR: Check two matrices for conformability and equality; print both if they differ.
To invoke this proc, copy the following line and amend the matrix arguments appropriately
call checkerr (x1,x2); /* CHECK FOR DIFFERENCE BETWEEN MATRICES*/
*/
PROC(0) = CHECKERR(x1,x2);
local errvec, err;
screen on ; format /rd 5,0;
print "PROC CHECKERR: Compare Matrices: X1 is " rows(x1) " by" cols(x1) ", X2 is " rows(x2) " by" cols(x2);
format /rd 7,3;
IF rows(x1)==rows(x2);
IF cols(x1)==cols(x2);
GOTO NEXTBIT;
ENDIF;
ENDIF;
print "Matrices not conformable";
retp;
NEXTBIT:
errvec= vec(x1-x2);
err=errvec'*errvec;
IF err > 10^(-16);
print "First Matrix is: " x1; print "Second Matrix is:" x2;
format /rd 20,14;
print "Problem: Matrices differ. Average absolute deviation is: " sqrt(err/ (rows(x1)*cols(x1)));
retp;
ENDIF;
print "Matrices are identical to within a sum of squared differences of 10^(-16)";
retp;
endp;
/* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX */
/* PROC CORRXY: Calculate h-by-1 vector of squared corr. coefficients between a given m-by-h
matrix x (which gives the values of h indexes for m countries) and a given m-by-1 vector y */
PROC (1) = CORRXY (x,y);
local h, result, jj, temp;
h = cols(x); /* # of indexes */
result=zeros(h,1);
jj=1;
DO WHILE jj <= h;
temp = corrx(x[.,jj]~y);
result[jj,1] = temp[1,2]^2; /* N.B. r squared */
jj=jj+1;
ENDO;
retp(result);
endp;
/* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX */
/* Proc GEARYINC: Solve M.theta = theta for Geary real incomes and calculate world prices.
Input: n-by-m matrices of budget shares W and quantities Q.
W is always the matrix of actual budget shares.
If Q is matrix of actual quantities, then output is the Geary real incomes and world prices;
If Q is matrix of virtual quantities, output is GAIA real incomes and Geary-Konus world prices; */
/**************************************************************************/
PROC(2) = Gearyinc(W,Q);
local n, m, qtot, qrel, MM, Timenow, ZGeary , ppi, qtothat ;
n=rows(Q);
m=cols(Q);
/**************************************************************************/
Timenow = hsec;
qtot = Q*ones(m,1); /* Form n-by-1 vector of total world consumption levels of each good */
qtothat = diagrv (zeros(n,n),qtot) ; /* Form diagonal matrix with vector qhat on princ. diag. */
MM = Q' (W/qtothat); /* m-by-m matrix whose largest eigenvector equals the Geary relative incomes */
/* Solve M.theta=theta for ZGeary, with element #m normalised to one */
ZGeary = (trimr( (trimr(MM',m-1,0))', 0,1) / (eye(m-1) - trimr( trimr(MM',0,1)', 0,1) ) ) | eye(1);
ppi = (W/qtothat) * ZGeary;
/**************************************************************************/
/* PRINT RESULTS */
screen Off; format /rd 5,4; print "Time taken (seconds): " (hsec - Timenow) ;
print "ppi: " ppi';
ppi = ppi ./ ppi[n,1]; /* Normalise world prices by setting that of good n equal to one: expenditures are not affected */
print "ppi: " ppi';
print "ZGeary: " ZGeary';
/**************************************************************************/
finish:
retp(ZGeary, ppi);
endp;
/* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX */
/* PROC GINICOEF: Given an m-by-h matrix x, which gives the values of h indexes for m countries,
calculate the h-by-1 vector of Gini coefficients for the columns of x */
PROC (1) = GINICOEF (x);
local h, m, result, k, jj, xtemp;
m=rows(x); /* # of countries */
h=cols(x); /* # of indexes */
result=zeros(h,1);
k = seqa(1,1,m);
jj=1;
DO WHILE jj <= h;
xtemp = rev( sortc(x[.,jj],1) ); /* Sort column jj in Descending order */
result[jj,1] = (m + 1 - (2*(k'*xtemp))/(m*meanc(x[.,jj]))) / (m-1);
jj=jj+1;
ENDO;
retp(result);
endp;
/* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX */
/* PROC LAMOFP: For QUAIDS: Given an n-by-m price matrix p, compute the m-by-1 vector of
lambda(p) functions [typical element: sigma_i { lambda_i * logp_ij }] in each country */
PROC LAMBOFP(p);
retp( logprel(p)'*lambda );
endp;
/* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
Proc LOGALFAP: Calculates QUAIDS "subsistence" price index log{ alpha(p) } */
PROC logalfap(p);
/* Input: p is an arbitrary n-by-mm price matrix; mm is calculated by this procedure (i.e., NOT
inputted), so when "p" equals pi, it is set equal to one.
Globals: n; m; alpha is an (n-1)-by-1 vector of coefficients;
ggamma is an (n-1)-by-(n-1) matrix of gamma coefficients.
Returns: logalfap: the mm-by-1 vector of alpha(p) functions, one for each country */
local alphan, logprn, jjj, gammabit, mm, answer;
mm = cols(p);
gammabit=zeros(mm,1); /* mm-by-1 vector of portion of alpha(p) consisting of quadratic form in
p'*gamma*p: this has to be calculated separately each time */
alphan = alpha | 1-sumc(alpha); /* Full n-by-1 vector of alpha's adds to one (not zero)
so it must be calculated; it would not be right to have "(logprn)'*alpha" in alphaofp below */
/* Compute logprn: (n-1)-by-mm matrix of log prices, relative to price of good n */
logprn=log(p[1:n-1,.])-log(p[n,.]);
jjj=1; /* Country index for next loop */
DO WHILE jjj<=mm;
gammabit[jjj,1]=(logprn[.,jjj])'*ggamma*logprn[.,jjj]; /* Quadratic form in ggamma for jjj */
jjj=jjj+1;
ENDO;
answer = (log(p))'*alphan + gammabit/2;
retp(answer);
endp;
/* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX */
/* PROC LOGPREL: Given an n-by-m price matrix p, compute the (n-1)-by-m matrix of log
prices, all relative to the price of good n */
PROC LOGPREL(p);
retp( log(p[1:n-1,.]) - log(p[n,.]) );
endp;
/* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX */
/* PROC LOGYSTAR: For QUAIDS: Given an n-by-m2 price matrix p, and an m1-by-1 vector of
utility levels u, compute the m1-by-m2 vector of u.beta/(1-u.lambda) functions.
This is mainly for use with pi as the (n-by-1) price vector (so m2=1).
With m2=m1=m, this is an input to the star matrix of Allen indexes */
PROC (1) = LOGYSTAR(p,u);
retp( (u*bbetaofp(p)')./(1 - u*lambofp(p)') );
endp;
/* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX */
/* PROC PRINTOUT: Given a h-by-m matrix of h indexes "results", calculate summary statistics,
and (optionally) print summary table and Figures 4-6. */
PROC (0) = printout(results, EKS, ZGeary, mrel, printifS,graphif);
Local m, countrynum, CV, CorrEKS, CorrGear, Gini, Trans ;
m = rows(results);
countrynum=seqa(1,1,m);
CV = (stdc(results)./meanc(results));
CorrEKS = (corrxy(results,EKS)) ;
CorrGear = (corrxy(results,ZGeary));
Gini = GiniCoef (results);
Trans = Transfer (results);
/************************************************************************/
IF PRINTIFS==1;
screen on; format /rd 3,0;
print " SUMMARY TABLE (Relative to country number " mrel ")";
print " Expend. EKS Geary GAIA";
format /rd 6,3; /* This format is best for table; for input to excel graph use next line */
format /rd 10,7;
print countrynum~results;
print " Expend. EKS Geary GAIA";
print "Average: " (meanc(results))';
print "Std. Dev.: " (stdc(results))';
print " CV: " CV';
format /rd 6,4; /* This format is best for table; for input to excel graph use next line */
format /rd 10,7;
print "Corr(EKS):" CorrEKS';
print "Corr(G): " CorrGear';
print "Gini: " Gini';
format /rd 6,3; /* This format is best for table; for input to excel graph use next line */
format /rd 10,7;
print "Transfer: " Trans';
ENDIF;
/***********************************************************************/
IF GRAPHIF==0; /* SET UP GRAPH (PROVIDED "GRAPHIF">0) */
retp;
ENDIF;
graphset;
fonts ("complex simgrma"); /* load fonts: \203 (standard, with serif), \202 (greek) */
_pdate = ""; /* "": don't print date; to print date, type " " */
call window4;
/***********************************************************************/
/* Print Fig. 4: Correlations with EKS and Geary Indexes */
title ("Fig. 4: Correlations with EKS and Geary Indexes");
xlabel ("Correlation with EKS Index");
ylabel ("Correlation with Geary Index");
xtics(.9935, 1.000, .0005, 1); ytics( .9935,1.000, .0005, 1); /* Axis tick marks: symmetric.
.9935 needed to show all values of r^2 (.9965 for r). */
axmargin (2.25,1,0,1 ) ; /* Centre graph in window */
_plotsiz = {5 5}; /* Ensure a square plot size */
_plctrl = -1; /* plot data points only; no connecting lines */
xy (CorrEKS, CorrGear);
/***********************************************************************/
/* Print Fig. 5: Gini plotted against CV */
nextwind;
graphset; /* Need to reset, since a square plot is not required this time */
fonts ("complex simgrma"); /* load fonts: \203 (standard, with serif), \202 (greek) */
title ("Fig. 5: Gini Coefficient and Coefficient of Variation");
xlabel ("Coefficient of Variation");
ylabel ("Gini Coefficient");
_plctrl = -1; /* plot data points only; no connecting lines */
xy (CV, Gini);
/************************************************************************/
/* Print Fig. 6: Gini plotted against Transfer */
nextwind; /* No need to reset with "graphset;" this time: same specs. as Fig. 5 */
title ("Fig. 6: Gini Coefficient and Implied Transfer");
xlabel ("Implied Transfer (in billions US$)");
ylabel ("Gini Coefficient");
_plctrl = -1; /* plot data points only; no connecting lines */
xy (Trans, Gini);
/***********************************************************************/
endwind;
retp;
endp;
/* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
Proc READDATA: Read n-by-m p (price) and q (quantity) matrices from files;
small files for TESTIF=1; general ones for 1970, 1975, 1980 or 1985 otherwise */
proc (3) = readdata(n,nyear,testif); /* Input: n (# commodities), nyear and testif (=1 for test only) */
local m, p, q; /* Returns: m (# countries); p; q */
LOAD PATH = c:\me\research\geary2\empirics\gaussin ; /* Read p and q from this directory */
screen on; format /rd 4,0;
IF TESTIF==1;
n=11; /* Do not alter this: Fixed for small sample data set */
m=6; /* do. */
nyear=1970; /* do. */
LOAD p[n,m] = 70p-6.PRN; /* p: 11-by-6 */
LOAD q[n,m] = 70q-6.PRN; /* q: 11-by-6 */
output file = C:\ME\RESEARCH\GEARY2\EMPIRICS\GAUSSOUT\TEST70-6.OUT reset;
goto finish;
ENDIF;
IF NYEAR==1970; /* Value of NYEAR determines where data are read and how many countries */
m=16; /* m = # countries */
LOAD p[n,m] = 70P-16.PRN; /* p: n-by-m */
LOAD q[n,m] = 70Q-16.PRN; /* q: n-by-m */
output file = C:\ME\RESEARCH\GEARY2\EMPIRICS\GAUSSOUT\1970-16.OUT reset;
goto finish;
ENDIF;
IF NYEAR==1975; /* Value of NYEAR determines where data are read and how many countries */
m=34; /* m = # countries */
LOAD p[n,m] = 75P-34.PRN; /* p: n-by-m */
LOAD q[n,m] = 75Q-34.PRN; /* q: n-by-m */
output file = C:\ME\RESEARCH\GEARY2\EMPIRICS\GAUSSOUT\1975-34.OUT reset;
goto finish;
ENDIF;
IF NYEAR==1980; /* Value of NYEAR determines where data are read and how many countries */
m=60; /* m = # countries */
LOAD p[n,m] = 80P-60.PRN; /* p: n-by-m */
LOAD q[n,m] = 80Q-60.PRN; /* q: n-by-m */
output file = C:\ME\RESEARCH\GEARY2\EMPIRICS\GAUSSOUT\1980-60.OUT reset;
goto finish;
ENDIF;
IF NYEAR==1985; /* Value of NYEAR determines where data are read and how many countries */
m=64; /* m = # countries */
LOAD p[n,m] = 85P-64.PRN; /* p: n-by-m */
LOAD q[n,m] = 85Q-64.PRN; /* q: n-by-m */
output file = C:\ME\RESEARCH\GEARY2\EMPIRICS\GAUSSOUT\1985-64.OUT reset;
goto finish;
ENDIF;
Print "Error"; /* End program if NYEAR not known */
finish:
print "\L RESULTS FOR " nyear " DATA WITH" n " GOODS AND" m " COUNTRIES \L";
retp(m,p,q);
endp;
/* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
Proc SORTDATA: Calculate expenditure vector z and budget share matrix w; sort p and q
matrices by expenditure; and scale p and z by sample mean.
Returns: p: n-by-m matrix, sorted by country expenditure and scaled appropriately ;
q: n-by-m matrix, sorted by country expenditure and scaled appropriately ;
pprimeq: m-by-m matrix: (k,j) element gives country j's q's valued at k's p's (all scaled);
z: m-by-1 vector of country expenditures, in descending order and scaled appropriately ;
zraw: m-by-1 vector equal to z but NOT scaled ;
w: n-by-m matrix of budget shares, sorted by country expenditure;
w_trim: (n-1)-by-m matrix of budget shares on all commodities except #n;
poverz: (For LES) n-by-m matrix of prices deflated by expenditure;
poz_trim: (For LES) (n-1)-by-m matrix of prices deflated by expenditure, except #n. */
PROC(9) = SORTDATA(p,q, printif);
local m, pprimeq, z, zraw, expij, w, w_trim, poverz, poz_trim, pscale, zscale, countrynum;
local znew;
m=cols(p); /* m = # countries = # columns of p (and also q) */
/***********************************************************************/
/* SORT P AND Q MATRICES BY COUNTRY EXPENDITURE */
pprimeq=p'*q; /* calculate m-by-m P'Q matrix: (k,j) element gives country j's q's valued at k's p's */
z=diag(pprimeq); /* pick diagonal of P'Q: m-by-1 vector of country expenditures */
p=trimr((rev(sortc(z~p',1)))',1,0); /* sort columns (N.B.) of p in descending (N.B.) order of z */
q=trimr((rev(sortc(z~q',1)))',1,0); /* sort columns (N.B.) of q in descending (N.B.) order of z */
pprimeq=p'*q; /* recalculate P'Q matrix for new ranking of countries */
z=diag(pprimeq); /* do. for z */
zraw = z ; /* Save unscaled expenditures */
/***********************************************************************/
/* CALCULATE EXPENDITURE AND AVERAGE BUDGET SHARE MATRICES */
expij=p.*q; /* Calculate n-by-m matrix of expenditure on each commodity */
w = expij./z'; /* Deflate expend. by country expend. to get n-by-m matrix of budget shares */
w_trim=trimr(w,0,1); /* (n-1)-by-m matrix of budget shares on all commodities except #n */
/* N.B. Budget shares are independent of normalisation */
/***********************************************************************/
/* For Semiflexible QUAIDS: Scale p and z, and then scale q compatibly */
pscale = meanc(p'); /* scaling factor for prices: set equal to one at sample mean */
zscale = meanc(z); /* scaling factor for expenditure: set equal to one at sample mean */
/* pscale = ones(n,1); /* Default scaling factor for prices: leave unchanged */
zscale = ones(1,1); /* Default scaling factor for expenditure: leave unchanged */ */
p = p./pscale; /* scale prices */
z = z./zscale; /* scale expenditure */
q = ( q.* pscale ) ./ zscale ; /* Scale quantities compatibly */
pprimeq=p'*q; /* recalculate P'Q matrix for scaled data */
/***********************************************************************/
/* For LES: Calculate full and trimmed matrices of prices deflated by expenditure */
poverz = p./z';
poz_trim = trimr(poverz,0,1);
/***********************************************************************/
/* Print p, q, z and w if PRINTIF=1 */
IF PRINTIF==1;
countrynum = seqa(1,1,m);
SCREEN On; format /rd 8,3; print "z and q:" countrynum~z~q';
print "p:" countrynum~p';
print "w:" countrynum~w';
SCREEN OFF; print poverz'; print "poz_trim prime :" poz_trim';
ENDIF;
/***********************************************************************/
retp(p, q, pprimeq, z, zraw, w, w_trim, poverz, poz_trim);
endp;
/* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX */
/* PROC TRANSFER: Calculate the h-by-1 vector of transfers implied by the UN benchmark
for foreign aid (0.7% of GNP) for the columns of a given m-by-h matrix x, which gives
the values of h indexes for m countries. */
PROC (1) = TRANSFER (x);
local h, m, result, k, jj, xtemp, share, mrelcons, UScons, mUS, OECDyes, OECDind, OECDpop,
Pop1980;
m=rows(x); /* # of countries */
h=cols(x); /* # of indexes */
result=zeros(h,1);
IF TESTIF==1; /* Proc does not work on test data set (TESTIF=1) */
retp(result);
ENDIF;
share = 0.007; /* fraction of income to be transferred */
mrelcons = .11256; /* per capita consumption (US$'000's) in country # mrel (for which z=1.0)*/
UScons=7.91; /* per capita consumption (US$'000's) in US (country # = 5)*/
mUS = 5; /* Country # of US, relative to which results are to be expressed */
/**************************************************************************/
/* Generate m-by-1 vector OECDyes, whose entries equal 1 if country is an OECD member */
OECDyes = zeros(m,1);
OECDind = {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 20 23}; /* Values of countrynum for OECD
members */
k=cols(OECDind); /* Calculate number of OECD members in the sample [=18 in 1980, m=60] */
jj=1;
DO WHILE jj<= k;
OECDyes[OECDind[1,jj],1] = 1;
jj = jj+1;
ENDO;
/**************************************************************************/
load path = c:\me\research\geary2\empirics\gaussin;
load pop1980; /* m-by-1 vector of population levels (in millions) in 1980 */
OECDpop = Pop1980 .* OECDyes; /* = 1980 Population for OECD members; = 0 otherwise */
/**************************************************************************/
xtemp = (x./x[mUS,.])*UScons; /* Express all indexes in terms of per capita cons. ('000 US$) */
xtemp= xtemp.*OECDpop; /* Express all indexes in terms of total cons. (billions US$) */
xtemp = xtemp*share; /* Final answer is a fraction "share" of this (billions US$) */
/* format /rd 8,3; print countrynum~x~xtemp; */ /* Decomment to check */
result = sumc(xtemp);
finish:
retp(result);
endp;
/* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX */
/* PROC WINDOW1: Procedure (with zero returns) which sets up 1-window graph */
proc(0) = window1;
BEGWIND; /* Set up 1 window */
MAKEWIND(11.69,8.27,0,0,0); /* Create a single window the full size of the page (allowing
overlapping graphs) */
_pdate = "\201JPN "; /* Prints date string if any characters are in quotes */
endp;
/* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX */
/* PROC WINDOW4: Procedure (with zero returns) which sets up 4-window graph */
proc(0) = window4;
BEGWIND; /* Set up 4 windows */
WINDOW(2,2,0);
_pdate = "" ; /*"\201JPN"; */ /* Prints date string if any characters are in quotes */
SETWIND(1); /* First Window */
endp;