function [ff,a0,a1,a2,a3,H,wbar,kbar,pmx,pmc,pmi,pc,pinvest,px,pkprime,Rk,z,Gamma,Gam_muG,ybar,we, ...
    E0,E1,E2,E3,F0,F1,F2,F3,intensity0,intensity1,intensity2,intensity3, ...
    G0,G1,G2,G3,M0,M1,M2,M3,Vw,v0,v1,v2,v3,v,L,m,BigQ,l0,l1,l2,l3,vtilde0,vtilde1,vtilde2,vtilde3, ...
    chi0,chi1,chi2,chi3,Om0,Om1,Om2,Om3,Om4,Jw,recruiting,kappa,w,Jzplus,bu,f,Vxzplus,Vzplus1,Vzplus2,Vzplus3,Vzplus0, ...
    eta,ystar, Jtildezplus0, fcallprime0, fcallprime1, fcallprime2, fcallprime3, ...
    gcallprime0, gcallprime1, gcallprime2, gcallprime3, ...
    sigma,sigma_a,omegabar,rkbar,AL,psizplus,Uzplus,delta,varphi] = ...
    ctw_steadystate_solvecore_10(xxx,pibar,mupsi,tauk,muzplus,R,alphaa,Rf,epsilon,lambdad,...
    etag,Rx,etax,omegax,etaf,omegai,etai,Upsilon,zetac,betaa,b,tauc,zetah,sigmaL,omegac,...
    etac,Rnustar,tauw,tauy,etaa,taud,s,phalo,pitildew,lambdamc,mc,phalomi,phalomc,...
    phalomx,Rstar,lambdami,lambdamx,phalox,lambdax,unemp,recruitshare,mu,weshare,...
    Fomegabar,taumx,taumc,taumi,mcmi,mcmx,mcmc,taux,mcx,hours_target,iy_target,nk_target,...
    rho,sigmam,sigma_match,iota,BigN,bshare,F,xy_target,SW,SE,delta)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The model code is based on Lawrence J. Christiano, Mathias Trabandt and Karl Walentin, 2011. Introducing financial frictions and unemployment into a small open economy model, Journal of Economic Dynamics and Control, Elsevier, vol. 35(12), pages 1999-2041.
% Adapted by Ginters Buss for Latvia as a small open economy in a monetary
% union.
% 
% If these codes are useful, cite: 
% Ginters Buss (2015) "Search-and-matching frictions and labor market dynamics in Latvia", Working Papers 2015/04, Latvijas Banka.
% 
% Copyright of this version: Latvijas Banka
% These codes are distributed by Latvijas Banka for academic purposes. 
% If in doubt, contact Latvijas Banka.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%backing out the guesses for the variables to be solved for
sigma_a=xxx(1)*.1708;
omegabar=xxx(2)*.3949;
rkbar=xxx(3)*.1102;
AL=xxx(4)*2.5e5;     % % old 1000
v0=xxx(5)/103;
psizplus=xxx(6)*3.1238;
Uzplus=xxx(7)*23.05;   
%delta=xxx(8)/33;
a0=xxx(8)*.6233;
sigma=xxx(9)*.4112;
varphi=xxx(10)*1.8257;

%due to no steady state wage dispersion:
v1=v0;
v2=v0;
v3=v0;
a1=a0;
a2=a0;
a3=a0;

%equations to be solved
ff=NaN*zeros(size(xxx));

%First, calculate steady states based on guesses. Then, adjust guesses to
%make list of equations below hold with equality
H=hours_target;

%Real Wage
wbar=(mc/(taud*(1/(1-alphaa))^(1-alphaa)*(1/alphaa)^alphaa*rkbar^alphaa*Rf^(1-alphaa)/epsilon))^(1/(1-alphaa));

% Capital labor ratio
kbar=H*muzplus*((taud*mupsi^alphaa*wbar*Rf)/((1-alphaa)*mc))^(1/alphaa);

%Relative import prices
pmx=taumx*Rnustar*varphi/mcmx;
pmc=taumc*Rnustar*varphi/mcmc;
pmi=taumi*Rnustar*varphi/mcmi;

