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

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, ... 

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

         fl_risk,it_zgridno, it_agridno,mt_trans_z,fl_beta, fl_mu, fl_sig,fl_rho, fl_tolvfi ... 

        ] = varargin{:}; 

    bl_print  = false; 

    bl_plot   = false; 

else

    close all;

    fl_ahi         = 50;

    fl_zhi         = 1.8;

    it_agridno     = 15;

    it_zgridno     = 5;

    ar_a           = linspace(0,fl_ahi,it_agridno);

    fl_phi         = 0.3;

    fl_risk        = 1;

    fl_alpha       = 0.4;

    fl_theta       = 0.79-fl_alpha;

    fl_delta       = 0.05;

    fl_kappa       = 0;

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

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

    fl_sig         = 0.2;     

    fl_lambda      = 3;      

    fl_beta       = 0.92; 

    [ar_z, mt_trans_z] = mytauchen_z(fl_mu,fl_rho,fl_sig,it_zgridno,fl_lambda);    

    ar_z           = exp(ar_z);

    [fl_r,fl_w] = ...

        deal(0.05,2);

    fl_tolvfi    = 10^-12;

    bl_print     = false;

    bl_plot      = false; 

end 

fl_R  = fl_r + fl_delta; 

%% Value function Iteration

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

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

mt_oploc = zeros(it_agridno,it_zgridno);             % 

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

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

         fl_risk,it_zgridno, it_agridno); 

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

while fl_crit > fl_tolvfi 

mt_evf = mt_vf*mt_trans_z'; 

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

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

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

    mt_diff= mt_Tvf-mt_vf; 

    fl_crit=norm(mt_diff); 

    mt_vf = mt_Tvf; 

end 

mt_pf  = ar_a(mt_oploc); 

mt_con = mt_coh - mt_pf;  

if(bl_print) 

    disp('Below is the Value function');

    disp(mt_vf);

    disp('Below is the Policy Function');

    disp(mt_pf);

    disp('Below is the Location of optimal choice of assets');

    disp(mt_oploc);

 end

%% Printing Outputs

%% Plot Value function, Policy Function and Consumption

% 3D Value function

if(bl_plot) 

figure(1)

chart=surf(ar_z,ar_a,mt_vf);

title('Value Function');

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_2OC/Valuefunction_3D.png')

end

% 2D Value function

figure(2)

%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

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

%clr = jet(numel(chart));

%for m = 1:numel(chart) 

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

%end

title('Value Function');

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_2OC/Valuefunction_2D.png')

end

figure(3)

for i = 1: it_zgridno

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

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

hold on

end

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

title('Policy Function');

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_2OC/Policyfunction.png')

end

figure(4)

for i = 1: it_zgridno

a(i)  = plot((mt_con(:,i,1)));

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

hold on

end

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

title('Consumption');

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_2OC/Consumption.png')

end

end

%% Profile off

if(bl_profile) 

profile off;

profile viewer;

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

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

end