function [mt_vf, mt_pf, mt_oploc, mt_o, mt_k, mt_l]=hhsave_VFI_3OC(varargin)

close all;

bl_profile = true;      % Switch off profile if running a tester/calling from another function

bl_saveimg = true;

if(bl_profile)

profile off;

profile on;

end 

addpath(genpath('/Users/sidhantkhanna/Documents/GitHub/BKS modified/')); 

if ~isempty(varargin) 

    [ar_a,ar_z,ar_n, ...

        fl_alpha,fl_theta,fl_delta,fl_kappa,fl_r,fl_w,fl_phi,fl_ahi,fl_zhi, ...

         fl_risk,it_zgridno, it_agridno,it_ngridno, mt_trans_z,mt_trans_n,fl_beta, fl_mu_z, fl_sig_z, ...

        fl_rho_z,fl_lambda_z,fl_mu_n, fl_sig_n,fl_rho_n,fl_lambda_n, fl_tolvfi ...

        ] = varargin{:};

    bl_print  = false;

    bl_plot   = false;

else 

    fl_ahi         = 50; 

    fl_alo         = 0; 

    fl_zhi         = 2.2; 

    it_agridno     = 100; 

    it_zgridno     = 7; 

    it_ngridno     = 2; 

 %  ar_a           = linspace(fl_alo ,fl_ahi,it_agridno);

 % mp_grid_control = containers.Map('KeyType','char', 'ValueType','any');

 %  mp_grid_control('grid_powerspace_power') = 2;

 %   ar_a            = ff_saveborr_grid(fl_alo, fl_ahi, it_agridno,  mp_grid_control)

    ar_a           = ff_saveborr_grid(fl_alo, fl_ahi, it_agridno, 'grid_powerspace'); 

    ar_a           = ar_a'; 

    fl_phi         = 0.5; 

    fl_risk        = 1.5; 

    fl_alpha       = 0.4; 

    fl_theta       = 0.79-fl_alpha; 

    fl_delta       = 0.05; 

    fl_kappa       = 0; 

    fl_mu_z        = 0;                   % mean of AR(1) entrepeneurial productivity process 

    fl_rho_z       = 0.9;                 % persistence parameter of the AR(1) entrepreneurial productivity process 

    fl_sig_z       = 0.02;  

    fl_lambda_z    = 3;  

    fl_beta        = 0.92;  

 %  [ar_z, mt_trans_z] = mytauchen_z(fl_mu_z,fl_rho_z,fl_sig_z,it_zgridno,fl_lambda_z); 

    [ar_z, mt_trans_z] = ffy_rouwenhorst(fl_rho_z,fl_sig_z,it_zgridno);     

    ar_z           = exp(ar_z); 

    P1=mt_trans_z^1000; 

    sd=P1(1,:); 

    el=sd*ar_z; 

    ar_z =ar_z/el; 

    ar_z         = ar_z';    

    ar_n           = [0,1]; 

    mt_trans_n     = [0.1,0.9;0.9,0.1]; 

    [fl_r,fl_w] = ... 

        deal(0.02,1); 

    fl_tolvfi    = 10^-12; 

    bl_print     = true; 

    bl_plot      = true; 

end 

fl_R  = fl_r + fl_delta; 

%% Value function Iteration

mt_vf = zeros(it_agridno,it_zgridno,it_ngridno);             % Value function 

mt_Tv = zeros(it_agridno,it_zgridno,it_ngridno);             % Next guess of value function 

mt_oploc = zeros(it_agridno,it_zgridno,it_ngridno);             % 

[mt_util,mt_coh, mt_k, mt_l, mt_o] = hhsave_util_3OC(ar_a,ar_z,ar_n, ... 

        fl_alpha,fl_theta,fl_delta,fl_kappa,fl_r,fl_w,fl_phi,fl_ahi,fl_zhi, ... 

         fl_risk,it_zgridno, it_agridno,it_ngridno); 

fl_crit = 1;                      % difference between current and next guess of value function 

mt_evf = zeros(it_agridno,it_zgridno,it_ngridno);  

while fl_crit > fl_tolvfi 

