/* 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;