Dynamic Discrete Programming Using the Compecon Toolbox
Don Ferguson March 14, 2005
1
The Solver: ddpsolve()
The Compecon Toolbox provides tools for solving nite and innite horizon discrete dynamic programming problems, both stochastic and deterministic. The default algorithms are backward recursion for nite horizon problems and policy iteration for innite horizon problems. Value function iteration can also be selected for innite horizon problems using optset. The Solver The solver is the ddpsolve() function [v,x,pstar]=ddpsolve(model,vinit) where v x pstar model vinit the optimal value function, an n-vector the optimal policy function, an n-vector the transition matrix for the optimal policy, an nxn matrix a structure variable that species the properties of the model an optional variable, the initial value function, an n-vector
The Model Structure The model structure variable has elds discount reward transfunc transprob horizon vterm , the discount factor, a scalar on (0,1] f , the reward function, an n m matrix the transition function for deterministic problems, an n m matrix of values on {1, . . . , m} P the stochastic transition matrices, an m n n array (a set of m n n probability matrices) T , the number of time periods in a nite horizon problem (omit or set to inf for innite horizon problems) vT +1 , the terminal value function, an n 1 vector, the default is a null vector
The Options The optset() function can be used to select the following options tol maxit prtiters algorithm the convergence tolerance the maximum number of iterations 0/1 print each iteration newton (policy iteration) or funcit (value iteration), for innite horizon problems, default is newton 1
2
An Example
Consider the example used in class in which you are concerned with the management of a forest subject to infestation. The degree of infestation describes the state of the forest (1 low, 2 medium, 3 high). The actions are either to spray or not (1 spray, 2 no spray). It is assumed that a xed number of trees are harvested each period, but that the amount of usable lumber varies inversely with the degree of infestation. The reward depends on the amount of usable lumber that is harvested and on the cost of spraying. The state transitions are affected not only by spraying, but also by weather which can affect the spread of the infestation (cold winter leading to winter-kill, warm and wet summer weather leading to rapid reproduction, etc.). The following program illustrates how such a model could be implemented.
%ddpsolvetest.m -- Implements the example used in class involving a forest %subject to infestation. clear all; %The transition matrix with spraying. P1=[.9 .1 0; .7 .2 .1; .3 .4 .3]; %The transition matrix without spraying. P2=[.6 .3 .1;.3 .6 .1; .2 .3 .5]; %Construct the mxnxn P array. P(1,:,:)=P1; P(2,:,:)=P2; %The reward nxm function. %f=[4 6;2 4;-1 1]; %if infinite horizon stochastic, then never spray %f=[5 6;3 4; 0 1] %if infinite horizon stochastic, then always spray f=[4 6;3 4; 0 1] %if infinite horizon stochastic, spray only if heavily infested %The discount factor. delta=.9;
%%%PART ONE: INFINITE HORIZON STOCHASTIC MODEL %Construct the model for an infinite horizon model. model.discount=delta; model.transprob=P; model.reward=f; %Call the solver. [v,x,pstar]=ddpsolve(model); %Report the results. disp(the optimal value function)
2
disp(v) disp( ) disp(the optimal policy) disp(x) disp( ) disp(the realized transition probabilities) disp(pstar) pause();
%%%PART TWO: FINITE HORIZON STOCHASTIC MODEL %The same model but extended to a case in which the agent has a fixed term %lease of T years and pays a variable penalty depending on the degree of %infestation. clear model; %clear the previous model from memory %The terminal values. vterm=[0; -20; -40]; %The horizon. T=20; %Construct the infinite horizon model. The first three fiels are the same %as for the infinite horizon model. model.discount=delta; model.transprob=P; model.reward=f; model.horizon=T; model.vterm=vterm; %Call the solver. [v,x,pstar]=ddpsolve(model); %Report the results. disp(the optimal value function) disp(v) disp( ) disp(the optimal policy) disp(x) disp( ) disp(the realized transition probabilities) disp(pstar) pause();
%%%PART THREE: INFINITE HORIZON DETERMINISTIC MODEL
3
%The mxnxn transition probability array is replaced by an nxm transition %function. clear model; %clear the previous model from memory %The transition function. transfunc=[1 2;1 3;2 3]; %Construct the infinite horizon deterministic model. model.discount=delta; model.transfunc=transfunc; model.reward=f; %Call the solver. [v,x]=ddpsolve(model); %Report the results. disp(the optimal value function) disp(v) disp( ) disp(the optimal policy) disp(x)
3
Two Utility Functions: maxval() and valpol()
The function [v,x] = valmax(v,f,P,delta)
returns the value function and policy functions found as the solution to V (s) = max{f (s, x) +
x s
P (s |s, x)V (s )}, P (s |s, x)V (s )},
s
s = 1, . . . , n s = 1, . . . , n.
x(s) = argmax{f (s, x) + The function [pstar,fstar]=valpol(x,f,P)
returns the n n transition matrix pstar and the n 1 reward vector fstar that result from the adoption of the policy x. of Both these functions are used by ddpsolve(). They can be used independently to explore the consequences of adopting any policy, not just the optimal policy. 4
4
Monte Carlo Simulation: ddpsimul() and markov()
In models in which the state transitions are probabilistic the consequences of adopting the optimal policy will differ depending on the initial state and on which state transitions are realized. The Compecon Toolbox contains a function spath = ddpsimul(pstar,sinit,N) which will draw k samples of length N + 1 (the number or periods) from a stochastic process in which the state transitions are governed by the n n matrix pstar. The variable sinit is a k 1 vector of the indices of the initial states for each of the k samples. The output variable spath is a k N + 1 matrix in which row i contains the indices visited during sample i. The Compecon Toolbox also contains a function pi = markov(pstar) which computes the long run probabilities of visting each of the states (assuming that such probabilities exist) where pstar is the n n matrix of transition probabilities and pi is the n 1 vector of probabilities. The following program illustrates the use of these functions.
%testddpsimul.m -- Implements the example used in class involving a forest %subject to infestation and illustrates the use of ddpsimul() and markov(). clear all; %The transition matrix with spraying. P1=[.9 .1 0; .7 .2 .1; .3 .4 .3]; %The transition matrix without spraying. P2=[.6 .3 .1;.3 .6 .1; .2 .3 .5]; %Construct the mxnxn P array. P(1,:,:)=P1; P(2,:,:)=P2; %The reward nxm function. %f=[4 6;2 4;-1 1]; %if infinite horizon stochastic, then never spray %f=[5 6;3 4; 0 1] %if infinite horizon stochastic, then always spray f=[4 6;3 4; 0 1] %if infinite horizon stochastic, spray only if heavily infested %The discount factor. delta=.9; %Construct the model for an infinite horizon model. model.discount=delta;
5
model.transprob=P; model.reward=f; %Call the solver. [v,x,pstar]=ddpsolve(model); %Report the results. disp(the optimal value function) disp(v) disp( ) disp(the optimal policy) disp(x) disp( ) disp(the realized transition probabilities) disp(pstar) pause(); %Find the long run probabilities. pi=markov(pstar); disp(The long run probabilities are.) disp(pi) pause(); %Initialize the simulation. N=10; numsim=30; sinit=zeros(numsim,1); sinit(1:10,1)=1; sinit(11:20,1)=2; sinit(21:30,1)=3; %Perform the simulation. spath=ddpsimul(pstar,sinit,N); %Report the results of the simulation. disp(Starting from state 1 the sample paths for the states are) disp(spath(1:10,:)) disp( ) pause(); disp(Starting from state 2 the sample paths for the states are) disp(spath(11:20,:)) disp( ) pause(); disp(Starting from state 3 the sample paths for the states are) disp(spath(21:30,:))
6
5
Multiple States and Actions: gridmake() and getindex()
Suppose the you have two state variables s1 and s2 which can take on n1 and n2 discrete values respectively. This gives rise to a total of n = n1 n2 possible combinations, each of which can be viewed as a distinct state. We can index these states by converting each combination into a distinct index in the interval 1, . . . , n. To do this, the Compecon Toolbox contains a function S=gridmake(array1,array2) which translates combinations of the entries in two one dimensional arrays into an array of their grid positions and allows each combination to be represented by the index of their position in S. The function i=getindex(combination) then returns the position of a particular combination in the grid. Additional information and examples are contained in the following help les for the two functions. gridmake
GRIDMAKE Forms grid points USAGE X=gridmake(x); X=gridmake(x1,x2,x3,...); [X1,X2,...]=gridmake(x1,x2,x3,...); X=gridmake({y11,y12},x2,{y21,y22,y23});
Expands matrices into the associated grid points. If N is the dx2 matrix that indexes the size of the inputs GRIDMAKE returns a prod(N(:,1)) by sum(N(:,2)) matrix. The output can also be returned as either d matrices or sum(N(:,2)) matices If any of the inputs are grids, they are expanded internally Thus X=gridmake({x1,x2,x3}) X=gridmake(x1,x2,x3) and x={x1,x2,x3}; X=gridmake(x{:}) all produce the same output. Note: the grid is expanded so the first variable change most quickly. Example: X=gridmake([1;2;3],[4;5]) produces
7
1 2 3 1 2 3
4 4 4 5 5 5
The function performs an action similar to NDGRID, the main difference is in the increased flexability in specifying the form of the inputs and outputs. Also the inputs need not be vectors. X=gridmake([1;2;3],[4 6;5 7]) produces 1 4 6 2 4 6 3 4 6 1 5 7 2 5 7 3 5 7 See also: ndgrid
getindex
GETINDEX Finds the index value of a point USAGE i = getindex(s,S); INPUTS s : a 1xd vector or pxd matrix S : an nxd matrix OUTPUT i : a p-vector of integers in {1,...,n} indicating the row of S that most closely matches each row in s Example: S=[0 0; 0 1; 1 0; 1 1]; s=[0 0; 0 0; 1 0; 0 0; 1 1]; getindex(s, S) returns [1; 1; 3; 1; 4] Coded as a MEX file in C.
5.1
Example
See the asset replacement model with maintenance that is discussed in Sections 7.2.3 and 7.6.3. The model is implemented in demddp03.m which you can nd in the CEdemos directory.
8