clear all
clc

rhon = 0.035;    
beta = 0.99; 
zeta = 0.5; 
zetav = zeta; %hosios condition
nss = 0.9;
alphajss = 0.5;
qss = 0.9;
njss = alphajss*nss;

fun = @(x) root2d(x, qss, nss, njss, alphajss, rhon, beta, zeta);
x0 = [0.5,0.1,0.5,1,1,0.5,1,0.5];
lb = zeros(size(x0)); %lower bound of zero
ub = Inf*ones(size(x0));
options = optimoptions('lsqnonlin','FunctionTolerance',1e-14,'OptimalityTolerance',1e-9);
[A,resnorm,residual,exitflag] = lsqnonlin(fun,x0,lb,ub,options)

sjss = A(2)/A(1);
hjss = qss*A(2); 
fjss = hjss/sjss; 
nu = hjss/(sjss^zetav*A(2)^(1-zetav));
sji = sjss/0.2;

B = [rhon, beta, zeta, qss, nss, njss, alphajss];
C = [sjss, hjss, fjss, nu, sji, zetav];
A = [A,B,C];

save A A

function F = root2d(x, qss, nss, njss, alphajss, rhon, beta, zeta)

% employment
F(1) = rhon*njss - qss*x(2);
% consumption j
F(2) = x(3) - alphajss*x(4);
% budget constraint h
F(3) = x(4) - x(5)*nss;
% production functions
F(4) = x(6) - x(7)*njss;
% vacancy posting
F(5) = x(7) - x(5) - x(8)*((1-beta*(1-rhon))/qss);
% Nash wages
F(6) = x(5) - (1-zeta)*(x(7) + (x(8)*beta*(x(1)*2)));
% equilibrium
F(7) = x(3) - 0.97*x(6);
F(8) = x(8)*x(2) - 0.03*x(6);
end


