switch experiment

    case 'sizedep'

        p.taus        = Mu - 1;           % size-dependent subsidy
    
        Xnew          = fsolve('findequilibrium', [D; Y*1.05; N*1.1], options, w, z, p, 'planner', 'new');
        %Xnew          = fsolve('findequilibrium', [D; Y; N], options, w, z, p, 'planner', 'new');

    case 'efficient'
    
        p.taus        = 0; 
        Xnew          = fsolve('findequilibrium', [D*1.2; Y*1.6; N*1.2], options, w, z, p, 'planner', 'new'); %extra mult makes no differnece
        %Xnew          = fsolve('findequilibrium', [D; Y; N], options, w,
        %z, p, 'planner', 'new');

end

% New Steady State Allocations

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

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

Gnew          = (Znew/Dnew/(1 + p.taus) - 1)/Munew;
Qnew          = 1/(1 - p.beta*(1 - p.varphi))*Ynew/Nnew*Gnew; 
Lnew          = Lpnew + p.F*p.varphi*Nnew; 
Rnew          = 1/p.beta - 1 + p.delta; 

ynew          = qnew*Ynew; 

% Note change in focus to steady state welfare
Wolds = 1/(1 - p.beta)*(log(C)  - p.psi/(1 + p.nu)*L^(1 + p.nu)); 

Wnews = 1/(1 - p.beta)*(log(Cnew) - p.psi/(1 + p.nu)*Lnew^(1 + p.nu)); 

dWs   = (exp((1-p.beta)*(Wnews - Wolds)) - 1)*100;


GDPnew   = p.phi^(1/(p.theta - 1))*(1 - (1 - p.phi)*Z^(p.theta - 1))*Znew/(1 - (1 - p.phi)*Znew^(p.theta - 1))^(p.theta/(p.theta - 1))*Znew;
GDP   = p.phi^(1/(p.theta - 1))*(1 - (1 - p.phi)*Z^(p.theta - 1)*Mu^(-p.theta))*Z/(1 - (1 - p.phi)*Z^(p.theta - 1)*Mu^(1 - p.theta))^(p.theta/(p.theta - 1))*Z;
GDP1  = p.phi^(1/(p.theta - 1))*(1 - (1 - p.phi)*Z^(p.theta - 1)              )*Z/(1 - (1 - p.phi)*Z^(p.theta - 1)                 )^(p.theta/(p.theta - 1))*Z;


% 
% % Now start doing transition dynamics
% 
% global cfun fspace; 
% 
% npts      = 25; 
% 
% fspace    = fundefn('cheb', npts, 0.9, 1.5); 
% 
% Ngrid     = funnode(fspace); 
% 
% Zgrid     = zeros(npts, 1); 
% Ogrid     = zeros(npts, 1); 
% Mugrid    = zeros(npts, 1); 
% Dgrid     = zeros(npts, 1);
% Pgrid     = zeros(npts, 1);
% 
% options.Display = 'none';
% 
% for i = 1 : npts
% 
%     Dp      = fsolve('findequilibrium', Dp, options, w, z, p, 'planner', 'old', Ngrid(i));
% 
%     [~, ~, ~, Ogrid(i), Zgrid(i), ~, ~, ~, ~, ~, ~, ~, Mugrid(i), Pgrid(i)]  = findequilibrium(Dp, w, z, p, 'planner', 'old', Ngrid(i));
% 
%     Dgrid(i) = Dp; 
% 
% end
% 
% Ggrid  = (Zgrid./Dgrid/(1 + p.taus) - 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; Onew; 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)]; 
% 
% T    = length(Ct) - 1; 
% 
% Wold = 1/(1 - beta)*(log(C)  - psi/(1 + nu)*L^(1 + nu)); 
% 
% Wnew = (beta.^(0 : 1 : T - 1))*(log(Ct(2 : end)) - psi/(1 + nu)*Lt(2:end).^(1 + nu)) + beta^T/(1 - beta)*(log(Ct(end)') - psi/(1 + nu)*Lt(end)'.^(1 + nu));
% 
% dW   = (exp((1 - beta)*(Wnew - Wold)) - 1)*100;

clc

switch experiment

    case 'sizedep'

        Experiment = 'Size-dependent Subsidy';

    case 'efficient'

        Experiment = 'Efficient Allocation';

end

fprintf('<strong> Table 4: Implications of Alternative Policies, Benchmark Model </strong> \n'); 

fprintf('\n');

fprintf(' Steady state comparisons in percent - aggregate markup at %1.2f \n', Mu);  

fprintf('\n');

fprintf([' Experiment --> ' Experiment, ' \n']);  

fprintf('\n');

disp(table(round((Ynew/Y   - 1)*100, 1), round((Cnew/C   - 1)*100, 1), round((Lnew/L   - 1)*100, 1), round((Nnew/N   - 1)*100, 1), round((Knew/K   - 1)*100, 1), round((Znew/Z   - 1)*100, 1), round(dWs, 2), ...
    'VariableNames', {'Y','C','L','N','K','Z','Welfare,%'}))
