octave_packages/m/statistics/base/ols.m

   1 ## Copyright (C) 1996-2012 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
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License 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, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {[@var{beta}, @var{sigma}, @var{r}] =} ols (@var{y}, @var{x})
  21 ## Ordinary least squares estimation for the multivariate model
  22 ## @tex
  23 ## $y = x b + e$
  24 ## with
  25 ## $\bar{e} = 0$, and cov(vec($e$)) = kron ($s, I$)
  26 ## @end tex
  27 ## @ifnottex
  28 ## @w{@math{y = x*b + e}} with
  29 ## @math{mean (e) = 0} and @math{cov (vec (e)) = kron (s, I)}.
  30 ## @end ifnottex
  31 ##  where
  32 ## @tex
  33 ## $y$ is a $t \times p$ matrix, $x$ is a $t \times k$ matrix,
  34 ## $b$ is a $k \times p$ matrix, and $e$ is a $t \times p$ matrix.
  35 ## @end tex
  36 ## @ifnottex
  37 ## @math{y} is a @math{t} by @math{p} matrix, @math{x} is a @math{t} by
  38 ## @math{k} matrix, @math{b} is a @math{k} by @math{p} matrix, and
  39 ## @math{e} is a @math{t} by @math{p} matrix.
  40 ## @end ifnottex
  41 ##
  42 ## Each row of @var{y} and @var{x} is an observation and each column a
  43 ## variable.
  44 ##
  45 ## The return values @var{beta}, @var{sigma}, and @var{r} are defined as
  46 ## follows.
  47 ##
  48 ## @table @var
  49 ## @item beta
  50 ## The OLS estimator for @math{b}.
  51 ## @tex
  52 ## $beta$ is calculated directly via $(x^Tx)^{-1} x^T y$ if the matrix $x^Tx$ is
  53 ## of full rank.
  54 ## @end tex
  55 ## @ifnottex
  56 ## @var{beta} is calculated directly via @code{inv (x'*x) * x' * y} if the
  57 ## matrix @code{x'*x} is of full rank.
  58 ## @end ifnottex
  59 ## Otherwise, @code{@var{beta} = pinv (@var{x}) * @var{y}} where
  60 ## @code{pinv (@var{x})} denotes the pseudoinverse of @var{x}.
  61 ##
  62 ## @item sigma
  63 ## The OLS estimator for the matrix @var{s},
  64 ##
  65 ## @example
  66 ## @group
  67 ## @var{sigma} = (@var{y}-@var{x}*@var{beta})'
  68 ##   * (@var{y}-@var{x}*@var{beta})
  69 ##   / (@var{t}-rank(@var{x}))
  70 ## @end group
  71 ## @end example
  72 ##
  73 ## @item r
  74 ## The matrix of OLS residuals, @code{@var{r} = @var{y} - @var{x}*@var{beta}}.
  75 ## @end table
  76 ## @seealso{gls, pinv}
  77 ## @end deftypefn
  78
  79 ## Author: Teresa Twaroch <twaroch@ci.tuwien.ac.at>
  80 ## Created: May 1993
  81 ## Adapted-By: jwe
  82
  83 function [beta, sigma, r] = ols (y, x)
  84
  85   if (nargin != 2)
  86     print_usage ();
  87   endif
  88
  89   if (! (isnumeric (x) && isnumeric (y)))
  90     error ("ols: X and Y must be numeric matrices or vectors");
  91   endif
  92
  93   if (ndims (x) != 2 || ndims (y) != 2)
  94     error ("ols: X and Y must be 2-D matrices or vectors");
  95   endif
  96
  97   [nr, nc] = size (x);
  98   [ry, cy] = size (y);
  99   if (nr != ry)
 100     error ("ols: number of rows of X and Y must be equal");
 101   endif
 102
 103   if (isinteger (x))
 104     x = double (x);
 105   endif
 106   if (isinteger (y))
 107     y = double (y);
 108   endif
 109
 110   ## Start of algorithm
 111   z = x' * x;
 112   [u, p] = chol (z);
 113
 114   if (p)
 115     beta = pinv (x) * y;
 116   else
 117     beta = u \ (u' \ (x' * y));
 118   endif
 119
 120   if (isargout (2) || isargout (3))
 121     r = y - x * beta;
 122   endif
 123   if (isargout (2))
 124
 125     ## z is of full rank, avoid the SVD in rnk
 126     if (p == 0)
 127       rnk = columns (z);
 128     else
 129       rnk = rank (z);
 130     endif
 131
 132     sigma = r' * r / (nr - rnk);
 133   endif
 134
 135 endfunction
 136
 137
 138 %!test
 139 %! x = [1:5]';
 140 %! y = 3*x + 2;
 141 %! x = [x, ones(5,1)];
 142 %! assert (ols(y,x), [3; 2], 50*eps)
 143
 144 %!test
 145 %! x = [1, 2; 3, 4];
 146 %! y = [1; 2];
 147 %! [b, s, r] = ols (x, y);
 148 %! assert (b, [1.4, 2], 2*eps);
 149 %! assert (s, [0.2, 0; 0, 0], 2*eps);
 150 %! assert (r, [-0.4, 0; 0.2, 0], 2*eps);
 151
 152 %!test
 153 %! x = [1, 2; 3, 4];
 154 %! y = [1; 2];
 155 %! [b, s] = ols (x, y);
 156 %! assert (b, [1.4, 2], 2*eps);
 157 %! assert (s, [0.2, 0; 0, 0], 2*eps);
 158
 159 %!test
 160 %! x = [1, 2; 3, 4];
 161 %! y = [1; 2];
 162 %! b = ols (x, y);
 163 %! assert (b, [1.4, 2], 2*eps);
 164
 165 %% Test input validation
 166 %!error ols ();
 167 %!error ols (1);
 168 %!error ols (1, 2, 3);
 169 %!error ols ([true, true], [1, 2]);
 170 %!error ols ([1, 2], [true, true]);
 171 %!error ols (ones (2,2,2), ones (2,2));
 172 %!error ols (ones (2,2), ones (2,2,2));
 173 %!error ols (ones(1,2), ones(2,2));