%Relative consumption price
pc=(1-omegac+omegac*pmc^(1-etac))^(1/(1-etac));

%Relative investment price
pinvest=(1-omegai+omegai*pmi^(1-etai))^(1/(1-etai));

%Relative export price
px=(taux*Rx/(mcx*varphi))*(omegax*pmx^(1-etax)+1-omegax)^(1/(1-etax));

%Relative price of capital
pkprime=pinvest/Upsilon;

%Rental rate on capital
Rk=(rkbar*(1-tauk)+(1-delta)*pkprime+tauk*delta*mupsi*pkprime/pibar)/(mupsi*pkprime/pibar);

%some abbreviations
z       =   (log(omegabar)+sigma^2/2)/sigma;
Gamma   =   normcdf(z-sigma)+omegabar*(1-normcdf(z));
Gam_muG =   (1-mu)*normcdf(z-sigma)+omegabar*(1-normcdf(z));

%output
ybar=(phalo^(lambdad/(lambdad-1))*epsilon*(kbar/H/(mupsi*muzplus))^alphaa)...
    /(1+phalo^(lambdad/(lambdad-1))*(1/mc-phalo^(lambdad/(1-lambdad))))*H;

%endowment of entrepreneurs
we=weshare*ybar;

%calligraphic epsilons
E0=1-normcdf( (log(a0)+0.5*sigma_a^2)/sigma_a - sigma_a);
E1=1-normcdf( (log(a1)+0.5*sigma_a^2)/sigma_a - sigma_a);
E2=1-normcdf( (log(a2)+0.5*sigma_a^2)/sigma_a - sigma_a);
E3=1-normcdf( (log(a3)+0.5*sigma_a^2)/sigma_a - sigma_a);

%calligraphic F
F0=normcdf( (log(a0)+0.5*sigma_a^2)/sigma_a );
F1=normcdf( (log(a1)+0.5*sigma_a^2)/sigma_a );
F2=normcdf( (log(a2)+0.5*sigma_a^2)/sigma_a );
F3=normcdf( (log(a3)+0.5*sigma_a^2)/sigma_a );

%intensities
intensity0=( psizplus*wbar*((1-tauy)/(1+tauw))*(E0/(1-F0))/(zetah*AL) )^(1/sigmaL);
intensity1=( psizplus*wbar*((1-tauy)/(1+tauw))*(E1/(1-F1))/(zetah*AL) )^(1/sigmaL);
intensity2=( psizplus*wbar*((1-tauy)/(1+tauw))*(E2/(1-F2))/(zetah*AL) )^(1/sigmaL);
intensity3=( psizplus*wbar*((1-tauy)/(1+tauw))*(E3/(1-F3))/(zetah*AL) )^(1/sigmaL);

%G updater
G0=1;
G1=pitildew/pibar/muzplus;
G2=pitildew^2/pibar^2/muzplus^2;
G3=pitildew^3/pibar^3/muzplus^3;

%calligraphic M updater
M0=(1-F0);
M1=M0*(1-F1);
M2=M1*(1-F2);
M3=M2*(1-F3);

%Vw - marginal value to worker of a change in the wage
Vw=(1-tauy)/(1+tauw)*(G0*M0*intensity0+betaa*rho*G1*M1*intensity1+(betaa*rho)^2*G2*M2*intensity2+(betaa*rho)^3*G3*M3*intensity3);

%total vacancies
v=v0+v1+v2+v3;

%total employment
L=1-unemp;

%matching function
m=sigmam*(1-L)^sigma_match*v^(1-sigma_match);

%probability of filling a vacancy
BigQ=m/v;

%calculate employment per cohort
l0=BigQ*v0/(1-rho*(1-F0));
l1=l0;
l2=l0;
l3=l0;

%vtilde's - vacancy rates
vtilde0=BigQ^iota*v0/((1-F0)*l0);
vtilde1=BigQ^iota*v1/((1-F1)*l1);
vtilde2=BigQ^iota*v2/((1-F2)*l2);
vtilde3=BigQ^iota*v3/((1-F3)*l3);

