Matlab中仿射強度過程特徵函式傅里葉反演的PDF計算

鑑於其特徵函式，我正在嘗試使用傅立葉反演公式來繪製仿射隨機強度簡化形式信用模型的 PDF。

仿射過程的特徵函式 $ \lambda(t) $ 通常給出為

$$ \phi_{\lambda(t)}(u) = \mathrm{E}[e^{iu\lambda(t)}] = \exp(A(t-s,iu)+B(t-s,iu)\lambda(s)) $$ PDF的傅里葉反演公式是

$$ f_{\lambda(t)}(x)=\frac{1}{\pi}\int_0^\infty \mathrm{\Re}[e^{-iux}\phi_{\lambda(t)}(u)]du $$ 採用 CIR 流程（我知道 CIR 有一個 $ \chi^2 $ 封閉式 PDF 和此處使用的 CIR 僅用於說明），其係數：

$$ 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} $$ 然後在matlab腳本中，使用積分求積，我（嘗試）計算PDF $ \lambda $ -pointsX = (0:0.005:0.1)使用T=1下面的程式碼。

顯然有一個問題（很可能在fcnPhi下面） - 非常感謝這裡的任何幫助

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

需要求解具有復雜邊界條件的 Char 函式的 Riccati 方程，沒有一般封閉形式的解（除了“基本”仿射模型、CIR 等）a la Duffie http://www.mit.edu/~junpan/dps.pdf

以下 CIR 的解決方案

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');

引用自：https://quant.stackexchange.com/questions/12614