function u = copularnd(type,varargin)
%COPULARND Random values from a copula.
%   U = COPULARND('Gaussian',RHO,N) returns a matrix U of random values
%   generated from a Gaussian copula with correlation parameter RHO.  Each
%   column of U comes from a Uniform(0,1) marginal distribution.  If RHO is
%   a scalar correlation coefficient, U is an N-by-2 matrix generated from
%   a bivariate Gaussian copula.  If RHO is a P-by-P correlation matrix, U
%   is an N-by-P matrix generated from a P-variate Gaussian copula.
%
%   U = COPULARND('t',RHO,NU,N) returns a matrix U of random values
%   generated from a t copula with correlation parameter RHO and degrees of
%   freedom NU.  Each column of U comes from a Uniform(0,1) marginal
%   distribution. If RHO is a scalar correlation coefficient, U is an
%   N-by-2 matrix generated from a bivariate t copula.  If RHO is a P-by-P
%   correlation matrix, U is an N-by-P matrix generated from a P-variate t
%   copula.
%
%   U = COPULARND(TYPE,ALPHA,N) returns an N-by-2 matrix of random values
%   generated from the bivariate Archimedean copula with parameter ALPHA.
%   TYPE is one of 'Clayton', 'Frank', or 'Gumbel'.  Each column of U comes
%   from a Uniform(0,1) marginal distribution.
%
%   Example:
%      % Determine the linear correlation coefficient corresponding to a
%      % bivariate Gaussian copula having a rank correlation of -0.5
%      tau = -0.5
%      rho = copulaparam('gaussian',tau)
%
%      % Generate dependent beta random values using that copula
%      u = copularnd('gaussian',rho,100)
%      b = betainv(u,2,2)
%
%      % Verify that those pairs have a sample rank correlation approximately
%      % equal to tau
%      tau_sample = kendall(b)

%   Written by Peter Perkins, The MathWorks, Inc.
%   Revision: 1.0  Date: 2003/09/05
%   This function is not supported by The MathWorks, Inc.
%
%   Requires MATLAB R13, including the Statistics Toolbox.


switch lower(type)

% Elliptical copulas
%
% Random vectors from these copulas can be generated by creating random
% vectors from the multivariate normal (or t) distribution, then
% transforming them to uniform marginals using the normal (or t) CDF.
case 'gaussian'
    if nargin < 3
        error('Requires three input arguments for the Gaussian copula.');
    end
    rho = varargin{1};
    n = varargin{2};
    if numel(rho) == 1
        if (rho < -1 | 1 < rho)
            error('RHO must be a correlation coefficient between -1 and 1, or a positive semidefinite correlation matrix.');
        end
        u = normcdf(mvnrnd([0 0],[1 rho; rho 1],n));
    else
        if ~iscor(rho)
            error('RHO must be a correlation coefficient between -1 and 1, or a positive semidefinite correlation matrix.');
        end
        p = size(rho,1);
        u = normcdf(mvnrnd(zeros(1,p),rho,n));
    end

case 't'
    if nargin < 4
        error('Requires four input arguments for the t copula.');
    end
    rho = varargin{1};
    nu = varargin{2};
    n = varargin{3};
    if ~(0 < nu)
        error('NU must be positive for the t copula.');
    end
    if numel(rho) == 1
        if (rho < -1 | 1 < rho)
            error('RHO must be a correlation coefficient between -1 and 1, or a positive semidefinite correlation matrix.');
        end
        u = tcdf(mvtrnd([1 rho; rho 1],nu,n),nu);
    else
        if ~iscor(rho)
            error('RHO must be a correlation coefficient between -1 and 1, or a positive semidefinite correlation matrix.');
        end
        u = tcdf(mvtrnd(rho,nu,n),nu);
    end

% Archimedean copulas
%
% Random pairs from these copulae can be generated sequentially: first
% generate u1 as a uniform r.v.  Then generate u2 from the conditional
% distribution F(u2 | u1; alpha) by generating uniform random values, then
% inverting the conditional CDF.
case {'clayton' 'frank' 'gumbel'}
    if nargin < 3
        error('Requires three input arguments for an Archimedean copula.');
    end
    alpha = varargin{1};
    if numel(alpha) ~= 1
        error('ALPHA must be a scalar.');
    end
    
    switch lower(type)
    case 'clayton'
        if alpha < 0
            error('ALPHA must be nonnegative for the Clayton copula.');
        end
        n = varargin{2};
        u1 = rand(n,1);
        % The inverse conditional CDF has a closed form for this copula.
        p = rand(n,1);
        if alpha > sqrt(eps)
            u2 = u1.*(p.^(-alpha./(1+alpha)) - 1 + u1.^alpha).^(-1./alpha);
        else
            u2 = p;
        end
        u = [u1 u2];

    case 'frank'
        n = varargin{2};
        u1 = rand(n,1);
        % The inverse conditional CDF has a closed form for this copula.
        p = rand(n,1);
        if abs(alpha) > log(realmax)
            u2 = (u1 < 0) + sign(alpha).*u1; % u1 or 1-u1
        elseif abs(alpha) > sqrt(eps)
            u2 = -log((exp(-alpha.*u1).*(1-p)./p + exp(-alpha))./(1 + exp(-alpha.*u1).*(1-p)./p))./alpha;
%             u2 = -log(1 + (1-exp(alpha))./(exp(alpha) + exp(alpha.*(1-u1)).*(1-p)./p))./alpha;
        else
            u2 = p;
        end
        u = [u1 u2];

    case 'gumbel'
        if alpha < 1
            error('ALPHA must be greater than or equal to 1 for the Gumbel copula.');
        end
        n = varargin{2};
        u1 = rand(n,1);
        % The inverse conditional CDF does not have a closed form for this
        % copula.  The inversion must be done numerically.
        p = rand(n,1);
        u2 = condCDFinv(@gumbelCondCDF,u1,p,alpha);
        u = [u1 u2];
    end

otherwise
    error('Unrecognized copula type: ''%s''',type);
end


function p = gumbelCondCDF(u1,u2,alpha)
%GUMBELCONDCDF Conditional distribution function for the Gumbel copula.
%   P = gumbelCondCDF(U,ALPHA) returns the conditional cumulative distibution
%   function of U2, given U1, for the Gumbel copula.
% logC = -((-log(u1)).^alpha + (-log(u2)).^alpha).^(1./alpha);
nlog1 = -log(u1);
nlog2 = -log(u2);
maxlog = max(nlog1,nlog2);
logC = -maxlog .* ((nlog1./maxlog).^alpha + (nlog2./maxlog).^alpha).^(1./alpha);
p = (-nlog1./logC).^(alpha-1) .* exp(logC)./u1;


function u2 = condCDFinv(condCDF,u1,p,alpha)
%CONDCDFINV Inverse conditional distribution function
%   U2 = CONDCDFINV(CONDCDF,U1,P,ALPHA) returns U2 such that
%
%      CONDCDF(U1,U2,ALPHA) = P,
%
%  where CONDCDF is a function handle to a function that computes the
%  conditional cumulative distribution function of U2 given U1, for an
%  archimedean copula with parameter ALPHA.
%
% CONDCDFINV uses a simple binary chop search.  Newton's method or the
% secant method would probably be faster.

lower = zeros(size(p));
upper = ones(size(p));
width = 1;
tol = 1e-12;
while width > tol
    u2 = .5 .* (lower + upper);
    lo = feval(condCDF,u1,u2,alpha) < p;
    lower(lo) = u2(lo);
    upper(~lo) = u2(~lo);
    width = .5 .* width;
end