Message Boards Message Boards

0
|
6547 Views
|
1 Reply
|
0 Total Likes
View groups...
Share
Share this post:

Fix recursive algorithm that doesn't call certain parts?

Posted 7 years ago

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
POSTED BY: k b
Posted 7 years ago

Perhaps there is a size limit or perhaps there was another reason that the end of your notebook appears to have been cut off and not shown. If you could edit your post to have all the code then it would be more likely that someone could suggest how to get it to work.

POSTED BY: Bill Simpson
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract