clear all
close all

%%%%%%%%%%%%%%%%%%%%%%%
% Simulations with Dynare

ff    = 5000;
beta  = 1/(1+.05/12);

s     = 0.01705;
pss   = 0.2842;
hh    = .295;
gam   = .6;
thets = (pss/hh)^(1/gam);


% Hall
gam   = .5;
s     = 0.031;
pss   = 0.285;
thess = .72;
b     = .7;

% Shimer
gam   = .72;
s     = 0.035;
pss   = 0.6;
thess = .72;
b     = .4;

% HM
% gam   = .5;
% s     = 0.031;
% pss   = 0.285;
% thess = .72;
% b     = .975;

% HM and Shimer
%gam   = .5;
%s     = 0.035;
%pss   = 0.6;
%thess = .72;
%b     = .975;

us    = s/(s+pss);
thets = thess;
vacs  = thets*us;

hh    = pss/(thess^gam);
kappa = (beta * (1 - gam - (1-gam)*b)/(1-beta*(1-s)))*pss/thets;

disp('kappa  hh   q')
disp([kappa  hh   hh*(thess^(gam-1))])

%break

%set_dynare_seed(31415);
dynare('piss_AKL2','noclearall')
dynare_AKL_2=oo_.endo_simul;
save dynare_AKL_2


dynare('piss_AKL1bis','noclearall')
dynare_AKL_1=oo_.endo_simul;
save dynare_AKL_1


unemp_dy2 = 1-dynare_AKL_2(1,:);
jfr_dy2   = dynare_AKL_2(4,:);
vac_dy2   = dynare_AKL_2(5,:);
unemp_dy1 = 1-dynare_AKL_1(1,:);
jfr_dy1   = dynare_AKL_1(4,:);
vac_dy1   = dynare_AKL_1(5,:);


figure
times = 1:1:ff;
plot(times(1:ff),us*ones(1,ff),times(1:ff),mean(unemp_dy2(1,1:ff))*ones(1,ff),...
    times(1:ff),unemp_dy1(1,1:ff),times(1:ff),unemp_dy2(1,1:ff))
xlabel('Time','FontSize',20);
ylabel('Unemployment rate','FontSize',20);
set(gca,'Fontsize',16);
legend('u(ss)','E[u]','Dynare oder 1','Dynare oder 2')

figure
times = 1:1:ff-1;
plot(times(1:ff-1),vacs*ones(1,ff-1),times(1:ff-1),mean(vac_dy2(1,1:ff-1))*ones(1,ff-1),...
    times(1:ff-1),vac_dy1(1,1:ff-1),times(1:ff-1),vac_dy2(1,1:ff-1))
xlabel('Time','FontSize',20);
ylabel('Vacancy rate','FontSize',20);
set(gca,'Fontsize',16);
legend('v(ss)','E[v]','Dynare oder 1','Dynare oder 2')

figure
times = 1:1:ff;
plot(times(1:ff),pss*ones(1,ff),times(1:ff),mean(jfr_dy2(1:ff))*ones(1,ff),...
    times(1:ff),jfr_dy1(1:ff),times(1:ff),jfr_dy2(1:ff))
xlabel('Time','FontSize',20);
ylabel('Job finding rate','FontSize',20);
set(gca,'Fontsize',16);
legend('p(\theta(ss))','E[p(\theta)]','Dynare oder 1','Dynare oder 2')


figure
subplot(221),plot(1:500,us*ones(1,500),1:500,mean(unemp_dy2(1,:))*ones(1,500),...
    1:500,unemp_dy1(1,1501:2000),1:500,unemp_dy2(1,1501:2000))
xlabel('Time','FontSize',16);
ylabel('x = u','FontSize',16);
title('Unemployment rate','FontSize',20)
set(gca,'Fontsize',16);
%legend('u(ss)','E[u]','Dynare oder 1','Dynare oder 2')

subplot(222), plot(1:500,vacs*ones(1,500),1:500,mean(vac_dy2(1,:))*ones(1,500),...
    1:500,vac_dy1(1,1501:2000),1:500,vac_dy2(1,1501:2000))
xlabel('Time','FontSize',16);
ylabel('x = v','FontSize',16);
title('Vacancy rate','FontSize',20)
set(gca,'Fontsize',16);
%legend('v(ss)','E[v]','Dynare oder 1','Dynare oder 2')

subplot(223),plot(1:500,pss*ones(1,500),1:500,mean(jfr_dy2(:))*ones(1,500),...
    1:500,jfr_dy1(1501:2000),1:500,jfr_dy2(1501:2000))
