octave_packages/statistics-1.1.3/copulacdf.m

   1 ## Copyright (C) 2008 Arno Onken <asnelt@asnelt.org>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {@var{p} =} copulacdf (@var{family}, @var{x}, @var{theta})
  18 ## @deftypefnx {Function File} {} copulacdf ('t', @var{x}, @var{theta}, @var{nu})
  19 ## Compute the cumulative distribution function of a copula family.
  20 ##
  21 ## @subheading Arguments
  22 ##
  23 ## @itemize @bullet
  24 ## @item
  25 ## @var{family} is the copula family name. Currently, @var{family} can
  26 ## be @code{'Gaussian'} for the Gaussian family, @code{'t'} for the
  27 ## Student's t family, @code{'Clayton'} for the Clayton family,
  28 ## @code{'Gumbel'} for the Gumbel-Hougaard family, @code{'Frank'} for
  29 ## the Frank family, @code{'AMH'} for the Ali-Mikhail-Haq family, or
  30 ## @code{'FGM'} for the Farlie-Gumbel-Morgenstern family.
  31 ##
  32 ## @item
  33 ## @var{x} is the support where each row corresponds to an observation.
  34 ##
  35 ## @item
  36 ## @var{theta} is the parameter of the copula. For the Gaussian and
  37 ## Student's t copula, @var{theta} must be a correlation matrix. For
  38 ## bivariate copulas @var{theta} can also be a correlation coefficient.
  39 ## For the Clayton family, the Gumbel-Hougaard family, the Frank family,
  40 ## and the Ali-Mikhail-Haq family, @var{theta} must be a vector with the
  41 ## same number of elements as observations in @var{x} or be scalar. For
  42 ## the Farlie-Gumbel-Morgenstern family, @var{theta} must be a matrix of
  43 ## coefficients for the Farlie-Gumbel-Morgenstern polynomial where each
  44 ## row corresponds to one set of coefficients for an observation in
  45 ## @var{x}. A single row is expanded. The coefficients are in binary
  46 ## order.
  47 ##
  48 ## @item
  49 ## @var{nu} is the degrees of freedom for the Student's t family.
  50 ## @var{nu} must be a vector with the same number of elements as
  51 ## observations in @var{x} or be scalar.
  52 ## @end itemize
  53 ##
  54 ## @subheading Return values
  55 ##
  56 ## @itemize @bullet
  57 ## @item
  58 ## @var{p} is the cumulative distribution of the copula at each row of
  59 ## @var{x} and corresponding parameter @var{theta}.
  60 ## @end itemize
  61 ##
  62 ## @subheading Examples
  63 ##
  64 ## @example
  65 ## @group
  66 ## x = [0.2:0.2:0.6; 0.2:0.2:0.6];
  67 ## theta = [1; 2];
  68 ## p = copulacdf ("Clayton", x, theta)
  69 ## @end group
  70 ##
  71 ## @group
  72 ## x = [0.2:0.2:0.6; 0.2:0.1:0.4];
  73 ## theta = [0.2, 0.1, 0.1, 0.05];
  74 ## p = copulacdf ("FGM", x, theta)
  75 ## @end group
  76 ## @end example
  77 ##
  78 ## @subheading References
  79 ##
  80 ## @enumerate
  81 ## @item
  82 ## Roger B. Nelsen. @cite{An Introduction to Copulas}. Springer,
  83 ## New York, second edition, 2006.
  84 ## @end enumerate
  85 ## @end deftypefn
  86
  87 ## Author: Arno Onken <asnelt@asnelt.org>
  88 ## Description: CDF of a copula family
  89
  90 function p = copulacdf (family, x, theta, nu)
  91
  92   # Check arguments
  93   if (nargin != 3 && (nargin != 4 || ! strcmpi (family, "t")))
  94     print_usage ();
  95   endif
  96
  97   if (! ischar (family))
  98     error ("copulacdf: family must be one of 'Gaussian', 't', 'Clayton', 'Gumbel', 'Frank', 'AMH', and 'FGM'");
  99   endif
 100
 101   if (! isempty (x) && ! ismatrix (x))
 102     error ("copulacdf: x must be a numeric matrix");
 103   endif
 104
 105   [n, d] = size (x);
 106
 107   lower_family = lower (family);
 108
 109   # Check family and copula parameters
 110   switch (lower_family)
 111
 112     case {"gaussian", "t"}
 113       # Family with a covariance matrix
 114       if (d == 2 && isscalar (theta))
 115         # Expand a scalar to a correlation matrix
 116         theta = [1, theta; theta, 1];
 117       endif
 118       if (any (size (theta) != [d, d]) || any (diag (theta) != 1) || any (any (theta != theta')) || min (eig (theta)) <= 0)
 119         error ("copulacdf: theta must be a correlation matrix");
 120       endif
 121       if (nargin == 4)
 122         # Student's t family
 123         if (! isscalar (nu) && (! isvector (nu) || length (nu) != n))
 124           error ("copulacdf: nu must be a vector with the same number of rows as x or be scalar");
 125         endif
 126         nu = nu(:);
 127       endif
 128
 129     case {"clayton", "gumbel", "frank", "amh"}
 130       # Archimedian one parameter family
 131       if (! isvector (theta) || (! isscalar (theta) && length (theta) != n))
 132         error ("copulacdf: theta must be a vector with the same number of rows as x or be scalar");
 133       endif
 134       theta = theta(:);
 135       if (n > 1 && isscalar (theta))
 136         theta = repmat (theta, n, 1);
 137       endif
 138
 139     case {"fgm"}
 140       # Exponential number of parameters
 141       if (! ismatrix (theta) || size (theta, 2) != (2 .^ d - d - 1) || (size (theta, 1) != 1 && size (theta, 1) != n))
 142         error ("copulacdf: theta must be a row vector of length 2^d-d-1 or a matrix of size n x (2^d-d-1)");
 143       endif
 144       if (n > 1 && size (theta, 1) == 1)
 145         theta = repmat (theta, n, 1);
 146       endif
 147
 148     otherwise
 149       error ("copulacdf: unknown copula family '%s'", family);
 150
 151   endswitch
 152
 153   if (n == 0)
 154     # Input is empty
 155     p = zeros (0, 1);
 156   else
 157     # Truncate input to unit hypercube
 158     x(x < 0) = 0;
 159     x(x > 1) = 1;
 160
 161     # Compute the cumulative distribution function according to family
 162     switch (lower_family)
 163
 164       case {"gaussian"}
 165         # The Gaussian family
 166         p = mvncdf (norminv (x), zeros (1, d), theta);
 167         # No parameter bounds check
 168         k = [];
 169
 170       case {"t"}
 171         # The Student's t family
 172         p = mvtcdf (tinv (x, nu), theta, nu);
 173         # No parameter bounds check
 174         k = [];
 175
 176       case {"clayton"}
 177         # The Clayton family
 178         p = exp (-log (max (sum (x .^ (repmat (-theta, 1, d)), 2) - d + 1, 0)) ./ theta);
 179         # Product copula at columns where theta == 0
 180         k = find (theta == 0);
 181         if (any (k))
 182           p(k) = prod (x(k, :), 2);
 183         endif
 184         # Check bounds
 185         if (d > 2)
 186           k = find (! (theta >= 0) | ! (theta < inf));
 187         else
 188           k = find (! (theta >= -1) | ! (theta < inf));
 189         endif
 190
 191       case {"gumbel"}
 192         # The Gumbel-Hougaard family
 193         p = exp (-(sum ((-log (x)) .^ repmat (theta, 1, d), 2)) .^ (1 ./ theta));
 194         # Check bounds
 195         k = find (! (theta >= 1) | ! (theta < inf));
 196
 197       case {"frank"}
 198         # The Frank family
 199         p = -log (1 + (prod (expm1 (repmat (-theta, 1, d) .* x), 2)) ./ (expm1 (-theta) .^ (d - 1))) ./ theta;
 200         # Product copula at columns where theta == 0
 201         k = find (theta == 0);
 202         if (any (k))
 203           p(k) = prod (x(k, :), 2);
 204         endif
 205         # Check bounds
 206         if (d > 2)
 207           k = find (! (theta > 0) | ! (theta < inf));
 208         else
 209           k = find (! (theta > -inf) | ! (theta < inf));
 210         endif
 211
 212       case {"amh"}
 213         # The Ali-Mikhail-Haq family
 214         p = (theta - 1) ./ (theta - prod ((1 + repmat (theta, 1, d) .* (x - 1)) ./ x, 2));
 215         # Check bounds
 216         if (d > 2)
 217           k = find (! (theta >= 0) | ! (theta < 1));
 218         else
 219           k = find (! (theta >= -1) | ! (theta < 1));
 220         endif
 221
 222       case {"fgm"}
 223         # The Farlie-Gumbel-Morgenstern family
 224         # All binary combinations
 225         bcomb = logical (floor (mod (((0:(2 .^ d - 1))' * 2 .^ ((1 - d):0)), 2)));
 226         ecomb = ones (size (bcomb));
 227         ecomb(bcomb) = -1;
 228         # Summation over all combinations of order >= 2
 229         bcomb = bcomb(sum (bcomb, 2) >= 2, end:-1:1);
 230         # Linear constraints matrix
 231         ac = zeros (size (ecomb, 1), size (bcomb, 1));
 232         # Matrix to compute p
 233         ap = zeros (size (x, 1), size (bcomb, 1));
 234         for i = 1:size (bcomb, 1)
 235           ac(:, i) = -prod (ecomb(:, bcomb(i, :)), 2);
 236           ap(:, i) = prod (1 - x(:, bcomb(i, :)), 2);
 237         endfor
 238         p = prod (x, 2) .* (1 + sum (ap .* theta, 2));
 239         # Check linear constraints
 240         k = false (n, 1);
 241         for i = 1:n
 242           k(i) = any (ac * theta(i, :)' > 1);
 243         endfor
 244
 245     endswitch
 246
 247     # Out of bounds parameters
 248     if (any (k))
 249       p(k) = NaN;
 250     endif
 251
 252   endif
 253
 254 endfunction
 255
 256 %!test
 257 %! x = [0.2:0.2:0.6; 0.2:0.2:0.6];
 258 %! theta = [1; 2];
 259 %! p = copulacdf ("Clayton", x, theta);
 260 %! expected_p = [0.1395; 0.1767];
 261 %! assert (p, expected_p, 0.001);
 262
 263 %!test
 264 %! x = [0.2:0.2:0.6; 0.2:0.2:0.6];
 265 %! p = copulacdf ("Gumbel", x, 2);
 266 %! expected_p = [0.1464; 0.1464];
 267 %! assert (p, expected_p, 0.001);
 268
 269 %!test
 270 %! x = [0.2:0.2:0.6; 0.2:0.2:0.6];
 271 %! theta = [1; 2];
 272 %! p = copulacdf ("Frank", x, theta);
 273 %! expected_p = [0.0699; 0.0930];
 274 %! assert (p, expected_p, 0.001);
 275
 276 %!test
 277 %! x = [0.2:0.2:0.6; 0.2:0.2:0.6];
 278 %! theta = [0.3; 0.7];
 279 %! p = copulacdf ("AMH", x, theta);
 280 %! expected_p = [0.0629; 0.0959];
 281 %! assert (p, expected_p, 0.001);
 282
 283 %!test
 284 %! x = [0.2:0.2:0.6; 0.2:0.1:0.4];
 285 %! theta = [0.2, 0.1, 0.1, 0.05];
 286 %! p = copulacdf ("FGM", x, theta);
 287 %! expected_p = [0.0558; 0.0293];
 288 %! assert (p, expected_p, 0.001);