Time

Hal Varian
April 1992

Symbolic dynamic programming

Dynamic programming is a powerful way to solve optimization problems.  We consider the problem described in the text here, but in order to keep things simple we will ignore the uncertainty.   Recall that the one-period objective function is Log[c] and the discount factor is alpha.  We use w to denote wealth and R is the total return on savings.

Two-period case

U[c_] := Log[c] + alpha Log[(w-c)R]
solution1=Solve[D[U[c],c]==0,c][[1]]

Substituting back into U[c] gives us the indirect utility function.  We clear c after we're done so it won't interfere with the calculations further on.

Simplify[U[c]/.solution1]
Clear[c]

Calculating the value function by recursion

Let's try to program Mathematica to do the dynamic programming recursion all by itself. As it turns out the only thing that we have to do is to make the consumption variable, c,  "local" to the routine where it is being calculated.   Otherwise, Mathematica gets confused by all the values of c that are floating around.  To do this, we use the Module construction. Here's how it works for a 3-period problem.

V3[w_,t_]:=
 Module[{c},
 Log[c] + alpha*V3[(w-c)*R,t+1]/.Solve[D[Log[c] 
        + alpha*V3[(w-c)*R,t+1],c]==0,c][[1]]]

V3[w_,3] := Log[w]

TableForm[{V3[w,3],V3[w,2],Simplify[V3[w,1]]}]

Solving for consumption

Once we have the value function, it is easy to solve for consumption in a given period. 
First let's check that we get the same values as before one period before the end.

V3[w,2]
Solve[D[Log[c] + alpha*V3[(w-c)*r,3],c]==0,c][[1]]

Yep, that checks.  Here's first period consumption in a 3-period problem.

Solve[D[Log[c] + alpha*V3[(w-c)*r,2],c]==0,c][[1]]

To solve for first-period consumption in, say, a 5-period problem, we have to go back and define a new value function for that problem.

V5[w_,t_]:= 
 Module[{c},
 Log[c] + alpha*V5[(w-c)*R,t+1]/.Solve[D[Log[c] 
        + alpha*V5[(w-c)*R,t+1],c]==0,c][[1]]]

V5[w_,5] := Log[w]

Solve[D[Log[c] + alpha*V5[(w-c)*r,2],c]==0,c][[1]]

Note that consumption is independent of the rate of return, as we argued in the text.

Speeding things up

Actually, we should define the value function a little bit differently if we want to use Mathematica efficiently.  When we calculate V[w,1], we have to calculate all the other values of V[w,t].  When we calculate V[w,2] we have to calculate those values all over again.  We can tell Mathematica to "remember" the values of the V[w,t] function by using the following construction.   (The "f" stands for "fast".)

V5f[w_,t_]:= V5f[w,t] =
 Module[{c},
 Log[c] + alpha*V5f[(w-c)*R,t+1]/.Solve[D[Log[c] 
        + alpha*V5f[(w-c)*R,t+1],c]==0,c][[1]]]

V5f[w_,5] := Log[w]

Timing[V5[w,1];]
Timing[V5[w,2];]

Timing[V5f[w,1];]
Timing[V5f[w,2];]

Note that the second way of calculating the value function is a lot faster.

Stochastic dynamic programming

One of the strengths of dynamic programming is that it is easy to handle stochastic problems.  Consider the savings problem discussed above.  Suppose that R takes value Ru with probability p and Rd with probability 1-p.  Then the expected value function is

V3s[w_,t_]:= V3s[w,t] =
 Module[{c},
 Log[c] + alpha*(p*V3s[(w-c)*Ru,t+1]+
             (1-p)*V3s[(w-c)*Rd,t+1])/.Solve[D[Log[c] 
        + alpha*(p*V3s[(w-c)*Ru,t+1]+
             (1-p)*V3s[(w-c)*Rd,t+1]),c]==0,c][[1]]]
V3s[w_,3] := Log[w]

We solve for optimal consumption in period 1:

