function f = APDes(alpha_bar,alpha,lda);
% APDES   Alpha_bar-expected shortfall from a standard APD random 
%   variable with parameters alpha and lambda.
%   X = APDES(ALPHA_BAR,ALPHA,LDA) returns the ALPHA_BAR-expected shorfall of a standard 
%   APD random variable with shape parameters ALPHA, 0 < ALPHA < 1, and LDA, LDA > 0. 
%
% Inputs:
%   ALPHA_BAR = Nx1 vector of the probabilities alpha_bar, 0 < ALPHA_BAR(n) < 1, for 1 <= n <= N;
%   ALPHA = Nx1 vector of the shape parameter alpha, 0 < ALPHA(n) < 1, for 1 <= n <= N;
%   LDA = Nx1 vector the shape parameter lda, LDA(n) > 0, for 1 <= n <= N;
%
% Outputs:
%   Y = Nx1 vector of ALPHA_BAR-expected shortfalls corresponding to the specified shape values ALPHA and LDA.
%
% 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).
%
%   APDES uses GAMINV 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_bar alpha lda] = distchck(3,alpha_bar,alpha,lda);

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

%   Initialize Y to zero.
y = zeros(size(alpha));
a = 2*(alpha.^lda).*((1-alpha).^lda)./(alpha.^lda+(1-alpha).^lda);

k = find(alpha_bar <= 0 | alpha_bar >= 1 | alpha >= 1 | alpha <= 0 | lda <= 0);
if any(k),
    tmp  = NaN;
    y(k) = tmp(ones(size(k)));
end

k = find(alpha_bar > 0  &  alpha_bar < 1 & alpha > 0 & alpha < 1 & lda > 0);
if (~any(k(:))), return; end
alpha_bark = alpha_bar(k); alphak = alpha(k); ldak = lda(k);

% Computes the alpha_bar-ES when alpha_bar < alpha. Uses GAMINV.
k1 = find(alpha_bar < alpha & alpha_bar > 0 & alpha < 1 & lda > 0);
if any(k1(:))
    alpha_bark1 = alpha_bar(k1); alphak1 = alpha(k1); ldak1 = lda(k1); ak1 = a(k1);
    yk1 = gaminv(1-alpha_bark1./alphak1,1./ldak1,1);  
    yk1 = 1 - gammainc(yk1,2./ldak1);
    yk1 = (alphak1./alpha_bark1).*(alphak1./(ak1.^(1./ldak1))).*(gamma(2./ldak1)./gamma(1./ldak1)).*yk1;
    yk1 = yk1 + AlphaQuantInv(alphak1,ldak1,alpha_bark1);
else
    yk1 = [];
end

% Computes the alpha_bar-ES when alpha_bar = alpha. Uses GAMINV.
k3 = find(alpha_bar == alpha & alpha > 0 & alpha < 1 & lda > 0);
if any(k3(:))
    alpha_bark3 = alpha_bar(k3); alphak3 = alpha(k3); ldak3 = lda(k3); ak3 = a(k3);
    yk3 = (alphak3./(ak3.^(1./ldak3))).*(gamma(2./ldak3)./gamma(1./ldak3));
else
    yk3 = [];
end

% Computes the alpha_bar-ES when alpha_bar > alpha. Uses GAMINV.
k2 = find(alpha_bar > alpha &  alpha_bar < 1 & alpha > 0 & lda > 0);
if any(k2(:))
    alpha_bark2 = alpha_bar(k2); alphak2 = alpha(k2); ldak2 = lda(k2); ak2 = a(k2);
    yk2 = gaminv(1-(1-alpha_bark2)./(1-alphak2),1./ldak2,1);
    yk2 = gammainc(yk2,2./ldak2);
    yk2 = ((1-alphak2)./alpha_bark2).*((1-alphak2)./(ak2.^(1./ldak2))).*(gamma(2./ldak2)./gamma(1./ldak2)).*yk2;
    yk2 = (alphak2./alpha_bark2).*(alphak2./(ak2.^(1./ldak2))).*(gamma(2./ldak2)./gamma(1./ldak2)) - yk2;
    yk2 = yk2 + AlphaQuantInv(alphak2,ldak2,alpha_bark2);
else
    yk2 = [];
end

yk = [yk1; yk3; yk2];

% Store the values in the correct place
y(k) = yk;