/***
This file contains GAUSS routines for several rank-based estimators.
The first few procedures are the objective functions for these
Estimators. Following the objective functions, there are procedures
for optimization. The procedures "mre" and "mrc" at the end of this
file are sample calls to the (Nelder-Mead) optimization routines.
LAST REVISION: June 1, 1999
Please e-mail any comments and/or questions to Jason Abrevaya at
jason.abrevaya@gsb.uchicago.edu.
***/
/* global variables for the dependent and independent variables */
clear dep, indep;
/* global variables for the Nelder-Mead simplex optimization method */
clear p, y, ndim, psum, ftol, lam1, lam2, nfunk;
/* global variables for the maximum rank correlation estimator */
clear mrcflag, mrctree, mrcroot;
/********
OBJECTIVE FUNCTIONS
********/
/*
mrcobj --- procedure that calculates the objective function for the maximum
rank correlation estimator of Han (1987); the algorithm used in
this procedure, which requires O(n log n) calculations, is
adapted from Abrevaya (1999, Economics Letters)
Argument: theta -- coefficient vector
Globals: dep -- vector of dependent variables (nobs x 1)
indep -- matrix of independent variables (nobs x k)
mrcflag -- flag (1 indicates tree structure has been created)
mrctree -- matrix which stores all the tree information
mrcroot -- node at the top of the tree
Output: the value of the maximum rank correlation objective function
Note: If there are no ties in the index values, columns 6-9 of
mrctree are not needed (and the code below could be modified
to compute the objective function faster).
*/
proc mrcobj(theta);
local nobs, temp, leftlim, rightlim, node, i, retval;
nobs = rows(dep);
if mrcflag==0; /* create tree if not created already */
mrcflag = 1;
/* form vector of unique (& sorted) y values */
mrctree = unique(sortc(dep,1),1);
/* form blank array for the tree data
2nd column --- left child
3rd column --- right child
4th column --- number of hits
5th column --- number of times down left subtree
6th column --- maximum index down left subtree
7th column --- number of times with max index down left subtree
8th column --- maximum index for hits
9th column --- number of times with max index for hits */
mrctree = mrctree~zeros(rows(mrctree),8);
mrcroot = trunc((1+rows(mrctree))/2);
temp = mrcroot~1~rows(mrctree);
do while temp/=0;
node = temp[1,1];
leftlim = temp[1,2];
rightlim = temp[1,3];
if node>leftlim;
mrctree[node,2] = trunc((leftlim+(node-1))/2);
temp = temp|( trunc((leftlim+(node-1))/2)~leftlim~(node-1) );
endif;
if nodemrctree[node,6]);
mrctree[node,6] = temp[i,2];
mrctree[node,7] = 1;
else;
mrctree[node,7] = mrctree[node,7]+1;
endif;
node = mrctree[node,2];
else;
retval = retval + mrctree[node,4] + mrctree[node,5] - (temp[i,2]==mrctree[node,6]).*(mrctree[node,7]) - (temp[i,2]==mrctree[node,8]).*(mrctree[node,9]);
node = mrctree[node,3];
endif;
endo;
if (mrctree[node,9]==0) or (temp[i,2]>mrctree[node,8]);
mrctree[node,8] = temp[i,2];
mrctree[node,9] = 1;
else;
mrctree[node,9] = mrctree[node,9]+1;
endif;
mrctree[node,4] = mrctree[node,4]+1;
retval = retval + mrctree[node,5] - (temp[i,2]==mrctree[node,6]).*(mrctree[node,7]);
i = i+1;
endo;
retp(-retval);
endp;
/*
mrcobj2 --- procedure that calculates the objective function for the maximum
rank correlation estimator of Han (1987); the algorithm used in
this procedure, which requires O(n^2) calculations, takes
advantage of Gauss's quick handling of vector calculations
(unlike mrcobj3 below)
Argument: theta -- coefficient vector
Globals: dep -- vector of dependent variables (nobs x 1)
indep -- matrix of independent variables (nobs x k)
Output: the value of the maximum rank correlation objective function
*/
proc mrcobj2(theta);
local i, j, nobs, retval, temp;
nobs = rows(dep);
/* sort based on index values to make the double-loop faster */
temp = sortc(dep~(indep*theta),2);
retval = 0;
i = 1;
do while i<=(nobs-1);
j=i+1;
/* skip over index values which tie the i-th index value */
do while temp[j,2]==temp[i,2];
j=j+1;
if j==(nobs+1);
break;
endif;
endo;
/* loop through all larger index values & count how many dep. variables are larger */
do while j<=nobs;
if temp[j,1]>temp[i,1];
retval = retval+1;
endif;
j=j+1;
endo;
i = i+1;
endo;
retp(-retval);
endp;
/*
mrcobj3 --- procedure that calculates the objective function for the maximum
rank correlation estimator of Han (1987); the algorithm used in
this procedure, which requires O(n^2) calculations, does not
take advantage of Gauss's quick handling of vector calculations
(i.e., slower than mrcobj2 above)
Argument: theta -- coefficient vector
Globals: dep -- vector of dependent variables (nobs x 1)
indep -- matrix of independent variables (nobs x k)
Output: the value of the maximum rank correlation objective function
*/
proc mrcobj3(theta);
local i, nobs, retval, temp;
nobs = rows(dep);
/* sort based on index values to make the double-loop faster */
temp = sortc(dep~(indep*theta),2);
/* loop through the data to calculate the objective function */
retval = 0;
i = 1;
do while i<=(nobs-1);
retval = retval + (temp[i+1:nobs,2].>temp[i,2])'*(temp[i+1:nobs,1].>temp[i,1]);
i = i+1;
endo;
retp(-retval);
endp;
/*
uniqrank --- procedure that returns the ranks of a numerical vector; unlike
the built-in Gauss procedure rankindx, this procedure gives tied
values the same rank (e.g., the vector 1|1|2|2|3 has associated
ranks 1|1|3|3|5)
Argument: vect -- numerical vector
Output: the ranks of the entries in vect (in the same order as the input)
*/
proc uniqrank(vect);
local nobs, temp, rankvect, i, numties, retval;
nobs = rows(vect);
temp = sortc(vect~seqa(1,1,nobs),1);
rankvect = zeros(nobs,1);
i = 1;
numties = 0;
do while itemp[j,1];
retval = retval+1;
elseif temp[i,2]>temp[j,2] and temp[i,1]=5000;
"Error: NMAX exceeded!";
endif;
nfunk = nfunk+2;
ytry = amotry(&f,ihi,-1);
if ytry=y[inhi];
ysave = y[ihi];
ytry = amotry(&f,ihi,0.5);
if ytry>=ysave;
i = 1;
do while i<=ndim+1;
if i/=ilo;
p[i,.] = norm(0.5*(p[i,.]'+p[ilo,.]'))';
y[i] = f(p[i,.]');
endif;
i = i+1;
endo;
psum = sumc(p);
/* psum = norm(psum); */
nfunk = nfunk+ndim;
endif;
else;
nfunk = nfunk-1;
endif;
endo;
retp(nfunk);
endp;
/*
* amotry -- Nelder-Mead "amotry" procedure for simplex method (converted
* to GAUSS from Numerical Recipes code)
*/
proc amotry(&f,ihi,fac);
local j, fac1, fac2, ytry, ptry;
local f:proc;
fac1 = (1-fac)/ndim;
fac2 = fac1-fac;
ptry = psum*fac1 - p[ihi,.]'fac2;
ptry = norm(ptry);
ytry = f(ptry);
if ytry