% 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