PDF Calculation by Fourier Inversion of Characteristic Function for Affine Intensity Process in Matlab

I'm trying to use the Fourier inversion formula to plot the PDF of an Affine Stochastic Intensity Reduced Form Credit Model, given its characteristic function.

The characteristic function of an affine process $\lambda(t)$ is commonly given as

$$\phi_{\lambda(t)}(u) = \mathrm{E}[e^{iu\lambda(t)}] = \exp(A(t-s,iu)+B(t-s,iu)\lambda(s))$$

The Fourier inversion formula for PDF is

$$f_{\lambda(t)}(x)=\frac{1}{\pi}\int_0^\infty \mathrm{\Re}[e^{-iux}\phi_{\lambda(t)}(u)]du$$

Taking a CIR process (I’m aware that CIR has a $\chi^2$ Closed-Form PDF and the use of CIR here is just for illustration) which has coefficients:

$$A(T)=\frac{2\kappa\theta}{\sigma^2}\log\left(\frac{2\gamma e^{\frac{1}{2}(\kappa+\gamma)T} }{(\kappa+\gamma)(e^{\gamma T}-1)+2\gamma}\right)$$

$$B(T)=\frac{2 (e^{\gamma T}-1) }{(\kappa+\gamma)(e^{\gamma T}-1)+2\gamma}$$

In matlab script then, using quadrature for the integral, I (try to) calculate the PDF at the $\lambda$-points X = (0:0.005:0.1) for T=1 with the code below.

Clearly there is a problem though (quite probably with fcnPhi below ) - Would greatly appreciate any help here

kappa = .07;
theta = .2;
sigma = .06;
lambda0 = .06;
T = 1;
gamma = 1;

A = ((2*kappa*theta)/(sigma^2))* log(2*gamma*exp(0.5*(kappa+gamma)*(T))./((kappa+gamma)*(exp(gamma*(T))-1)+2*gamma));
B = 2*(exp(gamma*(T))-1)/((kappa+gamma)*(exp(gamma*(T))-1)+2*gamma);

fcnPhi = @(u)( exp(u.*(A + B*lambda0)) );

X = (0:0.005:0.1)';
for i = 1:size(X,1)
x = X(i);
fcnPdfIntgrl = @(u)( real( exp(-1i.*u.*x) .* fcnPhi(u) ) );
pdf_X(i,1) = (1/pi) * integral(fcnPdfIntgrl,0,10000);
end
plot(X,pdf_X);

• Have you tried the matlab-function ifft? It seems to make more sense than doing the quadrature explicitly. Depending on the implementation, you might have to normalize your vector afterwards though. – vanguard2k Jun 11 '14 at 6:56
• Hi Vanguard, Ultimately I will use ifft but I stuck with the quadrature here to avoid an unnecessary complication to the issue at hand. There is a problem with my interpretation of how to use the char function fcnPhi I think - and so I will have the same problem with the ifft. I really need to solve this quadrature version before proceeding I think as ifft doesn't solve the problem at hand (but obviously ultimately is a more efficient computational implementation) – StudentT Jun 11 '14 at 11:38
• What exactly is the problem? Are you getting results you don't expect? Or, is it an error? Post either the error, or the unexpected results (and describe what the expected results would look like) and I can probably help you. – Fletch Jun 13 '14 at 20:59
• well, the output should be a pdf, if you can run the code you'll see its not. But actually, I've just solved the problem - it is indeed the characteristic function. Although there is another issue with matlab's integral function not being able to handle the small tolerance required. I'll post the solution when I get around to it – StudentT Jun 13 '14 at 21:06

Need to solve Riccati Equations with Complex Boundary Conditions for Char Function, no general closed form solution (apart from the "Basic" Affine model, CIR etc) a la Duffie http://www.mit.edu/~junpan/dps.pdf

Solution for CIR below

kappa = 2;
theta = .05;
sigma = .02;
lambda0 = .03;
T = 1;

X = (0.03:0.001:0.07)';

c = 2*kappa/(sigma^2*(1 - exp(-kappa*T))); % Scaling Factor
q = 2*theta*kappa/sigma^2 - 1; % Deg-of-Freedom
u = c*exp(-kappa*T)*lambda0;
v = c*X;

pdf_1 = 2*c*ncx2pdf(2*v,2*q+2,2*u);

fcnA = @(u)( -(2*kappa*theta/(sigma^2)) * log( 1 - ((sigma^2)/(2*kappa)) .* 1i.*u.*(1-exp(-kappa*T)) )  );
fcnB = @(u)( (1i.*u.*exp(-kappa*T)) ./ (1 - ((sigma^2)/(2*kappa)).*1i.*u.*(1-exp(-kappa*T)) ) );
fcnPhi = @(u)( exp( fcnA(u) + fcnB(u)*lambda0 ) );

pdf_2 = NaN(size(X));
for i = 1:size(X,1)
x = X(i);
fcnPdfIntgrl = @(u)( real( exp(-1i.*u.*x) .* fcnPhi(u) ) );
pdf_2(i,1) = (1/pi) * integral(fcnPdfIntgrl,0,Inf,'RelTol',1e-10,'AbsTol',1e-10);
end

subplot(6,1,1:2);
plot(X,pdf_1);
title('Probability Density Function by \chi^2');
subplot(6,1,3:4);
plot(X,pdf_2);
title('Probability Density Function by Fourier Inversion of CF');

U = (0:10:4000);
subplot(6,1,5);
plot(U,real(fcnPhi(U)));
title('Characteristic Function');
subplot(6,1,6);
plot(U,fcnPdfIntgrl(U));
title('Fourier Inversion Integrand');