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} =} copulapdf (@var{family}, @var{x}, @var{theta})
18 ## Compute the probability density function of a copula family.
20 ## @subheading Arguments
24 ## @var{family} is the copula family name. Currently, @var{family} can
25 ## be @code{'Clayton'} for the Clayton family, @code{'Gumbel'} for the
26 ## Gumbel-Hougaard family, @code{'Frank'} for the Frank family, or
27 ## @code{'AMH'} for the Ali-Mikhail-Haq family.
30 ## @var{x} is the support where each row corresponds to an observation.
33 ## @var{theta} is the parameter of the copula. The elements of
34 ## @var{theta} must be greater than or equal to @code{-1} for the
35 ## Clayton family, greater than or equal to @code{1} for the
36 ## Gumbel-Hougaard family, arbitrary for the Frank family, and greater
37 ## than or equal to @code{-1} and lower than @code{1} for the
38 ## Ali-Mikhail-Haq family. Moreover, @var{theta} must be non-negative
39 ## for dimensions greater than @code{2}. @var{theta} must be a column
40 ## vector with the same number of rows as @var{x} or be scalar.
43 ## @subheading Return values
47 ## @var{p} is the probability density of the copula at each row of
48 ## @var{x} and corresponding parameter @var{theta}.
51 ## @subheading Examples
55 ## x = [0.2:0.2:0.6; 0.2:0.2:0.6];
57 ## p = copulapdf ("Clayton", x, theta)
61 ## p = copulapdf ("Gumbel", x, 2)
65 ## @subheading References
69 ## Roger B. Nelsen. @cite{An Introduction to Copulas}. Springer,
70 ## New York, second edition, 2006.
74 ## Author: Arno Onken <asnelt@asnelt.org>
75 ## Description: PDF of a copula family
77 function p = copulapdf (family, x, theta)
84 if (! ischar (family))
85 error ("copulapdf: family must be one of 'Clayton', 'Gumbel', 'Frank', and 'AMH'");
88 if (! isempty (x) && ! ismatrix (x))
89 error ("copulapdf: x must be a numeric matrix");
94 if (! isvector (theta) || (! isscalar (theta) && size (theta, 1) != n))
95 error ("copulapdf: theta must be a column vector with the same number of rows as x or be scalar");
102 if (n > 1 && isscalar (theta))
103 theta = repmat (theta, n, 1);
106 # Truncate input to unit hypercube
110 # Compute the cumulative distribution function according to family
111 lowerarg = lower (family);
113 if (strcmp (lowerarg, "clayton"))
115 log_cdf = -log (max (sum (x .^ (repmat (-theta, 1, d)), 2) - d + 1, 0)) ./ theta;
116 p = prod (repmat (theta, 1, d) .* repmat (0:(d - 1), n, 1) + 1, 2) .* exp ((1 + theta .* d) .* log_cdf - (theta + 1) .* sum (log (x), 2));
117 # Product copula at columns where theta == 0
118 k = find (theta == 0);
124 k = find (! (theta >= 0) | ! (theta < inf));
126 k = find (! (theta >= -1) | ! (theta < inf));
128 elseif (strcmp (lowerarg, "gumbel"))
129 # The Gumbel-Hougaard family
130 g = sum ((-log (x)) .^ repmat (theta, 1, d), 2);
131 c = exp (-g .^ (1 ./ theta));
132 p = ((prod (-log (x), 2)) .^ (theta - 1)) ./ prod (x, 2) .* c .* (g .^ (2 ./ theta - 2) + (theta - 1) .* g .^ (1 ./ theta - 2));
134 k = find (! (theta >= 1) | ! (theta < inf));
135 elseif (strcmp (lowerarg, "frank"))
138 error ("copulapdf: Frank copula PDF implemented as bivariate only");
140 p = (theta .* exp (theta .* (1 + sum (x, 2))) .* (exp (theta) - 1))./ (exp (theta) - exp (theta + theta .* x(:, 1)) + exp (theta .* sum (x, 2)) - exp (theta + theta .* x(:, 2))) .^ 2;
141 # Product copula at columns where theta == 0
142 k = find (theta == 0);
147 k = find (! (theta > -inf) | ! (theta < inf));
148 elseif (strcmp (lowerarg, "amh"))
149 # The Ali-Mikhail-Haq family
151 error ("copulapdf: Ali-Mikhail-Haq copula PDF implemented as bivariate only");
153 z = theta .* prod (x - 1, 2) - 1;
154 p = (theta .* (1 - sum (x, 2) - prod (x, 2) - z) - 1) ./ (z .^ 3);
156 k = find (! (theta >= -1) | ! (theta < 1));
158 error ("copulapdf: unknown copula family '%s'", family);
170 %! x = [0.2:0.2:0.6; 0.2:0.2:0.6];
172 %! p = copulapdf ("Clayton", x, theta);
173 %! expected_p = [0.9872; 0.7295];
174 %! assert (p, expected_p, 0.001);
177 %! x = [0.2:0.2:0.6; 0.2:0.2:0.6];
178 %! p = copulapdf ("Gumbel", x, 2);
179 %! expected_p = [0.9468; 0.9468];
180 %! assert (p, expected_p, 0.001);
183 %! x = [0.2, 0.6; 0.2, 0.6];
185 %! p = copulapdf ("Frank", x, theta);
186 %! expected_p = [0.9378; 0.8678];
187 %! assert (p, expected_p, 0.001);
190 %! x = [0.2, 0.6; 0.2, 0.6];
191 %! theta = [0.3; 0.7];
192 %! p = copulapdf ("AMH", x, theta);
193 %! expected_p = [0.9540; 0.8577];
194 %! assert (p, expected_p, 0.001);