% Now start doing transition dynamics

global cfun fspace; 

npts       = 5; 

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

if Markuptarget < 1.1
   
    fspace     = fundefn('cheb', npts, 0.65*N, 1.20*N); 

else
    
    fspace     = fundefn('cheb', npts, 0.80*N, 1.20*N); 

end

Ngrid      = funnode(fspace); 
Zgrid      = zeros(npts, 5); 

nsgrid     = zeros(S, npts, 5); 

% add firms incrementally, so when we add new firms, old continue

for kt = 1 : 5

unif       = rand(Nmax, S);   % fix this once so firms incrementally enter each sector 
ns         = nssave; 

for it = 1 : npts 
    
    if Ngrid(it) > N 
        
    ms        = poissrnd(Ngrid(it) - Ngrid(it-1), S, 1);

    ns        = ns  + ms;                               % add extra firms

    z         = zsave; 

    for i = 1 : S
   
      z(ns(i) + 1 : end, i) = eps;                      % kill extra firms (set z = eps)
    
    end

    z         = sort(z, 1, 'descend');

    nsmax     = max(ns); 
    z         = z(1 : nsmax, :);                        % no reason to carry extra entries

    else
     
    ns        = nssave;     
    z         = zsave; 

    for i = 1 : S
   
      z(ns(i) + 1 : end, i) = eps;                      % kill extra firms (set z = eps)
    
    end
    
    die       = unif < (N - Ngrid(it))/N;
    
    z(die)    = eps; 
        
    z         = sort(z, 1, 'descend');
           
    ns        = sum(z > eps)';        
    
    nsmax     = max(ns); 
    z         = z(1 : nsmax, :);                        % no reason to carry extra entries

    end
    
    nsgrid(:, it, kt) = ns; 
    
    temp           = sum(z.^(p.gamma - 1)).^(1/(p.gamma - 1)); 
    Zgrid(it, kt)  = mean(temp.^(p.sigma - 1)).^(1/(p.sigma - 1));
    
end
 
end

cfun           = funfitxy(fspace, kron(ones(5, 1), Ngrid), Zgrid(:));

nn             = nodeunif(1000, Ngrid(1), Ngrid(end));

close all

switch experiment
    
    case 'efficient'

        Mugrid         = ones(size(Ngrid));           % remove aggregate markup distortion
        
        Ggrid          = funeval(cfun, fspace, Ngrid, 1)./funeval(cfun, fspace, Ngrid); 

    case 'sizedep'
        
        Mugrid         = ones(size(Ngrid))*Musave;    % keep aggregate markup on to isolate role of misallocation
        
        Ggrid          = funeval(cfun, fspace, Ngrid, 1)./funeval(cfun, fspace, Ngrid)./Mugrid; 

end

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

switch experiment
    
    case 'efficient'

  p.xis        = 0; 
  p.xie        = 0;                                    % optimal entry subsidy, leads to 0.00397% welfare gains

  Xnew         = fsolve('findequilibrium', [Y*1.5; N*1.1], options, p);

    case 'sizedep' 
    
  p.xis        = 0; 
  p.xie        = 0; 
  
  Xnew         = fsolve('findequilibrium', [Y*1.05; N*0.85], options, p);
     
end

% Find New Steady State

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

[err, Wnew, Lnew, Lpnew, Knew, Bnew, Cnew, Qnew, Znew, Munew, Gnew]  = findequilibrium(Xnew,   p);

Rnew           = p.r + p.delta; 

% Transition dynamics

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

x1_ss  = [K;    N];                                                                                     save x1_ss x1_ss; 
x2_ss  = [Ynew; Cnew; Bnew; Lnew; Lpnew; Knew; Nnew; Wnew; Rnew; Qnew; Znew; Munew; Gnew];              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; 
Lt(1)  = L; 
Lpt(1) = Lp; 
Kt(1)  = K; 
Nt(1)  = N; 
Wt(1)  = W; 
Rt(1)  = R; 
Zt(1)  = Z; 
Mut(1) = Mu; 

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 7: Implications of Alternative Policies, Oligopoly </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(dW, 2), ...
    'VariableNames', {'Y','C','L','N','K','Z','Welfare,%'}))
