%  entropy_options.m 
%  Entropy and cumulants for Poisson disaster example 
%  Notation 
%    log re = log r1 + N(mu,sigma^2) + j N(theta,delta^2)
%    j = # poisson jumps, intensity omega 
%  Backus, Chernov, and Martin, "Equity index options"
%  Written:  April 2009 and after 
format compact 
clear all
%close all

disp(' ')
disp('Option model (Merton poisson-normal mixture)') 
disp('--------------------------------------------')
disp(' ')

%  ---------------------------------------------------------------------------------------
%  1. Cumulants 
syms s w mu mus theta thetas real 
syms sigma sigmas delta deltas omega omegas phi positive 
ep = 0.0400;
sigre = 0.1800;

% true parameter values (BCJ) 
omega_num = 1.5120;  theta_num = -0.0259; delta_num = 0.0407;
mu_num = ep - omega_num*theta_num 
sigma_num = sqrt(sigre^2 - omega_num*(theta_num^2+delta_num^2))

% risk-neutral parameter values (BCJ) 
omegas_num = 1.5120; deltas_num = 0.0981; thetas_num = log(1-0.0482)-deltas_num^2/2
mus_num = - sigma_num^2/2 - omegas_num*(exp(thetas_num+deltas_num^2/2)-1) 
phi_num = (delta_num/deltas_num)^2

% generating functions -- several flavors 
%  p and p* (p* same as p with diff parameters) 
cgf_p = s*mu + (sigma*s).^2/2 + omega*(exp(theta*s+(s*delta)^2/2)-1);      

%  log p*/p 
%  h_pi and k from appendix 
cgf_w = s*(s-1)*(mus-mu)^2/(2*sigma^2);
mgf_z1 = phi^(s/2)*((1-s*(1-phi)))^(-1/2)*exp(s*(s-1)*(theta-thetas)^2/(2*deltas^2*(1-s*(1-phi)))); 
cgf_z = s*(omega-omegas) + omega*((omegas/omega)^s*mgf_z1-1);
cgf = cgf_w + cgf_z;
%return 

% some cumulants (another check) (take derivatives of cum gen function, evaluate at s=0) 
disp('Cumulants (1 to 5)') 
cgf1 = diff(cgf_p,s,1);
cgf2 = diff(cgf_p,s,2);
cgf3 = diff(cgf_p,s,3);
cgf4 = diff(cgf_p,s,4);
cgf5 = diff(cgf_p,s,5);

disp(' ')
disp('***for p') 
kappa1 = double(subs(cgf1,[s mu sigma omega theta delta],[0 mu_num sigma_num omega_num theta_num delta_num]))
kappa2 = double(subs(cgf2,[s mu sigma omega theta delta],[0 mu_num sigma_num omega_num theta_num delta_num]))
kappa3 = double(subs(cgf3,[s mu sigma omega theta delta],[0 mu_num sigma_num omega_num theta_num delta_num]))
kappa4 = double(subs(cgf4,[s mu sigma omega theta delta],[0 mu_num sigma_num omega_num theta_num delta_num]))
kappa5 = double(subs(cgf5,[s mu sigma omega theta delta],[0 mu_num sigma_num omega_num theta_num delta_num]))
gamma1 = kappa3/kappa2^1.5
gamma2 = kappa4/kappa2^2

disp(' ')
disp('***for ps') 
kappa1 = double(subs(cgf1,[s mu sigma omega theta delta],[0 mus_num sigma_num omegas_num thetas_num deltas_num]))
kappa2 = double(subs(cgf2,[s mu sigma omega theta delta],[0 mus_num sigma_num omegas_num thetas_num deltas_num]))
kappa3 = double(subs(cgf3,[s mu sigma omega theta delta],[0 mus_num sigma_num omegas_num thetas_num deltas_num]))
kappa4 = double(subs(cgf4,[s mu sigma omega theta delta],[0 mus_num sigma_num omegas_num thetas_num deltas_num]))
kappa5 = double(subs(cgf5,[s mu sigma omega theta delta],[0 mus_num sigma_num omegas_num thetas_num deltas_num]))
gamma1 = kappa3/kappa2^1.5
gamma2 = kappa4/kappa2^2

%return 
disp(' ')
disp('*** for log ps/p') 
cgf = subs(cgf,[mu sigma omega theta delta mus omegas thetas deltas phi],...
        [mu_num sigma_num omega_num theta_num delta_num mus_num omegas_num thetas_num deltas_num phi_num]);
cgf_odd  = (cgf - subs(cgf,s,-s))/2;
cgf_even = (cgf + subs(cgf,s,-s))/2;

cgf1 = diff(cgf,s,1);
cgf2 = diff(cgf,s,2);
cgf3 = diff(cgf,s,3);
cgf4 = diff(cgf,s,4);
cgf5 = diff(cgf,s,5);
kappa1 = double(subs(cgf1,s,0))
%return
kappa2 = subs(cgf2,s,0)
kappa3 = subs(cgf3,s,0)
kappa4 = subs(cgf4,s,0)
kappa5 = subs(cgf5,s,0)
gamma1 = kappa3/kappa2^1.5
gamma2 = kappa4/kappa2^2
entropy = subs(cgf,s,1) - kappa1 
varover2 = kappa2/2
odd  = subs(cgf_odd,s,1) - kappa1 
even = subs(cgf_even,s,1) - kappa2/2 

%return
    
%  Individual cumulants 
disp(' ')
disp('Individual cumulants for pricing kernel') 
ncum = 10;      % make this whatever we want 
icum=[1:ncum]';
kappa = 0.0*icum;
kappap = kappa;
kappaps = kappa;

cgf;
cgf_ps = subs(cgf_p,[mu sigma omega theta delta],[mus_num sigma_num omegas_num thetas_num deltas_num]);
cgf_p  = subs(cgf_p,[mu sigma omega theta delta],[mu_num sigma_num omega_num theta_num delta_num]);

t0 = cputime;
for j = 1:length(icum)
    jcum = icum(j);
    kappaj = diff(cgf,s,jcum);
    kappa(j) = subs(kappaj,[s],[0])/gamma(jcum+1); 
    kappapj = diff(cgf_p,s,jcum);
    kappap(j) = subs(kappapj,[s],[0])/gamma(jcum+1); 
    kappapsj = diff(cgf_ps,s,jcum);
    kappaps(j) = subs(kappapsj,[s],[0])/gamma(jcum+1); 
end 
cumulant_time_secs = cputime - t0

[kappa kappap kappaps]

jfactorial = gamma(icum+1);
%return 

%  figure parameters 
FontSize = 14;
FontName = 'Helvetica';  % or 'Times' 
LineWidth = 2;

clf
figure(1) 
subplot(2,1,1), bar(icum(2:ncum),kappap(2:ncum).*jfactorial(2:ncum),'b')
set(gca,'LineWidth',LineWidth,'FontSize',FontSize,'FontName',FontName)
title('Cumulants of equity return','FontSize',FontSize,'FontName',FontName) 
subplot(2,1,2), bar(icum(2:ncum),kappa(2:ncum),'b')
set(gca,'LineWidth',LineWidth,'FontSize',FontSize,'FontName',FontName)
xlabel('Order j','FontSize',FontSize,'FontName',FontName)

print -depsc cumulants_options_ppsm.eps

return