# Simulating discrete time stochastic dynamic systems

What is the canonical way of simulating discrete time stochastic dynamical systems in Mathematica using the new functionality of Random processes?

To take a concrete example, lets consider the optimal gambling problem. A gambler comes to a casino with an initial fortune $x_1$ and let $X_n$ denote his fortune at time $n$. At each time he bets a fraction $\alpha x_n$. At each time he wins with probability $p$ and loses with probability $(1-p)$. Thus, the dynamics of the system can be written as:

$$X_{n+1} = \begin{cases} (1+\alpha) X_n, & \text{with probability } p \\ (1-\alpha) X_n, & \text{with probability } 1-p \\ \end{cases}$$

I want to simulate this system for 20 time steps with $x_1 = 12$, $p=0.8$, and then plot 50 sample paths and the mean value of $X_{20}$.

-

## 1 Answer

I would proceed like the following. It will be natural to propose that the win-event occurs following BinomialDistribution with probability $p=0.8$ so that we can use the built-in BinomialProcess to simulate the win and losses in $20$ time steps for $50$ sample paths.

timstep = 20;
win = BinomialProcess[.8];
samplepaths=50;
process = RandomFunction[Evaluate[win], {0, timstep - 1}, samplepaths];
ListLinePlot[process, InterpolationOrder -> 0, Frame -> True,
PlotStyle -> {{Dashed, Opacity[.9]}}]


Then just use FoldList on the TemporalData to get your result simulation of $50$ sample paths and the mean in a simple functional style!

x1 = 12;
alpha = .15; (* assume some value for alpha *)
WinEvents = process["States"];
ST[WinEvents_?ListQ, x1_, alpha_, SimTime_] := Module[{},
FoldList[Max[If[#2 == 1, (1 + alpha) #1, (1 - alpha) #1], 0.] &,
x1,
Prepend[Differences@WinEvents, First@WinEvents]]
];
simdat = ST[#, x1, alpha, timstep] & /@ WinEvents;


For visualization

Show[ListLinePlot[Legended[Mean[simdat], "Mean"], PlotRange -> All,
PlotStyle -> {{AbsoluteThickness[3.5], Opacity[.9], Red}}],
ListLinePlot[simdat, PlotRange -> All,
PlotStyle -> {{Dashed, Opacity[.7]}}],
ListPlot[simdat, PlotRange -> All], Frame -> True]


-
Thanks. this is exactly what I was looking for. I did not know about process["States"] and now I can read up more about it. –  Aditya Apr 15 '13 at 9:57
@Aditya You are welcome! Mathematica doc has got lots of hidden gems..wish you a happy treasure hunting ;) –  PlatoManiac Apr 15 '13 at 10:01
An alternative option is to use BernoulliProcess instead of BinomialProcess. Then, one would not need to take the differences in ST function. –  Aditya Apr 15 '13 at 10:11
You are right! I randomly chose a distribution for this example.. –  PlatoManiac Apr 15 '13 at 11:10