Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

First timer here and have been finding these boards very useful in learning Mathematica.

I'm trying to implement a numerical procedure for the Hull-White trinomial tree in Mathematica. Despite using memoization, I'm finding it very slow when I get a few steps into the tree.

In the first step I need to solve an equation to get an alpha(x). Only then can I proceed with solving the range of Q(x+1,j). I was expecting the fact that since the x-steps are stored in memory, the x+1 steps would be done rather quickly.

EDIT:

From Q and alpha I can get the rates I am trying to simulate. Once I have my rates, I can calculate a payout at the end of the tree (nodevalue). I then go back through the tree to get the value of the nodevalue(0,0) i.e: The price of such a payout on day1.

This fails for a larger number of steps and is rather slow. Any tips on optimising? I get the following message for 50 steps which takes about 30 seconds:

In[436]:= Timing[nodevalue[0, 0]]

During evaluation of In[436]:= Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>

During evaluation of In[436]:= Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>

During evaluation of In[436]:= Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>

During evaluation of In[436]:= General::stop: Further output of Solve::ratnz will be suppressed during this calculation. >>

Out[436]= {3.719, 0.0100418}

I have included some minor changes from my initial post. Also highlighted in comments where I think the code may be slowing me down.

Unprotect[M]; 
Unprotect[Q]; 
Unprotect[P]; 

(*input here*)

tenor = 5; 
steps = 50; (*<- Code Fails for a number bigger than 51...*)
a = 0.1; 
sigma = 0.01; 
f[r_] := r; g[x_] := x; 
zerorate = 3; dfData = {{0., 1.}, {0.25, 0.992023586787659}, 
      {0.5, 0.984368170142745}, {0.75, 0.977171337237306}, {1., 
   0.970130543530332}, 
      {1.25, 0.962597668870528}, {1.5, 0.955366107894445}, {1.75, 
   0.948429429702776}, 
      {2., 0.941083365231486}, {2.25, 0.933654274758575}, {2.5, 
   0.926206583554808}, 
      {2.75, 0.918741332047699}, {3., 0.911259555496316}, {3.25, 
   0.903801785010502}, 
      {3.5, 0.896334959688613}, {3.75, 0.888859978997039}, {4., 
   0.881295802567575}, 
      {4.25, 0.873576189971488}, {4.5, 0.865919722663361}, {4.75, 
   0.858328277043492}, 
      {5., 0.850550741079475}, {5.25, 0.842846646720007}, {5.5, 
   0.835217448771243}, 
      {5.75, 0.827664434924383}, {6., 0.819937775990933}, {6.25, 
   0.812205585525199}, 
      {6.5, 0.804552663713649}, {6.75, 0.796980136312276}, {7., 
   0.789237583106047}, 
      {7.25, 0.781505055652225}, {7.5, 0.77402686027691}, {7.75, 
   0.766557832213359}, 
      {8., 0.759098675520764}, {8.25, 0.751650086883238}, {8.5, 
   0.744212755551351}, 
      {8.75, 0.736787363285735}, {9., 0.729374584302738}, {9.25, 
   0.721894203851475}, 
      {9.5, 0.714508798782454}, {9.75, 0.707218557903375}, {10., 
   0.699782433930927}, 
      {10.25, 0.690835814909544}, {10.5, 0.681956382812677}, 
      {10.75, 0.673147423947596}, {11., 0.664118563293887}, 
      {11.25, 0.655066776897885}, {11.5, 0.646092713575233}, 
      {11.75, 0.637100851957785}, {12., 0.627994883834108}, 
      {12.25, 0.618863849860742}, {12.5, 0.61001555866861}, 
      {12.75, 0.601452720972706}, {13., 0.592305728759494}, {13.25, 
   0.58344808791873}, 
      {13.5, 0.57459130321309}, {13.75, 0.565737328419695}, {14., 
   0.556888097001247}, 
      {14.25, 0.548045521596268}, {14.5, 0.539211493523719}, 
      {14.75, 0.530387882302119}, {15., 0.521576535183242}, 
      {15.25, 0.513266517651005}, {15.5, 0.505068785221317}, 
      {15.75, 0.496894639770595}, {16., 0.488656405115787}, 
      {16.25, 0.480444844598625}, {16.5, 0.472349543622504}, 
      {16.75, 0.464370791774091}, {17., 0.456245727354391}, 
      {17.25, 0.448152213475503}, {17.5, 0.440178388577009}, 
      {17.75, 0.432324309215909}, {18., 0.424331272572483}, 
      {18.25, 0.416202749238728}, {18.5, 0.408369203376619}, 
      {18.75, 0.400828120870626}, {19., 0.392815561804551}, 
      {19.25, 0.385097371741437}, {19.5, 0.377419607786398}, 
      {19.75, 0.369783221910831}, {20., 0.36218914499444}}; 
P[t_] := P[t] = Interpolation[dfData, t]; 
payout[r_, k_] := payout[r, k] = Max[r - k, 0]; 
strike = 3.19/100; 

(*definitions*)
M = Exp[(-a)*deltat] - 1; 
V = sigma^2*((1 - Exp[-2*a*deltat])/(2*a)); 
deltat = tenor/steps; 
jMax = IntegerPart[-0.184/M] + 1; 
jMin = -jMax; 
deltax = Sqrt[3*V]; 
probMiddle[j_] := {1/6 + (j^2*M^2 + j*M)/2, 2/3 - j^2*M^2, 
   1/6 + (j^2*M^2 - j*M)/2}; 
