% 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