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

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

close all;

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

bl_saveimg = false;

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

        ] = varargin{:}; 

    bl_print  = false; 

    bl_plot   = false; 

else

    fl_alo         = 0;

    fl_ahi         = 10;

    fl_zhi         = 2;

    it_agridno     = 100;

    it_zgridno     = 7;

    it_ngridno     = 2;

   %ar_a           = linspace(0,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', mp_grid_control);

    ar_a           = ar_a';

    fl_phi         = 0.2;

    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           = exp(ar_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';   

    mt_trans_n =[0.7,0.3;0.7,0.3];

    ar_n =[0,1];

    [fl_r,fl_w] = ...

        deal(0.02,2);

    fl_tolvfi    = 10^-12;

    bl_print       = true;

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

end 

fl_R  = fl_r + fl_delta; 

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; 

% Calling VFI

[mt_vf, mt_pf, mt_oploc,mt_coh, mt_o, mt_k, mt_l]=hhsave_VFI_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, mt_trans_z,mt_trans_n,fl_beta, fl_mu_z, fl_sig_z,fl_rho_z, ... 

     fl_lambda_z, 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,it_ngridno); 

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

while(fl_crit1>fl_toldis) 

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

    for i= 1:it_agridno 

        for j= 1:it_zgridno 

            for c= 1:it_ngridno 

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

                    for g=1:it_zgridno 

                        for h = 1: it_ngridno 

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

                end 

            end 

            end 

        end 

    end 

    fl_sum1(it_s+1)=sum(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)); 

    mt_diff1 = mt_diff1(:); 

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));

     st_ar_a = strsplit(num2str(ar_a));

    for i=1:it_agridno

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

    end

    disp('Table for Stationary Distribution at l=0');

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

    tb_dis.Properties.RowNames = st_a;

    disp(tb_dis); 

    disp('Table for Stationary Distribution at l=1');

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

    tb_dis.Properties.RowNames = st_a;

    disp(tb_dis);

 end 

if(bl_profile) 

profile off; 

profile viewer;

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

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

end