for i = 1:it_ngridno      

mt_evf(:,:,i) = mt_vf(:,:,i)*mt_trans_z'; 

end 

for j = 1:it_zgridno      

mt_evf1(:,j,:) = squeeze(mt_evf(:,j,:))*mt_trans_n'; 

end 

mt_ev1 = repmat(mt_evf1, [1 1 1 it_agridno]); 

mt_ev1 = permute(mt_ev1,[4 2 3 1]);     

[mt_Tvf,mt_oploc] = max((mt_util+fl_beta.*(mt_ev1)),[], 4); 

    mt_diff = mt_Tvf-mt_vf; 

    mt_diff = mt_diff(:); 

    fl_crit = norm(mt_diff); 

    mt_vf   = mt_Tvf; 

end 

mt_pf  = ar_a(mt_oploc); 

mt_con = mt_coh - mt_pf;  

mt_savefraccoh =  mt_pf./mt_coh; 

if(bl_print) 

    st_ar_a = strsplit(num2str(ar_a)); 

    for i=1:it_agridno 

        st_a(i) =strcat('a=', st_ar_a(i)); 

    end 

    st_ar_z = strsplit(num2str(ar_z)); 

    disp('Table for value function at l=0'); 

    tb_VF = array2table(mt_vf(:,:,1),'VariableNames',{'z1','z2','z3','z4','z5','z6','z7'}); 

    tb_VF.Properties.RowNames = st_a; 

    disp(tb_VF); 

    disp('Table for value function at l=1'); 

    tb_VF = array2table(mt_vf(:,:,2),'VariableNames',{'z1','z2','z3','z4','z5','z6','z7'}); 

    tb_VF.Properties.RowNames = st_a; 

    disp(tb_VF); 

    disp('Table for Policy function at l=0'); 

    tb_pf = array2table(mt_pf(:,:,1),'VariableNames',{'z1','z2','z3','z4','z5','z6','z7'}); 

    tb_pf.Properties.RowNames = st_a; 

    disp(tb_pf); 

    disp('Table for Policy function at l=1'); 

    tb_pf = array2table(mt_pf(:,:,2),'VariableNames',{'z1','z2','z3','z4','z5','z6','z7'}); 

    tb_pf.Properties.RowNames = st_a; 

    disp(tb_pf); 

    disp('Table for Consumption at l=0 '); 

    tb_con = array2table(mt_con(:,:,1),'VariableNames',{'z1','z2','z3','z4','z5','z6','z7'}); 

    tb_con.Properties.RowNames = st_a; 

    disp(tb_con); 

    disp('Table for Consumption at l=1 '); 

    tb_con = array2table(mt_con(:,:,2),'VariableNames',{'z1','z2','z3','z4','z5','z6','z7'}); 

    tb_con.Properties.RowNames = st_a; 

    disp(tb_con); 

    disp('Table for Savings/COH at l=0'); 

    tb_savefraccoh = array2table(mt_savefraccoh(:,:,1),'VariableNames',{'z1','z2','z3','z4','z5','z6','z7'}); 

    tb_savefraccoh.Properties.RowNames = st_a; 

    disp(tb_savefraccoh); 

    disp('Table for Savings/COH at l=1'); 

    tb_savefraccoh = array2table(mt_savefraccoh(:,:,2),'VariableNames',{'z1','z2','z3','z4','z5','z6','z7'}); 

    tb_savefraccoh.Properties.RowNames = st_a; 

    disp(tb_savefraccoh); 

end 

%% Plot Value function and Policy Function

% 3D Value function - l=0

if(bl_plot) 

figure(1) 

surf(ar_z,ar_a,mt_vf(:,:,1)); 

title('Value Function, l=0'); 

xlabel('Productivity') ; 

ylabel('Assets'); 

zlabel('Value'); 

xt = get(gca, 'XTick');                                 

set(gca, 'XTick', xt, 'XTickLabel', xt*fl_ahi/it_agridno); 

if(bl_saveimg) 

saveas(gcf, '/Users/sidhantkhanna/Documents/GitHub/BKS modified/code/Households/figures/hhsave_VFI_3OC/Valuefunction_3D_0.png') 

