%% scriptExp02.m - Compare DP and Fixed policies
clear all
scriptSetData;
% Out-of-sample scenarios for comparing policies
rng default
numOutSample = 500;
emplOutPaths = SampleEmplPaths(empl0,transMatrix,timeHorizon,numOutSample);
retOutPaths = lognrnd(mu, sigma, numOutSample, timeHorizon) - 1;
% use log utility
% utilFun = @(x) log(x);
% use power utility
delta = 3;
utilFun = @(x) x.^(1-delta)/(1-delta);

%% Learn and evaluate fixed policy
% In sample scenarios for learning fixed policy
numInSample = 1000;
emplInPaths = SampleEmplPaths(empl0,transMatrix,timeHorizon,numInSample);
retInPaths = lognrnd(mu, sigma, numInSample, timeHorizon) - 1;
% Learn fixed policy
[fractCons, fractRisky] = FindFixedPolicy(W0, empl0, utilFun, gamma, ...
        incomeValues, riskFree, retInPaths, emplInPaths);
% Run fixed policy
[aveUtilFix, utilValsFix, consPathsFix, wealthPathsFix] = ...
    RunFixedPolicy([fractCons, fractRisky], W0, empl0, utilFun, gamma, ...
        incomeValues, riskFree, retOutPaths, emplOutPaths);

%% Learn and evaluate DP policy
% set up discretizations for DP
numScenarios = 20;
[values, retProbs] = MakeScenariosQ(numScenarios,mu,sigma);
retScenarios = values - 1; 
Wmin = 1; Wmax = 200; numPoints = 30;
wealthValues = [linspace(Wmin,Wmax,numPoints), 300, 500];
wealthGrid = repmat({wealthValues(:)}, timeHorizon, 1);
% Learn DP policy
splineList = FindDPPolicy(utilFun, gamma, wealthGrid, incomeValues, ...
    transMatrix, retScenarios, retProbs, riskFree);
% Run DP policy
[utilValsDP, alphaPathsDP, consPathsDP, wealthPathsDP] = ...
    RunDPPolicy(splineList, W0, empl0, utilFun, gamma, ...
    incomeValues, transMatrix, retScenarios, retProbs, riskFree, ...
    retOutPaths, emplOutPaths);
aveUtilDP = mean(utilValsDP);
% Compare results 
aveUtilDP
aveUtilFix
%%
figure(1)
hold on
plot(mean(consPathsDP), 'k', 'LineWidth', 1);
plot(mean(consPathsFix), 'k--', 'LineWidth', 1);
title('average consumption path')
hold off

%%
figure(2)
hold on
plot(mean(wealthPathsDP), 'k', 'LineWidth', 1);
plot(mean(wealthPathsFix), 'k--', 'LineWidth', 1);
title('average wealth path')
hold off

%%
figure(3)
hold on
plot(mean(alphaPathsDP), 'k', 'LineWidth', 1);
yline(fractRisky, 'k--', 'LineWidth', 1);
title('average alpha path')
hold off