function [mt_k_star, mt_l_star]=optikl(varargin)

close all; 

bl_profile = false;      % 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, ... 

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

         it_zgridno, it_agridno, ... 

        ] = varargin{:}; 

    bl_print  = false; 

    bl_plot   = false; 

else

    close all;

    fl_ahi         = 5;

    fl_zhi         = 2.2;

    it_agridno     = 5;

    it_zgridno     = 6;

    ar_a = linspace(0,fl_ahi,it_agridno);

    ar_z = linspace(1,fl_zhi,it_zgridno);

    fl_phi         = 0.5;

    fl_alpha = 0.4;

    fl_theta = 0.79-fl_alpha;

    fl_delta = 0.05;

    fl_kappa = 0;

    [fl_r,fl_w,fl_ahi] = ...

        deal(0.05,3,100);

    bl_print       = true;

    bl_plot        = true;

end 

%% Set Parameters

fl_R = fl_r + fl_delta; 

%% Set Controls

bl_constrained = true; 

if( fl_phi == 1) 

    bl_constrained = false;

end

bl_feas = true; % 

%% Mesh Grids

[a_m,z_m] = meshgrid(ar_a,ar_z);    % creates 2 matrices of size it_zgridno x it_agridno  

                                    % a_m has it_zgridno no. of identical rows, each row is ar_a. So each column

                                    % has identical a element repeated it_zgridno times

                                    % z_m has it_agridno no of identical columns, each column is ar_z. So each

                                    % row has identical z element repeated it_agridno times

ar_a_v = a_m(:);                    % matrix unlayered one column under another 

ar_z_v = z_m(:);                    % matrix unlayered one column under another 

ar_k_unc = exp((log(fl_R./(ar_z_v.*fl_alpha))- fl_theta*log(fl_R./(ar_z_v.*fl_alpha))+ fl_theta.*log(fl_w./(ar_z_v.*fl_theta)))./(fl_alpha+fl_theta-1)); 

ar_l_unc = exp((log(fl_w./(ar_z_v.*fl_theta))+ fl_alpha*log(fl_R./(ar_z_v.*fl_alpha))- fl_alpha.*log(fl_w./(ar_z_v.*fl_theta)))./(fl_alpha+fl_theta-1)); 

mt_k_unc = reshape(ar_k_unc, [it_zgridno,it_agridno]); 

mt_l_unc = reshape(ar_l_unc, [it_zgridno,it_agridno]); 

mt_k_unc = mt_k_unc'; 

mt_l_unc = mt_l_unc'; 

%% Optimal Unconstrained 

if (~bl_constrained) 

    % in the constrained case, a can not matter, but z does matter

    if(bl_feas)

    feas  = (1+fl_r)*fl_kappa;

    feas_v = ar_z_v.*ar_k_unc.^fl_alpha.*ar_l_unc.^fl_theta;

    for i = 1:numel(ar_a_v)

       if (feas_v(i)<feas)

       ar_k_unc(i) = 0;

       ar_l_unc(i) = 0;

       end

    end

    end

    mt_k_unc = reshape(ar_k_unc, [it_zgridno,it_agridno]);

    mt_l_unc = reshape(ar_l_unc, [it_zgridno,it_agridno]);

    mt_k_unc = mt_k_unc';

    mt_l_unc = mt_l_unc';

    mt_k_star = mt_k_unc;

    mt_l_star = mt_l_unc;

if( bl_print)

    disp('Below is mt_k_star')

    disp(size(mt_k_star));

    disp(mt_k_star);

    disp('Below is mt_l_star')

    disp(size(mt_l_star));

    disp(mt_l_star);

end

