I'm doing a finance project for my differential equations class to explore stochastic processes. I'm exploring IRA retirement funds.

From the project write-up:

Standing assumption: at time of retirement, you have 500 thousand dollars invested in an IRA retirement account, from which you will withdraw 40 thousand dollars per year until the account is exhausted—call this length of time the lifetime of the account.

Question 1: Assume the IRA pays (continuously) at a risk-free annual interest rate of 4%. In this case, the balance B is described by the IVP:

B’ = .04*B – 40, B(0) = 500

Easy to solve but I've included it so that those who answer know the context of the project and are financially-savvy. The next part involves the stochastic process:

Adding randomness: For remaining questions, we will assume the interest rate is not constant over time, but has some randomness built in. The simplest approach is to say that the annual interest rate is given by:

r = .04 + v*dW(t)

where dW(t) is the (standard) Wiener process (i.e., Brownian motion) and v is a measure of the volatility. (In effect, this makes the interest rate follow a Brownian motion path, wandering around 4%. The larger the volatility, the faster it wanders.) Thus the IRA’s balance is described by a stochastic Ito equation saying

dB = (.04 + v*dW(t))*B – 40, B(0) = 500

In Mathematica, I've described this equation as proc1. The WienerProcess[0, 0.25] implies 0 drift and 25% volatility.

proc1 = ItoProcess[\[DifferentialD]B[t] == (.04*B[t] - 40) \[DifferentialD]t +
.04*B[t]*\[DifferentialD]w[t], B[t], {B, 500}, t, w \[Distributed] WienerProcess[0, 0.25]]

data1 = RandomFunction[proc1, {0., 20., 0.01}, 50]

data1 is the RandomFunction using proc1 as its process with from 0 to 20, applying the process at every 0.01 step. It does this to generate 50 "walks" and returns the TemporalData[] structure associated with this RandomFunction

I've never seen TemporalData[] before and what I've seen in the Function Navigator hasn't been helpful to me.

What I'm trying to do:

Determine the total number of "paths" with "state" > 0 at "time" == targetLife or, can also determine the total number of "paths" with "state" < 0 at "time" == targetLife

So far, I've been working on the single element case:

state = Table[Last[data1["Path"][[i]], {i, 1, Length[data1["Path"]-1]}]
lifetime = Length[TakeWhile[state, # > 0&]]-1

I subtract 1 from lifetime to account for the initial condition.

It seems to work, but I can't seem to extend it to cover all paths. Any help with this would be greatly appreciated.

Thanks in advance!


2 Answers 2


If you have version 10 or higher you can use a combination of TimeSeriesMap and TimeSeriesThread.

bool = TimeSeriesMap[Boole[# > 0] &, data1];
tot = TimeSeriesThread[Total, bool];

enter image description here

The total tot is a TemporalData object so we could use the "SliceData" property to obtain the number of paths at a time t that have positive state.

tot["SliceData", 17.]
(* {40.} *)

However, since there is only one path, it makes sense to convert to TimeSeries for easier processing...

tot = TimeSeries@tot
(* 40. *)
  • $\begingroup$ Could you explain a little more about what the x-axis is measuring? It looks like individual time steps. Also, I do not understand how this counts the number of paths in TemporalData who have a positive balance after n years - could you elaborate? $\endgroup$
    – ATAX
    Dec 10, 2014 at 2:56
  • $\begingroup$ Could you be more explicit with what you mean by "values"? I like the way you tested for the "state", but the graphs look as if they are just testing one "path" at each time step. Am I just interpreting the graph incorrectly? Could you possibly add, as label or some other object, the number of "paths" that exceed a target lifetime (the original D.E. gave a lifetime of ~17.33 years) on your graph? Sorry if I'm coming off as demanding, I just do not see why this way gives tells me the total number of paths that exceed a target lifetime. $\endgroup$
    – ATAX
    Dec 10, 2014 at 3:19
  • $\begingroup$ I'll also add what I have so far - its not as "nice" as a logical test and nowhere near as pretty, but the output is the important part - I'm finishing it right now. $\endgroup$
    – ATAX
    Dec 10, 2014 at 3:20
  • $\begingroup$ The graphs show the total number of paths at time t that have positive state. So, if we look at the last graph there are 50 paths with positive state until about t = 16.5 or so. At around t = 17.5 the number of paths with positive state is around 20. $\endgroup$
    – Andy Ross
    Dec 10, 2014 at 3:25
  • $\begingroup$ Ah, I see now! Thank you very much! I do like your way better than mine, and I am assuming I can use the Epilog (or similar ListPlot manipulations) to get numerical answers from the graph if I don't wish to graphically estimate? I will post what I have as an answer to my own question if you were curious. $\endgroup$
    – ATAX
    Dec 10, 2014 at 3:33

Sorry for the horrible coding practices, but I couldn't get it to work otherwise.

First I initialize a list:

LifetimeList1 = {}

I then use a For loop to append the total lifetime of each "path" to LifetimeList1

For[i = 1, i <= Length[data1["Paths"]], i++, 
      Table[Last[data1["Path", i][[j]]], {j, 1, Length[data1["Path", i]]}], # > 0 &]]]];

Now that I have the lifetime of each path in a list, its easier to tell which paths exceed a target lifetime using Select[]

EDIT: Corrected logical error in code


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.