xlabel('Time','FontSize',16);
ylabel('x = jfr','FontSize',16);
title('Job finding rate','FontSize',20)
set(gca,'Fontsize',16);
h = legend('$x(ss)$','$E[x]$','$\frac{x_t - x(ss)}{x(ss)}$ order 1','$\frac{x_t - x(ss)}{x(ss)}$ order 2');
set(h,'Interpreter','latex') 

%break

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Statistics: moments on the simulated data 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

qunemp_dy1 = max(QUARTER(unemp_dy1'),0.001);
qunemp_dy2 = max(QUARTER(unemp_dy2'),0.001);
qjfr_dy1   = max(QUARTER(jfr_dy1'),0);
qjfr_dy2   = max(QUARTER(jfr_dy2'),0);

lambdaHP=1600;%10^5;
[trendt,qunemp_dy1_hp] = HPFILTER(log(qunemp_dy1),lambdaHP);
[trendt,qunemp_dy2_hp] = HPFILTER(log(qunemp_dy2),lambdaHP);
[trendt,qjfr_dy1_hp]   = HPFILTER(log(jfr_dy1),lambdaHP);
[trendt,qjfr_dy2_hp]   = HPFILTER(log(jfr_dy2),lambdaHP);

disp('     u1      u2     jfr1       jfr2')
disp([std(qunemp_dy1_hp(100:end)) std(qunemp_dy2_hp(100:end)) std(qjfr_dy1_hp(100:end)) std(qjfr_dy2_hp(100:end))])

break


%%%%%%%%%%%%%%%%%%%%%%%%
% NL Model : 
% value function iteration 
% Discrete 

% AR 1 on Labour productivity (lambda=10^5)
lambda = .9;%.95;%
sigma  = .009;%.0065;%.015;%.015/2;%.02;%.009;
m      = 3;
N      = 100;%20;

n  = 1000;
ff = 5000;

[Tran,s,probst,alambda,asigmay]=markovappr(lambda,sigma,m,N);
% Tran is the transition matrix of the Markov chain
% s is the discretized state space of y_t
% alambda is the theoretical first order autoregression coefficient for Markov chain
% asigma is the theoretical standard deviation for Markov chain Y_t

Prob = Tran;
Mean = probst;
yref = 1;
grid = s;
y    = (1+grid)*yref;    

%%%%%%%%%%%%%%%%%%%%%%%

F    = 0;%-.3; % Firing costs
r    = .05/12;
beta = 1/(1+r);
gam  = 0.5;
z    = 0.975;%.98;%
eta  = 1-gam;
%s    = 0.0311;
s    = 0.01705;%.02575;%0.0255;

% si rho = .95
%s = .0165;


% wage à la Hall and Milgrom
w  =  eta*yref+(1-eta)*z;

% forward: firm value with the bargaining set
Rj  = y'-w;     
% guess
crit = 1;
J0   = .1*ones(N,1);
JJnew = J0;
while crit>1e-10
    JJ = Rj+max(-F,beta*(1-s)*Prob*JJnew);
    crit = max((JJ-JJnew)./JJnew);
    JJnew = JJ;
end
% forward: firm value without the bargaining set
JJfaux = inv(eye(N)-beta*(1-s)*Prob)*Rj;

figure,
plot(y,JJnew,'-',y,JJfaux)
title('Value functions')
xlabel('productivity: a')
ylabel('Value: J(a)')
legend('with max','no max')

% the expections
Theta  = beta*Prob*JJnew;
Thetaf = beta*Prob*JJfaux;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% An ad hoc calibration to find Ep
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%hh = .105;%.115;
%c  = .107;%.1

hh = .295;%.37;%.115;
c  = .15;%.1


for jj=1:N
    qthet(jj) = c/Theta(jj);
    qthetf(jj) = c/Thetaf(jj);
end

% Search of the grid point at which one of the 3 regimes occurs
N0 = min(find(F+Theta>0));
if N0==1
    N00 = min(find(Theta>0));
    if N00==1
        Ns = min(find(qthet<1));
    else
        Ns = N00+min(find(qthet(N00:N)<1));
    end
else
    N00 = N0+min(find(Theta(N0:N)>0));
    if N00==1
        Ns = min(find(qthet<1));
    else
        Ns = N00+min(find(qthet(N00:N)<1));
    end
end

% idem for the naive model
Nsf = min(find(qthetf<1 & qthetf>0));

% The FOC
if N0>1
for jj=1:N0-1
    thet(jj) = 0;
    sr(jj)   = 1;
    lamb(jj) = c - Theta(jj);
    nu(jj)   = 0;
    phi(jj)  = -F - Theta(jj);
end
for jj=N0:Ns-1
    thet(jj) = 0;
    sr(jj)   = s;
    lamb(jj) = c - Theta(jj);
    nu(jj)   = F + Theta(jj);
    phi(jj)  = 0;
end
for jj=Ns:N
    thet(jj) = ((hh/c)*Theta(jj)).^(1/(1-gam)); %((hh/c)./Theta(jj)).^(1/(1-gam));
    sr(jj)   = s;
    nu(jj)   = F + Theta(jj);
    phi(jj)  = 0;
    lamb(jj) = 0;
end
else
for jj=1:Ns-1
    thet(jj) = 0;
    sr(jj)   = s;
    lamb(jj) = c - Theta(jj);
    nu(jj)   = F + Theta(jj);
    phi(jj)  = 0;
end
for jj=Ns:N
    thet(jj) = ((hh/c)*Theta(jj)).^(1/(1-gam)); %((hh/c)./Theta(jj)).^(1/(1-gam));
    sr(jj)   = s;
    nu(jj)   = F + Theta(jj);
    phi(jj)  = 0;
    lamb(jj) = 0;
end
end

% The FOC in the naive model
for jj=1:N
    thetf(jj) = ((hh/c)*max(0,Thetaf(jj))).^(1/(1-gam)); %((hh/c)./Thetaf(jj)).^(1/(1-gam));
end


% Job find rate for each productivity level
pthet  = hh*thet.^gam;
pthetf = hh*thetf.^gam;
 

figure,
plot(y,max(0,lamb),'-o',y,max(0,nu),'-s',y,max(0,phi),'-*')
xlabel('productivity: a')
ylabel('Multipliers')
legend('\lambda (v \geq 0)','\nu (l \geq 0)','\phi (l \leq n)')


figure,
subplot(221),plot(y,pthet,'-o')
title('Job finding rate with max')
xlabel('productivity: a')
ylabel('JFR: p(\theta(a))')
subplot(222),plot(y,pthetf,'-o')
title('Job finding rate, no max')
xlabel('productivity: a')
ylabel('JFR: p(\theta(a))')
subplot(223),plot(y,sr,'-o')
title('Job separation rate with max')
xlabel('productivity: a')
ylabel('JSR: s(a)')
subplot(224),plot(y,s*ones(N,1),'-o')
title('Job finding rate, no max')
xlabel('productivity: a')
ylabel('JSR: s(a)')


% Etimated job finding rate using the model
Ep   = Mean'*pthet';
Es   = Mean'*sr';
uss  = s/(s+Ep); % implied unemployment rate

% The use of these moment in the detreministic model
pthets = Ep;
thets  = (pthets/hh)^(1/gam);
us     = uss;
vacs   = thets*us;

%break

%%%%%%%%%%%%%%%%%%%%%%%
% Continuous productivity of Dynare toward a discret process

prodchoc = dynare_AKL_2(2,:);
for ii=1:ff%5000
    if prodchoc(ii)<=y(1)
        ydis(ii)  = y(1);
        state(ii) = 1;
    elseif prodchoc(ii)>=y(N)
        ydis(ii)  = y(N);
        state(ii) = N;
    elseif (prodchoc(ii)>y(1)) & (prodchoc(ii)<y(N))
        for jj=2:N
            if (prodchoc(ii)>y(jj-1)) & (prodchoc(ii)<=y(jj))
                ydis(ii)  = y(jj);
                state(ii) = jj;
            end
        end
    end
end

nb = 0; 
for i=1:ff%5000   % iteration on the number of periods
% Model with occ binding constraint
uprime(1) = us;                 % initial level of unemployment
nprime(1) = 1-us;               % initial level of employment
yy(1)     = nprime(1)*mean(y);
% Model without occ binding constraint
uprimef(1) = us;                 % initial level of unemployment
nprimef(1) = 1-us;               % initial level of employment
yyf(1)     = nprimef(1)*mean(y);

  for ii=1:N
    if state(i) == ii
            probthetf(i) = pthetf(ii);
            destf(i)     = s*nprimef(i);
            destrf(i)    = s;
            thetaf(i)    = thetf(ii);
            uprimef(i+1) = s*nprimef(i) + (1-pthetf(ii))*uprimef(i);
            nprimef(i+1) = (1-s)*nprimef(i) + pthetf(ii)*uprimef(i);
            yyf(i+1)     = y(ii)*nprimef(i+1);
            if ii<=N0-1
                uprime(i)   = 1; 
                uprime(i+1) = (1-pthet(ii));
                nprime(i+1) = pthet(ii);
                if nprime(i)>0
                    dest(i)     = nprime(i);
                    destr(i)    = 1;
                else
                    dest(i)     = 0;
                    destr(i)    = 0;
                end 
                nprime(i)   = (1-s)*(nprime(i)-dest(i)); 
                probthet(i) = pthet(ii);
                theta(i)    = thet(ii);
                yy(i+1)     = y(ii)*nprime(i+1);
            else
                probthet(i) = pthet(ii);
                dest(i)     = s*nprime(i);
                destr(i)    = s;
                theta(i)    = thet(ii);
                uprime(i+1) = s*nprime(i) + (1-pthet(ii))*uprime(i);
                nprime(i+1) = (1-s)*nprime(i) + pthet(ii)*uprime(i);
                yy(i+1)     = y(ii)*nprime(i+1);
            end
            if nprime(i) == 0
                ysn(i) = 0;
            else
                ysn(i) = y(ii);
                nb = nb+1;
            end
        end
    end
    vacf(i)     = thetaf(i)*uprimef(i);
    welfcf(i)   = yyf(i) + uprimef(i)*z - c*thetaf(i)*uprimef(i);

    vac(i)     = theta(i)*uprime(i);
    welfc(i)   = yy(i) + uprime(i)*z - c*theta(i)*uprime(i);
    welfcs(i)  = (1-us)*yref + us*z - c*thets*us;        

       

 unemp      = mean(uprime(100:end));
 emp        = mean(nprime(100:end));
 prod       = mean(yy(100:end));
 welfarec   = mean(welfc(100:end));
 welfarecs  = mean(welfcs(100:end));
 vacm       = mean(vac(100:end));
 thetuv     = mean(theta(100:end));
 probthetL  = probthet(2:end);
 stdyy      = std(yy(:));
 stdprob    = std(probthet(:));
 stdprobl   = std(log(probthet(:)));
 stdu       = std(uprime(:));
 aprob      = mean(probthet(:));
 adestr     = mean(destr(:));
 
 unempf      = mean(uprimef(100:end));
 empf        = mean(nprimef(100:end));
 prodf       = mean(yyf(100:end));
 welfarecf   = mean(welfcf(100:end));
 vacmf       = mean(vacf(100:end));
 thetuvf     = mean(thetaf(100:end));
 aprobf      = mean(probthetf(:));
 adestrf     = mean(destrf(:));
 
end

 umean         = unemp;
 nmean         = emp;
 ymean         = prod;
 vmean         = vacm;
 welfarecmean  = welfarec;
 welfarecsmean = welfarecs;
 stdyymean     = stdyy;
 costc         = (welfarecsmean/welfarecmean)-1;
 stdprobmean   = stdprob;
 stdprobmeanl  = stdprobl;
 probmean      = aprob;
 destmean      = adestr;
 stduu         = stdu;
 thetamean     = thetuv;

 umeanf        = unempf;
 nmeanf        = empf;
 ymeanf        = prodf;
 welfarecmeanf = welfarecf;
 vmeanf        = vacmf;
 thetameanf    = thetuvf;
 probmeanf     = aprobf;
 destmeanf     = adestrf;


% distributions 
 
[ll,cc]=size(ysn);
snumysn = ll*cc;
for ii=1:N
    numysn(ii)  = length(find(ysn==y(ii)));
    freqysn(ii) = numysn(ii)/nb;%/snumysn;
    numyy(ii)   = length(find(ydis==y(ii)));
    freqy(ii)   = numyy(ii)/snumysn;
end

incp = .01;
Nump = 1/incp;
for j=1:ff
    p(1,j)     = length(find(uprime<incp));
    pf(1,j)    = length(find(uprimef<incp));
    p(Nump,j)  = length(find(uprime>=1-incp & uprime<=1));
    pf(Nump,j) = length(find(uprimef>=1-incp & uprimef<=1));
    for ii=2:Nump-1
        p(ii,j) = length(find(uprime>=(ii-1)*incp & uprime<ii*incp));
        pf(ii,j) = length(find(uprimef>=(ii-1)*incp & uprimef<ii*incp));
    end
end
[llu,ccu] = size(uprime);
snumu     = llu*ccu;
for i=1:Nump
    pm(i) = mean(p(i,j),2)/snumu;
end
[lluf,ccuf] = size(uprimef);
snumuf     = lluf*ccuf;
for i=1:Nump
    pmf(i) = mean(pf(i,j),2)/snumuf;
end


unemp_dy2 = 1-dynare_AKL_2(1,:);
jfr_dy2   = dynare_AKL_2(4,:);
vac_dy2   = dynare_AKL_2(5,:);
unemp_dy1 = 1-dynare_AKL_1(1,:);
jfr_dy1   = dynare_AKL_1(4,:);
vac_dy1   = dynare_AKL_1(5,:);

jfrss    = pthets;
destss   = s;
Ejfr     = probmean;
Edest    = destmean;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% figures


figure,bar(y,[Mean freqy' freqysn'])
title('Distributions of a')
xlabel('productivity: a')
ylabel('f(a)')
legend('Egodic distribution','empirical distribution','observed distribution')

figure,bar([incp:incp:1],[pm' pmf'])
title('Distributions of u')
xlabel('unemployment: u')
ylabel('f(u)')
legend('With occ. binding C','Without occ. binding C')


figure
times = 1:1:ff;
plot(times(1:ff),uprime(1:ff),times(1:ff),uprimef(1:ff),...
    times(1:ff),us*ones(1,ff),times(1:ff),umean*ones(1,ff),...
    times(1:ff),unemp_dy1(1,1:ff),times(1:ff),unemp_dy2(1,1:ff))
xlabel('Time','FontSize',20);
ylabel('Unemployment rate','FontSize',20);
set(gca,'Fontsize',16);
legend('Model with occ binding C','Model without occ binding C','u(ss)','E[u]','Dynare oder 1','Dynare oder 2')

figure
times = 1:1:ff-1;
plot(times(1:ff-1),vac(1:ff-1),times(1:ff-1),vacf(1:ff-1),...
    times(1:ff-1),vacs*ones(1,ff-1),times(1:ff-1),vmean*ones(1,ff-1),...
    times(1:ff-1),vac_dy1(1,1:ff-1),times(1:ff-1),vac_dy2(1,1:ff-1))
xlabel('Time','FontSize',20);
ylabel('Vacancy rate','FontSize',20);
set(gca,'Fontsize',16);
legend('Model with occ binding C','Model without occ binding C','v(ss)','E[v]','Dynare oder 1','Dynare oder 2')

figure,
[XX,HH1,HH2] = plotyy(times(1:ff-1),probthet(1:ff-1),times(1:ff-1),destr(1:ff-1));
set(get(XX(1),'Ylabel'),'String','JFR (left scale)','Fontsize',16) 
set(get(XX(2),'Ylabel'),'String','JSR (right scale)','Fontsize',16) 
xlabel('Time','FontSize',20);
legend('JFR','JSR')

times = 1:1:ff;
figure
subplot(121),plot(times(1:ff),probthet(1:ff),times(1:ff),probthetf(1:ff),...
    times(1:ff),jfrss*ones(1,ff),times(1:ff),Ejfr*ones(1,ff),...
    times(1:ff),jfr_dy1(1:ff),times(1:ff),jfr_dy2(1:ff))
xlabel('Time','FontSize',20);
ylabel('Job finding rate','FontSize',20);
set(gca,'Fontsize',16);
legend('Model with occ binding C','Model without occ binding C','p(\theta(ss))','E[p(\theta)]','Dynare oder 1','Dynare oder 2')

subplot(122),plot(times(1:ff),destr(1:ff),times(1:ff),destrf(1:ff),...
    times(1:ff),destss*ones(1,ff),times(1:ff),Edest*ones(1,ff),...
    times(1:ff),destss*ones(1,ff),times(1:ff),destss*ones(1,ff))
xlabel('Time','FontSize',20);
ylabel('Job separation rate','FontSize',20);
set(gca,'Fontsize',16);
legend('Model with occ binding C','Model without occ binding C','s','E[s]','Dynare oder 1','Dynare oder 2')


disp('------------------------------------------')
disp('Macro with business cycle with occ binding C')
disp('   emploi chomage  s moyen  p moyen theta moyen')
disp([nmean umean  destmean probmean thetamean]);
disp('Macro with business cycle without occ binding C')
disp('   emploi chomage  s moyen  p moyen theta moyen')
disp([nmeanf umeanf  destmeanf probmeanf thetameanf]);
disp('------------------------------------------')
disp('Macro without business cycle')
disp('   emploi chomage  s moyen  p moyen  theta moyen')
disp([1-us us  s pthets thets]);
disp('------------------------------------------')
disp('std(consumption)   std(u)')
disp([stdyymean  stduu])
disp('------------------------------------------')
disp('Cost of BC E[Prod]  Prodss ')
disp([costc ymean (1-us)*yref]);
disp('------------------------------------------')
disp('vacancies')
disp('volatile and stablised')
disp([vmean  vacs]);

