function xxx = lect01FFToptionnormal
clear;
format compact;
% parameters for the diffusion model
mu    = .0;
sigma = .3;
t = .25;
% parameters
a = 600;            % endpoint of char fun grid (0,+a)
N =  4096;          % number of grid points
eta = a/N;          % char fun grid size
b = pi/eta;         % +/- bounds for the log-strike grid (-b,+b)
lambda = 2*pi/a;    % log-strike grid size
sprintf('Parameter eta    (  char fun grid size): %2.4f',eta)
sprintf('Parameter lambda (log-strike grid size): %2.4f',lambda)
% dumping parameter as in Carr-Madan
aa = 1.25;
% create grids
u = (0:N-1) * eta;
x = -b + (0:N-1) * lambda;
% create char fun and invert
h = cfn(u, mu, sigma, t, aa);
h2 = exp(i*b*u) .* h * eta;
g   = fft(h2);
%plot(u,real(h2),u,imag(h2))
g2  = real( g .* exp(-aa*x) / pi);  % prices based on char fun
gbs = bs(exp(x), mu, sigma, t);     % Black-Scholes prices
% plot them together
plot(exp(x),g2,exp(x),1-exp(x), exp(x), gbs )
axis([.7 1.3 0 .3])

function y = cfn(th, mu, sig, t, aa)
% The characteristic function of the log-price
% Adjusted to incorporate the Simpson weighting scheme
th1 = th - (aa+1)*i;
f1 = exp(-mu*t) * exp( (mu-.5*sig^2)*t*i*th1 - .5*(sig^2)*t*(th1.^2) );
f2 = aa^2 + aa - (th.^2) + i*(2*aa+1)*th;
y  = f1 ./ f2;
% Create Simpson's weights
N  = size(th,2);
q1 = (-1).^(1:N);
q2 = eye(1,N);
S  = ( 3 + q1 - q2 )/3;
y  = y .* S;

function y = bs(x, r, sig, t)
% The simple Black-Scholes European call option Price
d1 = ( -log(x)+( r + .5*sig^2 )*t ) / sig/sqrt(t);
d2 = d1 - sig*sqrt(t);
n1 = normcdf(d1);
n2 = normcdf(d2);
y = n1 - x.*n2*exp(-r*t);