function [VC, CDF, VCdelta, VCgamma, VCvega] = digital_option_price(distr, S0, K, x, T, CP, sv1, sv2, sv3, rf, zeta, N, L)
% This function allows one to price digital options and resulting risk-neutral CDFs
% for a given known underlying distribution (see the function: moment_gen_function.m)
%
% The below is lightly adapted code from Nicola Fusari's Matlab Toolbox for
% Option Pricing (including comments). This code edits (with some
% notational hacks) his code to obtain the price of a digital option
% (should work for digital calls and puts, but I've tested it only for
% calls, so set CP=1) following Fang and Oosterlee (2008). The option pays
% K if S_T>K and 0 otherwise. It then transforms the resulting price to a
% risk-neutral CDF evaluated at the same strikes, and following Nicola's
% code, spits out Greeks for the digital option if desired.
%
%
% Nicola's code:
%
%
% The option price approximation is based on the paper by Fang & Oosterlee 
% (2008), SIAM J. Sci. Comput. Vol 31, No. 2, pp. 826-848.
%
% I've tried to keep the notation as close as possible to the one in the
% paper.
%
% OUTPUT:   OPTION PRICE
%           RN CDF
%           OPTION DELTA 
%           OPTION GAMMA
%           OPTION VEGA
% INPUT:    Distr ---> select the distribution of the underlying (see
%                      function: moment_gen_function.m) 
%           S0    ---> underlying initial value
%           K     ---> column vector of strikes
%           x     ---> vector that contains the parameters for the selected
%                      distribution (see function: moment_gen_function.m) 
%           T     ---> maturity (in YEARS)
%           CP    ---> Call/Put Flag  (call=1, put=-1)
%           sv1   ---> value of the first variance component (if needed)
%           sv2   ---> value of the first variance component (if needed)
%           rf    ---> risk-free rate used to discound the option price
%           d     ---> dividend yield
% REMARK: S0 and K must have the same length ---> COLUMN VECTORS
% Nicola Fusari (2011):
% You can distribute and use the code. Please cite the author if you use
% the code in in your research.

% Compute numerical derivatives of the MGF to computed MEAN and VARIANCE of
% the underlying distribution to decide WHERE to TRUNCATE the DENSITY
% approximation:

delta = 0.01; % 0.01

cf1 = log(moment_gen_function(distr, x, -1.*delta./2  , 0, 0, 0, sv1, sv2, sv3, T, rf, zeta));
cf0 = log(moment_gen_function(distr, x,            0  , 0, 0, 0, sv1, sv2, sv3, T, rf, zeta));
CF1 = log(moment_gen_function(distr, x,  1.*delta./2  , 0, 0, 0, sv1, sv2, sv3, T, rf, zeta));

%   http://mathworld.wolfram.com/FiniteDifference.html

d1 = (CF1 - cf1)./delta;
d2 = (CF1 + cf1 - 2.*cf0)./((delta/2)^2);  


a = real(d1 - L*sqrt(abs(d2)));
b = real(d1 + L*sqrt(abs(d2)));

if CP==1
    c = 0; d = b;
else
    c = a; d = 0;
end

kay = 0:N-1;

dt = T;

X = log(S0./K);
uu = (kay.*pi)./repmat(b-a,1,N);

% Psi_k
tmp3 = zeros(1,numel(uu));
tmp3(2:end) = (sin(uu(2:end).*repmat(d-a,1,N-1)) - sin(uu(2:end).*repmat(c-a,1,N-1))).*(1./uu(2:end));
tmp3(1) = d - c;
PSI = tmp3;

CHIPSIdiff = PSI;

% Characteristic function of log(S_T/S_t)
[cf,~, beta1, ~, ~]  = moment_gen_function(distr, x, 1i.*uu, zeros(1,length(uu)), zeros(1,length(uu)), zeros(1,length(uu)), sv1, sv2, sv3, dt, rf, zeta);

vk = repmat(CP,numel(K),N).*repmat(K,[1,numel(CHIPSIdiff)]).*repmat(CHIPSIdiff,[numel(K),1]).*repmat(2./(b-a),numel(K),N);
fk = real(repmat(cf,[numel(X),1]).*exp(repmat(1i,numel(K),N).*repmat(uu,[numel(X),1]).*(repmat(X,[1,size(cf,2)])-repmat(a,numel(K),N))));

fk(:,1) = 0.5.*fk(:,1);
VC = max(repmat(exp(-rf*dt),numel(K),1).*sum((fk.*vk),2),0.0);

CDF = 1-VC.*exp(rf*T).*(1./K);

% Delta
fk2 = real(repmat(cf,[numel(X),1]).*exp(repmat(1i,numel(K),N).*repmat(uu,[numel(X),1]).*(repmat(X,[1,size(cf,2)])-repmat(a,numel(K),N))).*repmat(1i,numel(K),N).*repmat(uu,[numel(X),1])); %for delta
fk2(:,1) = 0.5.*fk2(:,1);
VCdelta = repmat(exp(-rf*dt),numel(K),1).*sum((fk2.*vk),2).*(1./S0);

% Gamma
fk3 = real(repmat(cf,[numel(X),1]).*exp(repmat(1i,numel(K),N).*repmat(uu,[numel(X),1]).*(repmat(X,[1,size(cf,2)])-repmat(a,numel(K),N))).*(repmat(-1i,numel(K),N).*repmat(uu,[numel(X),1])+(1i.*repmat(uu,[numel(X),1])).^2 )); %for gamma
fk3(:,1) = 0.5.*fk3(:,1);
VCgamma = repmat(exp(-rf*dt),numel(K),1).*sum((fk3.*vk),2).*(1./(S0.^2));

% Vega
fk4 = real(repmat(cf,[numel(X),1]).*repmat(beta1(end,:),[numel(X),1]).*exp(repmat(1i,numel(K),N).*repmat(uu,[numel(X),1]).*(repmat(X,[1,size(cf,2)])-repmat(a,numel(K),N))));
fk4(:,1) = 0.5.*fk4(:,1);
VCvega = max(exp(-rf*dt).*sum((fk4.*vk),2),0).*2.*sqrt(sv1(1));
