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

close all; 

bl_profile = false;      % 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_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') = 1;

 %  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', mp_grid_control);

    ar_a           = ar_a';

    fl_phi         = 0;

    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.2; 

    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