/* Least Absolute Deviation
Author: Ron Schoenberg	rons@u.washington.edu
Date: May 29, 1995
Provided without guarantee for public non-commercial use.

Gauss code for least absolute deviation
estimation: given y_i and x_i, i = 1,...,n,
estimate b such that
                      \sum |y_i - b*x_i|
is minimized.
The following code solves this as a re-iterated weighted least squares
problem where the weights are the inverse of the absolute values of the
residuals.
*/
   rndseed 45456;
   nobs = 100;
   numx = 3;
   b = .5 * ones(numx+1,1);

   x = ones(nobs,1)~rndn(nobs,numx);
   y = x * b + rndn(nobs,1)^4; /* 4th power used for kurtotic error term */

   print "LS estimates   " (y / x)';
   print;

   b_old = 0;
   b_new = b;
   w = x;
   do until abs(b_new - b_old) < 1e-3;
       b_old = b_new;
       b_new = invpd(w'x)*(w'y);
       w = x./abs(y - x * b_new);
   endo;
   print "LAD estimates  "  b_new';



OUTPUT:

   LS estimates   3.6536  1.1954  0.0869  0.0421

   LAD estimates  0.8297  0.6028  0.3644  0.5249


These results support the contention that the LAD estimates are more stable
for fat-tailed errors.
