(*********************************************************************)
(*                                   APO                             *)
(*  A Mathematica code based on the homotopy analysis method (HAM)   *)
(*  to gain the optimal exercise boundary of American put option     *)
(*  without dividends.                                               *)
(*                                                                   *)
(*  This version is for Mathematica 5.2                              *)
(*                                                                   *)
(*                                                                   *)
(*                                                  Dr. Shijun Liao  *)
(*                                    Shanghai Jiao Tong University  *)
(*                                                            China  *)
(*                                                   December, 2010  *)
(*********************************************************************)
<<Calculus`Pade`;
<<Graphics`Graphics`;

(* Define approx[f] for Taylor expansion of f *)
approx[f_] := Module[{temp},
temp[0] = Series[f, {t, 0, OrderTaylor}]//Normal;
temp[1] = temp[0] /. t^(n_.)*Derivative[j_][DiracDelta][0] -> 0;
temp[2] = temp[1] /. t^(n_.)*DiracDelta[0] -> 0;
temp[3] = temp[2] /. DiracDelta[0] -> 0;
temp[4] = temp[3] /. Derivative[j_][DiracDelta][0] -> 0;
temp[5] = N[temp[4],60]//Expand;
If[KeyCutOff == 1, temp[5] = temp[5]//Chop];
temp[5]
];

(* Define approx2[f] for Taylor expansion of f *)
approx2[f_] := Module[{temp},
temp[0] = Expand[f];
temp[1] = temp[0] /. Derivative[n_][DiracDelta][t] -> dd[n];
temp[2] = temp[1] /. DiracDelta[t] -> dd[0];
temp[3] = Series[temp[2],{t, 0, OrderTaylor}]//Normal;
temp[4] = temp[3] /. dd[0] -> DiracDelta[t];
temp[5] = temp[4] /. dd[n_] -> Derivative[n][DiracDelta][t];
temp[6] = N[temp[5],60]//Expand;
If[KeyCutOff == 1, temp[6] = temp[6]//Chop];
temp[6]
];

(*  Define GetLK[n]  *)  
lamda      := (1 - gamma)/2 - 1/2*Sqrt[(4 p + (1 + gamma)^2)];
kernel[s_] := -s^lamda/lamda;
lK0[0]      = -1/lamda;
lK0[i_]    := D[kernel[s], {s, i}] /. s -> 1 // Expand;

GetLK[m0_,m1_,Nappr_]:= Module[{temp,K1,K2,lK1,lK2},
For[i = Max[m0,0], i <= m1, i++, 
   K[i]   = invl[lK0[i]]; 
   K1[i]  = K[i]/.{Derivative[_][DiracDelta][t_]->0,DiracDelta[t_]->0};
   K2[i]  = Collect[K[i]-K1[i],{DiracDelta[t],Derivative[Blank[]][DiracDelta][t]}];
   temp   = Series[K1[i],{t,0,Nappr}]//Normal;
   lK1[i] = LaplaceTransform[temp, t, p];
   lK2[i] = LaplaceTransform[K2[i],t, p];
   LK[i]  = Collect[lK1[i] + lK2[i], p];
  ];
];

(*  define Getf[n] and Getg[n]  *)
mu[m_,n_]    := If[m == 1, b[n], Sum[mu[m-1,i]*b[n-i],{i,m-1,n-1}]];
psi[n_,m_]   := dV[n,m]/m!;
alpha[n_,i_] := Sum[psi[n,m]*mu[m,i],{m,1,i}];
beta[n_,i_]  := Sum[(m+1)*psi[n,m+1]*mu[m,i],{m,1,i}];
f[0]      := 0;
g[0]      := 1;
Getf[n_]  := Sum[alpha[j,n-j],{j,0,n-1}];
Getg[n_]  := Sum[beta[j,n-j] ,{j,0,n-1}];

(*  Define Getb[n]  *)
b[0]  := 1;
BB[0] := 1;
B[0]  := X;
Getb[n_] := Module[{temp},
If[n == 1,
   b[1] = c0*dV[0, 0] // Expand,
   temp = b[n - 1] + c0*(b[n - 1] + dV[n - 1, 0] + f[n - 1] )//Expand;
   b[n] = approx[temp];
   ];
];

(*  Define GetDV[m,n]  *)
GetDV[m_, n_] := Module[{temp},
If[n == 1, DV[m, 1] = -g[m],        
   temp[1]  = Expand[LK[n]*Lg[m]];
   DV[m, n] = invl[temp[1]];
   ];    
DV[m, n] = approx[DV[m, n]];
];

(*   Define dV[m,n]   *)
dV[m_,n_] := Module[{temp},
If[NumberQ[flag[m,n]],
   Goto[100],
   GetDV[m,n];
   flag[m,n] = 1 
   ];
Label[100];
DV[m,n]
];

(*   Define hp[f_,m_,n_]        *)
hp[f_,m_,n_]:= Module[{k,i,df,res,q},
df[0] = f[0];
For[k = 1, k <= m+n, k++, df[k] = f[k] - f[k-1]//Expand ];
res = df[0] + Sum[df[i]*q^i,{i,1,m+n}];
Pade[res,{q,0,m,n}]/.q->1
];

(*  Get [m,n] Pade approximant of B  *)  
pade[order_]:= Module[{temp,s,i,j},
temp[0] = BB[order] /. t^i_. -> s^(2*i);
temp[1] = Pade[temp[0],{s,0,OrderTaylor,OrderTaylor}];
If[KeyCutOff == 1, temp[1] = temp[1]//Chop];
BBpade[order] = temp[1] /. s^j_. -> t^(j/2);
Bpade[order]  = X*BBpade[order]/. t -> (sigma^2*t/2);
];

(*define the inverse Laplace transformation*)
invl[Sqrt[p]] := -1/(2*Sqrt[Pi]*t^(3/2));

invl[p^n_] := Module[{temp, nInt},
nInt = IntegerPart[n];
If[n > 1/2 && n > nInt,
   Goto[100], 
   temp[2] = InverseLaplaceTransform[p^n, p, t];
   Goto[200];
   ];
Label[100];
temp[1] = -1/2/Sqrt[Pi]/t^(3/2);
temp[2] = D[temp[1], {t, nInt}];
Label[200];
temp[2]//Expand
];

invl[d_./(c_. + a_.*Sqrt[4p + b_.])] := Module[{temp},
temp[1] = d/(4a)*Exp[-b*t/4];
temp[2] = 2/Sqrt[Pi*t];
temp[3] = c/a*Exp[c^2*t/(4a^2)]*Erfc[c*Sqrt[t]/(2a)];
temp[1]*(temp[2]-temp[3])//Expand
];

invl[d_./(p*(c_. + a_.*Sqrt[4p + b_.]))]:= Module[{temp},
temp[1] = Sqrt[b]*Erf[Sqrt[b*t]/2];
temp[2] = c/a*Exp[-(b-(c/a)^2)*t/4]*Erfc[c*Sqrt[t]/(2a)];
temp[3] = -1/(b - (c/a)^2)*d/a*(c/a-temp[1]-temp[2]);
temp[3]//Expand
];

invl[p^i_.*Sqrt[c_.*p + a_.]] := Module[{temp}, 
temp = D[-Exp[-a*t/c]/(2*c*Sqrt[Pi]*(t/c)^(3/2)),{t, i}];
temp//Expand
];

invl[Sqrt[c_.*p + a_.]] := -Exp[-a*t/c]/(2*c*Sqrt[Pi]*(t/c)^(3/2));
invl[f_]     := InverseLaplaceTransform[f, p, t] // Expand;
invl[p_Plus] := Map[invl, p];
invl[c_*f_]  := c*invl[f] /; FreeQ[c, p];

(*  Main code   *)
ham[m0_, m1_] := Module[{temp, k, n},
If[m0 == 1,
Print[" Strike price = ?"];
temp[0] = Input[];
If[!NumberQ[temp[0]],Goto[100]];
X = IntegerPart[temp[0]*10^10]/10^10;
Print[" Risk-free interest rate = ?"];
temp[0] = Input[];
If[!NumberQ[temp[0]],Goto[100]];
r = IntegerPart[temp[0]*10^10]/10^10;
Print[" Volatility = ?"];
temp[0] = Input[];
If[!NumberQ[temp[0]],Goto[100]];
sigma = IntegerPart[temp[0]*10^10]/10^10;
Print[" Time to expiry = ?"];
temp[0] = Input[];
If[!NumberQ[temp[0]],Goto[100]];
T = IntegerPart[temp[0]*10^10]/10^10;
gamma   = 2*r/sigma^2;
texp    = sigma^2*T/2;
Bp      = X*gamma/(1 + gamma);
Label[100];
If[!NumberQ[gamma], 
   X     = .;
   r     = .;
   sigma = .;
   gamma = .;
   T     = .;
   ];
Print["--------------------------------------------------------------"];
Print[" INPUT PARAMETERS: "];
Print["   Strike price            (X) =  ",X," ($) "];
Print["   Risk-free interest rate (r) =  ",r];
Print["   Volatility          (sigma) =  ",sigma];
Print["   Time to expiry          (T) =  ",T," (year)"];
Print["--------------------------------------------------------------"];
Print[" CORRESPONDING PARAMETERS: "];
Print["   gamma  =  ",gamma];
Print["   dimensionless time to expiry   (texp)  =  ",texp//N];
Print["   perpetual optimal exercise price (Bp)  =  ",Bp//N,"($)"];
Print["--------------------------------------------------------------"];
Print[" CONTROL PARAMETERS: "];
Print["   OrderTaylor  =  ",OrderTaylor];
Print["   c0           =  ",c0];
Print["--------------------------------------------------------------"];
KeyCutOff = If[OrderTaylor < 80  &&  NumberQ[gamma], 1, 0];
If[KeyCutOff == 1, 
   Print["Command Chop is used to simplify the result"],
   Print["Command Chop is NOT used "]
   ];
If[NumberQ[gamma],
   Print["Pade technique is used"],
   Print["Pade technique is NOT used"]
   ];
Clear[flag,DV];
];
For[k = Max[1, m0], k <= m1, k++,
    Print[" k  =  ", k];
    If[k == 1, GetKK[]; GetBJ[]; GetKn[]];
    If[k == 1, Lg[0] = LaplaceTransform[g[0], t, p]];
    If[k == 1, GetLK[0,2,OrderTaylor], GetLK[k+1,k+1,OrderTaylor]];
    Getb[k];
    BB[k] = Collect[BB[k - 1] + b[k], t];
    temp[0] = X*BB[k] /. t-> (sigma^2*t/2)//Expand;
    B[k] = Collect[temp[0],t];
    If[NumberQ[gamma],pade[k]];
    temp[1] = Getg[k];
    temp[2] = Getf[k];
    g[k] = approx2[temp[1]];
    f[k] = approx2[temp[2]];
    Lg[k] = LaplaceTransform[g[k], t, p];
    If[NumberQ[gamma] && NumberQ[sigma] && NumberQ[X],
       Print[" Optimal exercise price at the time to expiration = ", B[k]/.t->T//N];
       Print["          Modified result given by Pade technique = ",Bpade[k]/.t->T//N];
       ];
    ];
Print["  Well done !"];
If[NumberQ[gamma] && NumberQ[sigma] && NumberQ[X],
   Plot[{Bp, B[m1], Bpade[m1]}, {t, 0, 1.25*T}, PlotRange -> {0.8*Bp, X}, 
         PlotStyle ->{RGBColor[1, 0, 0], RGBColor[0, 1, 0], RGBColor[0, 0,1]}];
   Print["  Order of homotopy-approximation : ",m1];
   Print["  Green line : optimal exercise boundary B in polynomial "];
   Print["  Blue line  : optimal exercise boundary B by Pade method "];
   Print["  Red line   : perpetual optimal exercise price "];
   ];
];

(*  Dimensionless analytic formula given by Kuske and Keller  *)
GetKK[] := Module[{alpha},
alpha = -Log[9*Pi*gamma^2*t]/2;
KK0 = Exp[-2*Sqrt[alpha*t]]; 
KK  = X*KK0 /. t->sigma^2/2*t;
];

(*  Dimensionless analytic formula given by Kuske and Keller  *)
GetBJ[] := Module[{alpha},
Bp0   = gamma/(1+gamma);
alpha = -Log[4*E*gamma^2*t/(2 - Bp0^2)]/2;
BJ0 = Exp[-2*Sqrt[alpha*t]]; 
BJ  = X*BJ0 /. t->sigma^2/2*t;
];

(*  Dimensionless analytic formula given by Knessl  *)
GetKn[] := Module[{z},
z   = Abs[Log[4*Pi*gamma^2*t]];
Kn0 = Exp[-Sqrt[2*t*z]*(1+1/z^2)];
Kn  = X*Kn0 /. t->sigma^2/2*t;
];

(* Define the order of Taylor's series expansion *)
OrderTaylor   = 8;

(* Assign the convergence-control parameter *)
c0 = -1;

(*   Get 8th-order homotopy-approximation of B  *)
ham[1,8];

(*   Get 10th-order homotopy-approximation of B  *)
ham[8,10]