%12/05_ SS running the .mod file at the end. 
clear
clc

% population shares
gammaba = 0.12;
gammabn = 0.29;
gammasa = 0.12;
gammasn = 0.47;

%discount factors
betas = 0.99;
betab = 0.98;

zeta = 0.74; %sector 1 consumption share
eta = 1; %inverse of Frisch elasticity

rss = 1/betas; %gross interest rate 
bettass = betas; %discount factor
kappa = 0.9;   %SS loan to value ratio 

%labour shares targets 
alpha1ba = 0.0175;
alpha1bn = 0.3744;
alpha1sa = 0.0175;
alpha1sn = 0.5906;
alpha2ba = 0.1475;
alpha2bn = 0.2735;
alpha2sa = 0.1475;
alpha2sn = 0.4315;

%Fsolve 
fun = @(x)root2d(x,alpha1ba, alpha1bn, alpha1sa, alpha1sn, alpha2ba, alpha2bn, alpha2sa, alpha2sn, gammaba, gammabn, gammasa, gammasn, betas, zeta);
x0 = [0.6,1.1,1,0.5,0.5,0.5,0.5,1,1,1,1,1,1,1,1,1,1];
A = fsolve(fun,x0);

wass = A(1); 
wnss = A(2);
yss = A(3);
zetaba = A(4);
zetabn = A(5);
zetasa = A(6);
zetasn = A(7);
cbass = A(8);
cbnss = A(9);
csass = A(10);
csnss = A(11);
nbass = A(12);
nbnss = A(13);
nsass = A(14);
nsnss = A(15);
bbass = A(16);
bsass = A(17);

%output and consumption
y1ss = zeta*yss;    
y2ss = (1-zeta)*yss;
c1ss = y1ss; 
c2ss = y2ss;

% borrowing and housing
bbnss = bbass;
bsnss = bsass;
qss = 1; 
hss = bbass/(kappa*betas*qss);

% consumption
c1bass = zeta*cbass; 
c2bass = (1-zeta)*cbass; 
c1sass = zeta*csass; 
c2sass = (1-zeta)*csass;
c1bnss = zeta*cbnss; 
c2bnss = (1-zeta)*cbnss; 
c1snss = zeta*csnss; 
c2snss = (1-zeta)*csnss;

% check equilibrium in good markets
y1_test = gammaba*c1bass + gammabn*c1bnss + gammasa*c1sass + gammasn*c1snss;
y2_test = gammaba*c2bass + gammabn*c2bnss + gammasa*c2sass + gammasn*c2snss;

% Lagranfe multipliers on borrowing constraints
deltabass = (1-betab/betas)*betas/cbass; 
deltabnss = (1-betab/betas)*betas/cbnss; 

% sectoral labour
n1bass = zetaba*nbass;
n2bass = (1-zetaba)*nbass; 
n1bnss = zetabn*nbnss;
n2bnss = (1-zetabn)*nbnss; 
n1sass = zetasa*nsass;
n2sass = (1-zetasa)*nsass; 
n1snss = zetasn*nsnss;
n2snss = (1-zetasn)*nsnss; 

%productivities
A1ss = y1ss/((n1bass^alpha1ba)*(n1sass^alpha1sa)*(n1bnss^alpha1bn)*(n1snss^alpha1sn));
A2ss = y2ss/((n2bass^alpha2ba)*(n2sass^alpha2sa)*(n2bnss^alpha2bn)*(n2snss^alpha2sn));

%weights on housing
omegahsa = ((1-betas)*hss)/csass; 
omegahsn = ((1-betas)*hss)/csnss; 
omegahba = ((1-betab)/cbass - kappa*deltabass)*hss; 
omegahbn = ((1-betab)/cbnss - kappa*deltabnss)*hss; 

%weights on labour
omegansa = wass/(csass*(nsass^eta));
omeganba = wass/(cbass*(nbass^eta));
omegansn = wnss/(csnss*(nsnss^eta));
omeganbn = wnss/(cbnss*(nbnss^eta));

A = [omegahba, omegahbn, omegahsa, omegahsn, omeganba, omeganbn, ...
    omegansa, omegansn, zetaba, zetabn, zetasa, zetasn, alpha1ba, ...
    alpha1bn, alpha1sa, alpha1sn, alpha2ba, alpha2bn, alpha2sa, alpha2sn, ...
    gammaba, gammabn, gammasa, gammasn, betas, betab, zeta, eta,...
    kappa, yss, nbass, nsass, nbnss, nsnss, wass, wnss, qss, hss];

save A A

dynare MM_120521

function F = root2d(x,alpha1ba, alpha1bn, alpha1sa, alpha1sn, alpha2ba, alpha2bn, alpha2sa, alpha2sn, gammaba, gammabn, gammasa, gammasn, betas, zeta)
% labour demand functions
F(1) = alpha1ba - (gammaba*x(1)*x(4)*x(12))/(zeta*x(3));  
F(2) = alpha1bn - (gammabn*x(2)*x(5)*x(13))/(zeta*x(3)); 
F(3) = alpha1sa - (gammasa*x(1)*x(6)*x(14))/(zeta*x(3));
F(4) = alpha1sn - (gammasn*x(2)*x(7)*x(15))/(zeta*x(3)); 
F(5) = alpha2ba - (gammaba*x(1)*(1-x(4))*x(12))/((1-zeta)*x(3));  
F(6) = alpha2bn - (gammabn*x(2)*(1-x(5))*x(13))/((1-zeta)*x(3)); 
F(7) = alpha2sa - (gammasa*x(1)*(1-x(6))*x(14))/((1-zeta)*x(3));
F(8) = alpha2sn - (gammasn*x(2)*(1-x(7))*x(15))/((1-zeta)*x(3)); 
% total output
F(9) = x(3) - gammaba*x(8) - gammabn*x(9) - gammasa*x(10) - gammasn*x(11);
% consumption
F(10) = x(8)  - x(1)*x(12) - (1 - 1/betas)*x(16);  
F(11) = x(9)  - x(2)*x(13) - (1 - 1/betas)*x(16);  
F(12) = x(10) - x(1)*x(14) - (1 - 1/betas)*x(17);  
F(13) = x(11) - x(2)*x(15) - (1 - 1/betas)*x(17);
% profits
F(14) = gammaba*x(1)*x(4)*x(12) + gammabn*x(2)*x(5)*x(13) + gammasa*x(1)*x(6)*x(14) + gammasn*x(2)*x(7)*x(15) - zeta*x(3);
F(15) = gammaba*x(1)*(1-x(4))*x(12) + gammabn*x(2)*(1-x(5))*x(13) + gammasa*x(1)*(1-x(6))*x(14) + gammasn*x(2)*(1-x(7))*x(15) - (1-zeta)*x(3);
% relative wages target
F(16) = x(2) - x(1)*1.723;
% bond equilibrium
F(17) = x(16)*(gammaba+gammabn) + x(17)*(gammasa+gammasn);
end





