% Computes the call option price and risk-neutral CDF using numerical integration and Fourier inversion
% for the Christoffersen, Heston, and Jacobs (2013) model
% Based on code from Peter Christoffersen, Kris Jacobs, and Yoontae Jeon (GARCH Options Toolbox)
% S_t = d vector of spot stock prices
% K = d vector of Strike Prices
% tau = d vector of time to maturity
% h_t = d vector of variance processes
% compute_prices: set to 1 if want option prices, 0 if just rn probs

function [rn,c] = RN_CHJ(S_t,K,h_t,rF,tau,param,compute_prices)

if compute_prices==1
    t_Hk = @(u)imag(CF_CHJ(u-1i,S_t,h_t,rF,tau,param).*exp(-1i.*u.*log(K))./(1i*u+1))./u;
end
t_Hk2 = @(u)real(CF_CHJ(u,S_t,h_t,rF,tau,param).*exp(-1i.*u.*log(K))./(1i*u));

if compute_prices==1
    c = 0.5*S_t + exp(-rF*tau)/pi.*integral(t_Hk,0,1000,'ArrayValued',true);
else
    c = NaN;
end
rn = 0.5 - 1/pi.*integral(t_Hk2,0,1000,'ArrayValued',true);

% Characteristic function of CHJ Model
% u is an N-column vector (always 1)
% tau is a d-column vector of maturities
% Output is d vector
function [Psi] = CF_CHJ(u,S_t,h_t,rF,tau,param)

u1 = 1i*u;
N = numel(u); % always 1 when taking array-valued integral but leaving as is from GARCH toolbox code

T = max(tau);
AMat1 = zeros(N,T);
BMat1 = zeros(N,T);

alpha = param(1);   beta = param(2);   gamma = param(3);    omega = param(4);

% A and B and time T-1
AMat1(:,1) = u1*rF;
BMat1(:,1) = -0.5*u1 + 0.5*u1.^2;

% Recursion back to time t
for t = 2:T
    AMat1(:,t) = AMat1(:,t-1) + u1*rF + BMat1(:,t-1)*omega - 0.5*log(1-2*alpha*BMat1(:,t-1));
    BMat1(:,t) = -0.5*u1 + beta*BMat1(:,t-1) + (alpha*gamma^2)*BMat1(:,t-1) + (0.5*(u1).^2 + 2*alpha*gamma*BMat1(:,t-1).*(BMat1(:,t-1)*alpha*gamma - u1))./(1-2*alpha*BMat1(:,t-1));
end

Psi = transpose(exp(log(S_t').*u1 + AMat1(:,tau) + BMat1(:,tau).*repmat(h_t',N,1)));