function y=gam_dens(Z,c,v,rho,eps);
% Y=GAM_DENS(Z,C,V,RHO,EPS);
% calculates the one-step density: Y(t) = f(Z(t+1)|Z_t) where Z follows a
% autoregressive gamma process characterized by three variables C, V
% and RHO. In short, if we let:
%   + c     = 0.5*sigma^2*dt
%   + v     = 2*kappa*theta/Sigma^2
%   + rho   = 1-kappa*dt
% then in the time limit when dt approaches 0, z will approach the CIR
% process:
%   dz = kappa*(theta - z)*dt + sigma*z^0.5 dBt
% 
% EPS is the error tolerance level (recommended to be 1e-5);
%
% Y is a COLUMN vector one element less than Z (regardless of whether Z is
% a column or row vector.
%
% For more detailed applications of this code, see Discrete-time AffineQ
% Term Structure Models with Generalized Market Prices of Risk by Anh Le,
% Ken Singleton and Qiang Dai (the Review of Financial Studies)
%
% Coded by Anh Le (2009)

if nargin<5;
    eps = 1e-5;
end;
if size(Z,1)==1; Z=Z';end;

T = length(Z);
y = zeros(T-1,1);

Z2 = Z(2:end);
Z1 = Z(1:end-1);

v = ones(T-1,1)*v;

lambda = rho*Z1/c ;   % compute rho*Z_t/c and call it lambda
x = Z2/c;               % compute Z_{t+1}/c and call it x
num = x.*lambda;        % compute (rho*Z_t/c).*(Z_{t+1}/c) and call it num

% The gamma density is an infinite series with a bell shape. As such
% integration is best performed around the peak of the bell. Here, we
% identify the location of the peaks.
peak = floor(-0.5*(v+1)+(num+0.25*(v-1).^2).^0.5)+1;
a_peak = poisspdf(peak,lambda).*gampdf(x,v+peak,1)/c; % this will be needed for the integration later

% Find a point (n2) to the right of the peak that has the property a_n2/(1 -
% a_{n2+1}/a_n2) < eps
count = 1;
n2 = peak;
big_err = [1:(T-1)]';
while ~isempty(big_err);
    count = count+1;
    n2(big_err) = n2(big_err) + 2^count;
    big_err = find((n2(big_err)<1e10)&(1/c)*poisspdf(n2(big_err),lambda(big_err)).*gampdf(x(big_err),v(big_err)+n2(big_err))./(1- num(big_err)./(n2(big_err)+1)./(n2(big_err)+v(big_err))) > eps);
end; 
n1 = floor((peak + n2 - 1)/2);
n2 = min(n2,1e10); n1 = min(n1,1e10);

exact = find(n2-n1>1);
while ~isempty(exact);
    n_mid = floor(0.5*(n1(exact)+n2(exact)));
    an_mid = (1/c)*poisspdf(n_mid,lambda(exact)).*gampdf(x(exact),v(exact)+n_mid)./(1- num(exact)./(n_mid+1)./(n_mid+v(exact)));
    n1(exact(an_mid>eps)) = n_mid(an_mid>eps);
    n2(exact(an_mid<=eps)) = n_mid(an_mid<=eps);
    exact = find(n2-n1>1);
end;

% Find a point (n1) to the left of the peak such that similar tolerance
% level is observed. Since it can be shown that a_n decays faster as we
% walk to the left of a_k than as we walk to the right of a_k, we will
% use a simple way of picking n1: n1 is the reflection point of n2 through
% k, cut off at 0.
n1 = max(0,2*peak-n2);

ub_N = 30000; % if there are more than ub_N, we will use numerical technique. Otherwise, we just sum up all the terms.

sm = find(n2-n1+1<ub_N);
if ~isempty(sm);
    % For a reasonably small number of terms: I break the summation into two parts, both
    % starting from the peak of the graph, but one goes to the left and one
    % goes to the right. Why? while the straight forward way is to do the
    % sum from n1 to n2. Remember that the multiplier a_{n+1}/a_n = num/((n+1)(n+v))
    % therefore a sum from a_n1 to a_n2 can be done by creating a vector of
    % multipliers from n1 to n2 and then use some form of sum(cumprod(.)).
    % The problem with this method though is numerical error since our
    % starting point an_1 is likely to be a very small value and could
    % involve errors which may be compounded through the sum(cmprod(.)) all
    % the way until a_n2. As such, it is better to have a good starting
    % point such as the peak of the series (ak) and then scale down as we
    % walk to the left and the right. Since we are scaling down, errors can
    % only be smaller, not bigger.
    for j = 1:length(sm);
        s = 0;
        jj = sm(j);
        b = [peak(jj):n2(jj)-1];
        s = sum(cumprod([a_peak(jj), num(jj)./((b+1).*(b+v(jj)))]));
        b = [peak(jj)-1:-1:n1(jj)];
        s = s + a_peak(jj)*sum(cumprod((b+1).*(b+v(jj))/num(jj)));
        if s == 0; y(jj) = -1e8; else y(jj) = y(jj) + log(s); end;
    end;
end;

bg = find(n2-n1+1>=ub_N);
if ~isempty(bg);
    r1 = num(bg)./(n1(bg)+1)./(n1(bg)+v(bg)); % value of a_{n+1}/a_n when n=n1;
    r2 = num(bg)./(n2(bg)+1)./(n2(bg)+v(bg)); % value of a_{n+1}/a_n when n=n2;

    M = 100; % do the integration at M+1 different points
    % the idea here is to break the sum into smaller parts and
    % approximate those parts by numerical integration, the
    % precision of which depends on whether the terms don't draw
    % out too much of a nonlinear curve. To make sure that we have
    % reasonable curvature, the break points are determined as
    % points when the slopes (a_{n+1}/a_n) change by (r2-r1)/M as
    % determined by matrix r below:
    r = cumsum([r1 repmat((r2-r1)/M,1,M)],2);
    int_num = repmat(num(bg),1,M+1);
    int_v = repmat(v(bg),1,M+1);
    int_lambda = repmat(lambda(bg),1,M+1);

    % the indices corresponding to the slope matrix r
    int_k = floor(-0.5*(int_v+1)+(int_num./r+0.25*(int_v-1).^2).^0.5)+1;

    % map the indices to the unit interval
    int_x = (int_k - repmat(int_k(:,1),1,M+1))./repmat(int_k(:,end) - int_k(:,1),1,M+1);

    % computing the corresponding y values
    int_y = (1/c)*poisspdf(int_k,int_lambda).*gampdf(repmat(x(bg),1,M+1),int_v+int_k,1);

    % computing the corresponding first-order derivatives
    int_r_plus = int_num./(int_k+1)./(int_k+int_v); % a_{k+1}/a_k
    int_r_minus = int_k.*(int_k+int_v-1)./int_num; % a_{k-1}/a_k
    int_ydash = 0.5*(int_r_plus - int_r_minus).*int_y;

    % numerical integration
    s = (int_k(:,end) - int_k(:,1)).*sum(0.5*(int_x(:,2:end) - int_x(:,1:end-1)).*(int_y(:,2:end) + int_y(:,1:end-1)) ...
        + (1/12)*(int_x(:,2:end) - int_x(:,1:end-1)).^2.*(int_ydash(:,2:end) - int_ydash(:,1:end-1)),2);

    % log of density
    y(bg) = y(bg) + log(s);
end;