else 

    %% Optimal Constrained

    % calls function con_kl to evaluate the constraint on capital

    [ar_k_bar,ar_l_bar] = con_kl(ar_a_v, ar_z_v,ar_a,ar_z, fl_ahi,fl_zhi,fl_phi,fl_theta,fl_alpha,fl_w,fl_r,fl_delta,fl_kappa,it_agridno,it_zgridno); 

    mt_k_bar=reshape(ar_k_bar, [it_zgridno,it_agridno]); 

    mt_l_bar=reshape(ar_l_bar, [it_zgridno,it_agridno]); 

    mt_k_bar = mt_k_bar'; 

    mt_l_bar = mt_l_bar'; 

    k3        = [ar_k_unc,ar_k_bar];       

    [k4,id]   = min(abs(k3'));            % K star is min of k_unc and k_con. k3 transposed for the correct order for id 

    ar_k_star = k4'; 

    l3        =[ar_l_unc,ar_l_bar]; 

    ar_l_star = zeros(numel(ar_a_v),1);    % assigning l star in constrained case 

    for i = 1:numel(ar_a_v) 

        ar_l_star(i,1)=l3(i,id(i)); 

    end 

    if(bl_feas) 

    feas  = (1+fl_r)*fl_kappa; 

    feas_v = ar_z_v.*ar_k_star.^fl_alpha.*ar_l_star.^fl_theta; 

    for i = 1:numel(ar_a_v) 

       if (feas_v(i)<feas) 

       ar_k_star(i) = 0;

       ar_l_star(i) = 0;

       end

    end 

    end 

mt_k_star=reshape(ar_k_star, [it_zgridno,it_agridno]); 

mt_l_star=reshape(ar_l_star, [it_zgridno,it_agridno]); 

 mt_k_star=mt_k_star'; 

 mt_l_star=mt_l_star'; 

if( bl_print) 

    disp('Below is mt_k_unc')

    disp(size(mt_k_unc));

    disp(mt_k_unc);

    disp('Below is mt_k_con')

    disp(mt_k_bar);

    disp('Below is mt_k_star')

    disp(size(mt_k_star));

    disp(mt_k_star);

    disp('Below is mt_l_unc')

    disp(size(mt_l_unc));

    disp(mt_l_unc);

    disp('Below is mt_l_con')

    disp(mt_l_bar);

    disp('Below is mt_l_star')

    disp(mt_l_star);

end

if(bl_plot) 

%% 3D Plots for k* and l* 

figure (1)

a1 = surf(ar_z, ar_a, mt_k_bar , 'FaceColor','r', 'FaceAlpha',0.5, 'EdgeColor','r'); M1 = 'k bar';

hold on

a2 = surf(ar_z, ar_a, mt_k_unc , 'FaceColor','b', 'FaceAlpha',0.5, 'EdgeColor','b'); M2 = 'k unc';

%hold on

%a3 = surf(ar_z, ar_a, mt_k_star,'FaceColor','g', 'FaceAlpha',0.5, 'EdgeColor','g'); M3 = 'k star' ;

hold off

xlabel('z');

ylabel('a');

zlabel('k(unconstrained/constrained/optimal)');

LEG = legend([a1;a2], M1, M2, 'Location','north')

LEG.FontSize = 20;

view(150,15);

set(gca,'FontSize',17)

if(bl_saveimg)

saveas(gcf, '/Users/sidhantkhanna/Documents/GitHub/BKS modified/code/Firms/figures/optikl/Optimum_k.png')

end

figure (2)

b1 = surf(ar_z, ar_a, mt_l_bar, 'FaceColor','r', 'FaceAlpha',0.5, 'EdgeColor','r'); N1 = 'l bar';

hold on

b2 = surf(ar_z, ar_a, mt_l_unc, 'FaceColor','b', 'FaceAlpha',0.5, 'EdgeColor','b'); N2 = 'l unc';

%hold on

%b3 = surf(ar_z, ar_a, mt_l_star,'FaceColor','g', 'FaceAlpha',0.5, 'EdgeColor','g'); N3 = 'l star' ;

hold off

xlabel('z');

ylabel('a');

zlabel('l(unconstrained/constrained/optimal)');

LEG = legend([b1;b2], N1, N2, 'Location','north')

LEG.FontSize = 20;

view(150,15);

set(gca,'FontSize',17)

if(bl_saveimg)

saveas(gcf, '/Users/sidhantkhanna/Documents/GitHub/BKS modified/code/Firms/figures/optikl/Optimum_l.png')

end

%% 2D plot for k*

% Jet color Graph All

figure(3);

hold on;

% one line, unconstrained choice

pl_unc = plot(ar_z, mt_k_unc(1,:));

pl_unc.HandleVisibility = 'on';

pl_unc.LineStyle = '-';

pl_unc.Color = 'black';

pl_unc.LineWidth = 4;

% one constrained choice line, low a

pl_unc = plot(ar_z, mt_k_bar(2,:));

pl_unc.HandleVisibility = 'on';

pl_unc.LineStyle = '--';

pl_unc.Color = 'red';

pl_unc.LineWidth = 3;

% one constrained choice line, high a

pl_unc = plot(ar_z, mt_k_bar(4,:));

pl_unc.HandleVisibility = 'on';

pl_unc.LineStyle = '-.';

pl_unc.Color = 'blue';

pl_unc.LineWidth = 3;

% additional options

grid on;

grid minor;

xlim([1, 2.2]);

ylim([0, 10]);

title('Unconstrained Choices and Endogenous Borrowing Constraints');

ylabel('Capital Levels');

xlabel({'Productivity Levels'});

legend({'Unconstrained Optimal Choice',...

        'Borrowing Bound low A', ...

        'Borrowing Bound high A'}, 'Location', 'best', 'Interpreter','latex');

if(bl_saveimg)

saveas(gcf, '/Users/sidhantkhanna/Documents/GitHub/BKS modified/code/Firms/figures/optikl/Optimum_2D.png')

end

hold off;

end

end 

if(bl_profile) 

profile off;

profile viewer;

st_file_name = '/Users/sidhantkhanna/Documents/GitHub/BKS modified/code/Profile/Firms/optikl';

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

end