%  entropy_poisson.m 
%  Entropy and cumulants for Poisson disaster example 
%  Two components, normal and poisson 
%  Backus, Chernov, and Martin, "Equity index options"
%  Written:  April 2009 and after 
format compact 
clear all
close all

disp(' ')
disp('Entropy and cumulants for poisson-normal consumption growth') 
disp('---------------------------------------------------------------')
disp(' ')

%  ---------------------------------------------------------------------------------------
%  1. Symbolic math calculation of cgf etc for log g
syms s omega theta sigma alpha delta real 

cgf_w = (sigma*s).^2/2;
cgf_z = omega*(exp(theta*s+(s*delta)^2/2)-1) - omega*theta*s;       % demeaned 
cgf = cgf_w + cgf_z;

% some cumulants (another check) (take derivatives of cum gen function, evaluate at s=0) 
disp('Cumulant formulas (1 to 5)') 
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 = factor(subs(cgf1,s,0))
kappa2 = factor(subs(cgf2,s,0))
kappa3 = factor(subs(cgf3,s,0))
kappa4 = factor(subs(cgf4,s,0))
kappa5 = factor(subs(cgf5,s,0))

gamma1 = kappa3/kappa2^1.5
gamma2 = kappa4/kappa2^2

%return 

%  ---------------------------------------------------------------------------------------
%  2. Summary statistics for entropy and components 
%     computation:  entropy + odd/even components   
%
%  parameter inputs
sigma_g = 0.035; mu_g = 0.020;                  % Table 1 ballpark numbers for log cons growth    
kappa1g = mu_g; kappa2g = sigma_g^2;            % same numbers, new notation 
ep = 0.040;  logr1 = 0.020; lambda = 18/3.5;    % returns and leverage 

% Barro (QJE, 2006, Fig IB) [mean jump = -0.35, std dev = 0.13] 
%dy = -[0.17 0.22 0.27 0.32 0.37 0.42 0.47 0.52 0.57 0.67]'; n = [2 12 14 7 9 5 1 2 4 4]';
%prob = n/sum(n); Edy = prob'*dy Stdy = sqrt(prob'*(dy.*dy) - Edy^2)

omega_num = 0.01; theta_num = -0.3; delta_num = 0.15;      % adapted from Barro (QJE, 2006, table 1) 
mu_num = kappa1g - omega_num*theta_num
sigma_num = sqrt(sigma_g^2-omega_num*(theta_num^2+delta_num^2))

cgf_odd  = (cgf - subs(cgf,s,-s))/2;
cgf_even = (cgf + subs(cgf,s,-s))/2;

disp(' ')
disp('Components of entropy (alpha, entropy, variance, odd, even)')  
alphasome = [2 5 10 5.1893];
for i = 1:length(alphasome)
%    disp(' ')
    alpha = alphasome(i);
    entropy = subs(cgf,[s sigma omega theta delta],[-alpha sigma_num omega_num theta_num delta_num]);
    odd = double(subs(cgf_odd,[s sigma omega theta delta],[-alpha sigma_num omega_num theta_num delta_num]));
    even = subs(cgf_even,[s sigma omega theta delta],[-alpha sigma_num omega_num theta_num delta_num]);
    var = double(subs(kappa2*alpha^2/2,[s sigma omega theta delta],[-alpha sigma_num omega_num theta_num delta_num]));
    [alpha entropy var odd even-var] 
end     

varpercent = var./entropy
oddpercent = odd./entropy
evenpercent = (even-var)./entropy

%return 

%  ---------------------------------------------------------------------------------------
%  3. Choosing alpha to match equity premium 

disp('Value of alpha that hits ep for lognormal case')
alpha_ln = (2*ep/(lambda*kappa2g) + lambda)/2  
ent_ln = alpha_ln^2*kappa2g/2
logr1_ln = alpha_ln*kappa1g - ent_ln 

disp(' ')
disp('Value of alpha that hits ep for poisson case') 

alphalast = 3;
alphanow = 7;
eplast = subs(cgf,[s sigma omega theta delta],[-alphalast sigma_num omega_num theta_num delta_num]) ...
                - subs(cgf,[s sigma omega theta delta],[lambda-alphalast sigma_num omega_num theta_num delta_num]); 
epnow = subs(cgf,[s sigma omega theta delta],[-alphanow sigma_num omega_num theta_num delta_num]) ...
                - subs(cgf,[s sigma omega theta delta],[lambda-alphanow sigma_num omega_num theta_num delta_num]);   
tol = 1.e-8;
maxit = 50;
for it = 1:maxit        %  secant method, function is equity premium - ep (second is data)
    alphanext = alphanow - (epnow-ep)*(alphanow-alphalast)/(epnow-eplast);
    epnext = subs(cgf,[s sigma omega theta delta],[-alphanext sigma_num omega_num theta_num delta_num]) ...
                - subs(cgf,[s sigma omega theta delta],[lambda-alphanext sigma_num omega_num theta_num delta_num]); 
    if max(abs(epnext-ep)) < tol, break, it, end 
    alphalast = alphanow;
    alphanow = alphanext;
    eplast = epnow;
    epnow = epnext;
