price=EuropeanCompound(s,k1,k2,r,T1,T2,vol,d,cp)

function price=EuropeanCompound(s,k1,k2,r,T1,T2,vol,d,cp)

% calculate european compound option;
% cp=1 call on call; cp=2 call on put;
% cp=3 put on call; cp=4 put on put;
% k2: compound option exercise price; k1:underlying option exercise price;
% T1:compound option exercise date: T2: underlying option exercise date;
%decide whether it s a call or put;
if cp==1 | cp==3
    callput=1; %call option
elseif cp==2 | cp==4
    callput=2; %put option
end
%calculate the critical price making option price=k
sstar=critical(s,k1,k2,r,T2-T1,vol,d,callput);
a1=(log(s/sstar)+(r-d+0.5*vol^2)*T1)/(vol*sqrt(T1));
a2=a1-vol*sqrt(T1);
b1=(log(s/k2)+(r-d+0.5*vol^2)*T2)/(vol*sqrt(T2));
b2=b1-vol*sqrt(T2);
Tratio=sqrt(T1/T2);
if cp==1
    price=s*exp(-d*T2)*cbnd(a1,b1,Tratio)-k2*exp(-r*T2)*cbnd(a2,b2,Tratio)-exp(-r*T1)*k1*normcdf(a2);
elseif cp==4
    price=s*exp(-d*T2)*cbnd(a1,-b1,-Tratio)-k2*exp(-r*T2)*cbnd(a2,-b2,-Tratio)+exp(-r*T1)*k1*normcdf(a2);
elseif cp==2
    price=k2*exp(-r*T2)*cbnd(-a2,-b2,Tratio)-exp(-r*T1)*k1*normcdf(-a2)-s*exp(-d*T2)*cbnd(-a1,-b1,Tratio);
elseif cp==3
    price=k2*exp(-r*T2)*cbnd(-a2,b2,-Tratio)+exp(-r*T1)*k1*normcdf(-a2)-s*exp(-d*T2)*cbnd(-a1,b1,-Tratio);
end

    

function price1=critical(s,k1,k2,r,T,vol,d,callput)
maxiter = 5000;
tol = 1e-6;
o = optimset('MaxIter', maxiter, 'TolFun', tol);
[ss, fval, exitFlag] = fzero(@corp1,[0.00001, 5*s],o,k1,k2,r,T,vol,d,callput);
price1=ss;

function cd=corp1(ss,k1,k2,r,T,vol,d,callput)
[callprice,putprice]=blsprice(ss,k2,r,T,vol,d);
if callput==1
    cd=k1-callprice;
else
    cd=k1-putprice;
end