probUpward[j_] := {1/6 + (j^2*M^2 - j*M)/2, -(1/3) - j^2*M^2 + 2*j*M, 
       7/6 + (j^2*M^2 - 3*j*M)/2}; 
probDownward[j_] := {7/6 + (j^2*M^2 + 3*j*M)/2, -(1/3) - j^2*M^2 - 2*j*M, 
       1/6 + (j^2*M^2 + j*M)/2}; 

q[k_, j_] := Which[jMin < k < jMax, probMiddle[k], jMax <= k, probDownward[k], 
         k <= jMin, 
    probUpward[k]][[Which[jMin < k < jMax, k - j + 2, jMax <= k, 
          k - j + 1, k <= jMin, j - k + 1]]]; 

(*Tree Geometry, suspect sb here could be improved or more elegant*)
sb[i_Integer /; Inequality[0, LessEqual, i, Less, steps], 
       j_Integer /; jMin <= j <= jMax] := 
  Which[i - Abs[j] == 0, If[j > 0, {-1}, {1}], 
         i - Abs[j] == 1 || Abs[j] == jMax, 
    If[j == 0, {0}, If[j > 0, {-1, 0}, {0, 1}]], 
         j > jMax || j < jMin, If[j > 0, {-1, 0}, {0, 1}], 
    jMax - 2 == jMin + 2 && 
           i > jMax && j == 0, {-2, -1, 0, 1, 2}, j == jMax - 2 && i > jMax, 
         {-1, 0, 1, 2}, j == jMin + 2 && i > jMax, {-2, -1, 0, 1}, 
    0 < 1, {-1, 0, 1}] + 
       j; 

n[m_] := Min[jMax, m]


(* =====      THINK BELOW IS WHERE THINGS ARE SLOWING DOWN    =====*)


bigAlpha[0] = zerorate; 
Q[0, 0] = 1; 

bigAlpha[m_Integer /; Inequality[0, LessEqual, m, Less, steps]] := 
     bigAlpha[m] = alpha[m] /. NSolve[P[(m + 1)*deltat] == 

      Sum[Q[m, j]*Exp[(-g[alpha[m] + j*deltax])*deltat], {j, -n[m], n[m]}], 
     alpha[m], 
           Reals]; 


Q[m_Integer /; Inequality[0, LessEqual, m, Less, steps], 
       j_Integer /; jMin <= j <= jMax] := 
     Q[m, j] = Sum[Q[m - 1, sb[m, j][[k]]]*q[sb[m, j][[k]], j]*
           Exp[(-g[bigAlpha[m - 1] + sb[m, j][[k]]*deltax])*deltat], {k, 1, 
     Length[sb[m, j]]}]; 
rate[i_Integer /; Inequality[0, LessEqual, i, Less, steps], 
       j_Integer /; jMin <= j <= jMax] := 
  rate[i, j] = g[bigAlpha[i] + j*deltax]; 
jEnd = Min[jMax, steps]; 


nodevalue[i_Integer /; Inequality[0, LessEqual, i, Less, steps], 
       j_Integer /; jMin <= j <= jMax] := nodevalue[i, j] = 
       Which[i == steps - 1, payout[rate[i, j], strike], jMin < j < jMax, 
         probMiddle[j] . {nodevalue[i + 1, j + 1], nodevalue[i + 1, j], 
             nodevalue[i + 1, j - 1]}, jMax <= j, probDownward[j] . 
           {nodevalue[i + 1, j], nodevalue[i + 1, j - 1], 
      nodevalue[i + 1, j - 2]}, 
         j <= jMin, 
    probUpward[j] . {nodevalue[i + 1, j + 2], nodevalue[i + 1, j + 1], 
             nodevalue[i + 1, j]}]; 

More info on Hull-White tree. Trying to implement pg13-14 in the journal: (PDF)

share|improve this question

migrated from stackoverflow.com Jun 6 '12 at 10:25

This question came from our site for professional and enthusiast programmers.

    
Thanks Heike. Have done that, hopefully its more readable. Is there a way for me to move this post to StackExchange or should I just go ahead and post there? –  jwhatthe Jun 6 '12 at 10:16
    
Thanks Mr.Wizard for sending it over. Much appreciated. I've included the rest of my code as suggested. –  jwhatthe Jun 6 '12 at 10:58
2  
Why the Unprotect? User defined functions/variables aren't normally protected. BTW be aware that using names starting with uppercase letters may cause naming conflicts. In the case of single letter variables, the letters "C", "D", "E", "I", "K", "N", "O" are reserved words, so M, Q and P are OK, but better safe than sorry. –  Sjoerd C. de Vries Jun 6 '12 at 11:06
    
Yes, I realised that it's bad practice to use capital letters so I wanted to Unprotect to be safe. This is to match the notation used in the paper as clsoe as possible. –  jwhatthe Jun 6 '12 at 11:14
5  
"Unprotect to be safe". This is a contradictio in terminis. What if you had unprotected N and redefined it? If you insist on using capitals couldn't you use \[ScriptCapitalM] or so? –  Sjoerd C. de Vries Jun 6 '12 at 11:21

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.