Solve[D[Log[c] + alpha*(p*V3s[(w-c)*Ru,2]+
     (1-p)*V3s[(w-c)*Rd,2]),c]==0,c][[1]]

Note that even in the stochastic case, consumption doesn't depend on p.  That is, optimal consumption is independent of the probability distribution of the rate of return.

Numeric dynamic programming

By the numbers

Of course only a few dynamic programming problems have explicit closed-form expressions.  However, we can use exactly the same techniques to solve problems numerically.  First we'll try the problem we've already investigate since we can check the numerical formula against the symbolic one.   The Mathematica function FindMinimum[f[x],{x,x0}] returns a list {fmax,{x->xmax}} where fmax  is a (local) minimum value of f[x] starting at the valued x0.

Here's the expression for the optimal consumption with wealth w at time t in a 3-period problem.

c3n[w_,t_] := c3n[w,t] =
       c/.FindMinimum[-Log[c] - alpha*V3n[(w-c)*R,t+1],{c,1}][[2]]
c3n[w_,3] := w

Here's the expression for the value function:

V3n[w_,t_]:= V3n[w,t] =
        Log[c3n[w,t]] + alpha*V3n[(w-c3n[w,t])*R,t+1]
V3n[w_,3] := N[Log[w]]

Set the parameters so that we can do this numerically.

R=1.5
alpha=.9

Now let's calculate

{V3n[100,1],V3n[100,2],V3n[100,3]}

Compare this to the values found by the symbolic calculation:

{V3[100,1],V3[100,2],N[V3[100,3]]}

Pretty close!  Now let's try the same calculation for a 5-period problem:

c5n[w_,t_] := c5n[w,t] =
       c/.FindMinimum[-Log[c] - alpha*V5n[(w-c)*R,t+1],{c,1}][[2]]
c5n[w_,5] := w
V5n[w_,t_]:= V5n[w,t] =
        Log[c5n[w,t]] + alpha*V5n[(w-c5n[w,t])*R,t+1]
V5n[w_,5] := N[Log[w]]

{V5n[100,3],V5n[100,4],V5n[100,5]}

[I tried to compute V5n[100,2] but I got tired of waiting for it to finish.  I don't think that it should take as long as it does; perhaps there is a bug somewhere...]

Question:  Why are these the same as  V3n[100,1], V3n[100,2] and V3n[100,1]?

Answer

Because each one measures the value of having wealth of 100 0, 1, and 2 periods before the termination date.

Stochastic numeric dynamic programming

Of course you can do the stochastic programming numerically too.  Here's a little example for a 3-period problem.

c3ns[w_,t_] := c3ns[w,t] =
       c/.FindMinimum[-Log[c] - alpha*(p*V3ns[(w-c)*Ru,t+1]
             +(1-p)*V3ns[(w-c)*Rd,t+1]),{c,1}][[2]]
c3ns[w_,3] := w
V3ns[w_,t_]:= V3ns[w,t] =
        Log[c3ns[w,t]] + alpha*(p*V3ns[(w-c3ns[w,t])*Ru,t+1]
                           +(1-p)*V3ns[(w-c3ns[w,t])*Rd,t+1])
V3ns[w_,3] := N[Log[w]]

We compare the symbolic and the numeric calculations.  Of course we have to provide numeric values for all parameters.  (The calculation of V3ns[100,1] is rather slow.)

Ru=1.5;
Rd=.5;
p=.5;
alpha=.9;

{N[V3s[100,3]],V3s[100,2],V3s[100,1]}

{V3ns[100,3],V3ns[100,2],V3ns[100,1]}

This is getting to be a pretty slow calculation.  In practice the value function for a numerical stochastic dynamic programming function is typically approximated by a step function so that interpolation methods can be used to find the value function at some points.  The  best method of interpolation will vary from problem to problem. 

We've only investigated problems with one state variable, wealth.  If there are several state variables the same principles apply, but the calculations become unmanageable if there are more than 4 or 5 state variables.