%chi's - new arrivals
chi0=BigQ^(1-iota)*vtilde0;
chi1=BigQ^(1-iota)*vtilde1;
chi2=BigQ^(1-iota)*vtilde2;
chi3=BigQ^(1-iota)*vtilde3;

%Omega updater
Om0=1-F0;
Om1=(1-F1)*(chi0+rho)*(1-F0);   % equivalently Om1=(1-F1)*Om0*(chi0+rho)
Om2=(1-F2)*Om1*(chi1+rho);
Om3=(1-F3)*Om2*(chi2+rho);
Om4=(1-F0)*Om3*(chi3+rho);
% note that the above implies that all Omega are equal, and the value is
% 1-F0

%in the case of union bargaining, set MONOP=1
MONOP=0;
if MONOP == 0
    OOm0=1-F0;
    OOm1=(1-F1)*OOm0*rho;
    OOm2=(1-F2)*OOm1*rho;
    OOm3=(1-F3)*OOm2*rho;
    OOm4=(1-F0)*OOm3*rho;
else
    OOm0=Om0;
    OOm1=Om1;
    OOm2=Om2;
    OOm3=Om3;
end

Jw=-(G0*OOm0*intensity0+betaa*G1*OOm1*intensity1+betaa^2*G2*OOm2*intensity2+betaa^3*G3*OOm3*intensity3);

%calculate kappa for given recruitmentshare
recruiting=recruitshare*ybar;
kappa=recruiting*2/(vtilde0^2*l0*(1-F0)+vtilde1^2*l1*(1-F1)+vtilde2^2*l2*(1-F2)+vtilde3^2*l3*(1-F3));


%Nash wage in steady state (FOC for vacancies for j=N-1 used here)
w=(wbar*E0*intensity0-kappa*vtilde3/(BigQ^(1-iota))/betaa+kappa*(1-F0)*(vtilde0^2/2+vtilde0*rho/(BigQ^(1-iota))))/(intensity0*wbar*(1-F0));


%surplus of agency
Jzplus=1/(1-betaa^BigN*Om4/(1-F0))*(((wbar*(E0)/(1-F0)-G0*w*wbar)*intensity0 - kappa*vtilde0^2/2)*Om0 ...
    +((wbar*(E1)/(1-F1)-G1*w*wbar)*intensity1 - kappa*vtilde1^2/2)*Om1*betaa ...
    +((wbar*(E2)/(1-F2)-G2*w*wbar)*intensity2 - kappa*vtilde2^2/2)*Om2*betaa*betaa...
    +((wbar*(E3)/(1-F3)-G3*w*wbar)*intensity3 - kappa*vtilde3^2/2)*Om3*betaa^3);

