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

    fl_alo         = 0;
    fl_ahi         = 50;
    fl_zhi         = 2.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.5;
 %  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.5;
    fl_risk        = 1.5;
    fl_alpha       = 0.4;
    fl_theta       = 0.79-fl_alpha;
    fl_delta       = 0.08;
    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';
    mt_trans_n =[0.3,0.7;0.3,0.7];
    ar_n =[0,1];
    ar_z           = exp(ar_z)';
    [fl_r,fl_w] = ...
        deal(0.04,3);
    fl_tolvfi    = 10^-12;
    fl_toldis    = 10^-12;              % Tolerance level for convergence stationary distribution
    bl_print  = true;
    bl_plot   = true;

[a_m,z_m] = meshgrid(ar_a,ar_z);
a_m       = a_m';
z_m       = z_m';
a_m_3OC   = repmat (a_m, [1 1 it_ngridno]);
z_m_3OC   = repmat (z_m, [1 1 it_ngridno]);

ar_phi = [0:0.25:1];
ar_y = zeros(1,length(ar_phi));
ar_y_form = zeros(1,length(ar_phi));
ar_y_inf = zeros(1,length(ar_phi));
ar_kdemand = zeros(1,length(ar_phi));
ar_ldemand = zeros(1,length(ar_phi));
ar_ksupply = zeros(1,length(ar_phi));
ar_lsupply = zeros(1,length(ar_phi));
ar_percent_ent  = zeros(1,length(ar_phi));
ar_percent_se = zeros(1,length(ar_phi));
ar_percent_worker = zeros(1,length(ar_phi));
ar_mean_tal_ent = zeros(1,length(ar_phi));
ar_sd_tal_ent = zeros(1,length(ar_phi));
ar_mean_tal_se = zeros(1,length(ar_phi));
ar_sd_tal_se = zeros(1,length(ar_phi));
ar_meanl = zeros(1,length(ar_phi));
ar_sd_l = zeros(1,length(ar_phi));



for g=1:length(ar_phi)
    fl_phi=ar_phi(g);

[fl_y,mt_dis, mt_vf,mt_pf,mt_wage,mt_profit,mt_income,mt_coh, mt_o, mt_k, mt_l] = PE_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, fl_toldis );
ar_y(g)                 = fl_y;
mt_y                 = z_m_3OC.*mt_dis.*mt_k.^fl_alpha.*mt_l.^fl_theta;
mt_y_inf             = zeros(it_agridno,it_zgridno,it_ngridno);
mt_y_form            = zeros(it_agridno,it_zgridno,it_ngridno);
mt_y_inf(mt_o==1)    = mt_y(mt_o==1);
mt_y_form(mt_o==2)   = mt_y(mt_o==2);
fl_y_inf             = sum(sum(sum(mt_y_inf)));
fl_y_form            = sum(sum(sum(mt_y_form)));
ar_y_inf(g)          = fl_y_inf;
ar_y_form(g)         = fl_y_form;


mt_ent               = zeros(it_agridno,it_zgridno,it_ngridno);
mt_se                = zeros(it_agridno,it_zgridno,it_ngridno);
mt_worker            = zeros(it_agridno,it_zgridno,it_ngridno);

mt_ent(mt_o==2)      = 1;
mt_se (mt_o==1)      = 1;
mt_worker(mt_o==0)   = 1;

ar_kdemand(g)   = sum(sum(sum(mt_dis.*mt_k)));
ar_ldemand(g)   = sum(sum(sum(mt_dis.*mt_l)));

% Capital and labor supplied

ar_ksupply(g)   = sum(sum(sum(mt_dis.*mt_pf)));
ar_lsupply(g)   = sum(sum(sum(mt_dis.*mt_worker)));


% Percent employed in each occupations

ar_percent_ent(g)       = sum(sum(sum(mt_dis.*mt_ent)));
ar_percent_se(g)        = sum(sum(sum(mt_dis.*mt_se)));
ar_percent_worker(g)    = sum(sum(sum(mt_dis.*mt_worker)));

% Talent of entrepreneurs

mt_dis_ent          = mt_dis;
mt_dis_ent(mt_o~=2) = 0;
mt_z_m_3OC_ent      = z_m_3OC;
mt_z_m_3OC_ent(mt_o~=2)= 0;
ar_z_m_3OC_ent      = mt_z_m_3OC_ent(:);
ar_dis_ent          = mt_dis_ent(:);
ar_dis_ent(ar_z_m_3OC_ent==0)   = [];
ar_z_m_3OC_ent(ar_z_m_3OC_ent==0) = [];
fl_length_ent = length(ar_dis_ent);

%fl_mean_tal_ent = (sum(sum(sum(z_m_3OC.*mt_dis_ent))))/(sum(sum(sum(mt_dis_ent))));
fl_mean_tal_ent = sum(ar_z_m_3OC_ent.*ar_dis_ent)/(sum(sum(sum(mt_dis_ent))));
ar_mean_tal_ent(g) = fl_mean_tal_ent;
W = ar_dis_ent;
X = ar_z_m_3OC_ent;
Wmean = fl_mean_tal_ent;
N = fl_length_ent;
ar_sd_tal_ent(g) = sqrt((sum(W.*(X - Wmean).^2)/N) / ((N-1)*sum(W)/N));

%Talent of Informal firm owners
mt_dis_se          = mt_dis;
mt_dis_se(mt_o~=1) = 0;

mt_z_m_3OC_se      = z_m_3OC;
mt_z_m_3OC_se(mt_o~=1)= 0;
ar_z_m_3OC_se      = mt_z_m_3OC_se(:);
ar_dis_se          = mt_dis_se(:);
ar_dis_se(ar_z_m_3OC_se==0)   = [];
ar_z_m_3OC_se(ar_z_m_3OC_se==0) = [];
fl_length_se = length(ar_dis_se);

