varargout=panjer(F,a,b,dx,n,str)

function varargout=panjer(F,a,b,dx,n,str)
% PANJER - approximate compound distribution using panjer recursion
%
% [gl,gu,x]=panjer(F,a,b,dx,n) uses Panjer recursion to calculate bounds
% on the distribution of the compound random variable
%   S = \sum_{i=1}^N X_i
% where N is a random variable of the Panjer class (see below) 
% and the X_i are independent identically distributed non-negative random
% variables with cumulative distribution function F(x).
%
% [gl,x]=panjer(F,a,b,dx,n,'lo') returns only the "lower" distribution
% [gu,x]=panjer(F,a,b,dx,n,'up') returns only the "upper" distribution
% [gl,gu,x]=panjer(F,a,b,dx,n,'both') returns both
% [gl,gu,x,GL,GU]=panjer(F,a,b,dx,n,'both') returns the cumulative distributions
% [gl,x,GL]=panjer(F,a,b,dx,n,'lo') returns the cumulative distribution
% [gu,x,GU]=panjer(F,a,b,dx,n,'up') returns the cumulative distribution
%
% Inputs:
%   F is a function handle to a cumulative distribution
%     where F(x)=0 for x<0.
%   a and b are scalar constants satisfying a<=1 and a+b>=0.
%   dx is the distance between evaluation points
%   n is the number of points, i.e. x=dx*(0:n-1);
%   str (optional) is one of 'lo', 'up', 'both'
%
% Outputs:
%   gl is the probability mass function g(x) of the "lower" estimate of the 
%      distribution of S, which assumes that the probability distribution
%      of X between x(k-1) and x(k) is concentrated at x(k).
%      By "lower" we mean GL<=G(x) where G is the cumulative distribution.
%   gu is the probability mass function g(x) of the "upper" estimate of the 
%      distribution of S, which assumes that the probability distribution
%      of X between x(k-1) and x(k) is concentrated at x(k-1).
%      By "upper" we mean G(x)<=GU where G is the cumulative distribution.
%   x  is the set of points at which gl and gu are evaluated, x=dx*(0:n-1)
%   GL is the lower estimate of the cumulative distribution
%      GL = cumsum(gl)
%   GU is the upper estimate of the cumulative distribution
%      GU = cumsum(gu)
%
% The Panjer class of random variables are all counting variables that
% that satisfy P[N=k]=p_k, p_k=(a+b/k)*p_{k-1} with a+b>=0
%
%   Poisson: a=0, b=L, p_0=exp(-L)
%   Binomial: a=-p/(1-p), b=p*(n+1)/(1-p), p_0=(1-p)^n
%   Negative Binomial: a=1-p, b=(1-p)*(r-1), p_0=p^r.
%
%
% Example (plot distribution of compound Poisson/Lognormal):
%   n = 1e4;
%   xmax = 100;
%   lam = 10;
%   mu = 0;
%   sig = 1;
%   [gl,gu,xp,GL,GU]=panjer(@(x)logncdf(x,mu,sig),0,lam,xmax/n,n);
%   plot(xp,GL,'b-',xp,GU,'r-');
%
% References:
%   Panjer, Harry H. (1981) "Recursive evaluation of a family of compound
%    distributions.", ASTIN Bulletin (International Actuarial Association)
%    12 (1): 2226. http://www.casact.org/library/astin/vol12no1/22.pdf.  
%

% Author: Ben Petschel 17/3/2009
%
% Change history:
%  17/3/2009 - created
%  18/3/2009 - minor comment updates
%


% see which discretisation to calculate:
if nargin>5,
  switch str
    case 'lo'
      dolo = true;
      doup = false;
    case 'up'
      dolo = false;
      doup = true;
    case 'both'
      dolo = true;
      doup = true;
  end;
elseif nargin==5
  dolo = true;
  doup = true;
else
  error('panjer: requires 5 input arguments');
end;

% discretization points:
x=dx*(0:n-1);

% calculate p0
if a==0,
  % Poisson
  p0 = exp(-b);
elseif a<0,
  % binomial: p=a/(a-1), n=-b/a-1
  p0 = (1-a/(a-1))^(-b/a-1);
elseif a>0,
  % negative binomial: p=1-a, r=1+b/a
  p0 = (1-a)^(1+b/a);
end;

% calculate prob(X=0) and prob[x(i)<X<=x(i+1)]
F0=F(0);
dF=diff([F(x),1]);

% panjer recursion is:
%   g(x_j)=(1/(1-a*f(x_0)))*sum_{i=1}^j(a+b*i/j)*f(x_i)*g(x_{j-i})
%   g(x_0)= p_0/(1-a*f(x_0))^(1+b/a) if a~=0 or p0*exp(b*f(x_0)) if a=0
% however index starts at zero so use:
%   g(j+1)=(1/(1-a*f(1)))*sum_{i=1}^j(a+b*i/j)*f(i+1)*g(j+1-i)
% also rewrite as:
%   g(j+1:end)=g(j+1:end)+c*(a+b*(1:n-j)./(j:n-1)).*f(2:end-j+1)*g(j)
%   = g(j+1:end)+c*(a-b*(j-1)+b./(j:n-1)).*f(2:end-j+1)*g(j)
if dolo,
  % for lower bound on panjer recursion
  fl=[F0,dF(1:end-1)];
  gl=zeros(1,n);
  cl=1/(1-a*fl(1));
  if a==0,
    gl(1)=p0*exp(b*fl(1));
  else
    gl(1)=p0/(1-a*fl(1))^(1+b/a);
  end;
  for j=1:n-1,
    % vectorized formula (can speed up results for larger arrays):
    gl(j+1:end)=gl(j+1:end)+cl*(a+b-b*(j-1)./(j:n-1)).*fl(2:end-j+1)*gl(j);
    % % non-vectorized formula:
    % gl(j+1)=cl*sum( (a+b*(1:j)/j).*fl(2:j+1).*gl(j:-1:1) );
  end;
end; % if dolo,
if doup,
  % for upper bound on panjer recursion
  fu=dF;
  fu(1)=fu(1)+F0;
  gu=zeros(1,n);
  cu=1/(1-a*fu(1));
  if a==0,
    gu(1)=p0*exp(b*fu(1));
  else
    gu(1)=p0/(1-a*fu(1))^(1+b/a);
  end;
  for j=1:n-1,
    % vectorized formula (can speed up results for larger arrays):
    gu(j+1:end)=gu(j+1:end)+cu*(a+b-b*(j-1)./(j:n-1)).*fu(2:end-j+1)*gu(j);
    % % non-vectorized formula:
    % gu(j+1)=cu*sum( (a+b*(1:j)/j).*fu(2:j+1).*gu(j:-1:1) );
  end;
end; % if doup,

% assign outputs
if dolo && doup,
  % assign 
  varargout{1} = gl;
  if nargout > 1,
    varargout{2} = gu;
  end;
  if nargout > 2,
      varargout{3} = x;
  end;
  if nargout > 3,
    varargout{4} = cumsum(gl);
  end;
  if nargout > 4,
    varargout{5} = cumsum(gu);
  end;
elseif dolo && ~doup,
  varargout{1} = gl;
  if nargout > 1,
    varargout{2} = x;
  end;
  if nargout > 2,
    varargout{3} = cumsum(gl);
  end;
elseif doup && ~dolo,
  varargout{1} = gu;
  if nargout > 1,
    varargout{2} = x;
  end;
  if nargout > 2,
    varargout{3} = cumsum(gu);
  end;
end;