end 

figure(2) 

surf(ar_z,ar_a,mt_vf(:,:,2)); 

title('Value Function, l=1'); 

xlabel('Productivity') ; 

ylabel('Assets'); 

zlabel('Value'); 

xt = get(gca, 'XTick');                                 

set(gca, 'XTick', xt, 'XTickLabel', xt*fl_ahi/it_agridno); 

if(bl_saveimg) 

saveas(gcf, '/Users/sidhantkhanna/Documents/GitHub/BKS modified/code/Households/figures/hhsave_VFI_3OC/Valuefunction_3D_1.png') 

end 

% 2D Value function

figure(3) 

%chart = plot(mt_vf);

for i = 1: it_zgridno 

a(i)  = plot((mt_vf(:,i,1)), '.'); 

b{i}  = sprintf(['Z shock =', num2str(ar_z(i))]); 

hold on 

end 

clr = jet(numel(a)); 

for m = 1:numel(a) 

   set(a(m),'Color',clr(m,:)) 

end 

legend(b,'Location','South') 

%clr = jet(numel(chart));

%for m = 1:numel(chart) 

 %  set(chart(m),'Color',clr(m,:))

%end

title('Value Function, l=0'); 

xlabel('Asset') ; 

ylabel('Value'); 

xt = get(gca, 'XTick');                                 

set(gca, 'XTick', xt, 'XTickLabel', xt*fl_ahi/it_agridno); 

if(bl_saveimg) 

saveas(gcf, '/Users/sidhantkhanna/Documents/GitHub/BKS modified/code/Households/figures/hhsave_VFI_3OC/Valuefunction_2D_0.png') 

end 

figure(4) 

%chart = plot(mt_vf);

for i = 1: it_zgridno 

a(i)  = plot((mt_vf(:,i,2)),'.'); 

b{i}  = sprintf(['Z shock =', num2str(ar_z(i))]); 

hold on 

end 

clr = jet(numel(a)); 

for m = 1:numel(a) 

   set(a(m),'Color',clr(m,:)) 

end 

legend(b,'Location','South') 

%clr = jet(numel(chart));

%for m = 1:numel(chart) 

 %  set(chart(m),'Color',clr(m,:))

%end

title('Value Function, l=1'); 

xlabel('Asset') ; 

ylabel('Value'); 

xt = get(gca, 'XTick');                                 

set(gca, 'XTick', xt, 'XTickLabel', xt*fl_ahi/it_agridno); 

if(bl_saveimg) 

saveas(gcf, '/Users/sidhantkhanna/Documents/GitHub/BKS modified/code/Households/figures/hhsave_VFI_3OC/Valuefunction_2D_1.png') 

end 

figure(5) 

for i = 1: it_zgridno 

a(i)  = plot((mt_pf(:,i,1)),'.'); 

b{i}  = sprintf(['Z shock =', num2str(ar_z(i))]); 

hold on 

end 

clr = jet(numel(a)); 

for m = 1:numel(a) 

   set(a(m),'Color',clr(m,:)) 

end 

legend(b,'Location','north') 

title('Policy Function,l=0'); 

xlabel('Asset') ; 

ylabel('Asset Choice'); 

xt = get(gca, 'XTick');                                 

set(gca, 'XTick', xt, 'XTickLabel', xt*fl_ahi/it_agridno); 

if(bl_saveimg) 

saveas(gcf, '/Users/sidhantkhanna/Documents/GitHub/BKS modified/code/Households/figures/hhsave_VFI_3OC/Policyfunction_0.png') 

end 

figure(6) 

for i = 1: it_zgridno 

a(i)  = plot((mt_pf(:,i,2)),'.'); 

b{i}  = sprintf(['Z shock =', num2str(ar_z(i))]); 

hold on 

end 

clr = jet(numel(a)); 

for m = 1:numel(a) 

   set(a(m),'Color',clr(m,:)) 

end 

legend(b,'Location','north') 

title('Policy Function, l=1'); 

xlabel('Asset') ; 

ylabel('Asset choice'); 

xt = get(gca, 'XTick');                                 

