I've been reading about functions in Mathematica and playing with this for hours, but I'm clearly missing something. The following code correctly calculates "fail", the number of failed scenarios.

I want to convert the code between the dashes to a function of the variable "spend" that returns the value of fail as the value of that function, in other words f[spend]:=fail. Can someone please show me how to do this while defining the code as a Module?


n = 10000; (* Number of samples to generate  *)
years = 30; (* Number  of years in each simulated scenario *)
mu = .056; (* Average portfolio rate of return *)
sigma = .11; (* Portfolio standard deviation *)
portfolio = 1000000; (* Initial portfolio value *)

withdrawal = .042; (* Annual spending percentage of initial portfolio value *)
returns = Table[1, {n}]; (* one scenario of market returns *)
spend = withdrawal*portfolio; 

(* --------------------------- *)

fail = 0;

For[k = 1, k <= n, k++,
 returns = RandomVariate[LogNormalDistribution[mu, sigma], {years}];
 x = portfolio;

 For[i = 1, i <= years, i++,
  x = returns[[i]]*x - spend
 fail = fail + If[x <= 0, 1, 0];

(* --------------------------- *)

Print["Failed: ", fail]

Print["Percent failed: ", N[fail/n, 3]]
  • $\begingroup$ Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. $\endgroup$
    – bbgodfrey
    May 3, 2015 at 15:50

1 Answer 1


You just need fail as the last expression in the module, so that it gets returned:

f[spend_] := Module[{fail = 0, returns, x},
  For[k = 1, k <= n, k++, 
   returns = RandomVariate[LogNormalDistribution[mu, sigma], {years}];
   x = portfolio;
   For[i = 1, i <= years, i++, x = returns[[i]]*x - spend];
   fail = fail + If[x <= 0, 1, 0];];

A faster approach will be to generate all the returns at once and use Fold to do the iterative calculation and Total[UnitStep[...]] to count the negative results:

f[spend_] := Total[UnitStep[-Fold[Times[##] - spend &, portfolio, #] & /@ 
    RandomVariate[LogNormalDistribution[mu, sigma], {n, years}]]]
  • $\begingroup$ Thanks! Still cant see what I was doing wrong, but I think I was messing up the syntax around that last "fail" in the function. Your alternative is elegant. I'll need to study that for a while. Out of curiosity, should mu, sigma and years also be localized in the first approach? (I'm switching to your function, regardless.) $\endgroup$ May 3, 2015 at 17:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.