Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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
share|improve this question
1  
try assigning k[0],y[0] after you Clear[k,y] .. –  george2079 Mar 6 at 20:27
    
thanks for the reply! I do not get the mistake $RecursionLimit::reclim: Recursion depth of 1024 exceeded. >> and y[1], r[1] are calculated correctly. however in the next iteration y[2]=y[1], which is not the case when manually solving the steps. –  user12806 Mar 6 at 21:14

2 Answers 2

up vote 1 down vote accepted

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.

share|improve this answer
    
yes :-) sweet, thank you. –  user12806 Mar 6 at 21:39

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]]];
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.