% Apply the optimal DP consumption-saving policy
function [utilVals, alphaPaths, consPaths, wealthPaths] = ...
    RunDPPolicy(valFunsList, W0, empl0, utilFun, gamma, ...
    incomeValues, transMatrix, inRetScenarios, inRetProbs, riskFree, ...
    outRetScenarios, outEmplScenarios)
[numOutPaths, timeHorizon] = size(outRetScenarios);
numEmplStates = length(incomeValues);
numInScenarios = length(inRetScenarios);
alphaPaths = zeros(numOutPaths, timeHorizon);
consPaths = zeros(numOutPaths, timeHorizon+1);
wealthPaths = zeros(numOutPaths, timeHorizon+1);
utilVals = zeros(numOutPaths, 1);
options = optimoptions(@fmincon,'Display','off');
% start simulation scenarios
for k = 1: numOutPaths
    fprintf(1, 'Working on sample path %d of %d\n', k, numOutPaths);
    finWealth = W0;
    emplNow = empl0;
    wealthNow = finWealth + incomeValues(emplNow);
    % Now carry out forward simulation
    for t = 0:(timeHorizon-1)
        wealthPaths(k,t+1) = wealthNow;
        x0 = [0.5*wealthNow; 0.5]; % starting point for fmincon
        xOut = fmincon(@(x) -objfun(x, wealthNow, emplNow),x0,[],[],[],...
            [],zeros(2,1),[wealthNow;1],[],options);
        % collect path info
        cons = xOut(1); alpha = xOut(2);
        consPaths(k, t+1) = cons; alphaPaths(k,t+1) = alpha;
        utilVals(k) = utilVals(k) + gamma^(t)*feval(utilFun, cons);
        wealthNow = incomeValues(outEmplScenarios(k, t+1))+...
            (wealthNow-cons)*...
            (1+riskFree+alpha*(outRetScenarios(k,t+1)-riskFree));        
    end % for time
    % consume all available wealth at time horizon
    wealthPaths(k,timeHorizon+1) = wealthNow;
    consPaths(k,timeHorizon+1) = wealthNow;
    utilVals(k) = utilVals(k)+gamma^(timeHorizon)*feval(utilFun,wealthNow);
end % for k in simulated scenarios
% nested subfunction
    function outVal = objfun(x, totalWealth, emplState)
        C = x(1);
        alpha = x(2);
        outVal = 0;
        for hh=1:numInScenarios
            futWealth = (totalWealth-C) * ...
                (1+riskFree+alpha*(inRetScenarios(hh)-riskFree));
            for jj=1:numEmplStates
                outVal = outVal + inRetProbs(hh)*transMatrix(emplState,jj) * ...
                    ppval(valFunsList{t+1,jj}, futWealth);
            end % for jj
        end % for hh
        outVal = feval(utilFun,C)+gamma*outVal;
    end % nested function
end