end    

it
alpha = alphanext
epnext
error = epnext-ep

%return 

disp(' ')
disp('Risk-neutral parameters') 
mus_num = mu_num - alpha*sigma_num^2
omegas_num = omega_num*exp(-alpha*theta_num+(alpha*delta_num)^2/2)
thetas_num = theta_num - alpha*delta_num^2

disp(' ')
disp('Skewness and kurtosis (true, risk-neutral, log m)') 
gamma1_p = subs(gamma1,[sigma omega theta delta],[sigma_num omega_num theta_num delta_num])
gamma2_p = subs(gamma2,[sigma omega theta delta],[sigma_num omega_num theta_num delta_num])
gamma1_ps = subs(gamma1,[sigma omega theta delta],[sigma_num omegas_num thetas_num delta_num])
gamma2_ps = subs(gamma2,[sigma omega theta delta],[sigma_num omegas_num thetas_num delta_num])
gamma1_m = -subs(gamma1,[sigma omega theta delta],[sigma_num omega_num theta_num delta_num])
gamma2_m = subs(gamma2,[sigma omega theta delta],[sigma_num omega_num theta_num delta_num])

%return 

%  ---------------------------------------------------------------------------------------
%  4. Plots of entropy, logr1, and ep v. alpha 
%  figure parameters 
FontSize = 14;
FontName = 'Helvetica';  % or 'Times' 
LineWidth = 2;

alphagrid = [0.1:0.1:12]';
zerogrid = 0*alphagrid;
entropy_ln = zerogrid;
entropy_ds = zerogrid;
entropy_bm = zerogrid;
logr1_data = zerogrid + logr1;
ep_data = zerogrid + ep;
ep_ds = zerogrid;

%  vectorize?? [would this work if we added .*s to cgf?] 
for it = 1:length(alphagrid)
    alpha = alphagrid(it);
    entropy_ds(it) = subs(cgf,[s sigma omega theta delta],[-alpha sigma_num omega_num theta_num delta_num]);        
    entropy_bm(it) = subs(cgf,[s sigma omega theta delta],[-alpha sigma_num omega_num -theta_num delta_num]);        
    ep_ds(it) = subs(cgf,[s sigma omega theta delta],[-alpha sigma_num omega_num theta_num delta_num]) ...
                - subs(cgf,[s sigma omega theta delta],[lambda-alpha sigma_num omega_num theta_num delta_num]);   
end    

entropy_ln = alphagrid.^2*kappa2g/2; 
logr1_ln = alphagrid*kappa1g - entropy_ln; 
logr1_ds = alphagrid*kappa1g - entropy_ds;
ep_ln = (alphagrid.^2-(lambda-alphagrid).^2)*kappa2g/2;

%return 

figure(1)
plot(alphagrid,ep_data,'LineWidth',LineWidth,'Color','k','LineStyle','-')
set(gca,'LineWidth',LineWidth,'FontSize',FontSize,'FontName',FontName)
hold on
axis([0 12 0 0.4]) 
xlabel('Risk Aversion \alpha','FontSize',FontSize,'FontName',FontName)
ylabel('Entropy of Pricing Kernel L(m)','FontSize',FontSize,'FontName',FontName)
text(0.30,0.06,'equity premium','FontSize',FontSize,'FontName',FontName)
plot(alphagrid,entropy_ln,'LineWidth',LineWidth,'Color','r','LineStyle','-')
text(8.5,0.07,'normal','FontSize',FontSize,'FontName',FontName)
print -depsc entropy_pn_log.eps
%
plot(alphagrid,entropy_ds,'LineWidth',LineWidth,'Color','b','LineStyle','-')
text(8.5,0.2,'disasters','FontSize',FontSize,'FontName',FontName)
print -depsc entropy_pn_ds.eps
%
plot(alphagrid,entropy_bm,'LineWidth',LineWidth,'Color','m','LineStyle','-')
text(8.5,0.01,'booms','FontSize',FontSize,'FontName',FontName)
print -depsc entropy_pn_bm.eps
print -depsc entropy_pn_all.eps

return 

ratio = [alphagrid entropy_ds./entropy_ln];
i2 = find(ratio==2);
i8 = find(ratio==8)
disp('Selected risk aversion and entropy ratios (total to lognormal)') 
ratio(i2,:)
ratio(i8,:)

