1 ## Copyright (C) 2008 Arno Onken <asnelt@asnelt.org>
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
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
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/>.
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.
21 ## @subheading Arguments
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.
33 ## @var{x} is the support where each row corresponds to an observation.
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
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.
54 ## @subheading Return values
58 ## @var{p} is the cumulative distribution of the copula at each row of
59 ## @var{x} and corresponding parameter @var{theta}.
62 ## @subheading Examples
66 ## x = [0.2:0.2:0.6; 0.2:0.2:0.6];
68 ## p = copulacdf ("Clayton", x, theta)
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)
78 ## @subheading References
82 ## Roger B. Nelsen. @cite{An Introduction to Copulas}. Springer,
83 ## New York, second edition, 2006.
87 ## Author: Arno Onken <asnelt@asnelt.org>
88 ## Description: CDF of a copula family
90 function p = copulacdf (family, x, theta, nu)
93 if (nargin != 3 && (nargin != 4 || ! strcmpi (family, "t")))
97 if (! ischar (family))
98 error ("copulacdf: family must be one of 'Gaussian', 't', 'Clayton', 'Gumbel', 'Frank', 'AMH', and 'FGM'");
101 if (! isempty (x) && ! ismatrix (x))
102 error ("copulacdf: x must be a numeric matrix");
107 lower_family = lower (family);
109 # Check family and copula parameters
110 switch (lower_family)
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];
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");
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");
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");
135 if (n > 1 && isscalar (theta))
136 theta = repmat (theta, n, 1);
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)");
144 if (n > 1 && size (theta, 1) == 1)
145 theta = repmat (theta, n, 1);
149 error ("copulacdf: unknown copula family '%s'", family);
157 # Truncate input to unit hypercube
161 # Compute the cumulative distribution function according to family
162 switch (lower_family)
165 # The Gaussian family
166 p = mvncdf (norminv (x), zeros (1, d), theta);
167 # No parameter bounds check
171 # The Student's t family
172 p = mvtcdf (tinv (x, nu), theta, nu);
173 # No parameter bounds check
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);
182 p(k) = prod (x(k, :), 2);
186 k = find (! (theta >= 0) | ! (theta < inf));
188 k = find (! (theta >= -1) | ! (theta < inf));
192 # The Gumbel-Hougaard family
193 p = exp (-(sum ((-log (x)) .^ repmat (theta, 1, d), 2)) .^ (1 ./ theta));
195 k = find (! (theta >= 1) | ! (theta < inf));
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);
203 p(k) = prod (x(k, :), 2);
207 k = find (! (theta > 0) | ! (theta < inf));
209 k = find (! (theta > -inf) | ! (theta < inf));
213 # The Ali-Mikhail-Haq family
214 p = (theta - 1) ./ (theta - prod ((1 + repmat (theta, 1, d) .* (x - 1)) ./ x, 2));
217 k = find (! (theta >= 0) | ! (theta < 1));
219 k = find (! (theta >= -1) | ! (theta < 1));
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));
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);
238 p = prod (x, 2) .* (1 + sum (ap .* theta, 2));
239 # Check linear constraints
242 k(i) = any (ac * theta(i, :)' > 1);
247 # Out of bounds parameters
257 %! x = [0.2:0.2:0.6; 0.2:0.2:0.6];
259 %! p = copulacdf ("Clayton", x, theta);
260 %! expected_p = [0.1395; 0.1767];
261 %! assert (p, expected_p, 0.001);
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);
270 %! x = [0.2:0.2:0.6; 0.2:0.2:0.6];
272 %! p = copulacdf ("Frank", x, theta);
273 %! expected_p = [0.0699; 0.0930];
274 %! assert (p, expected_p, 0.001);
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);
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);