Increasing recursion speed in Hull-White trinomial tree calculation

Question

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)