So I have the following code (fairly long), everything works so far as it asks me for values on certain inputs X,r,sigma,T (as it should) . Before that you to manually input the order of Taylor series expansion, and it then uses that to run the script for K=1,2,3..., how many you asked it to run,. and plots the following graph. However, at the end of the script you have (Dimensionless analytic formula given by Kuske and Keller), (Dimensionless analytic formula given by Bunch and Johnson), (Dimensionless analytic formula given by Knessl). I like the script to call Kuske and Keller,Bunch and Johnson, Knessl - formulas and possibly print them (?) and also to call them separately. Any help will be appreciated.
(* Using the new command of Mathematica 8.0. *)
Pade := PadeApproximant;
(* 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 Ta