figure(2)
plot(alphagrid,logr1_ln,'LineWidth',LineWidth,'Color','r','LineStyle','-')
set(gca,'LineWidth',LineWidth,'FontSize',FontSize,'FontName',FontName)
hold on
plot(alphagrid,logr1_ds,'LineWidth',LineWidth,'Color','b','LineStyle','-')
plot(alphagrid,logr1_data,'LineWidth',LineWidth,'Color','k','LineStyle','-')
xlabel('Risk Aversion \alpha','FontSize',FontSize,'FontName',FontName)
ylabel('One-Period Interest Rate log r^1','FontSize',FontSize,'FontName',FontName)
text(1,-0.02,'sample mean','FontSize',FontSize,'FontName',FontName)
text(8,0.16,'normal','FontSize',FontSize,'FontName',FontName)
text(8,0.06,'disasters','FontSize',FontSize,'FontName',FontName)

print -depsc logr1_pn.eps

%return 

figure(3)
plot(alphagrid,ep_ln,'LineWidth',LineWidth,'Color','r','LineStyle','-')
set(gca,'LineWidth',LineWidth,'FontSize',FontSize,'FontName',FontName)
hold on
plot(alphagrid,ep_ds,'LineWidth',LineWidth,'Color','b','LineStyle','-')
plot(alphagrid,ep_data,'LineWidth',LineWidth,'Color','k','LineStyle','-')
xlabel('Risk Aversion \alpha','FontSize',FontSize,'FontName',FontName)
ylabel('Equity Premium','FontSize',FontSize,'FontName',FontName)
text(0.30,0.08,'sample mean','FontSize',FontSize,'FontName',FontName)
text(8,0.00,'normal','FontSize',FontSize,'FontName',FontName)
text(8,0.30,'disasters','FontSize',FontSize,'FontName',FontName)

print -depsc ep_pn.eps

%return 

figure(4)
plot(alphagrid,entropy_ds,'LineWidth',LineWidth,'Color','r','LineStyle','-')
set(gca,'LineWidth',LineWidth,'FontSize',FontSize,'FontName',FontName)
hold on
plot(alphagrid,ep_ds,'LineWidth',LineWidth,'Color','b','LineStyle','-')
plot(alphagrid,ep_data,'LineWidth',LineWidth,'Color','k','LineStyle','-')
xlabel('Risk Aversion \alpha','FontSize',FontSize,'FontName',FontName)
ylabel('Entropy and Equity Premium','FontSize',FontSize,'FontName',FontName)
text(0.30,0.12,'sample mean = lower bound','FontSize',FontSize,'FontName',FontName)
text(8,0.12,'equity premium','FontSize',FontSize,'FontName',FontName)
text(8,0.5,'entropy','FontSize',FontSize,'FontName',FontName)

print -depsc entropy_ep_pn.eps

%return 

%  ---------------------------------------------------------------------------------------
%  5. Individual cumulants 
%     quick and dirty calculations and plots 
%     computation:  take cumulants of log g and multiply by (-alpha)^j 
disp(' ')
ncum = 8;      % make this whatever we want 
icum=[1:ncum]';
kappa = 0.0*icum;

t0 = cputime;
%kappaj = cgf;   %  why is sequential version slower?  
for j = 1:length(icum)
    jcum = icum(j);
    kappaj = diff(cgf,s,jcum);
%    kappaj = diff(kappaj,s,1);
    kappa(j) = subs(kappaj,[s sigma omega theta delta],[0 sigma_num omega_num theta_num delta_num]);
end 
derivtime = cputime - t0

kappa;
disp('Standardized skewness and kurtosis') 
gamma1 = kappa(3)/kappa(2)^(3/2)
gamma2 = kappa(4)/kappa(2)^2 

kappaall = kappa(:,ones(length(alphasome),1));
icumall = icum(:,ones(length(alphasome),1));
alphaall = -alphasome(ones(size(icum)),:);
alphaall = alphaall./icumall;

alphacum = cumprod(alphaall);

% scale cumulants by (-alpha)^j/j
kappam = alphacum.*kappaall;
disp(' ')
disp('Individual cumulants [logg logm for alpha = (2,5,10)]') 
[kappa kappam] 

%return 

figure(5)
subplot(3,1,1), bar(icum(2:ncum),kappa(2:ncum),'b')
set(gca,'LineWidth',LineWidth,'FontSize',FontSize,'FontName',FontName)
ylabel('Cumulant','FontSize',FontSize,'FontName',FontName)

subplot(3,1,2), bar(icum(2:ncum),kappam(2:ncum,1),'b')
set(gca,'LineWidth',LineWidth,'FontSize',FontSize,'FontName',FontName)
ylabel('Contribution','FontSize',FontSize,'FontName',FontName)
text(6.7,0.00225,'risk aversion \alpha=2','FontSize',FontSize,'FontName',FontName)

subplot(3,1,3), bar(icum(2:ncum),kappam(2:ncum,3),'b')
set(gca,'LineWidth',LineWidth,'FontSize',FontSize,'FontName',FontName)
ylabel('Contribution','FontSize',FontSize,'FontName',FontName)
xlabel('Order j','FontSize',FontSize,'FontName',FontName)
text(6.7,0.08,'risk aversion \alpha=10','FontSize',FontSize,'FontName',FontName)

print -depsc cumulants_pn.eps

return