%average bshare (note that Gj and intensity_j may be different across
%cohorts. That's why we calc avg. bshare
bu=bshare/(1-tauy)*( (1-tauy)/(1+tauw)*w*wbar*intensity0*G0-zetah*AL*intensity0^(1+sigmaL)/((1+sigmaL)*psizplus)+...
    (1-tauy)/(1+tauw)*w*wbar*intensity1*G1-zetah*AL*intensity1^(1+sigmaL)/((1+sigmaL)*psizplus)+...
    (1-tauy)/(1+tauw)*w*wbar*intensity2*G2-zetah*AL*intensity2^(1+sigmaL)/((1+sigmaL)*psizplus)+...
    (1-tauy)/(1+tauw)*w*wbar*intensity3*G3-zetah*AL*intensity3^(1+sigmaL)/((1+sigmaL)*psizplus))/BigN;

%Prob. of finding a job
f=m/(1-L);

% Rewritten equation that relates value of being unemployed, Uzplus to
% value of beeing newly employed, Vxzplus
Vxzplus=(Uzplus*(1-betaa*(1-f))-bu*(1-tauy))/(betaa*f);

%Vzplus equations;
Vzplus0=(G3*(1-tauy)/(1+tauw)*w*wbar*intensity3...
    -zetah*AL*intensity3^(1+sigmaL)/((1+sigmaL)*psizplus)...
    +betaa*((1-rho+rho*F0)*Uzplus))/(1-betaa*rho*(1-F0));
Vzplus1=Vzplus0;
Vzplus2=Vzplus0;
Vzplus3=Vzplus0;

%bargaining power
eta=-(Vzplus0-Uzplus)*Jw/((Vw*Jzplus-(Vzplus0-Uzplus)*Jw));

%equations for the endogenous cutoff level
%first some helpful definitions
% Note that symmetry betw cohorts (based on full indexation) is
% assumed
fcallprime0=1/(a0*sigma_a*sqrt(2*3.141592653589793))*exp(-(log(a0)+0.5*sigma_a^2)^2/(2*sigma_a^2));
fcallprime1=fcallprime0;
fcallprime2=fcallprime0;
fcallprime3=fcallprime0;
gcallprime0=( -a0*fcallprime0*(1-F0)+E0*fcallprime0 ) / (1-F0)^2;
gcallprime1=gcallprime0;
gcallprime2=gcallprime0;
gcallprime3=gcallprime0;

intensity0_prime=1/sigmaL*intensity0^(1-sigmaL)*w*wbar*psizplus*(1-tauy)/(1+tauw)*gcallprime0/AL;  % ss version of eq "diffzeta_scaled"
% NOTE: Above WITH taxes

% implemented for symmetric cohorts only
Jtildezplus0=Jzplus/(1-F0);

% equivalently_1: Jtildezplus0=Jzplus*(chi0+rho);

% % equivalently_2:
% Jtildezplus0=wbar*(E0/(1-F0)-w)*intensity0-kappa/2*vtilde0^2 ...
% +betaa/(1-F0)*(...
%         (wbar*(E1/(1-F1)-G1*w)*intensity1-kappa/2*vtilde1^2)*Om1...
% +betaa*  (wbar*(E2/(1-F2)-G2*w)*intensity2-kappa/2*vtilde2^2)*Om2...
% +betaa^2*(wbar*(E3/(1-F3)-G3*w)*intensity3-kappa/2*vtilde3^2)*Om3...
% +betaa^3*Jzplus*Om4*(chi0+rho));

% equivalently_3: Jtildezplus0=((wbar*(E0)/(1-F0)-G0*w*wbar)*intensity0 - kappa*vtilde0^2/2)+betaa*Jzplus*(chi0+rho);  % from eq. on appendix page 29, obtained by plugging in eq. 4.19 into eq. B.60


% *** replace w new endog. separation eq's ****
% D1zplus=betaa*((wbar*E2-G2*wbar*w*(1-F2))*intensity2-kappa/2*vtilde2^2*(1-F2))*(chi1+rho)*(1-F1)+...
%     betaa^2*((wbar*E3-G3*wbar*w*(1-F3))*intensity3-kappa/2*vtilde3^2*(1-F3))*(chi1+rho)*(chi2+rho)*(1-F1)*(1-F2)+...
%     betaa^3*Jzplus*(chi1+rho)*(chi2+rho)*(chi3+rho)*(1-F1)*(1-F2)*(1-F3);
% 
% D2zplus=betaa*((wbar*E3-G3*wbar*w*(1-F3))*intensity3-kappa/2*vtilde3^2*(1-F3))*(chi2+rho)*(1-F2)+...
%     betaa^2*Jzplus*(chi2+rho)*(chi3+rho)*(1-F2)*(1-F3);
% 
% D3zplus=betaa*Jzplus*(chi3+rho)*(1-F3);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%LIST OF EQUATIONS TO BE SOLVED FOR%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%zero profits of banks
ff(1)=-nk_target+(1-Rk/R*Gam_muG);

%Feasibility
ystar=((1-etag)*ybar...    %(1-etag)*y
    -(1-omegai)*pinvest^(etai)*kbar/H*H*(1-(1-delta)/(muzplus*mupsi))/Upsilon...%investment
    -(1-omegac)*(pc^etac)...
    *(zetac*(muzplus-betaa*b)/(pc*(1+tauc)*(muzplus-b))...
    /psizplus)...
    -recruitshare*ybar...
    -(Gamma-Gam_muG)*Rk*pkprime*kbar/(pibar*muzplus))/((omegax*pmx^(1-etax)+1-omegax)^(etax/(1-etax))*(1-omegax)*phalox^(-lambdax/(lambdax-1))*px^(-etaf));%financial frictions monitoring costs

%Trade balance
ff(2)=(1-Rstar*s/pibar/muzplus)*etaa*phalo^(lambdad/(lambdad-1))*epsilon*(kbar/H/(mupsi*muzplus))^alphaa...%assets
    /(1+phalo^(lambdad/(lambdad-1))*(1/mc-phalo^(lambdad/(1-lambdad))))*H...
    +varphi*Rnustar*(... %imports
    omegac*(pc/pmc)^(etac)*phalomc^(lambdamc/(1-lambdamc))...%import consumption
    *(zetac*(muzplus-betaa*b)/(pc*(1+tauc)*(muzplus-b))...
    /psizplus)...
    +phalomi^(lambdami/(1-lambdami))*omegai*(pinvest/pmi)^(etai)*kbar/H*H*(1-(1-delta)/(muzplus*mupsi))/Upsilon...%import investment
    +phalomx^(lambdamx/(1-lambdamx))*omegax*(1/pmx*(omegax*pmx^(1-etax)+1-omegax)^(1/(1-etax)))^etax...%import export
    *phalox^(-lambdax/(lambdax-1))*px^(-etaf)*ystar )...
    -varphi*px*px^(-etaf)*ystar; %exports

%FOC for loan contract
ff(3)=(1-Gamma)*Rk/R+(1-Fomegabar)/(1-Fomegabar-mu*omegabar*lognpdf(omegabar,-sigma^2/2,sigma))*(Rk/R*Gam_muG-1);

%total employment
ff(4)=L-l0*(1-F0)-l1*(1-F1)-l2*(1-F2)-l3*(1-F3);

%total effective hours
ff(5)=H-intensity0*E0*l0-intensity1*E1*l1-intensity2*E2*l2-intensity3*E3*l3;

%cutoff determination
% old version: ff(6)=Jzplus-wbar*(E0-a0*(1-F0))*intensity0; % +(1-F0)*(Vzplus0-Uzplus)

% new: eeq90_ss, FOC for cutoff a0
% Below equation only hold without ss wage dispersion, symmetric cohorts
ff(6)= SW*(w*wbar*(1-tauy)/(1+tauw)-AL*intensity0^sigmaL/psizplus)*intensity0_prime +...
SE*wbar*((E0/(1-F0)-w)*intensity0_prime + gcallprime0*intensity0)...        
- (SW*(Vzplus0 - Uzplus) + SE*Jtildezplus0 )*fcallprime0 / (1-F0);

%Vxzplus equation
% old pre-October 28th correction:
% ff(7)=-Vxzplus+chi0*l0*(1-F0)/m*Vzplus1+chi1*l1*(1-F1)/m*Vzplus2+chi2*l2*(1-F2)/m*Vzplus3+chi3*l3*(1-F3)/m*Vzplus0;
% Below equation only hold without ss wage dispersion, so that F0=F1=F2=F3
ff(7)=-Vxzplus+(1-F0)*(chi0*l0*(1-F0)/m*Vzplus1+chi1*l1*(1-F1)/m*Vzplus2+chi2*l2*(1-F2)/m*Vzplus3+chi3*l3*(1-F3)/m*Vzplus0)+F0*Uzplus;

%match i/y to desired target
%i=kbar*(1-(1-delta)/(muzplus*mupsi))/Upsilon;
%ff(8)=i*pinvest/ybar-iy_target;

%calc a0 subjet to normcdf, cutoff and standard deviation
ff(8)=normcdf((log(a0)+.5*sigma_a^2)/sigma_a)-F;

%Bancruptcy rate requirement for omega distribution
ff(9)=logncdf(omegabar,-sigma^2/2,sigma)-Fomegabar;

%match xy to desired target
% old incorrect version:
%ff(11)=px^(-etaf)*ystar*px/ybar-xy_target; 
ff(10)=px^(-etaf)*ystar*px*varphi/ybar-xy_target; 

