function [mt_dis, mt_o, mt_k, mt_l]= hhsave_dis_2OC(varargin)

addpath(genpath('/Users/sidhantkhanna/Desktop/Research/Matlab Code/My codes'));

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

        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_sigma,fl_rho, fl_tolvfi, fl_toldis ...

        ] = varargin{:};

    bl_print  = true;

    bl_plot   = false;

else 

    close all; 

    fl_ahi         = 50; 

    fl_zhi         = 2.2; 

    it_agridno     = 5; 

    it_zgridno     = 6; 

    ar_a           = linspace(0,fl_ahi,it_agridno); 

    fl_phi         = 0.5; 

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

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

    ar_z           = exp(ar_z); 

    fl_r           = 0.02; 

    fl_w           = 1.5; 

    fl_tolvfi    = 10^-12; 

    fl_toldis    = 10^-12;              % Tolerance level for convergence stationary distribution 

    bl_print       = true; 

end 

fl_R  = fl_r + fl_delta; 

% Calling VFI

[mt_vf, mt_pf, mt_oploc, mt_o, mt_k, mt_l]=hhsave_VFI_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,mt_trans_z,fl_beta, fl_mu, fl_sig,fl_rho, fl_tolvfi); 

%% Evaluating Stationary Distribution

it_s=1;     % Counting iteration for Convergence to Stationary Distribution 

fl_sum1(1)=0; % Variable to check if all elements of Stationary distribution Sum to 1 

fl_crit1=1; 

P(:,:,it_s+1)=zeros(it_agridno,it_zgridno); 

P(:,:,1)=ones(it_agridno,it_zgridno).*(1/(it_agridno*it_zgridno));  

while(fl_crit1>fl_toldis) 

    P(:,:,it_s+1)=zeros(it_agridno,it_zgridno); 

    for i= 1:it_agridno 

        for j= 1:it_zgridno 

            f=find((ar_a==mt_pf(i,j)));                    

                    for g=1:it_zgridno 

                    P(f,g,it_s+1)= P(f,g,it_s+1)+ mt_trans_z(j,g)*P(i,j,it_s); %evaluates Probability in a particular row and column in the next period 

                end 

            end 

    end 

    fl_sum1(it_s+1)=sum(sum(P(:,:,it_s+1))); % sum of all elements of the matrix P(:,:,s) in a particular iteration, which must be equal to 1 for all s 

    % we observe that the vector sum1 has value 1 in all the cells on running the code

    mt_diff1=(P(:,:,it_s+1)- P(:,:,it_s)); 

fl_crit1=norm(mt_diff1); 

it_s=it_s+1; 

end  

fl_Ea=sum(sum(P(:,:,it_s),2).*ar_a'); 

mt_dis=P(:,:,it_s); 

%% Printing outputs 

if(bl_print) 

    disp('Below is the Stationary Distribution'); 

    disp((mt_dis)); 

 end 

if(bl_profile) 

profile off; 

profile viewer;

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

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

end