Too slow computation in mathematica (while R runs the same code very fast)

Question

I have made the following code for simulating an evolution of a set of variables in a discreet-time framework. But the mathematica computation is way too slow. I did the same simulation using R and the result comes out right away. But with mathematica it takes forever. I need to simulate the model for at least 100 periods, and there is no problem with R. But with mathematica, I cannot go over time 12. From time 13, it takes forever. Is there any way I can run the simulation with mathematica as fast as R? The mathematica code looks like the following:

Time period:

period = Range[3, 13]

The paramter value:

p = 0.7
q = 0.2
k = 0.3
id = 0.05
il = 0.1
TP = 3
TP = 3
TC = 1
TF = 3

Value for some initial periods:

Q[0] = Q[1] = Q[2] = 1
X[0] = X[1] = X[2] = 1
FK[0] = FK[1] = FK[2] = 1
FH[0] = FH[1] = FH[2] = 1
S[0] = S[1] = S[2] = 0.5
C[0] = C[1] = C[2] = 0.5
P[0] = P[1] = P[2] = 0.5
BK[0] = BK[1] = BK[2] = 0.1
BH[0] = BH[1] = BH[2] = 0.1
LK[0] = LK[1] = LK[2] = 1
LH[0] = LH[1] = LH[2] = 1
r[0] = r[1] = r[2] = 0.1

A System Equations (or functions):

C[t_] := ((1 + p*q)/(1 + q))*S[t - TF] + p*id*FK[t - TF] - p*id*LK[t - TF] + BK[t]

S[t_] := (1 - k)*C[t] + k*C[t - TC] + (((1 - p)*q)/(1 + q))*S[t - TC] + id*(1 - p)*FK[t - TC] - il*(1 - p)*LK[t - TC] + id*FH[t - TC] - il*LH[t - TC] + BH[t]

FK[t_] := FK[t - 1] + S[t - 1]/(1 + q) + p*((q/(1 + q))*S[t - 1] + id*FK[t - 1] - il*LK[t - 1]) - C[t - 1] + BK[t - 1]

FH[t_] := FH[t - 1] + k*C[t - 1] + (1 - p)*((q/(1 + q))*S[t - 1] + id*FK[t - 1] - il*LK[t - 1]) + id*FH[t - 1] - il*LH[t - 1] - (k*C[t - 1 - TC] + (1 - p)*((q/(1 + q))*S[t - 1 - TC] + id*FK[t - 1 - TC] - il*LK[t - 1 - TC]) + id*FH[t - 1 - TC] - il*LH[t - 1 - TC]) + BH[t - 1]

P[t_] := C[t - TP]

Q[t_] := Q[t - 1] + C[t - 1] - P[t - 1]

X[t_] := X[t - 1] + P[t - 1] - 1/(1 + q)*S[t - 1]

gBK[t_] := 0.03 + 0.5*r[t] - 0.2*il

gBH[t_] := 0.05 - 0.05*q - 0.01*il

BK[t_] := (1 + gBK[t - 1])*BK[t - 1]

BH[t_] := (1 + gBH[t - 1])*BH[t - 1]

LK[t_] := LK[t - 1] + BK[t - 1]

LH[t_] := LH[t - 1] + BH[t - 1]

r[t_] := (q/(1 + q)*S[t])/(Q[t] + X[t])

Then I run the following:

ListPlot[Map[C, period]]

And I never get the result once the period goes beyond 12. Thank you!

Answer 1

Two observations about your code:

1) Don't start your symbols' names with uppercase. The convention is that only System defined symbols start with uppercase (and you fell into the trap with C[n])
2) Learn about memoization. It's needed for recursive definitions going faster (that is your case).

Here you have your (modified) code running instantaneously:

p = 0.7;
q = 0.2;
k = 0.3;
id = 0.05;
il = 0.1;
TP = 3;
TP = 3;
TC = 1;
TF = 3;
Q[0] = Q[1] = Q[2] = 1;
X[0] = X[1] = X[2] = 1;
FK[0] = FK[1] = FK[2] = 1;
FH[0] = FH[1] = FH[2] = 1;
S[0] = S[1] = S[2] = 0.5;
c[0] = c[1] = c[2] = 0.5;
P[0] = P[1] = P[2] = 0.5;
BK[0] = BK[1] = BK[2] = 0.1;
BH[0] = BH[1] = BH[2] = 0.1;
LK[0] = LK[1] = LK[2] = 1;
LH[0] = LH[1] = LH[2] = 1;
r[0] = r[1] = r[2] = 0.1;
c[t_] := c[t] = ((1 + p*q)/(1 + q))*S[t - TF] + p*id*FK[t - TF] - 
   p*id*LK[t - TF] + BK[t]
S[t_] := S[t] = (1 - k)*c[t] + 
   k*c[t - TC] + (((1 - p)*q)/(1 + q))*S[t - TC] + 
   id*(1 - p)*FK[t - TC] - il*(1 - p)*LK[t - TC] + id*FH[t - TC] - 
   il*LH[t - TC] + BH[t]
 FK[t_] := FK[t] = FK[t - 1] + S[t - 1]/(1 + q) + 
   p*((q/(1 + q))*S[t - 1] + id*FK[t - 1] - il*LK[t - 1]) - c[t - 1] + BK[t - 1]
 FH[t_] :=  FH[t] = FH[t - 1] + 
   k*c[t - 1] + (1 - p)*((q/(1 + q))*S[t - 1] + id*FK[t - 1] - 
      il*LK[t - 1]) + id*FH[t - 1] - il*LH[t - 1] - (k*
      c[t - 1 - TC] + (1 - p)*((q/(1 + q))*S[t - 1 - TC] + 
        id*FK[t - 1 - TC] - il*LK[t - 1 - TC]) + id*FH[t - 1 - TC] - 
     il*LH[t - 1 - TC]) + BH[t - 1]
P[t_] := P[t] = c[t - TP]
Q[t_] := Q[t] = Q[t - 1] + c[t - 1] - P[t - 1]
X[t_] := X[t] = X[t - 1] + P[t - 1] - 1/(1 + q)*S[t - 1]
gBK[t_] := 0.03 + 0.5*r[t] - 0.2*il
gBH[t_] := 0.05 - 0.05*q - 0.01*il
BK[t_] := BK[t] = (1 + gBK[t - 1])*BK[t - 1]
BH[t_] := BH[t] = (1 + gBH[t - 1])*BH[t - 1]
LK[t_] := LK[t] = LK[t - 1] + BK[t - 1]
LH[t_] := LH[t] = LH[t - 1] + BH[t - 1]
 r[t_] := (q/(1 + q)*S[t])/(Q[t] + X[t])

Then

ListPlot[c /@ Range[3, 60]]