experiment = 'entry';

switch experiment

    
case 'entry'
        
    p.xis        = 0; 
    
    p.xie        = 3;        % approximately triples N if Mu=1.15
    
    
    %p.xie        = 0.6;     % optimal entry subsidy, leads to 0.55668% welfare gains
    
    %optset('golden', 'tol', 1e-5); 
    
    %p.xie        = golden('findxie', 0, 1, D, Y, N, K, C, L, Lp, options, w, z, p); 
    
    Xnew         = fsolve('findequilibrium', [D; Y; N], options, w, z, p, 'market', 'new');

    [~, qnew, ~, ~, Znew, Wnew, Lpnew, Knew, Bnew, Unew, Uqnew, Cnew, Munew] = findequilibrium(Xnew, w, z, p, 'market', 'new');

    
case 'uniform'               % leads to 5.89706 welfare gains
    
    p.xis        = Mu - 1; 
    Xnew         = fsolve('findequilibrium', [D*1.2; Y*1.6; N*1.2], options, w, z, p, 'market', 'new');

    [~, qnew, ~, ~, Znew, Wnew, Lpnew, Knew, Bnew, Unew, Uqnew, Cnew, Munew] = findequilibrium(Xnew, w, z, p, 'market', 'new');


end

% New Steady State Allocations

Dnew          = Xnew(1); 
Ynew          = Xnew(2); 
Nnew          = Xnew(3); 


Gnew          = (1 + p.xie)*(1 + p.xis - 1/Munew);
Qnew          = 1/(1 - p.beta*(1 - p.varphi))*Gnew*Ynew/Nnew; 

Lnew          = Lpnew + p.F*p.varphi*Nnew; 
Rnew          = 1/p.beta - 1 + p.delta; 

ynew          = qnew*Ynew; 


