1 ## Copyright (C) 2012 Arno Onken
3 ## This program is free software: you can redistribute it and/or modify
4 ## it under the terms of the GNU General Public License as published by
5 ## the Free Software Foundation, either version 3 of the License, or
6 ## (at your option) any later version.
8 ## This program is distributed in the hope that it will be useful,
9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 ## GNU General Public License for more details.
13 ## You should have received a copy of the GNU General Public License
14 ## along with this program. If not, see <http://www.gnu.org/licenses/>.
17 ## @deftypefn {Function File} {@var{x} =} copularnd (@var{family}, @var{theta}, @var{n})
18 ## @deftypefnx {Function File} {} copularnd (@var{family}, @var{theta}, @var{n}, @var{d})
19 ## @deftypefnx {Function File} {} copularnd ('t', @var{theta}, @var{nu}, @var{n})
20 ## Generate random samples from a copula family.
22 ## @subheading Arguments
26 ## @var{family} is the copula family name. Currently, @var{family} can be
27 ## @code{'Gaussian'} for the Gaussian family, @code{'t'} for the Student's t
28 ## family, or @code{'Clayton'} for the Clayton family.
31 ## @var{theta} is the parameter of the copula. For the Gaussian and Student's t
32 ## copula, @var{theta} must be a correlation matrix. For bivariate copulas
33 ## @var{theta} can also be a correlation coefficient. For the Clayton family,
34 ## @var{theta} must be a vector with the same number of elements as samples to
35 ## be generated or be scalar.
38 ## @var{nu} is the degrees of freedom for the Student's t family. @var{nu} must
39 ## be a vector with the same number of elements as samples to be generated or
43 ## @var{n} is the number of rows of the matrix to be generated. @var{n} must be
44 ## a non-negative integer and corresponds to the number of samples to be
48 ## @var{d} is the number of columns of the matrix to be generated. @var{d} must
49 ## be a positive integer and corresponds to the dimension of the copula.
52 ## @subheading Return values
56 ## @var{x} is a matrix of random samples from the copula with @var{n} samples
57 ## of distribution dimension @var{d}.
60 ## @subheading Examples
65 ## x = copularnd ("Gaussian", theta);
71 ## x = copularnd ("t", theta, nu);
77 ## x = copularnd ("Clayton", theta, n);
81 ## @subheading References
85 ## Roger B. Nelsen. @cite{An Introduction to Copulas}. Springer, New York,
86 ## second edition, 2006.
90 ## Author: Arno Onken <asnelt@asnelt.org>
91 ## Description: Random samples from a copula family
93 function x = copularnd (family, theta, nu, n)
100 if (! ischar (family))
101 error ("copularnd: family must be one of 'Gaussian', 't', and 'Clayton'");
104 lower_family = lower (family);
106 # Check family and copula parameters
107 switch (lower_family)
111 if (isscalar (theta))
112 # Expand a scalar to a correlation matrix
113 theta = [1, theta; theta, 1];
115 if (! ismatrix (theta) || any (diag (theta) != 1) || any (any (theta != theta')) || min (eig (theta)) <= 0)
116 error ("copularnd: theta must be a correlation matrix");
120 if (! isscalar (d) || d != size (theta, 1))
121 error ("copularnd: d must correspond to dimension of theta");
130 if (! isscalar (n) || (n < 0) || round (n) != n)
131 error ("copularnd: n must be a non-negative integer");
140 if (isscalar (theta))
141 # Expand a scalar to a correlation matrix
142 theta = [1, theta; theta, 1];
144 if (! ismatrix (theta) || any (diag (theta) != 1) || any (any (theta != theta')) || min (eig (theta)) <= 0)
145 error ("copularnd: theta must be a correlation matrix");
147 if (! isscalar (nu) && (! isvector (nu) || length (nu) != n))
148 error ("copularnd: nu must be a vector with the same number of rows as x or be scalar");
154 if (! isscalar (n) || (n < 0) || round (n) != n)
155 error ("copularnd: n must be a non-negative integer");
160 # Archimedian one parameter family
162 # Default is bivariate
166 if (! isscalar (d) || (d < 2) || round (d) != d)
167 error ("copularnd: d must be an integer greater than 1");
171 # Default is one sample
175 if (! isscalar (n) || (n < 0) || round (n) != n)
176 error ("copularnd: n must be a non-negative integer");
179 if (! isvector (theta) || (! isscalar (theta) && size (theta, 1) != n))
180 error ("copularnd: theta must be a column vector with the number of rows equal to n or be scalar");
182 if (n > 1 && isscalar (theta))
183 theta = repmat (theta, n, 1);
187 error ("copularnd: unknown copula family '%s'", family);
196 # Draw random samples according to family
197 switch (lower_family)
200 # The Gaussian family
201 x = normcdf (mvnrnd (zeros (1, d), theta, n), 0, 1);
202 # No parameter bounds check
206 # The Student's t family
207 x = tcdf (mvtrnd (theta, nu, n), nu);
208 # No parameter bounds check
216 # Conditional distribution method for the bivariate case which also
217 # works for theta < 0
219 x(:, 2) = (1 + u(:, 1) .^ (-theta) .* (u(:, 2) .^ (-theta ./ (1 + theta)) - 1)) .^ (-1 ./ theta);
221 # Apply the algorithm by Marshall and Olkin:
222 # Frailty distribution for Clayton copula is gamma
223 y = randg (1 ./ theta, n, 1);
224 x = (1 - log (u) ./ repmat (y, 1, d)) .^ (-1 ./ repmat (theta, 1, d));
226 k = find (theta == 0);
228 # Produkt copula at columns k
231 # Continue argument check
233 k = find (! (theta >= -1) | ! (theta < inf));
235 k = find (! (theta >= 0) | ! (theta < inf));
240 # Out of bounds parameters
251 %! x = copularnd ("Gaussian", theta);
252 %! assert (size (x), [1, 2]);
253 %! assert (all ((x >= 0) & (x <= 1)));
258 %! x = copularnd ("t", theta, nu);
259 %! assert (size (x), [1, 2]);
260 %! assert (all ((x >= 0) & (x <= 1)));
264 %! x = copularnd ("Clayton", theta);
265 %! assert (size (x), [1, 2]);
266 %! assert (all ((x >= 0) & (x <= 1)));
271 %! x = copularnd ("Clayton", theta, n);
272 %! assert (size (x), [n, 2]);
273 %! assert (all ((x >= 0) & (x <= 1)));
279 %! x = copularnd ("Clayton", theta, n, d);
280 %! assert (size (x), [n, d]);
281 %! assert (all ((x >= 0) & (x <= 1)));