% DEMDDP10 Stochastic Cow Replacement Model
  fprintf('\nDEMDDP10 STOCHASTIC COW REPLACEMENT MODEL\n')
  close all

% Enter model parameters
  delta = 0.9;                  % discount factor  
  cost  = 500;                  % replacement cost  
  price = 150;                  % milk price

% Construct state space
  s = (1:10)';                 % lactation states
  cost = [400;500;600];           % productivity states
  n1 = length(s);
  n2 = length(cost);
  [S,COST] = gridmake(s,cost);    % combined state grid
  n = n1*n2;                    % number of states
  PTRANS=[0.5,0.4,0.1;0.2,0.6,0.2;0.1,0.4,0.5];

% Construct action space
  X = ['K','R'];                % keep or replace

% Construct reward function (actions keep=1, replace=2)
  y = (-0.25*S.^2+2*S+8);  % yield per lactation  
  f = [price*y price*y-COST];   % net revenue by action
  f(S==10,1) = -inf;           % force replace at lactation 10

% Construct state transition probability matrix
  P = zeros(2,n1,n2,n1,n2);
  for i=1:n1
  for j=1:n2
     if i<10
        P(1,i,j,i+1,j) = 1;      % Raise lactation number by 1, if keep
     else
        P(1,i,j,1,1) = PTRANS(j,1);      % Forced replacement after lactation 10
        P(1,i,j,1,2) = PTRANS(j,2);
        P(1,i,j,1,3) = PTRANS(j,3);
     end
     P(2,i,j,1,1) = PTRANS(j,1);         % Optional replacement
     P(2,i,j,1,2) = PTRANS(j,2);
     P(2,i,j,1,3) = PTRANS(j,3);
  end
  end
  P = reshape(P,2,n,n);

% Pach model structure
  clear model
  model.reward     = f;
  model.transprob  = P;
  model.discount   = delta;

% Solve infinite-horizon model using policy iteration
  [v,x,pstar] = ddpsolve(model);

%
