function m = APDstat(alpha,lda);
% APDSTAT.M     Statistical moments of Asymmetric Power Distribution (APD) random variables 
%   M = APDSTAT(ALPHA,LAD) computes the first four moments of an APD random
%   variable X with parameters ALPHA, 0 < ALPHA < 1, and LDA, LDA > 0.
%
% Inputs:
%   ALPHA = Nx1-vector of shape parameter alpha, 0 < ALPHA(n) < 1, for 1 <= n <= N; 
%   LDA = Nx1-vector of shape parameter lda, LDA(n) > 0 for 1 <= n <= N;
%
% Outputs:
%   M = Nx4 matrix of the first 4 moments of X.
%
% Comments: 
%   a standard APD(alpha,lda) random variable X has probability density:
%   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).
%
% 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 < 2, 
    error('Requires two input arguments'); 
end

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

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

% Computes A
a = (alpha.^lda)./(alpha.^lda+(1-alpha).^lda);

% First compute the matrix of non-centered moments e. Initialize e to zero.
e = zeros(size(a,1),4);

%   Return NaN if the arguments a and lda are outside their limits.
e(alpha <= 0 | alpha >= 1 | lda<=0,:) = NaN;     

k = find(alpha > 0 & alpha < 1 & lda>0);
ak = a(k); ldak = lda(k); r4 = [1,2,3,4]; one4 = [1,1,1,1]; onek = ones(size(ak));
if any(k)
    ek = gamma((1./ldak)*(1+r4))./gamma((1./ldak)*one4).*...
        ((2*ak*one4).^(-1*(1./ldak)*(1+r4)) + (-1).^(onek*r4).*(2*(1-ak*one4)).^(-1*(1./ldak)*(1+r4)))./...
        ((2*ak*one4).^(-1*(1./ldak)*one4) + (2*(1-ak*one4)).^(-1*(1./ldak)*one4));
end

e(k,:) = ek;

% compute m
m = e;
m1k = ek(:,1); 
m2k = ek(:,2)-ek(:,1).^2; 
m3k = (ek(:,3)-3*ek(:,1).*ek(:,2)+2*ek(:,1).^3)./(m2k.^1.5); 
m4k = (ek(:,4)-4*ek(:,1).*ek(:,3)+6*(ek(:,1).^2).*ek(:,2)-3*ek(:,1).^4)./(m2k.^2);
m(k,:) = [m1k, m2k, m3k, m4k];