Monte Carlo simulation of a financial market

Question

I'm fairly new to Mathematica and i want to review the average prices of this model. In order to do that i need to perform a Monte Carlo simulation on my financial markets model:

P[t_] := P[t] = P[t - 1] + μ*ED[t - 1] + RandomVariate[NormalDistribution[0, 0.03]];
ED[t_] := ED[t] = WC[t]*DC[t] + WF[t]*DF[t];
DC[t_] := DC[t] = b (P[t] - P[t - 1]);
DF[t_] := DF[t] = c (F - P[t]);
WF[t_] := WF[t] = 1/(1 + Exp[-β*a[t]]);
WC[t_] := WC[t] = 1 - WF[t];
a[t_] := a[t] = α0 + αn*(WF[t - 1] - WC[t - 1]) + αp (P[t - 1] - F)^2;

μ = 0.01;
β = 1;
b = 0.01;
c = 0.01;
α0 = 0;
αn = -0.4;
αp = 1;

F = 1;
P[0] = F;
P[1] = F;
P[2] = F;
P[3] = F + 0.01;

WF[1] = 0.2;
WC[1] = 0.8;

logprices = Table[P[t], {t, 1000, 6000}];
logreturns = Differences[logprices];
abslogreturns = Abs[logreturns];
Dis[t_] := Dis[t] = Abs[P[t] - F];
SumDis = Table[Dis[t], {t, 1000, 6000}];
Mean[SumDis]
Mean[logprices]
Kurtosis[logreturns]

P = Price, ED = excess demand, DC/DF = demand, WC/WF = weights, a = influences of weights, F = fundamental value, dis = distortion

A single run is supposed to yield the average price, the average distortion and the kurtosis of the logreturns. P, DC/DF and WF/WC will be plotted using ListLineplot.

I want to run this about 2000-5000 times and compile the resulting 2000-5000 average prices and distortions in a list or something along those lines so i can average them to see what a change of a variable does to those.

I read the Mathematica documentation on how to perform a Monte Carlo simulation, but I couldn't see how to transfer it to my model.

I tried the table function but the only result I achieved is {Null}.

My question is: how do I set this up correctly?

Answer 1

You use memoization in your code. In some cases that may be quite useful. Personally, I use it rarely, because it mixes data and operations on that data in a rather unfavorable way: After a single run of your simulation, you have to delete almost all you definitions. So the life cycle of the memoized data is pretty short; in the end it will be tiny fractions of a second. Moreover, you plan to produce huge loads of numerical data; these can be most efficiently stored and accessed in a simple array. Since you want to perform lots of simulations with many different parameters, you might also want decent performance of the code. Thus, I compiled your code into the following library function:

cf = Compile[{{noise, _Real, 1}, {F, _Real}, {b, _Real}, {c, _Real},
   {α0, _Real}, {αn, _Real}, {αp, _Real}, {β, _Real}, {μ, _Real}, {WF1, _Real}, {WC1, _Real}},
  Block[{P, ED, DC, DF, WF, WC, a, n},
   n = Length[noise];
   P = Table[0., {n}];
   ED = Table[0., {n}];
   DC = Table[0., {n}];
   DF = Table[0., {n}];
   WF = Table[0., {n}];
   WC = Table[0., {n}];
   a = Table[0., {n}];
   P[[1 ;; 3]] = {F, F, F + 0.01};
   WF[[1]] = WF1;
   WC[[1]] = WC1;
   Do[
    P[[i]] = P[[i - 1]] + μ ED[[i - 1]] + noise[[i]];
    DC[[i]] = b (P[[i]] - P[[i - 1]]);
    DF[[i]] = c (F - P[[i]]);
    a[[i]] = α0 + αn (WF[[i - 1]] - WC[[i - 1]]) + αp (P[[i - 1]] - F)^2;
    WF[[i]] = 1./(1. + Exp[-β  a[[i]]]);
    WC[[i]] = 1. - WF[[i]];
    ED[[i]] = WC[[i]]  DC[[i]] + WF[[i]] DF[[i]];
    , {i, 2, n}];
   {P, DC, DF, a, WF, WC, ED}
   ],
  CompilationTarget -> "C",
  RuntimeAttributes -> {Listable},
  Parallelization -> True,
  RuntimeOptions -> "Speed"
  ]

Now you can run 5000 simulations with 6000 iterations each within a second with

n = 6000;
m = 5000;
μ = 0.01;
β = 1;
b = 0.01;
c = 0.01;
α0 = 0;
αn = -0.4;
αp = 1;
F = 1.;
WF1 = 0.2;
WC1 = 0.8;
noise = RandomVariate[NormalDistribution[0, 0.03], {m, n}];
data = cf[noise, F, b, c, α0, αn, αp, β, μ, WF1, WC1];

In order to interpret the output: data has dimension {m,7,n}. You access {P, DC, DF, a, WF, WC, ED} (in that order) of the j-th iteration in the i-th simulation with data[[i,All,j]].