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(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;