fprintf('Various Checks of the Code\n');
fprintf('\n');
fprintf('Error in Kimball aggregator             = %9.2e \n',  Nnew*w'*Unew - 1);
fprintf('Error in Consumption-Leisure Choice     = %9.2e \n',  p.psi - Wnew/(Cnew*Lnew^p.nu));
fprintf('Error in Agg Prod Function              = %9.2e \n',  Ynew - Znew*(p.phi^(1/p.theta)*(Knew^p.alpha*Lpnew^(1 - p.alpha))^(1 - 1/p.theta) + (1 - p.phi)^(1/p.theta)*Bnew^(1 - 1/p.theta))^(p.theta/(p.theta - 1)));

fprintf('\n');
fprintf('\n');

fprintf('\n');
fprintf('\n');
fprintf('Compare to original allocations: 1. Old, 2. New 3. Change \n');
fprintf('\n');
fprintf('\n');
fprintf('Measure firms, N             = %9.3f %9.3f %9.3f\n',  [N,       Nnew,        Nnew/N   - 1]);
fprintf('Output, Y                    = %9.3f %9.3f %9.3f\n',  [Y,       Ynew,        Ynew/Y   - 1]);
fprintf('Consumption, C               = %9.3f %9.3f %9.3f\n',  [C,       Cnew,        Cnew/C   - 1]);
fprintf('Capital, K                   = %9.3f %9.3f %9.3f\n',  [K,       Knew,        Knew/K   - 1]);
fprintf('Employment, L                = %9.3f %9.3f %9.3f\n',  [L,       Lnew,        Lnew/L   - 1]);
fprintf('Employment (production), Lp  = %9.3f %9.3f %9.3f\n',  [Lp,      Lpnew,       Lpnew/Lp - 1]);
fprintf('TFP, Z                       = %9.3f %9.3f %9.3f\n',  [Z,       Znew,        Znew/Z   - 1]);
fprintf('True Markup, Mu              = %9.3f %9.3f %9.3f\n',  [Mu,      Munew,       Munew/Mu   - 1]);
    
fprintf('\n');
fprintf('\n');


% Now start doing transition dynamics

% See if we can use dynare: parameterize Z(N) with chebyshev polynomials

global cfun fspace; 

npts      = 11; 

optset('chebnode', 'nodetype', 2);              % Lobatto nodes (include endpoints)

fspace    = fundefn('cheb', npts, 1, 2); 

Ngrid     = funnode(fspace); 


Zgrid     = zeros(npts, 1); 
Mugrid    = zeros(npts, 1); 

options.Display = 'none';

for i = 1 : npts

    D      = fsolve('findequilibrium', D, options, w, z, p, 'market', 'old', Ngrid(i));

    [~, ~, ~, ~, Zgrid(i), ~, ~, ~, ~, ~, ~, ~, Mugrid(i)]  = findequilibrium(D, w, z, p, 'market', 'old', Ngrid(i));
    
end


Ggrid  = (1 + p.xie)*(1 + p.xis - 1./Mugrid);               % flow gains from varieties

cfun   = funfitxy(fspace, Ngrid, [Zgrid, Ggrid, Mugrid]);

xpar   = [p.beta; p.varphi; p.delta; p.psi; p.F; p.nu; p.alpha; p.theta; p.phi];                              save xpar xpar; 

x1_ss  = [K;    N];                                                                                           save x1_ss x1_ss; 
x2_ss  = [Ynew; Cnew; Bnew; Lpnew; Knew; Nnew; Wnew; Rnew; Qnew; Znew; Gnew; Znew/Munew; Munew];              save x2_ss x2_ss; 


dynare transition noclearall;

% Cleanup

filename = 'transition';

eval(['delete ' filename '*.m']);
eval(['delete ' filename '*.mat']);
eval(['delete ' filename '*.swp']);
eval(['delete ' filename '*.log']);
eval(['delete ' filename '*.eps']);
eval(['delete ' filename '*.pdf']);
eval(['delete ' filename '*.fig']);
eval(['delete ' filename '*.jnl']);


% Populate first period values from original steady state

Yt(1)  = Y; 
Ct(1)  = C; 
Bt(1)  = B; 
Lpt(1) = Lp; 
Kt(1)  = K; 
Nt(1)  = N; 
Wt(1)  = W; 
Rt(1)  = R; 
Zt(1)  = Z; 
Mut(1) = Mu; 

Lt     = [L; Lpt(2:end) + p.F*(Nt(2:end) - (1 - p.varphi)*Nt(1 : end-1))]; 
Xt     = [p.delta*K; Kt(2:end) - (1 - p.delta)*Kt(1 : end - 1)]; 

figure(50)

subplot(2,2,1)
hold on
plot([Yt(1 : 60)./Yt(1) - 1, Kt(1 : 60)./Kt(1) - 1])
title('Output, Capital','Interpreter','latex'); 
xlabel('years','Interpreter','latex');
h = legend('output', 'capital');
set(h,'Interpreter','latex')
xlim([-2, 62])

set(gca, 'ygrid', 'on');

subplot(2, 2, 2)
hold on
plot([Ct(1 : 60)./Ct(1) - 1, Lt(1 : 60)./Lt(1) - 1])
title('Consumption, Employment','Interpreter','latex'); 
xlabel('years','Interpreter','latex');
h = legend('consumption', 'employment');
set(h,'Interpreter','latex')
xlim([-2, 62])

set(gca, 'ygrid', 'on');


subplot(2, 2, 3)
hold on
plot(Zt(1 : 60)./Zt(1) - 1)
title('Productivity','Interpreter','latex'); 
xlabel('years','Interpreter','latex');
set(h, 'Interpreter', 'latex')
xlim([-2, 62])

set(gca, 'ygrid', 'on');

subplot(2, 2, 4)
hold on
plot(Nt(1 : 60)./Nt(1) - 1)
title('Mass of Producers','Interpreter','latex'); 
xlabel('years','Interpreter','latex');
set(h,'Interpreter','latex')
xlim([-2, 62])

set(gca, 'ygrid', 'on');


% check errors in Euler equations: 

figure(51)

subplot(2, 2, 1)
hold on
plot(p.beta*Ct(2:end-1)./Ct(3:end).*(Rt(3:end) + 1 - p.delta) - 1)
title('Euler Equation','Interpreter','latex'); 
xlabel('years','Interpreter','latex');

subplot(2, 2, 2)
hold on
plot(p.beta*Ct(2:end-1)./Ct(3:end).*Qt(3:end) - p.F*Wt(2:end - 1))
title('Free Entry','Interpreter','latex'); 
xlabel('years','Interpreter','latex');

subplot(2, 2, 3)
hold on
plot(Zt(2:end).*(p.phi^(1/p.theta)*(Kt(1:end-1).^p.alpha.*Lpt(2:end).^(1 - p.alpha)).^(1 - 1/p.theta) + (1 - p.phi)^(1/p.theta)*Bt(2:end).^(1 - 1/p.theta)).^(p.theta/(p.theta - 1)) - Ct(2:end) - Bt(2:end) - Kt(2:end) + (1 - p.delta)*Kt(1:end-1));
title('Resource Constraints','Interpreter','latex'); 
xlabel('years','Interpreter','latex');

subplot(2, 2, 4)
hold on
plot(Wt(2:end) - p.psi*Ct(2:end).*Lt(2:end).^p.nu);
title('Labor Choice','Interpreter','latex'); 
xlabel('years','Interpreter','latex');


T    = length(Ct) - 1; 

Wold = 1/(1 - p.beta)*(log(C)  - p.psi/(1 + p.nu)*L^(1 + p.nu)); 

Wnew = (p.beta.^(0 : 1 : T - 1))*(log(Ct(2 : end)) - p.psi/(1 + p.nu)*Lt(2:end).^(1 + p.nu)) + p.beta^T/(1 - p.beta)*(log(Ct(end)') - p.psi/(1 + p.nu)*Lt(end)'.^(1 + p.nu));

dW   = (exp((1 - p.beta)*(Wnew - Wold)) - 1)*100;

fprintf('\n')
fprintf('Welfare Gains                       = %9.5f \n',  dW);
fprintf('\n')

p.xie
