Solving a dynamic program via recursion

Question

I'm trying to use recursion to solve a joint inventory/ dynamic pricing problem as in monahan, petruzzi and zhao 2004.

I tried to solve for y[t] and k[t] via recursion with the following code:

m = 0.5;
λ = 0.1;
k[0] = 0;
y[0] = 0;
Clear[y];
Clear[k];
r[t_Integer, z_] := 
  1/z^m (-(1/λ)*E^(-λ*z) + 1/λ + k[t - 1]*λ*Integrate[(z - w)^m*E^(-λ*w), {w, 0, z}]) // N;
y[t_Integer] := y[t] = ArgMax[{r[t, j], j > y[t - 1]}, j] // N;
k[t_Integer] := k[t] = r[t, y[t]];
y[1]

I get the following message:

$RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>

When I solve the problem manually, it works fine:

b = 2;
m = 1 - 1/b;
λ = 1/187;
r1[z_] := 1/z^m (-(1/λ)*E^(-λ*z) + 1/λ ) // N;
y1 = ArgMax[{r1[g], g >= 0}, g];
y1
r1[y1]

234.953
8.72688

r2[z_] := 
  N[(1/z^m)*(-(1/λ)/E^(λ*z) + 1/λ + r1[y1]*λ*Integrate[(z - w)^m/E^(λ*w), {w, 0, z}])];
y2 = ArgMax[{r2[h], h >= y1}, h];
y2
r2[y2]

486.281
14.432

r3[z_] := 
  N[(1/z^m)*(-(1/λ)/E^(λ*z) + 1/λ + r2[y2]*λ*Integrate[(z - w)^m/E^(λ*w), {w, 0, z}])];
y3 = ArgMax[{r3[i], i >= y2}, i];
y3
r3[y3]

687.782
18.982

...

Could anybody tell me where's my mistake?

This is the updated code, after editing the suggestions:

m = 0.5;
\[Lambda] = 1/187;
r[0, z_] := 0;
Clear[k, y];
k[0] = 0;
y[0] = 0;
r[t_Integer, z_] := 
  1/z^m (-(1/\[Lambda])  Exp[-\[Lambda] z] + 1/\[Lambda] + 
      k[t - 1] \[Lambda] Integrate[(z - w)^m Exp[-\[Lambda] z], {w, 0,
          z}]) // N;
y[t_Integer] := y[t] = ArgMax[{r[t, j], j > y[t - 1]}, j] // N;
k[t_Integer] := k[t] = r[t, y[t]];
y[1]
k[1]
y[3]
k[3]

234.953

8.72688

234.953

11.3039

Answer 1

I modified your code to eliminate some numerically induced imaginary fuzz.

ClearAll[y, k];

m = 0.5;
λ = 1/187.;
k[0] = 0;
y[0] = 0;
r[t_Integer, z_] := 
  Chop @ N[1/z^m (-(1/λ)*E^(-λ*z) + 1/λ + k[t - 1]*λ*
    Integrate[(z - w)^m*E^(-λ*w), {w, 0, z}])];
y[t_Integer] := y[t] = Chop @ N @ ArgMax[{Re@r[t, j], j > y[t - 1]}, j];
k[t_Integer] := k[t] = Chop @ N @ r[t, y[t]];

I was able to get

y[1]

234.953

y[2]

486.281

y[3] 

687.782

So, perhaps my modifications are what you are looking for.

Answer 2

Thanks again for the help. Here's the whole commented code I wrote to produce my desired output in a nice form, after updating the parameters. It also simulates random demand und uses another recursion to calculate the inventory level and the optimal prices for the given inventory level:

(*Recursive Formulation of Joint Inventory/ Dynamic-Pricing Problem \
according to Monahan, Petruzzi and Zhao 2004 for price-dependent \
stationary periodic demand Subscript[D, t](p,A)=Subscript[Ap, t]^-b, \
where A is exponentially distributed and the decision maker has a one \
time ordering decision at the beginning of the horizon*)

ClearAll[y, k, q, d];

(*The parameters between begin and end can be changed*)

(*begin*)
(*constant price-elasticity*)
b = 2;

(*Intensity of A*)
\[Lambda] = 1/187;

(*No. of periods*)
T = 15;

(*Cost per unit*)
c = 2;
(*end*)

m = 1 - 1/b;

k[0] = 0;
y[0] = 0;

r[t_Integer, z_] := 
  Chop@N[1/z^m (-(1/\[Lambda])*Exp[-\[Lambda] z] + 1/\[Lambda] + 
       k[t - 1] \[Lambda] Integrate[(z - w)^m Exp[-\[Lambda] w], {w, 
          0, z}])];
y[t_Integer] := y[t] = Chop@N@ArgMax[{Re@r[t, j], j > y[t - 1]}, j];
k[t_Integer] := k[t] = Chop@N@Re@r[t, y[t]];

(*Optimal Inventory to stock for T periods*)
S = ((k[T] m)/c)^b;

(*Price depending on Inventory in stock at beginning of period*)
p[t_Integer] := If[q[t] == 0, 0, (y[t]/q[t])^(1 - m)];

(*Random price-dependent demand*)
d[t_] := d[t] = 
   If[t == 0, 0, -(1/\[Lambda]) Log[1 - RandomReal[]] p[t]^-b];

(*quantity in stock*)
q[T] = S;
q[t_Integer] := 
  q[t] = If[q[t + 1] - d[t + 1] > 0, q[t + 1] - d[t + 1], 0];

(*Optimal Revenue-to-go Function*)
R[t_] := k[t] q[t]^m;

v = Table[{t, y[t], k[t], R[t], q[t], d[t], p[t]}, {t, 0, T}];
Tablevh = 
  Prepend[v, {"t", 
    "\!\(\*SuperscriptBox[SubscriptBox[\(z\), \(t\)], \(*\)]\)", 
    "\!\(\*SuperscriptBox[SubscriptBox[\(r\), \(t\)], \(*\)]\)", 
    "\!\(\*SubscriptBox[\(R\), \(t\)]\)(\!\(\*SubscriptBox[\(I\), \(t\
\)]\))", "\!\(\*SubscriptBox[\(I\), \(t\)]\)", 
    "\!\(\*SubscriptBox[\(D\), \(t\)]\)(p,W)", 
    "\!\(\*SuperscriptBox[SubscriptBox[\(p\), \(t\)], \(*\)]\)"}];
Print@Grid[Tablevh, Frame -> All];
Print[StringForm["The maximum expected profit is ``.", 
   AccountingForm[R[T] - c q[T]^m, 
    IntegerLength@Round[R[T] - c q[T]^m] + 2]]];