| spectralrecovery(X,bid,ask,F,r,tau,nopt)
|
function out = spectralrecovery(X,bid,ask,F,r,tau,nopt)
%SPECTRALRECOVERY Recover option implied risk-neutral density. (***)
% OUT = SPECTRALRECOVERY(X,bid,ask,F,r,tau)
% OUT = SPECTRALRECOVERY(X,bid,ask,F,r,tau,nopt)
%
% Inputs:
% X [nx1] Strike prices
% bid [nx1] Put option bid quotes , corresponding to X
% ask [nx1] Put option ask quotes, corresponding to X
% F [1x1] Futures price
% r [1x1] Risk-free rate of interest (continuously compounded)
% tau [1x1] Time to maturity of the option set
% nopt Optimal: User provided spectral cutoff value.
%
% If no predetermined spectral cutoff value is supplied, the program
% searches for the smallest value N for which the algorithm yields a
% feasible result.
%
% OUT contains
% .N Smallest spectral cutoff for which all constraints hold
% .X 'Dense' strike grid [0,1,2,...,floor(B),B], B=2*F.
% .densityX Implied risk-neutral density for each point in X
% .returns 'Dense' return grid [-log(2):0.001:log(2)]
% .densityR Implied risk-neutral density for each point in RETURNS
% .prices 'Dense' put option prices for each point in X
% .POINTS Index vector of points with bid/ask provided.
%
% The spectral recovery method of Monnier (2013) takes available put option
% bid/ask quotes and estimates the smoothest RND compatible with the quotes
%
% (***)Monnier (2013)"Risk-neutral density recovery via spectral analysis"
%
% Matthias Held, matthias.held@whu.edu
% $Date: 2013/03/18$
%
% prepare solver
options = optimset( 'Algorithm','interior-point-convex', ...
'MaxIter',175,'TolFun',1E-8,'TolX',1E-8 , ...
'TolCon',1E-8, 'Display','none');
% reshape and prepare input data, sort by strike price
if nargin==6; nopt = [];end %
if size(X,2)>size(X,1); X=X'; end
if size(bid,2)>size(bid,1); bid=bid'; end
if size(ask,2)>size(ask,1); ask=ask'; end
temp = sortrows([X bid ask],1);
X = temp(:,1);
bid = temp(:,2);
ask = temp(:,3);
clear temp;
% choice of B
% set the cutoff strike B. "The higher the B, the higher we will need to go
% into the spectrum of gamma* [the put operator] since the smoothest RND
% that fits the dataw ill be more and more concentrated around the center
% of the intervall I=[0,B]." Monnier (2003), p. 21.
B = ceil(2*F);
xi = [0:1:B]';
% Set the 'dense' strike grid xi
% find those 'dense' strikes xi for which put quotes are observed,
% augment ask quotes for strike zero (P(0)=0) and for strike B.
n = length(xi);
xiBid = find(ismember(xi,X));
xiAsk = [1 ; xiBid ; n];
ask = [0 ; ask ; ask(end)+exp(-r*tau)*(xi(end)-X(end))];
% Iterate on N via Bisection method. Under 15 or so steps should suffice
if isempty(nopt)
NLeft = 4; %
NRight = 200; % maximum 200 (see beta function below)
Nwork = NRight; %
omHatR = 0; exR=-5;
else
NLeft = nopt;
NRight = nopt;
end
brSwitch = 0;
for kLoop = 1:20 % 15 steps should be enough, I guess
N = round(0.5*(NLeft+NRight));
% Bisection method applies only without nopt supplied by user.
if isempty(nopt)&(NRight==N)&(N~=200)
% reuse results from last step to save one iteration.
omHat = omHatR;exitflag= exR;break;
end
if N==200
brSwitch=1;
end
% Setup linear inequality and equality matrices
% Ineqal: bid <= P <= ask; P(0) <= 0; P(B) <= ask(end)+df*(B-X(end));
% P(.) is convex and incraesing.
% Eqal: q(0)=0, sum(q.*xi)=F
Omat = omegaMatrix(N,B); % matrix of inverted lambdas
phi2Temp= phi2(N,xi,B); % pricing matrix
% phi2Temp*omega_star = P(xi)
% Aieq<=bieq
Aieq = [ phi2Temp(xiAsk,:); % P(.) <=
-phi2Temp(xiBid,:); % -P(.)<=
A(xi)*phi2Temp]; % convex, increasing
bieq = [ask ; -bid ; b(xi,r,tau,F)]; % (ask,-bid,convex/increas.)'
% Aeq=beq
Aeq = psi([0:N],0,B)*Omat; % q(0) =
beq = 0; % 0
Aeq = [ Aeq ; % sum(q.*xi) =
exp(r*tau)*xi'*(psi2(N,xi,B)*Omat)];
beq = [beq ; F]; % F
% solve quadratic program
[omHat, ~ , exitflag] = quadprog( Omat^4,zeros(N+1,1),...
Aieq,bieq,Aeq,beq,...
[],[],[],options);
% convergence / choice of next N value
% if N=Nopt is supplied, stop the loop.
if ~isempty(nopt); break; end
% if convergence at N, enter left bracket, else right bracket
if brSwitch ==1
break;
end
if exitflag==1;
NRight = N;
omHatR = omHat;
exR = exitflag;
else;
NLeft=N;
end
end
% normalized return density estimation, cutoff at +/- 50 % return
StDiv = exp(-r*tau)*F;
dRet = 0.001;
RetHi = floor(log(2)/dRet)*dRet;
RetLo = -RetHi;
Ret = [RetLo:dRet:RetHi]';
qRet = q(omHat,r,tau,StDiv.*exp(Ret),B).*StDiv.*exp(Ret)*dRet;
% submit output structure fields
out.exitflag=exitflag;
out.N = N;
out.X = xi;
out.densityX= q(omHat,r,tau,xi,B);
out.returns = Ret;
out.densityR= qRet;
out.prices = P(omHat,xi,B);
out.points = xiBid;
end
% Linear inequality Aieq*m<=bieq of quadratic program
function out = A(xi)
% A(xi)*model_prices<=b(xi,r,tau,F)
% ensures that the resulting put option prices are convex and increasing
% in the strike xi.
N = length(xi);
T = zeros(2*N,N);
T(1:N-2,:) = T(1:N-2,:) + ...
[-diag(1./(xi(2:end-1)-xi(1:end-2))) zeros(N-2,2)];
T(1:N-2,:) = T(1:N-2,:) + ...
[zeros(N-2,1) ...
diag( 1./(xi(2:end-1)-xi(1:end-2)) + ...
1./(xi(3:end)-xi(2:end-1)) ) ...
zeros(N-2,1)];
T(1:N-2,:) = T(1:N-2,:) + ...
[zeros(N-2,2) diag(-1./(xi(3:end)-xi(2:end-1)))];
T(N-1:2*N-2,:) = diag(-ones(N,1));
T(2*N-1,N-1) = -1/(xi(N)-xi(N-1));
T(2*N-1,N) = 1/(xi(N)-xi(N-1));
T(2*N,1) = 1;
T(2*N,2) = -1;
out=T;
end
function out = b(xi,r,tau,F)
% A(xi)*model_prices<=b(xi,r,tau,F)
% ensures that the resulting put option prices are convex and increasing
% in the strike xi.
N =length(xi);
c = zeros(2*N,1);
c(1:N-2) = 0; % convexity
c(N-1:2*N-2) = -max(exp(-r*tau)*(xi-F),0); %df*[xi-F]+ < P,forall xi
c(2*N-1) = exp(-r*tau); % P(B)-P(B-1) < df
c(2*N) = 0; % P(0)-P(1) < 0
out = c;
end
% Pricing and density
function out = P(omega,xi,B)
% Pricing function. Using some input vector omega which is of length N,
% this function computes the model implied prices of options with strikes
% in xi using the equations in section 7.4 of Monnier's paper, p. 19.
N = length(omega)-1;
lx = length(xi);
out=zeros(lx,1);
for l=1:lx
out(l) = sum(omega.*phi([0:N]',xi(l),B));
end
end
function out = q(omega,r,tau,xi,B)
% Density function. Using some input vector omega of length N, q yields
% the option implied risk-neutral probability density at each strike
% in xi, using the equations in section 7.4 of Monnier's paper, p. 19.
N=length(omega)-1;
lx = length(xi);
out=zeros(lx,1);
for l=1:lx
out(l) =exp(r*tau)*sum(lambda([0:N]',B).^(-1).*omega.*psi([0:N]',xi(l),B));
end
end
% smoothness matrix
function out = omegaMatrix(N,B)
% The matrix Omega is defined in proposition 7.1 on p. 20. in his paper.
% Its fourth power reflects the smoothness of the resulting density.
out = (diag(lambda([0:N],B).^(-1)));
end
%% Singular Value Decomposition functions
function out = phi(k,xi,B)
out = a1(k,B).*exp(rho(k).*xi./B) + a2(k,B).*exp(-rho(k).*xi./B) + ...
a3(k,B).*cos(rho(k).*xi./B) + a4(k,B).*sin(rho(k).*xi./B);
end
function out = phi2(N,xi,B)
k = [0:N];
nr = length(xi);
nc = length(k);
out = zeros(nr,nc);
for l=1:nr
out(l,:) = a1(k,B).*exp(rho(k).*xi(l)./B) + ...
a2(k,B).*exp(-rho(k).*xi(l)./B) + ...
a3(k,B).*cos(rho(k).*xi(l)./B) + ...
a4(k,B).*sin(rho(k).*xi(l)./B);
end
end
function out = psi(k,xi,B)
out = a1(k,B).*exp(rho(k).*xi./B) + a2(k,B).*exp(-rho(k).*xi./B) - ...
a3(k,B).*cos(rho(k).*xi./B) - a4(k,B).*sin(rho(k).*xi./B);
end
function out = psi2(N,xi,B)
k = [0:N];
nr = length(xi);
nc = length(k);
out = zeros(nr,nc);
for l=1:nr
out(l,:) = a1(k,B).*exp(rho(k).*xi(l)./B) + ...
a2(k,B).*exp(-rho(k).*xi(l)./B) - ...
a3(k,B).*cos(rho(k).*xi(l)./B) - ...
a4(k,B).*sin(rho(k).*xi(l)./B);
end
end
function out = rho(k)
out=pi/2+k.*pi+(-1).^k.*beta(k);
end
% singular value function
function out = lambda(k,B)
out = (B./rho(k)).^2;
end
% coefficients in spectral decomposition
function out = a1(k,B)
out = 1./sqrt(B) .* (-1).^k ./ (exp(rho(k)) + (-1).^k);
end
function out = a2(k,B)
out = (-1).^k .* exp(rho(k)) .* a1(k,B);
end
function out = a3(k,B)
[r,c]=size(k);
out = -1./sqrt(B)*ones(r,c);
end
function out = a4(k,B)
out = 1./sqrt(B).*(1-(-1).^k .*exp(-rho(k)))./(1+(-1).^k.*exp(-rho(k)));
end
function out = beta(k)
% hardcoded values for solution u to fixed point equation:
% exp(pi/2+k*pi+(-1)^k*u)=(1+cos(u))/sin(u), p.4 in Monnier (2013)
% %b(1) = beta(0)
b = [ 0.304296874999999;0.0182979375620693;0.000775804281234741;
3.35526888492169E-05;1.44989240095019E-06;6.26556264750743E-08;
2.70759494442027E-09;1.17005786949584E-10;5.05627851300998E-12;
2.18501606381235E-13;9.44231055865944E-15;4.08039236703066E-16;
1.76329742232983E-17;7.61989906832813E-19;3.29285695516927E-20;
1.4229725131498E-21;6.14922178748496E-23;2.65731967710182E-24;
1.14833195326992E-25;4.96239231682844E-27;2.14444415972196E-28;
9.26698346394484E-30;4.00462665962646E-31;1.7305560914602E-32;
7.47841094872342E-34;3.23171439481034E-35;1.39655041709189E-36;
6.03504155754466E-38;2.60797792586203E-39;1.12700944921924E-40;
4.87024942210596E-42;2.10462560451071E-43;9.09491188491674E-45;
3.93026778811002E-46;1.69842270950123E-47;7.33955001457214E-49;
3.17170714422589E-50;1.37061893287199E-51;5.92298145359096E-53;
2.55955236413318E-54;1.10608286655496E-55;4.77981745882654E-57;
2.06554640981478E-58;8.92603537237651E-60;3.85728963001422E-61;
1.66688598791128E-62;7.20326748366261E-64;3.11281412271079E-65;
1.34516895069141E-66;5.81300210861437E-68;2.51202597988793E-69;
1.0855448537135E-70;4.69106465800331E-72;2.02719284701028E-73;
8.7602945995607E-75;3.78566654791975E-76;1.63593484775703E-77;
7.06951547958296E-79;3.05501464099165E-80;1.32019152990442E-81;
5.70506488658148E-83;2.46538207698245E-84;1.0653881956369E-85;
4.60395984054396E-87;1.9895514423896E-88;8.59763134129952E-90;
3.71537338044957E-91;1.60555841582145E-92;6.93824701490185E-94;
2.99828839395849E-95;1.29567789587734E-96;5.59913186886186E-98;
2.4196042692908E-99;1.0456058112381E-100;4.5184724062967E-102;
1.9526089724282E-103;8.4379884535621E-105;3.6463854334288E-106;
1.575746020784E-107;6.8094159746625E-109;2.9426154535309E-110;
1.271619436906E-111;5.4951658409081E-113;2.3746764749486E-114;
1.0261907509133E-115;4.4345723233009E-117;1.9163524591391E-118;
8.2813098533818E-120;3.5786784712102E-121;1.5464871894719E-122;
6.682977100181E-124;2.8879762616589E-125;1.2480077011904E-126;
5.3931302791802E-128;2.3305829107036E-129;1.0071361941009E-130;
4.3522301172363E-132;1.8807691655138E-133;8.1275404991538E-135;
3.5122287081591E-136;1.5177716431807E-137;6.5588859731481E-139;
2.8343516234503E-140;1.2248343938657E-141;5.2929893383167E-143;
2.2873080863704E-144;9.8843544688477E-146;4.2714168610694E-147;
1.8458465910475E-148;7.9766263712992E-150;3.4470128002996E-151;
1.4895892940632E-152;6.437099000024E-154;2.7817227004285E-155;
1.202091374088E-156;5.1947078385421E-158;2.2448367993894E-159;
9.7008193964322E-161;4.1921041648904E-162;1.811572467348E-163;
7.8285144532878E-165;3.3830078371113E-166;1.4619302415863E-167;
6.3175733967248E-169;2.7300710039146E-170;1.1797706521743E-171;
5.0982512533071E-173;2.2031541294862E-174;9.5206922474069E-176;
4.1142641659416E-177;1.7779347538261E-178;7.683152713008E-180;
3.3201913334816E-181;1.4347847690531E-182;6.2002671736177E-184;
2.6793783885376E-185;1.1578643867918E-186;5.0035856971286E-188;
2.1622454334165E-189;9.343909742688E-191;4.0378695187891E-192;
1.7449216334592E-193;7.5404900843292E-195;3.2585412218713E-196;
1.4081433401327E-197;6.0851391200938E-199;2.6296270456143E-200;
1.1363648823966E-201;4.910677914114E-203;2.1220963395445E-204;
9.1704097767052E-206;3.962893384877E-207;1.7125215077821E-208;
7.4004764433878E-210;3.1980358437824E-211;1.3819965958997E-212;
5.9721488046185E-214;2.5807995015207E-215;1.1152645859538E-216;
4.8194952748242E-218;2.0826927439022E-219;9.000131388405E-221;
3.8893094263271E-222;1.6807229829651E-223;7.2630627550915E-225;
3.1386538548829E-226;1.3563354088323E-227;5.8612566458022E-229;
2.5328784736287E-230;1.0945560461339E-231;4.7300060289896E-233;
2.0440210748331E-234;8.833015873825E-236;3.8170929267781E-237;
1.6495147726001E-238;7.1282063999879E-240;3.0803756760086E-241;
1.3311523095098E-242;5.7524142252015E-244;2.4858450957909E-245;
1.074235073086E-246;4.6421979697622E-248;2.006061281923E-249;
8.6690451234472E-251;3.7462533308605E-252;1.6188589162288E-253;
6.9957137077231E-255;3.0231243854192E-256;1.3063954638211E-257;
5.645435401422E-259;2.4397634197098E-260;1.0541809751804E-261;
4.5561391826726E-263;1.9688239427666E-264;8.510837770087E-266;
3.6776866659057E-267;1.5885285269736E-268;6.8660345654539E-270;
2.9691218129735E-271;1.2814558943987E-272;5.5320445505565E-274];
if size(k,2)>size(k,1)
b=b';
end
out=b(k+1);
end |
|
