clear all
close all
clc

%% Parameters

a = 1;
t = 0.15;
alpha = 0.30;
b = 0.43;
kappa = -0.05; 
rho = 0.95;

%% Variable steady state

u = 0.15;
fraction_public_emp=10/100;

lg=(1-u)*fraction_public_emp;
lp=1-u-lg;



etap=0.5; % Petrongolo & Pissarides (1999, 2000 ? looking into the black box...)
etag=0.5; % Petrongolo & Pissarides (1999, 2000 ? looking into the black box...)
beta= 0.985; % Steady state discount factor quarterly
stochdis=beta;
r=(1 - stochdis) / (stochdis);

lambdap=0.02; % Trouver une cible pour défendre cette valeur
lambdag=0.01; % Trouver une cible pour défendre cette valeur

mp=lp*lambdap;
mg=lg*lambdag;


s=0.5; % Trouver une cible pour défendre cette valeur
usg=s*u;
usp=(1-s)*u;


declared_vacancies=2164;  % number 2012
job_seekers_public=30756;  % number 2012
tightness_g = declared_vacancies/job_seekers_public;
vg=tightness_g*usg;


fraction_public_vac=0.37; % maybe too high (must be checked)
vp=vg/fraction_public_vac-vg;

meffp=mp/(usp^etap*vp^(1-etap));
meffg=mg/(usg^etag*vg^(1-etag));

qp=mp/vp;
qg=mg/vg;
pg=mg/usg;
pp=mp/usp;

fp=mp/u;
fg=mg/u;


z=0.62;
thetap=qp*beta*((1-alpha)*lp^(-alpha)-z/(1-t))*(1-b)/(1-beta*(1-lambdap) +beta*pp*b);
wp=(a*(1-alpha)*(lp)^(-alpha)+thetap*pp/qp)*b+(1-b)*z/(1-t);

Up =(z+thetap/qp*b/(1-b)*(1-t)*pp)/(1-beta);
Ug = Up;
Wp=((1-t)*wp+beta*lambdap*Up)/(1-beta*(1-lambdap));
Psi=(((Ug*(1-beta*(1-pg))-z)/(beta*pg)*(1-beta*(1-lambdag))-beta*lambdag*Ug)/(wp*(1-t)))-1;
wg=wp*(1+Psi);
Wg=((1-t)*wg+beta*lambdag*Ug)/(1-beta*(1-lambdag));

s_new=mg*(Wg-Up)/(mp*(Wp-Up)+mg*(Wg-Up));
thetag=thetap;
y = a*(lp^(1-alpha)- thetap*vp);
g=  a*(lg - thetag*vg); 
c = y+g; 


e_a = 0; 
e_w = 0;
e_v = 0;

%% Model

XX=[-1 + lp + lg + u;    % ok
% 2
-lp + (1-lambdap)*lp + mp;    % ok
% 3
-lg + (1-lambdag)*lg + mg ;   % ok
% 4
-mp + meffp*(usp^(etap))*(vp^(1-etap));    % ok
% 5
-mg + meffg*(usg^(etag))*(vg^(1-etag));    % ok
% 6
-usg + s*u ; % otherwise, s + usg/u, as in the research paper  of Gomes 
% 7
-qp + mp/vp ;   % ok
%8
-qg  + mg/vg ;   % ok
%9
-pp + mp/usp;
% 10
-pg + mg/usg;
% 11
-fp + mp / u ;    % ok the unconditional job-finding rates
%  12 
-fg + mg / u  ;   % ok the unconditional job-finding rates
% 13 
-r + (1 - stochdis) / (stochdis) ;  % ok The Euler Equation could also be written as 1 + beta*(1+r)*(c/c) 
%14
-stochdis + beta*((c/c));
%15
-usp + u - usg ;
%16
-thetap/qp + (stochdis*(a/a))*((1-alpha)*((lp)^(-alpha))- (wp/a) + (1-lambdap)*(thetap/qp)); % Firm's problem solution
%17
-wp + (a*(1-alpha)*(lp)^(-alpha)+thetap*pp/qp)*b+(1-b)*z/(1-t);
% 18 
-Wp + (1-t)*wp + stochdis*((1-lambdap)*Wp + lambdap*Up);    % ok The value of being employed
% 19
-Wg + (1-t)*wg + stochdis*((1-lambdag)*Wg + lambdag*Ug);    % ok The value of being employed
% 20
-Up + z + stochdis*(pp*Wp + (1-pp)*Up);    % ok The value of unemployment
% 21 
-Ug + z + stochdis*(pg*Wg + (1-pg)*Ug);
%22
-(mp*(Wp - Ug)) / (1 - s) +   (mg*(Wg - Ug)) / (s);    % ok The value of unemployment
% 23
-wg/wg + (wp/wp) + kappa*((wg/wp) - 1 - Psi) + e_w; % ok 
% 24
-y + a*(lp^(1-alpha)- (thetap*vp));
%25 
-g + a*(lg - (thetag*vg));
% 26 
-c+y+g;
% 27
-log(a) + rho*log(a) + e_a;
% 28
-vg + vg + e_v];

[[1:length(XX)]' XX]

lpss=lp;
lgss=lg;
uss=u;
mgss=mg;
mpss=mp;
vgss=vg;
vpss=vp;
uspss=usp;
usgss=usg;
qpss=qp;
qgss=qg;
fpss=fp;
fgss=fg;
rss=r;
ppss=pp;
pgss=pg;
wpss=wp;
wgss=wg;
Wpss=Wp;
Wgss=Wg;
Upss=Up;
Ugss=Ug;
Uss=Ug;
yss=y;
gss=g;
css=c;
ass=a;
betass=stochdis;
sss=s;

lambdapss = lambdap;
lambdagss = lambdag;
etapss = etap;
etagss = etag;
meffpss = meffp;
meffgss = meffg;
thetapss = thetap;
thetagss = thetag;
alphass = alpha;
bss = b;
Psiss = Psi;
kappass = kappa; 
betass = beta;
rhoss = rho;
zss=z;
tss=t;


% vecteur_tax=[0.05:0.01:0.2];
% for i=1:length(vecteur_tax)
%     tss=vecteur_tax(i);
%     
dynare newinteraction.mod noclearall
%  u_policy_tax(i)=oo_.steady_state(3);
% end


% plot(vecteur_tax,u_policy_tax)
