home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-bin / lib / octave / 1.1.1 / m / statistics / gls.m < prev    next >
Encoding:
Text File  |  1996-10-12  |  1.9 KB  |  68 lines
  1. # Copyright (C) 1994, 1995 John W. Eaton
  2. # 
  3. # This file is part of Octave.
  4. # 
  5. # Octave is free software; you can redistribute it and/or modify it
  6. # under the terms of the GNU General Public License as published by the
  7. # Free Software Foundation; either version 2, or (at your option) any
  8. # later version.
  9. # 
  10. # Octave is distributed in the hope that it will be useful, but WITHOUT
  11. # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12. # FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  13. # for more details.
  14. # 
  15. # You should have received a copy of the GNU General Public License
  16. # along with Octave; see the file COPYING.  If not, write to the Free
  17. # Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  
  19. function [BETA, v, R] = gls (Y, X, O)
  20.   
  21. # usage: [BETA, v [,R]] = gls (Y, X, O)
  22. #
  23. # Generalized Least Squares (GLS) estimation for the multivariate model
  24. #
  25. #   Y = X*B + E,  mean(E) = 0,  cov(vec(E)) = (s^2)*O
  26. #
  27. # with Y ...  T x p      As usual, each row of Y and X is an observation
  28. #      X ...  T x k      and each column a variable.
  29. #      B ...  k x p
  30. #      E ...  T x p
  31. #      O ... Tp x Tp.
  32. #
  33. # BETA is the GLS estimator for B.
  34. # v is the GLS estimator for s^2.
  35. # R = Y - X*BETA is the matrix of GLS residuals.
  36.  
  37. # Written by Teresa Twaroch (twaroch@ci.tuwien.ac.at) May 1993.
  38. # Dept of Probability Theory and Statistics TU Wien, Austria.
  39.  
  40.   if (nargin != 3)
  41.     usage ("[BETA, v [, R]] = gls (Y, X, O)");
  42.   endif
  43.  
  44.   [rx, cx] = size (X);
  45.   [ry, cy] = size (Y);
  46.   if (rx != ry)
  47.     error ("gls: incorrect matrix dimensions");  
  48.   endif
  49.  
  50.   O = O^(-1/2);
  51.   Z = kron (eye (cy), X);
  52.   Z = O * Z;
  53.   Y1 = O * reshape (Y, ry*cy, 1);
  54.   U = Z' * Z;
  55.   r = rank (U);
  56.  
  57.   if (r == cx*cy)
  58.     B = inv (U) * Z' * Y1;
  59.   else
  60.     B = pinv (Z) * Y1;
  61.   endif
  62.  
  63.   BETA = reshape (B, cx, cy);
  64.   R = Y - X * BETA;
  65.   v = (reshape (R, ry*cy, 1))' * (O^2) * reshape (R, ry*cy, 1) / (rx*cy - r);
  66.  
  67. endfunction
  68.