Message Boards Message Boards

0
|
5772 Views
|
6 Replies
|
4 Total Likes
View groups...
Share
Share this post:

Optimize code with Do Loop?

Posted 8 years ago

If anyone is willing to look.. I'd like to know how can I speed up the attached piece of code?

In R, a similar implementation runs 25% faster. Shouldn't mention perhaps, but C++ runs it ~40x faster.

Thank you.

Attachments:
POSTED BY: Sandu Ursu
6 Replies
Posted 8 years ago

Sander, I always thought that Block is faster (true that only marginally) than Module. Do you use it for any particular reason which makes Block look unfit?

POSTED BY: Sandu Ursu

Somehow, I would have expected more specialized functions such as ListConvolve or MovingAverage to perform better than qU Most[#] + qD Rest[#] &, but it is not true:

Do[Nest[qU Most[#] + qD Rest[#] &, call, n], {1000}] // Timing
Do[Nest[ListConvolve[{qD, qU}, #] &, call, n], {1000}] // Timing
Do[Nest[MovingAverage[#, {qU, qD}] &, call, n], {1000}] // Timing
POSTED BY: Gianluca Gorni
Posted 8 years ago

This is fantastic!

Sander, please see the attached file for my C++ code. Compiled with g++ it actually runs in 116 ms. [1000 iterations]


If compiled in VS 2015 it employs all the arsenal of SIMD and performs the iterations in ~85 ms.

Attachments:
POSTED BY: Sandu Ursu

Could you post your c code?

POSTED BY: Sander Huisman

Several ways are possible:

cf = Compile[{{call, _Real, 1}, {qU, _Real, 0}, {qD, _Real, 0}, {n, _Integer, 0}},
   Module[{c = call},
    Do[
     c = qU Most[c] + qD Rest[c];
     ,
     n
     ];
    c
    ]
   ];
ClearAll[DoItOrig, DoIt, DoIt2, DoItC]
DoItOrig[] := Module[{S, K, \[Alpha], \[Sigma], T, r, n, dt, R, u, d, qU, qD, call},
  S = 100;
  K = 100;
  \[Alpha] = 0.07;
  \[Sigma] = 0.2;
  T = 1;
  r = 0.03;
  n = 252;
  dt = T/n;
  R = E^(r dt);
  u = E^(\[Alpha] dt + \[Sigma] Sqrt[dt]);
  d = E^(\[Alpha] dt - \[Sigma] Sqrt[dt]);
  qU = (R - d)/(R (u - d));
  qD = (1 - R qU)/R;
  call = Clip[S u^Range[n, 0, -1] d^Range[0, n] - K, {0, \[Infinity]}];
  Do[call = qU call[[;; -2]] + qD call[[2 ;;]], {n}];
  call
  ]
DoIt[] := Module[{S, K, \[Alpha], \[Sigma], T, r, n, dt, R, u, d, qU, qD, call},
  S = 100;
  K = 100;
  \[Alpha] = 0.07;
  \[Sigma] = 0.2;
  T = 1;
  r = 0.03;
  n = 252;
  dt = T/n;
  R = E^(r dt);
  u = E^(\[Alpha] dt + \[Sigma] Sqrt[dt]);
  d = E^(\[Alpha] dt - \[Sigma] Sqrt[dt]);
  qU = (R - d)/(R (u - d));
  qD = (1 - R qU)/R;
  call = Ramp[S u^Range[n, 0, -1] d^Range[0, n] - K];
  Nest[qU Most[#] + qD Rest[#] &, call, n]
  ]
DoItC[] := Module[{S, K, \[Alpha], \[Sigma], T, r, n, dt, R, u, d, qU, qD, call},
  S = 100;
  K = 100;
  \[Alpha] = 0.07;
  \[Sigma] = 0.2;
  T = 1;
  r = 0.03;
  n = 252;
  dt = T/n;
  R = E^(r dt);
  u = E^(\[Alpha] dt + \[Sigma] Sqrt[dt]);
  d = E^(\[Alpha] dt - \[Sigma] Sqrt[dt]);
  qU = (R - d)/(R (u - d));
  qD = (1 - R qU)/R;
  call = Ramp[S u^Range[n, 0, -1] d^Range[0, n] - K];
  cf[call, qU, qD, n]
  ]
DoIt2[] := Module[{S, K, \[Alpha], \[Sigma], T, r, n, dt, R, u, d, qU, qD, call, pr},
  S = 100;
  K = 100;
  \[Alpha] = 0.07;
  \[Sigma] = 0.2;
  T = 1;
  r = 0.03;
  n = 252;
  dt = T/n;
  R = E^(r dt);
  u = E^(\[Alpha] dt + \[Sigma] Sqrt[dt]);
  d = E^(\[Alpha] dt - \[Sigma] Sqrt[dt]);
  qU = (R - d)/(R (u - d));
  qD = (1 - R qU)/R;
  call = Ramp[S u^Range[n, 0, -1] d^Range[0, n] - K];
  pr = qU^Range[n, 0, -1] qD^Range[0, n] Binomial[n, Range[0, n]];
  pr.call
  ]

With the following timings:

RepeatedTiming[a = DoItOrig[]]
RepeatedTiming[b = DoIt[]]
RepeatedTiming[c = DoItC[]]
RepeatedTiming[d = DoIt2[]]

{0.00187, {9.41536}}
{0.00072, {9.41536}}
{0.00044, {9.41536}}
{0.00048, 9.41536}

I don't think the Compile works as fast as it can because the array is resizing all the time. If you would keep the array fixed-length somehow it would be much faster I expect.

POSTED BY: Sander Huisman

This version runs faster:

S = 100;
K = 100;
\[Alpha] = 0.07;
\[Sigma] = 0.2;
T = 1;
r = 0.03;
n = 252;
dt = T/n;
R = E^(r dt);
u = E^(\[Alpha] dt + \[Sigma] Sqrt[dt]);
d = E^(\[Alpha] dt - \[Sigma] Sqrt[dt]);
qU = (R - d)/(R (u - d));
qD = (1 - R qU)/R;
call = Clip[S u^Range[n, 0, -1] d^Range[0, n] - K, {0, \[Infinity]}];
Do[Nest[qU Most[#] + qD Rest[#] &, call, n], {1000}] // Timing
POSTED BY: Gianluca Gorni
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