function r = APDrnd(alpha,lda,n);
% APDRND.M   Random matrices from a standard Asymmetric Power Distribution (APD).
%   R = APDRND(ALPHA,LDA,N) returns a matrix of random numbers chosen from 
%   a standard APD random variable with parameters ALPHA, 0 < ALPHA < 1, and LDA, LDA > 0.
%   If the parameters ALPHA and LDA are scalar, then R is an Nx1 vector. Alternatively,
%   if ALPHA and LAD are vectors of size Mx1 then R = APDRND(ALPHA,LDA,N)
%   returns an N by M matrix.
%
% Inputs:
%   ALPHA = Mx1-vector of shape parameter alpha, 0 < ALPHA(m) < 1 for 1 <= m <= M;
%   LDA = Mx1-vector of shape parameter lda, LDA(m) > 0 for 1 <= m <= M;
%   N = number of rows of the random vector/matrix;
%   
% Outputs:
%   R = NxM-matrix of standard APD random variables;
%
% Comments: 
%   a standard (scalar) APD random variable X with shape parameters alpha, 0 < alpha < 1,
%   and lda, lda > 0, has a pdf :
%   p(alpha,lda,x) = [A^(1/lda)/gamma(1+1/lda)]*...
%                   *exp[-(A/alpha^lda)*abs(x)^lda],        if x <= 0,
%                   *exp[-(A/(1-alpha)^lda)*abs(x)^lda],    if x > 0,
%   where A = (2*alpha^lda*(1-alpha)^lda)/(alpha^lda+(1-alpha)^lda).
%
%   APDRND uses GAMRND contained in Matlab Statistics Toolbox.
%
% References:
%   [1] Komunjer, I. (2006): "Asymmetric Power Distribution: Theory and 
%   Applications to Risk Measurement"
%
% Author: Ivana Komunjer - UCSD (komunjer@ucsd.edu)
% Version: 2.1        Date: 06/07/2006

if nargin < 3, 
    error('Requires three input arguments'); 
end

[errorcode alpha lda] = distchck(2,alpha,lda);

if errorcode > 0
    error('Requires non-scalar arguments to match in size.');
end

rows = n; 
columns = size(alpha,1);

% Initialize r to zero.
r = zeros(rows,columns); onen = ones(n,1);

%   Return NaN if the arguments ALPHA or LDA are outside their limit.
kn = find(alpha <= 0 | alpha >= 1 | lda <= 0);
r(:,kn) = NaN;     

k = find(alpha > 0 & alpha < 1 & lda > 0);
alphak = alpha(k); ldak = lda(k); onek = ones(size(alphak));
% Computes A
ak = 2*(alphak.^ldak).*((1-alphak).^ldak)./(alphak.^ldak + (1-alphak).^ldak);

% Step 1 : generates an NxM gamma random variable Z with shape parameter (1/lda,1)
zk = gamrnd(onen*(1./ldak)',onen*(onek)');
% Step 2: divides Z by A and raise to 1/lda power. This is W.
wk = (zk./(onen*ak')).^(onen*(1./ldak)');
% Step 3: generates a binomial random variable +1 (proba 1-alpha) and -1 (proba alpha).
rndsgn = 2*binornd(1,onen*(1 - alphak)') - 1;  
% Step 4: if rndsgn==-1 then R = -alpha*W and if rndsgn==+1 then R = (1-alpha)*W
rk = (rndsgn == -1).*((-1)*onen*alphak').*wk + (rndsgn == 1).*(onen*(1-alphak)').*wk;

r(:,k) = rk;