%fl_mean_tal_ent = (sum(sum(sum(z_m_3OC.*mt_dis_ent))))/(sum(sum(sum(mt_dis_ent))));
fl_mean_tal_se = sum(ar_z_m_3OC_se.*ar_dis_se)/(sum(sum(sum(mt_dis_se))));
ar_mean_tal_se(g) = fl_mean_tal_se;
W = ar_dis_se;
X = ar_z_m_3OC_se;
Wmean = fl_mean_tal_se;
N = fl_length_se;
ar_sd_tal_se(g) = sqrt((sum(W.*(X - Wmean).^2)/N) / ((N-1)*sum(W)/N));
%

% Firm size distributions
ar_l   = mt_l(:);
ar_dis = mt_dis(:);
ar_dis(ar_l==0) = [];
ar_l(ar_l==0)   = [];
fl_meanl        = sum(ar_dis.*ar_l)/(sum(ar_dis));
fl_length_l = length(ar_l);

ar_meanl(g) = fl_meanl;
W = ar_dis;
X = ar_l;
Wmean = fl_meanl;
N = fl_length_l;
ar_sd_l(g) = sqrt((sum(W.*(X - Wmean).^2)/N) / ((N-1)*sum(W)/N));
%

end

disp ('Degree of Contract enforcement friction is:')
    disp(ar_phi);
    disp('Output per capita is:')
    disp(ar_y);
    disp('Output per capita - informal')
    disp(ar_y_inf);
    disp('Output per capita is - formal')
    disp(ar_y_form);
    disp('Total capital per capita demanded is:')
    disp(ar_kdemand);
    disp('Total labor per capita demanded is:')
    disp(ar_ldemand);
    disp('Total capital per capita supplied is:')
    disp(ar_ksupply);
    disp('Total labor per capita supplied is:')
    disp(ar_lsupply);
    disp('Percentage entrepreneurs')
    disp(ar_percent_ent);
    disp('Percentage self-employed other than entrepreneurs')
    disp(ar_percent_se);
    disp('Percentage workers');
    disp(ar_percent_worker);
    disp('Mean talent of Entrepreneurs');
    disp(ar_mean_tal_ent);
    disp('Mean talent of Self Employed');
    disp(ar_mean_tal_se);
    disp('SD talent of Entrepreneurs');
    disp(ar_sd_tal_ent);
    disp('SD talent of Self Employed');
    disp(ar_sd_tal_se);
    disp('Mean of Firm size');
    disp(ar_meanl);
    disp('SD of firm size');
    disp(ar_sd_l);
    dlmwrite('test.csv',ar_phi,'delimiter',',','-append');
    dlmwrite('test.csv',ar_y,'delimiter',',','-append');
    dlmwrite('test.csv',ar_y_form,'delimiter',',','-append');
    dlmwrite('test.csv',ar_y_inf,'delimiter',',','-append');
    dlmwrite('test.csv',ar_kdemand,'delimiter',',','-append');
    dlmwrite('test.csv',ar_ldemand,'delimiter',',','-append');
    dlmwrite('test.csv',ar_ksupply,'delimiter',',','-append');
    dlmwrite('test.csv',ar_lsupply,'delimiter',',','-append');
    dlmwrite('test.csv',ar_percent_ent,'delimiter',',','-append');
    dlmwrite('test.csv',ar_percent_se,'delimiter',',','-append');
    dlmwrite('test.csv',ar_percent_worker,'delimiter',',','-append');
    dlmwrite('test.csv',ar_mean_tal_ent,'delimiter',',','-append');
    dlmwrite('test.csv',ar_sd_tal_ent,'delimiter',',','-append');
    dlmwrite('test.csv',ar_mean_tal_se,'delimiter',',','-append');
    dlmwrite('test.csv',ar_sd_tal_se,'delimiter',',','-append');
    dlmwrite('test.csv',ar_meanl,'delimiter',',','-append');
    dlmwrite('test.csv',ar_sd_l,'delimiter',',','-append');
Degree of Contract enforcement friction is:
         0    0.2500    0.5000    0.7500    1.0000

Output per capita is:
   1.0e+03 *

    0.0195    0.0341    0.0918    0.3878    2.0995

Output per capita - informal
    0.6435    0.6050    0.5314    0.4683    0.3527

Output per capita is - formal
   1.0e+03 *

    0.0188    0.0335    0.0913    0.3874    2.0991

Total capital per capita demanded is:
   1.0e+03 *

    0.0191    0.0316    0.0826    0.5577    6.9983

Total labor per capita demanded is:
    2.5305    4.4393   11.9396   50.4174  272.9332

Total capital per capita supplied is:
   25.0275   25.2314   25.5210   24.7930   24.2706

Total labor per capita supplied is:
    0.3791    0.3512    0.3137    0.2788    0.2406

Percentage entrepreneurs
    0.4585    0.4983    0.5518    0.6017    0.6562

Percentage self-employed other than entrepreneurs
    0.1625    0.1505    0.1345    0.1195    0.1031

Percentage workers
    0.3791    0.3512    0.3137    0.2788    0.2406

Mean talent of Entrepreneurs
    4.4264    4.2713    4.0963    3.9608    3.8361

Mean talent of Self Employed
    2.0088    1.9708    1.9118    1.8430    1.7449

SD talent of Entrepreneurs
    0.1148    0.1042    0.0975    0.0904    0.0838

SD talent of Self Employed
    0.0180    0.0182    0.0175    0.0159    0.0098

Mean of Firm size
    4.0752    6.8419   17.3980   69.9114  359.4182

SD of firm size
    0.2230    0.5740    2.3309   11.8576   67.0103