set(gca, 'XTick', xt, 'XTickLabel', xt*fl_ahi/it_agridno); 

if(bl_saveimg) 

saveas(gcf, '/Users/sidhantkhanna/Documents/GitHub/BKS modified/code/Households/figures/hhsave_VFI_3OC/Policyfunction_1.png') 

end 

figure(7) 

for i = 1: it_zgridno 

a(i)  = plot((mt_con(:,i,1)),'-o'); 

b{i}  = sprintf(['Z shock =', num2str(ar_z(i))]); 

hold on 

end 

clr = jet(numel(a)); 

for m = 1:numel(a) 

   set(a(m),'Color',clr(m,:)) 

end 

legend(b,'Location','north') 

title('Consumption,l=0'); 

xlabel('Asset') ; 

ylabel('Consumption'); 

xt = get(gca, 'XTick');                                 

set(gca, 'XTick', xt, 'XTickLabel', xt*fl_ahi/it_agridno); 

if(bl_saveimg) 

saveas(gcf, '/Users/sidhantkhanna/Documents/GitHub/BKS modified/code/Households/figures/hhsave_VFI_3OC/Consumption_0.png') 

end 

figure(8) 

for i = 1: it_zgridno 

a(i)  = plot((mt_con(:,i,2)),'-o'); 

b{i}  = sprintf(['Z shock =', num2str(ar_z(i))]); 

hold on 

end 

clr = jet(numel(a)); 

for m = 1:numel(a) 

   set(a(m),'Color',clr(m,:)) 

end 

legend(b,'Location','north') 

title('Consumption, l=1'); 

xlabel('Asset') ; 

ylabel('Consumption'); 

xt = get(gca, 'XTick');                                 

set(gca, 'XTick', xt, 'XTickLabel', xt*fl_ahi/it_agridno); 

if(bl_saveimg) 

saveas(gcf, '/Users/sidhantkhanna/Documents/GitHub/BKS modified/code/Households/figures/hhsave_VFI_3OC/Consumption_1.png') 

end 

figure(9) 

for i = 1: it_zgridno 

a(i)  = plot(mt_savefraccoh(:,i,1),'-o'); 

b{i}  = sprintf(['Z shock =', num2str(ar_z(i))]); 

hold on 

end 

clr = jet(numel(a)); 

for m = 1:numel(a) 

   set(a(m),'Color',clr(m,:)) 

end 

legend(b,'Location','South') 

title('SavingsfracCOH,l=0'); 

xlabel('Asset') ; 

ylabel('SavingsfracCOH'); 

ylim([0 1]) 

xt = get(gca, 'XTick');                                 

set(gca, 'XTick', xt, 'XTickLabel', xt*fl_ahi/it_agridno); 

if(bl_saveimg) 

saveas(gcf, '/Users/sidhantkhanna/Documents/GitHub/BKS modified/code/Households/figures/hhsave_VFI_3OC/Savings_COH_0.png') 

end 

figure(10) 

for i = 1: it_zgridno 

a(i)  = plot(mt_savefraccoh(:,i,2),'-o'); 

b{i}  = sprintf(['Z shock =', num2str(ar_z(i))]); 

hold on 

end 

clr = jet(numel(a)); 

for m = 1:numel(a) 

   set(a(m),'Color',clr(m,:)) 

end 

legend(b,'Location','South') 

title('SavingsfracCOH, l=1'); 

xlabel('Asset') ; 

ylabel('SavingsfracCOH'); 

ylim([0 1]) 

xt = get(gca, 'XTick');                                 

set(gca, 'XTick', xt, 'XTickLabel', xt*fl_ahi/it_agridno); 

if(bl_saveimg) 

saveas(gcf, '/Users/sidhantkhanna/Documents/GitHub/BKS modified/code/Households/figures/hhsave_VFI_3OC/Savings_COH_1.png') 

end 

end 

%% Plot Value function and Policy Function

if(bl_profile) 

profile off; 

profile viewer;

st_file_name = '/Users/sidhantkhanna/Documents/GitHub/BKS modified/code/Profile/Households/hhsave_VFI_3OC';

profsave(profile('info'), st_file_name);

end