%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%computation of the steady state
function [ys,check] = full_with_vacancy_steadystate(ys,junk)

global M_ check;
persistent LastSuccessFsolve;

% this ALMOST works, odd error msg "Too many output arguments"
% nrparams=size(M_.params);
% for cnt=1:nrparams
%     eval([M_.param_names(cnt,1:end) '= M_.params(' int2str(cnt) ')'])
% end

%read parameter values from M_
alphaa              	=	M_.params(	1	)	;
kappad             	=	M_.params(	2	)	;
betaa               	=	M_.params(	3	)	;
rhopi              	=	M_.params(	4	)	;
xid                	=	M_.params(	5	)	;
xix                	=	M_.params(	6	)	;
kappax             	=	M_.params(	7	)	;
ximc               	=	M_.params(	8	)	;
etaf               	=	M_.params(	9	)	;
ximx               	=	M_.params(	10	)	;
kappamc            	=	M_.params(	11	)	;
kappamx            	=	M_.params(	12	)	;
kappami            	=	M_.params(	13	)	;
ximi               	=	M_.params(	14	)	;
omegac             	=	M_.params(	15	)	;
etac               	=	M_.params(	16	)	;
omegai             	=	M_.params(	17	)	;
etai               	=	M_.params(	18	)	;
Spp                	=	M_.params(	19	)	;
muzplus            	=	M_.params(	20	)	;
mupsi              	=	M_.params(	21	)	;
delta              	=	M_.params(	22	)	;
b                  	=	M_.params(	23	)	;
sigmab             	=	M_.params(	24	)	;
sigmaa             	=	M_.params(	25	)	;
phitildea          	=	M_.params(	26	)	;
a                  	=	M_.params(	27	)	;
phitildes          	=	M_.params(	28	)	;
rhoR               	=	M_.params(	29	)	;
rpi                	=	M_.params(	30	)	;
ry                 	=	M_.params(	31	)	;
rq                 	=	M_.params(	32	)	;
rdeltapi           	=	M_.params(	33	)	;
rdeltay            	=	M_.params(	34	)	;
etax               	=	M_.params(	35	)	;
omegax             	=	M_.params(	36	)	;
phi                	=	M_.params(	37	)	;
rhomupsi           	=	M_.params(	38	)	;
epsilon            	=	M_.params(	39	)	;
pibar              	=	M_.params(	40	)	;
Upsilon            	=	M_.params(	41	)	;
zetac              	=	M_.params(	42	)	;
tauc               	=	M_.params(	43	)	;
tauk               	=	M_.params(	44	)	;
R                  	=	M_.params(	45	)	;
tauy               	=	M_.params(	46	)	;
tauw               	=	M_.params(	47	)	;
zetah              	=	M_.params(	48	)	;
nustar             	=	M_.params(	49	)	;
nuf                	=	M_.params(	50	)	;
nux                	=	M_.params(	51	)	;
etag               	=	M_.params(	52	)	;
lambdad            	=	M_.params(	53	)	;
lambdax            	=	M_.params(	54	)	;
lambdamc           	=	M_.params(	55	)	;
lambdami           	=	M_.params(	56	)	;
lambdamx           	=	M_.params(	57	)	;
wbar               	=	M_.params(	58	)	;
pic                	=	M_.params(	59	)	;
psizplus           	=	M_.params(	60	)	;
H                  	=	M_.params(	61	)	;
q                  	=	M_.params(	62	)	;
pimi               	=	M_.params(	63	)	;
mcmi               	=	M_.params(	64	)	;
mc                 	=	M_.params(	65	)	;
pix                	=	M_.params(	66	)	;
mcmx               	=	M_.params(	67	)	;
mcx                	=	M_.params(	68	)	;
pimx               	=	M_.params(	69	)	;
muz                	=	M_.params(	70	)	;
ybar               	=	M_.params(	71	)	;
u                  	=	M_.params(	72	)	;
kappaw             	=	M_.params(	73	)	;
AL                 	=	M_.params(	74	)	;
sigmaL             	=	M_.params(	75	)	;
Rk                 	=	M_.params(	76	)	;
varphi             	=	M_.params(	77	)	;
ystar              	=	M_.params(	78	)	;
mcmc               	=	M_.params(	79	)	;
px                 	=	M_.params(	80	)	;
pimc               	=	M_.params(	81	)	;
pmc                	=	M_.params(	82	)	;
pmi                	=	M_.params(	83	)	;
kbar               	=	M_.params(	84	)	;
i                  	=	M_.params(	85	)	;
c                  	=	M_.params(	86	)	;
pc                 	=	M_.params(	87	)	;
pkprime            	=	M_.params(	88	)	;
pinvest            	=	M_.params(	89	)	;
pmx                	=	M_.params(	90	)	;
pii                	=	M_.params(	91	)	;
pistar             	=	M_.params(	92	)	;
Rstar              	=	M_.params(	93	)	;
rhomuz             	=	M_.params(	94	)	;
rhoepsilon         	=	M_.params(	95	)	;
rhoUpsilon         	=	M_.params(	96	)	;
rhozetac           	=	M_.params(	97	)	;
rhozetah           	=	M_.params(	98	)	;
rhog               	=	M_.params(	99	)	;
rhotauy            	=	M_.params(	100	)	;
rhophitilde        	=	M_.params(	101	)	;
phitilde           	=	M_.params(	102	)	;
sig_muz            	=	M_.params(	103	)	;
sig_epsilon        	=	M_.params(	104	)	;
sig_upsilon        	=	M_.params(	105	)	;
sig_zetac          	=	M_.params(	106	)	;
sig_zetah          	=	M_.params(	107	)	;
sig_phitilde       	=	M_.params(	108	)	;
sig_epsR           	=	M_.params(	109	)	;
sig_pitarget       	=	M_.params(	110	)	;
sig_g              	=	M_.params(	111	)	;
sig_tauy           	=	M_.params(	112	)	;
s                  	=	M_.params(	113	)	;
Stilde             	=	M_.params(	114	)	;
Sprimetilde        	=	M_.params(	115	)	;
rkbar              	=	M_.params(	116	)	;
Rf                 	=	M_.params(	117	)	;
Rnustar            	=	M_.params(	118	)	;
x                  	=	M_.params(	119	)	;
Rx                 	=	M_.params(	120	)	;
varkappad          	=	M_.params(	121	)	;
varkappax          	=	M_.params(	122	)	;
varkappamx         	=	M_.params(	123	)	;
varkappami         	=	M_.params(	124	)	;
varkappamc         	=	M_.params(	125	)	;
varkappaw          	=	M_.params(	126	)	;
thetaw             	=	M_.params(	127	)	;
pibreve            	=	M_.params(	128	)	;
pitildedpi         	=	M_.params(	129	)	;
Kd                 	=	M_.params(	130	)	;
Fd                 	=	M_.params(	131	)	;
phalo              	=	M_.params(	132	)	;
pitildexpi         	=	M_.params(	133	)	;
Kx                 	=	M_.params(	134	)	;
Fx                 	=	M_.params(	135	)	;
phalox             	=	M_.params(	136	)	;
pitildemcpi        	=	M_.params(	137	)	;
Kmc                	=	M_.params(	138	)	;
Fmc                	=	M_.params(	139	)	;
phalomc            	=	M_.params(	140	)	;
cm                 	=	M_.params(	141	)	;
pitildemipi        	=	M_.params(	142	)	;
Kmi                	=	M_.params(	143	)	;
Fmi                	=	M_.params(	144	)	;
phalomi            	=	M_.params(	145	)	;
im                 	=	M_.params(	146	)	;
pitildemxpi        	=	M_.params(	147	)	;
Kmx                	=	M_.params(	148	)	;
Fmx                	=	M_.params(	149	)	;
phalomx            	=	M_.params(	150	)	;
xm                 	=	M_.params(	151	)	;
aofu               	=	M_.params(	152	)	;
taud               	=	M_.params(	153	)	;
taux               	=	M_.params(	154	)	;
taumc              	=	M_.params(	155	)	;
taumi              	=	M_.params(	156	)	;
taumx              	=	M_.params(	157	)	;
rhotaud            	=	M_.params(	158	)	;
rhotaux            	=	M_.params(	159	)	;
rhotaumc           	=	M_.params(	160	)	;
rhotaumi           	=	M_.params(	161	)	;
rhotaumx           	=	M_.params(	162	)	;
piw                	=	M_.params(	163	)	;
taub               	=	M_.params(	164	)	;
pitarget           	=	M_.params(	165	)	;
etaa               	=	M_.params(	166	)	;
a11                	=	M_.params(	167	)	;
a12                	=	M_.params(	168	)	;
a13                	=	M_.params(	169	)	;
a21                	=	M_.params(	170	)	;
a22                	=	M_.params(	171	)	;
a23                	=	M_.params(	172	)	;
a24                	=	M_.params(	173	)	;
a31                	=	M_.params(	174	)	;
a32                	=	M_.params(	175	)	;
a33                	=	M_.params(	176	)	;
a34                	=	M_.params(	177	)	;
sig_ystar          	=	M_.params(	178	)	;
c21                	=	M_.params(	179	)	;
sig_pistar         	=	M_.params(	180	)	;
c24                	=	M_.params(	181	)	;
c31                	=	M_.params(	182	)	;
c32                	=	M_.params(	183	)	;
sig_Rstar          	=	M_.params(	184	)	;
c34                	=	M_.params(	185	)	;
we                 	=	M_.params(	186	)	;
mu                 	=	M_.params(	187	)	;
Fomegabar          	=	M_.params(	188	)	;
n                  	=	M_.params(	189	)	;
omegabar           	=	M_.params(	190	)	;
sigma              	=	M_.params(	191	)	;
gammaa              	=	M_.params(	192	)	;
spread             	=	M_.params(	193	)	;
rhosigma           	=	M_.params(	194	)	;
rhogamma           	=	M_.params(	195	)	;
recruitshare       	=	M_.params(	196	)	;
bshare             	=	M_.params(	197	)	;
rhosigmam          	=	M_.params(	198	)	;
rhoeta             	=	M_.params(	199	)	;
BigN               	=	M_.params(	200	)	;
rho                	=	M_.params(	201	)	;
iota               	=	M_.params(	202	)	;
sigma_match        	=	M_.params(	203	)	;
kappa              	=	M_.params(	204	)	;
eta                	=	M_.params(	205	)	;
bu                 	=	M_.params(	206	)	;
Jw                 	=	M_.params(	207	)	;
Vw                 	=	M_.params(	208	)	;
BigQ               	=	M_.params(	209	)	;
v                  	=	M_.params(	210	)	;
m                  	=	M_.params(	211	)	;
sigmam             	=	M_.params(	212	)	;
f                  	=	M_.params(	213	)	;
w                  	=	M_.params(	214	)	;
Jzplus             	=	M_.params(	215	)	;
Vxzplus            	=	M_.params(	216	)	;
Uzplus             	=	M_.params(	217	)	;
recruiting         	=	M_.params(	218	)	;
L                  	=	M_.params(	219	)	;
pitildew           	=	M_.params(	220	)	;
G1                 	=	M_.params(	221	)	;
G2                 	=	M_.params(	222	)	;
G3                 	=	M_.params(	223	)	;
chi0               	=	M_.params(	224	)	;
chi1               	=	M_.params(	225	)	;
chi2               	=	M_.params(	226	)	;
chi3               	=	M_.params(	227	)	;
Om1                	=	M_.params(	228	)	;
Om2                	=	M_.params(	229	)	;
Om3                	=	M_.params(	230	)	;
Om4                	=	M_.params(	231	)	;
vtilde0            	=	M_.params(	232	)	;
vtilde1            	=	M_.params(	233	)	;
vtilde2            	=	M_.params(	234	)	;
vtilde3            	=	M_.params(	235	)	;
Vzplus0            	=	M_.params(	236	)	;
Vzplus1            	=	M_.params(	237	)	;
Vzplus2            	=	M_.params(	238	)	;
Vzplus3            	=	M_.params(	239	)	;
l0                 	=	M_.params(	240	)	;
l1                 	=	M_.params(	241	)	;
l2                 	=	M_.params(	242	)	;
l3                 	=	M_.params(	243	)	;
recruitsharePercent	=	M_.params(	244	)	;
SppScaled          	=	M_.params(	245	)	;
muzPercent         	=	M_.params(	246	)	;
mupsiPercent       	=	M_.params(	247	)	;
work_cap_para      	=	M_.params(	248	)	;
wesharePercent     	=	M_.params(	249	)	;
hours_target       	=	M_.params(	250	)	;
iy_target          	=	M_.params(	251	)	;
xy_target          	=	M_.params(	252	)	;
nk_target          	=	M_.params(	253	)	;
avg_Nash           	=	M_.params(	254	)	;
ylessd             	=	M_.params(	255	)	;
sigma_a            	=	M_.params(	256	)	;
epscall0           	=	M_.params(	257	)	;
epscall1           	=	M_.params(	258	)	;
epscall2           	=	M_.params(	259	)	;
epscall3           	=	M_.params(	260	)	;
fcall0             	=	M_.params(	261	)	;
fcall1             	=	M_.params(	262	)	;
fcall2             	=	M_.params(	263	)	;
fcall3             	=	M_.params(	264	)	;
intensity0         	=	M_.params(	265	)	;
intensity1         	=	M_.params(	266	)	;
intensity2         	=	M_.params(	267	)	;
intensity3         	=	M_.params(	268	)	;
mcall0             	=	M_.params(	269	)	;
mcall1             	=	M_.params(	270	)	;
mcall2             	=	M_.params(	271	)	;
mcall3             	=	M_.params(	272	)	;
a0                 	=	M_.params(	273	)	;
a1                 	=	M_.params(	274	)	;
a2                 	=	M_.params(	275	)	;
a3                 	=	M_.params(	276	)	;
SW                 	=	M_.params(	277	)	;
SE                 	=	M_.params(	278	)	;
gcallprime0        	=	M_.params(	279	)	;
gcallprime1        	=	M_.params(	280	)	;
gcallprime2        	=	M_.params(	281	)	;
gcallprime3        	=	M_.params(	282	)	;
fcallprime0        	=	M_.params(	283	)	;
fcallprime1        	=	M_.params(	284	)	;
fcallprime2        	=	M_.params(	285	)	;
fcallprime3        	=	M_.params(	286	)	;
Jtildezplus0       	=	M_.params(	287	)	;
Jtildezplus1       	=	M_.params(	288	)	;
Jtildezplus2       	=	M_.params(	289	)	;
Jtildezplus3       	=	M_.params(	290	)	;
rhosigma_a         	=	M_.params(	291	)	;
sig_sigma_a        	=	M_.params(	292	)	;
BigF               	=	M_.params(	293	)	;
BigFPercent        	=	M_.params(	294	)	;
ab                  =	M_.params(	295	)	;
sigmaLScaled        =	M_.params(	296	)	;
rRstar        =	M_.params(	297	)	;

%a check, since steady state wage dispersion not yet implemented
if abs(BigF) > .1e-10;
    ab=1;
    if abs(pibreve-pitarget) > .1e-10;
        error('fatal (ctw) must have pibreve equal to pitarget in this case. Steady state wage dispersion not yet implemented.');
    end;
    if abs(thetaw-1) > .1e-10;
        error('fatal (ctw) must have thetaw = 1 in this case. Steady state wage dispersion not yet implemented.');
    end
end

% set all working capital parameters to the one common estimated value
nuf=work_cap_para;
nustar=work_cap_para;
nux=work_cap_para;

%calculate no-indexation parameters
varkappad=1-kappad;
varkappax=1-kappax;
varkappamx=1-kappamx;
varkappami=1-kappami;
varkappamc=1-kappamc;
varkappaw=1-kappaw;

%some simple transformations
muz=1+muzPercent/100;
%mupsi=1+mupsiPercent/100;
Spp=SppScaled*10;
recruitshare=recruitsharePercent/100;
weshare=wesharePercent/100;
BigF= BigFPercent/100;

sigmaL=sigmaLScaled*10;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Steady State Calculations%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get mupsi residually, so that muzplus is always the same (and must be given in .mod-file
mupsi=(muzplus/muz)^((1-alphaa)/alphaa);
% Balanced growth factor
% *** check that following holds: muzplus=mupsi^(alphaa/(1-alphaa))*muz;

%Inflation rates
pic=pitarget;
pibar=pic;
pimx=pibar;
pimc=pibar;
pimi=pibar;
pii=pibar/mupsi;
piw=muzplus*pibar;
pix=pistar;

%Exchange rate growth and return on capital
s=pic/pistar;

% Nominal Interest Rate
R=(muzplus/betaa-taub)/(1-taub)*pibar;

%World Nominal Interest Rate
Rstar=(muzplus/betaa-taub)/(1-taub)*pistar;

%Working Capital Interest Rates
Rnustar=nustar*Rstar+1-nustar;
Rf=nuf*R+1-nuf;
Rx=nux*R+1-nux;

%Nonlinear pricing domestic intermediate
pitilded=pibar^kappad*pitarget^(1-kappad-varkappad)*pibreve^varkappad;

mc=1/lambdad*(1-betaa*xid*(pitilded/pibar)^(lambdad/(1-lambdad)))...
    /(1-betaa*xid*(pitilded/pibar)^(1/(1-lambdad)))...
    *((1-xid*(pitilded/pibar)^(1/(1-lambdad)))/(1-xid))^(1-lambdad);

phalo=((1-xid*(pitilded/pibar)^(1/(1-lambdad)))/(1-xid))^(1-lambdad)...
    /((1-xid*(pitilded/pibar)^(lambdad/(1-lambdad)))/(1-xid))^((1-lambdad)/lambdad);


%Nonlinear pricing exporters
pitildex=pix^kappax*pix^(1-kappax-varkappax)*pibreve^varkappax;

mcx=1/lambdax*((1-xix*(pitildex/pix)^(1/(1-lambdax)))/(1-xix))^((1-lambdax))...
    *(1-betaa*xix*(pitildex/pix)^(lambdax/(1-lambdax)))...
    /(1-betaa*xix*(pitildex/pix)^(1/(1-lambdax)));


phalox=((1-xix*(pitildex/pix)^(1/(1-lambdax)))/(1-xix))^(1-lambdax)...
    /((1-xix*(pitildex/pix)^(lambdax/(1-lambdax)))/(1-xix))^((1-lambdax)/lambdax);


%Nonlinear pricing consumption importer
pitildemc=pimc^kappamc*pitarget^(1-kappamc-varkappamc)*pibreve^varkappamc;

mcmc=1/lambdamc*((1-ximc*(pitildemc/pimc)^(1/(1-lambdamc)))/(1-ximc))^((1-lambdamc))...
    *(1-betaa*ximc*(pitildemc/pimc)^(lambdamc/(1-lambdamc)))...
    /(1-betaa*ximc*(pitildemc/pimc)^(1/(1-lambdamc)));

phalomc=((1-ximc*(pitildemc/pimc)^(1/(1-lambdamc)))/(1-ximc))^(1-lambdamc)...
    /((1-ximc*(pitildemc/pimc)^(lambdamc/(1-lambdamc)))/(1-ximc))^((1-lambdamc)/lambdamc);

%Nonlinear pricing investment importer
pitildemi=pimi^kappami*pitarget^(1-kappami-varkappami)*pibreve^varkappami;

mcmi=1/lambdami*((1-ximi*(pitildemi/pimi)^(1/(1-lambdami)))/(1-ximi))^((1-lambdami))...
    *(1-betaa*ximi*(pitildemi/pimi)^(lambdami/(1-lambdami)))...
    /(1-betaa*ximi*(pitildemi/pimi)^(1/(1-lambdami)));

phalomi=((1-ximi*(pitildemi/pimi)^(1/(1-lambdami)))/(1-ximi))^(1-lambdami)...
    /((1-ximi*(pitildemi/pimi)^(lambdami/(1-lambdami)))/(1-ximi))^((1-lambdami)/lambdami);

%Nonlinear pricing export importer
pitildemx=pimx^kappamx*pitarget^(1-kappamx-varkappamx)*pibreve^varkappamx;

mcmx=1/lambdamx*((1-ximx*(pitildemx/pimx)^(1/(1-lambdamx)))/(1-ximx))^((1-lambdamx))...
    *(1-betaa*ximx*(pitildemx/pimx)^(lambdamx/(1-lambdamx)))...
    /(1-betaa*ximx*(pitildemx/pimx)^(1/(1-lambdamx)));

phalomx=((1-ximx*(pitildemx/pimx)^(1/(1-lambdamx)))/(1-ximx))^(1-lambdamx)...
    /((1-ximx*(pitildemx/pimx)^(lambdamx/(1-lambdamx)))/(1-ximx))^((1-lambdamx)/lambdamx);

%Wage setting distortions and indexation
pitildew=pic^kappaw*pitarget^(1-kappaw-varkappaw)*pibreve^varkappaw*muzplus^thetaw;

%unemployment rate
unemp=1-L;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Solving for the core steady state equations using fsolve:%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% xxx=[0.1045
%     0.4930
%     0.1734
%     0.2836
%     0.0676
%     0.1065
%     0.2561
%     0.1849
%     0.9687
%     0.2864
%     0.1218];

% at ctw_employer_surplus_mode_0323a, an endo prior mode
% xxx=[1.2034
%     0.4881
%     0.4028
%     0.0416
%     0.0704
%     0.1783
%     0.1514
%     0.3278
%     0.6988
%     0.2902
%     0.2114];
xxx=[ 1 % sigma_a
    1   % omegabar
    1   % rkbar
    1   % AL
    1   % v0
    1   % psizplus
    1   % Uzplus
%    1   % delta
    1   % a0
    1   % sigma
    1]; % varphi 

   
%max # of fsolve tries
max_fsolve_tries=10;
fsolve_tries=1;
check=1; %Dynare steady state indicator (0 good, 1 bad)

%use Fsolve to solve the core
while fsolve_tries<=max_fsolve_tries && check==1,
    
    if fsolve_tries == 1
        if ~isempty(LastSuccessFsolve)
            xxx=LastSuccessFsolve;
        end
        
    % if equal 2: use mode guess from fixed vector above
        
    elseif fsolve_tries==3    % first mode, at calib values
        xxx=[ 1 % sigma_a
              1   % omegabar
              1   % rkbar
              1   % AL
              1   % v0
              1   % psizplus
              1   % Uzplus
             % 1   % delta
              1   % a0
              1   % sigma
              1]; % varphi
    elseif fsolve_tries>3
        %Generate random guesses from the interval [a, b].
        ten=8; lowerbound=0.005;
        % yields implicit upperbound of ten/2+lowerbound
        %xxx=ones(size(xxx))*(ten*rand*rand/2+lowerbound);
        xxx(1)=1*(ten*rand*rand/2+lowerbound); % sigma_a
        xxx(2)=1*(ten*rand*rand/2+lowerbound); % omegabar 
        xxx(3)=1*(ten*rand*rand/2+lowerbound); % rkbar
        xxx(4)=1*(ten*rand*rand/2+lowerbound); % AL 
        xxx(5)=1*(ten*rand*rand/2+lowerbound); % v0 
        xxx(6)=1*(ten*rand*rand/2+lowerbound); % psizplus 
        xxx(7)=1*(ten*rand*rand/2+lowerbound); % Uzplus 
        %xxx(8)=1*(ten*rand*rand/2+lowerbound); % delta
        xxx(8)=1*(ten*rand*rand/2+lowerbound); % a0 
        xxx(9)=1*(ten*rand*rand/2+lowerbound);% sigma                
        xxx(10)=1*(ten*rand*rand/2+lowerbound);% varphi         
    end

    try

        %solve core model equs with fsolve
        opts=optimset('Display','off','MaxIter',10000,'MaxFunEvals',200000,'TolFun',1e-10,'TolX',1e-10);
        [yyy,fval,exitflag]=fsolve(@ctw_steadystate_solvecore_10,xxx,opts,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,BigF,xy_target,SW,SE,delta);
        
        fsolve_tries=fsolve_tries+1;
        
        if   exitflag<=0
            %disp('Cannot find a steady state');
            ys=ones(M_.endo_nbr)*1e+2;
            check=1;
        elseif  exitflag>=1
%             [ff,a0,a1,a2,a3,H,wbar,kbar,pmx,pmc,pmi,pc,pinvest,px,pkprime,Rk,z,Gamma,Gam_muG,ybar,we, ...
%                 epscall0,epscall1,epscall2,epscall3,fcall0,fcall1,fcall2,fcall3,intensity0,intensity1,intensity2,intensity3, ...
%                 G0,G1,G2,G3,mcall0,mcall1,mcall2,mcall3,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,D1zplus,D2zplus,D3zplus,sigma,sigma_a,omegabar,rkbar,AL,psizplus,Uzplus,delta,varphi] = ...
%                 ctw_steadystate_solvecore(yyy,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,BigF,xy_target);
%             
            [ff,a0,a1,a2,a3,H,wbar,kbar,pmx,pmc,pmi,pc,pinvest,px,pkprime,Rk,z,Gamma,Gam_muG,ybar,we, ...
                epscall0,epscall1,epscall2,epscall3,fcall0,fcall1,fcall2,fcall3,intensity0,intensity1,intensity2,intensity3, ...
                G0,G1,G2,G3,mcall0,mcall1,mcall2,mcall3,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(yyy,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,BigF,xy_target,SW,SE,delta);
            
            LastSuccessFsolve=yyy;
                        
            %further steady states
            Jtildezplus1=Jtildezplus0;
            Jtildezplus2=Jtildezplus0;
            Jtildezplus3=Jtildezplus0;
            
            % nk
            nk=(pkprime*(1-Rk/R*Gam_muG));
            
            %use law of motion for net worth to calc gammaa
            gammaa=(1-we/nk/kbar)/((Rk-R-(Gamma-Gam_muG)*Rk)*pkprime/nk+R)*(pibar*muzplus);
            n=nk*kbar;
            spread=omegabar*Rk/(1-n/(pkprime*kbar))-R;
            
            %further steady states
            g=etag*ybar;
            q=varphi/pc;
            i=kbar*(1-(1-delta)/(muzplus*mupsi))/Upsilon;
            x=px^(-etaf)*ystar;
            c=zetac*(muzplus-betaa*b)/(pc*(1+tauc)*(muzplus-b))...
                /psizplus;
            phi=ybar*(1/mc-phalo^(lambdad/(1-lambdad)));
            a=etaa*ybar;
            k=kbar*u;
            im=omegai*(pinvest/pmi)^(etai)*i;
            xm=omegax*(1/pmx*(omegax*pmx^(1-etax)+1-omegax)^(1/(1-etax)))^etax...
                *phalox^(-lambdax/(lambdax-1))*px^(-etaf)*ystar;
            cm=omegac*(pc/pmc)^(etac)*c;
            
            %nonlinear pricing variables
            Kd=lambdad*psizplus*ybar*mc/(1-betaa*xid*(pitilded/pibar)^(lambdad/(1-lambdad)));
            Fd=psizplus*ybar/(1-betaa*xid*(pitilded/pibar)^(1/(1-lambdad)));
            
            Kx=lambdax*psizplus*varphi*px*x*mcx/(1-betaa*xix*(pitildex/pix)^(lambdax/(1-lambdax)));
            Fx=psizplus*varphi*px*x/(1-betaa*xix*(pitildex/pix)^(1/(1-lambdax)));
            
            Kmc=lambdamc*psizplus*pmc*mcmc*cm/(1-betaa*ximc*(pitildemc/pimc)^(lambdamc/(1-lambdamc)));
            Fmc=psizplus*pmc*cm/(1-betaa*ximc*(pitildemc/pimc)^(1/(1-lambdamc)));
            
            Kmi=lambdami*psizplus*pmi*mcmi*im/(1-betaa*ximi*(pitildemi/pimi)^(lambdami/(1-lambdami)));
            Fmi=psizplus*pmi*im/(1-betaa*ximi*(pitildemi/pimi)^(1/(1-lambdami)));
            
            Kmx=lambdamx*psizplus*pmx*mcmx*xm/(1-betaa*ximx*(pitildemx/pimx)^(lambdamx/(1-lambdamx)));
            Fmx=psizplus*pmx*xm/(1-betaa*ximx*(pitildemx/pimx)^(1/(1-lambdamx)));
            
            %capacity utilization parameter
            sigmab=rkbar/pinvest;
            aofu=0.5*sigmab*sigmaa+sigmab*(1-sigmaa)+sigmab*(sigmaa/2-1);
            
            %Creating Pitilde/Pi ratio variables
            pitildedpi=pitilded/pibar;
            pitildexpi=pitildex/pix;
            pitildemcpi=pitildemc/pimc;
            pitildemipi=pitildemi/pimi;
            pitildemxpi=pitildemx/pimx;
            
            %Data series steady states
            data_pid=0; %400*log(pibar);
            data_R=400*(R-1);
            data_wdiff=0 ;
            data_cdiff=0;
            data_idiff=0;
            data_qdiff=0;
            data_Hdiff=0;
            data_ydiff=0;
            data_xdiff=0;
            data_impdiff=0;
            data_pic=0; %400*log(pic);
            data_pii=0; %400*log(pii);
            data_ystardiff=0;
            data_pistar=400*log(pistar);
            data_Rstar=400*(Rstar-1);
            data_ndiff=0;
            data_unempdiff=0;
            %data_spreaddiff=0;
            data_spread=400*spread;
            data_gdiff=0;
            
            data_Hhat=0;
            data_unemphat=0;
            
            nxdivGDP = 0;
            
            avg_Nash=w*wbar*(l0*(1-fcall0)+ l1*G1*(1-fcall1) + l2*G2*(1-fcall2) + l3*G3*(1-fcall3))/L;
            ylessd  =ybar - mu*normcdf(z-sigma,0,1)*Rk*pkprime*kbar/(pibar*muzplus) ...
                -((kappa/2)*(vtilde0^2*l0*(1-fcall0)+vtilde1^2*l1*(1-fcall1)+vtilde2^2*l2*(1-fcall2)+vtilde3^2*l3*(1-fcall3)));
            
            Slev=.7028;
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %constructing the vector ys for Dynare%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            ys =[ log(mc); log(rkbar); log(wbar);  log(Rf); log(px); log(pix); log(k);
                log(Rnustar); log(H); log(pibar); log(mcx); log(pimc); log(mcmc);
                log(pimx); log(mcmx); log(pimi); log(mcmi); log(pic); log(pmc); log(pii); log(pmi);
                log(kbar);  log(i);  log(psizplus);  log(c);log(pkprime);  log(u);  log(pinvest); %log(s); 
                a;  log(q); log(x); log(pc); log(R);
                log(ybar); log(Rx); log(pmx); log(Rk); Stilde; Sprimetilde;
                log(pitildedpi); log(Kd); log(Fd); log(phalo);
                log(pitildexpi); log(Kx); log(Fx); log(phalox);
                log(pitildemcpi); log(Kmc); log(Fmc); log(phalomc); log(cm);
                log(pitildemipi); log(Kmi); log(Fmi); log(phalomi); log(im);
                log(pitildemxpi); log(Kmx); log(Fmx); log(phalomx); log(xm);
                aofu; log(piw);  log(Jzplus);log(w); log(Uzplus); log(Vxzplus);
                log(v);  log(G1); log(G2); log(G3); log(Om1);
                log(Om2); log(Om3); log(Om4);log(vtilde0); log(vtilde1); log(vtilde2); log(vtilde3);
                log(Vzplus0);log(Vzplus1); log(Vzplus2); log(Vzplus3); log(l0); log(l1);
                log(l2); log(l3); log(L); log(m); log(BigQ); log(pitildew); log(chi0);
                log(chi1); log(chi2); log(chi3); log(f); log(Vw); Jw; log(unemp);
                log(omegabar); log(n); log(spread); log(Fomegabar);
                log(muzplus); log(epsilon);  log(mupsi);  log(pitarget);
                log(Upsilon);  log(zetac);  log(Rstar); phitilde;  log(tauy); log(zetah);
                log(ystar); log(pistar);  log(g);  log(muz);
                log(taud); log(taux); log(taumc); log(taumi); log(taumx);
                log(sigmam);log(eta); log(sigma_a); log(sigma); log(gammaa);
                data_pid; data_R; data_wdiff; data_cdiff; data_idiff;
                data_qdiff; data_Hdiff; data_ydiff; data_xdiff; data_impdiff;
                data_pic; data_pii; data_ystardiff; data_pistar; data_Rstar;
                data_ndiff;data_unempdiff; data_spread; data_gdiff;
                log(avg_Nash); log(ylessd); log(epscall0); log(epscall1);
                log(epscall2); log(epscall3); log(fcall0); log(fcall1) ;log(fcall2) ;log(fcall3);
                log(intensity0); log(intensity1); log(intensity2); log(intensity3);
                log(mcall0); log(mcall1); log(mcall2); log(mcall3);
                log(a0) ;log(a1); log(a2) ;log(a3);
                % removed old: log(D1zplus); log(D2zplus); log(D3zplus);
                log(gcallprime0); log(gcallprime1); log(gcallprime2); log(gcallprime3);
                log(fcallprime0); log(fcallprime1); log(fcallprime2); log(fcallprime3);
                log(Jtildezplus0); log(Jtildezplus1); log(Jtildezplus2); log(Jtildezplus3);
                data_Hhat; data_unemphat; nxdivGDP; log(Slev); log(pkprime*kbar-n); a/ylessd];
            
            %Dynare Switch indicating good steady state
            check=0;
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %Checking the steady states%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            if check==0,
                vv=[rkbar,wbar,Rnustar,Rf,Rx,R,mc,mcx,mcmc,mcmi,mcmx,pibar,pix,pic,pii, ...
                    pimc,pimi,pc,px,pinvest,pmx,pmc,pmi,pkprime,kbar,sigmab,H,q,i,c,x, ...
                    psizplus,ybar,ystar,cm,xm,im,n,omegabar,spread];
                
                if min(vv) < 0
                    %disp('fatal (steadystate) steady state not admissible, negative values where they are not allowed');
                    ys=ys*1e+2; %passing a stupid steady state to dynare
                    check=1; %indicating this is a stupid steady state; Results in a low likelihood.
                end
                
                if BigQ>1 | BigQ<0 | Jzplus<0 | eta>1 | eta<0 |  w<0 |  Jw > 0 | Vw < 0 | H < 0 ...
                        | f<0 | f>1 | Vzplus0<0 | Vzplus1<0 | Vzplus2<0 | Vzplus3<0 ...
                        | v0<0 | v1<0| v2<0 | v3<0 | Uzplus<0,
                    %disp([BigQ;Jzplus;eta;w;Jw;Vw;H;intensity;f;Vzplus0;Vzplus1;Vzplus2;Vzplus3;Uzplus]);
                    disp('fatal: sign restrictions of labor market steady state violated')
                    ys=ys*1e+2;
                    check=1;
                    %keyboard
                end
            end
        end
        
    catch
        ys=ones(size(M_.endo_names,1))*1e+2;
        check=1;
        disp('Panic: Fsolve has a problem to terminate or other steady state calcs fail.');
        fsolve_tries=fsolve_tries+1;
    end
    
end

if fsolve_tries-1==max_fsolve_tries && check==1,
    %disp(['I tried ',num2str(max_fsolve_tries),' random guesses for fsolve.']);
    disp('I give up now and return no steady state.');
end


if fsolve_tries>2 && check==0,
    disp('Yeah, I found a steady state based on random initial guess.');
end

% saving s.s. values that are used as params in static file

M_.params(	1	)	=	alphaa              	;
M_.params(	2	)	=	kappad             	;
M_.params(	3	)	=	betaa               	;
M_.params(	4	)	=	rhopi              	;
M_.params(	5	)	=	xid                	;
M_.params(	6	)	=	xix                	;
M_.params(	7	)	=	kappax             	;
M_.params(	8	)	=	ximc               	;
M_.params(	9	)	=	etaf               	;
M_.params(	10	)	=	ximx               	;
M_.params(	11	)	=	kappamc            	;
M_.params(	12	)	=	kappamx            	;
M_.params(	13	)	=	kappami            	;
M_.params(	14	)	=	ximi               	;
M_.params(	15	)	=	omegac             	;
M_.params(	16	)	=	etac               	;
M_.params(	17	)	=	omegai             	;
M_.params(	18	)	=	etai               	;
M_.params(	19	)	=	Spp                	;
M_.params(	20	)	=	muzplus            	;
M_.params(	21	)	=	mupsi              	;
M_.params(	22	)	=	delta              	;
M_.params(	23	)	=	b                  	;
M_.params(	24	)	=	sigmab             	;
M_.params(	25	)	=	sigmaa             	;
M_.params(	26	)	=	phitildea          	;
M_.params(	27	)	=	a                  	;
M_.params(	28	)	=	phitildes          	;
M_.params(	29	)	=	rhoR               	;
M_.params(	30	)	=	rpi                	;
M_.params(	31	)	=	ry                 	;
M_.params(	32	)	=	rq                 	;
M_.params(	33	)	=	rdeltapi           	;
M_.params(	34	)	=	rdeltay            	;
M_.params(	35	)	=	etax               	;
M_.params(	36	)	=	omegax             	;
M_.params(	37	)	=	phi                	;
M_.params(	38	)	=	rhomupsi           	;
M_.params(	39	)	=	epsilon            	;
M_.params(	40	)	=	pibar              	;
M_.params(	41	)	=	Upsilon            	;
M_.params(	42	)	=	zetac              	;
M_.params(	43	)	=	tauc               	;
M_.params(	44	)	=	tauk               	;
M_.params(	45	)	=	R                  	;
M_.params(	46	)	=	tauy               	;
M_.params(	47	)	=	tauw               	;
M_.params(	48	)	=	zetah              	;
M_.params(	49	)	=	nustar             	;
M_.params(	50	)	=	nuf                	;
M_.params(	51	)	=	nux                	;
M_.params(	52	)	=	etag               	;
M_.params(	53	)	=	lambdad            	;
M_.params(	54	)	=	lambdax            	;
M_.params(	55	)	=	lambdamc           	;
M_.params(	56	)	=	lambdami           	;
M_.params(	57	)	=	lambdamx           	;
M_.params(	58	)	=	wbar               	;
M_.params(	59	)	=	pic                	;
M_.params(	60	)	=	psizplus           	;
M_.params(	61	)	=	H                  	;
M_.params(	62	)	=	q                  	;
M_.params(	63	)	=	pimi               	;
M_.params(	64	)	=	mcmi               	;
M_.params(	65	)	=	mc                 	;
M_.params(	66	)	=	pix                	;
M_.params(	67	)	=	mcmx               	;
M_.params(	68	)	=	mcx                	;
M_.params(	69	)	=	pimx               	;
M_.params(	70	)	=	muz                	;
M_.params(	71	)	=	ybar               	;
M_.params(	72	)	=	u                  	;
M_.params(	73	)	=	kappaw             	;
M_.params(	74	)	=	AL                 	;
M_.params(	75	)	=	sigmaL             	;
M_.params(	76	)	=	Rk                 	;
M_.params(	77	)	=	varphi             	;
M_.params(	78	)	=	ystar              	;
M_.params(	79	)	=	mcmc               	;
M_.params(	80	)	=	px                 	;
M_.params(	81	)	=	pimc               	;
M_.params(	82	)	=	pmc                	;
M_.params(	83	)	=	pmi                	;
M_.params(	84	)	=	kbar               	;
M_.params(	85	)	=	i                  	;
M_.params(	86	)	=	c                  	;
M_.params(	87	)	=	pc                 	;
M_.params(	88	)	=	pkprime            	;
M_.params(	89	)	=	pinvest            	;
M_.params(	90	)	=	pmx                	;
M_.params(	91	)	=	pii                	;
M_.params(	92	)	=	pistar             	;
M_.params(	93	)	=	Rstar              	;
M_.params(	94	)	=	rhomuz             	;
M_.params(	95	)	=	rhoepsilon         	;
M_.params(	96	)	=	rhoUpsilon         	;
M_.params(	97	)	=	rhozetac           	;
M_.params(	98	)	=	rhozetah           	;
M_.params(	99	)	=	rhog               	;
M_.params(	100	)	=	rhotauy            	;
M_.params(	101	)	=	rhophitilde        	;
M_.params(	102	)	=	phitilde           	;
M_.params(	103	)	=	sig_muz            	;
M_.params(	104	)	=	sig_epsilon        	;
M_.params(	105	)	=	sig_upsilon        	;
M_.params(	106	)	=	sig_zetac          	;
M_.params(	107	)	=	sig_zetah          	;
M_.params(	108	)	=	sig_phitilde       	;
M_.params(	109	)	=	sig_epsR           	;
M_.params(	110	)	=	sig_pitarget       	;
M_.params(	111	)	=	sig_g              	;
M_.params(	112	)	=	sig_tauy           	;
M_.params(	113	)	=	s                  	;
M_.params(	114	)	=	Stilde             	;
M_.params(	115	)	=	Sprimetilde        	;
M_.params(	116	)	=	rkbar              	;
M_.params(	117	)	=	Rf                 	;
M_.params(	118	)	=	Rnustar            	;
M_.params(	119	)	=	x                  	;
M_.params(	120	)	=	Rx                 	;
M_.params(	121	)	=	varkappad          	;
M_.params(	122	)	=	varkappax          	;
M_.params(	123	)	=	varkappamx         	;
M_.params(	124	)	=	varkappami         	;
M_.params(	125	)	=	varkappamc         	;
M_.params(	126	)	=	varkappaw          	;
M_.params(	127	)	=	thetaw             	;
M_.params(	128	)	=	pibreve            	;
M_.params(	129	)	=	pitildedpi         	;
M_.params(	130	)	=	Kd                 	;
M_.params(	131	)	=	Fd                 	;
M_.params(	132	)	=	phalo              	;
M_.params(	133	)	=	pitildexpi         	;
M_.params(	134	)	=	Kx                 	;
M_.params(	135	)	=	Fx                 	;
M_.params(	136	)	=	phalox             	;
M_.params(	137	)	=	pitildemcpi        	;
M_.params(	138	)	=	Kmc                	;
M_.params(	139	)	=	Fmc                	;
M_.params(	140	)	=	phalomc            	;
M_.params(	141	)	=	cm                 	;
M_.params(	142	)	=	pitildemipi        	;
M_.params(	143	)	=	Kmi                	;
M_.params(	144	)	=	Fmi                	;
M_.params(	145	)	=	phalomi            	;
M_.params(	146	)	=	im                 	;
M_.params(	147	)	=	pitildemxpi        	;
M_.params(	148	)	=	Kmx                	;
M_.params(	149	)	=	Fmx                	;
M_.params(	150	)	=	phalomx            	;
M_.params(	151	)	=	xm                 	;
M_.params(	152	)	=	aofu               	;
M_.params(	153	)	=	taud               	;
M_.params(	154	)	=	taux               	;
M_.params(	155	)	=	taumc              	;
M_.params(	156	)	=	taumi              	;
M_.params(	157	)	=	taumx              	;
M_.params(	158	)	=	rhotaud            	;
M_.params(	159	)	=	rhotaux            	;
M_.params(	160	)	=	rhotaumc           	;
M_.params(	161	)	=	rhotaumi           	;
M_.params(	162	)	=	rhotaumx           	;
M_.params(	163	)	=	piw                	;
M_.params(	164	)	=	taub               	;
M_.params(	165	)	=	pitarget           	;
M_.params(	166	)	=	etaa               	;
M_.params(	167	)	=	a11                	;
M_.params(	168	)	=	a12                	;
M_.params(	169	)	=	a13                	;
M_.params(	170	)	=	a21                	;
M_.params(	171	)	=	a22                	;
M_.params(	172	)	=	a23                	;
M_.params(	173	)	=	a24                	;
M_.params(	174	)	=	a31                	;
M_.params(	175	)	=	a32                	;
M_.params(	176	)	=	a33                	;
M_.params(	177	)	=	a34                	;
M_.params(	178	)	=	sig_ystar          	;
M_.params(	179	)	=	c21                	;
M_.params(	180	)	=	sig_pistar         	;
M_.params(	181	)	=	c24                	;
M_.params(	182	)	=	c31                	;
M_.params(	183	)	=	c32                	;
M_.params(	184	)	=	sig_Rstar          	;
M_.params(	185	)	=	c34                	;
M_.params(	186	)	=	we                 	;
M_.params(	187	)	=	mu                 	;
M_.params(	188	)	=	Fomegabar          	;
M_.params(	189	)	=	n                  	;
M_.params(	190	)	=	omegabar           	;
M_.params(	191	)	=	sigma              	;
M_.params(	192	)	=	gammaa              	;
M_.params(	193	)	=	spread             	;
M_.params(	194	)	=	rhosigma           	;
M_.params(	195	)	=	rhogamma           	;
M_.params(	196	)	=	recruitshare       	;
M_.params(	197	)	=	bshare             	;
M_.params(	198	)	=	rhosigmam          	;
M_.params(	199	)	=	rhoeta             	;
M_.params(	200	)	=	BigN               	;
M_.params(	201	)	=	rho                	;
M_.params(	202	)	=	iota               	;
M_.params(	203	)	=	sigma_match        	;
M_.params(	204	)	=	kappa              	;
M_.params(	205	)	=	eta                	;
M_.params(	206	)	=	bu                 	;
M_.params(	207	)	=	Jw                 	;
M_.params(	208	)	=	Vw                 	;
M_.params(	209	)	=	BigQ               	;
M_.params(	210	)	=	v                  	;
M_.params(	211	)	=	m                  	;
M_.params(	212	)	=	sigmam             	;
M_.params(	213	)	=	f                  	;
M_.params(	214	)	=	w                  	;
M_.params(	215	)	=	Jzplus             	;
M_.params(	216	)	=	Vxzplus            	;
M_.params(	217	)	=	Uzplus             	;
M_.params(	218	)	=	recruiting         	;
M_.params(	219	)	=	L                  	;
M_.params(	220	)	=	pitildew           	;
M_.params(	221	)	=	G1                 	;
M_.params(	222	)	=	G2                 	;
M_.params(	223	)	=	G3                 	;
M_.params(	224	)	=	chi0               	;
M_.params(	225	)	=	chi1               	;
M_.params(	226	)	=	chi2               	;
M_.params(	227	)	=	chi3               	;
M_.params(	228	)	=	Om1                	;
M_.params(	229	)	=	Om2                	;
M_.params(	230	)	=	Om3                	;
M_.params(	231	)	=	Om4                	;
M_.params(	232	)	=	vtilde0            	;
M_.params(	233	)	=	vtilde1            	;
M_.params(	234	)	=	vtilde2            	;
M_.params(	235	)	=	vtilde3            	;
M_.params(	236	)	=	Vzplus0            	;
M_.params(	237	)	=	Vzplus1            	;
M_.params(	238	)	=	Vzplus2            	;
M_.params(	239	)	=	Vzplus3            	;
M_.params(	240	)	=	l0                 	;
M_.params(	241	)	=	l1                 	;
M_.params(	242	)	=	l2                 	;
M_.params(	243	)	=	l3                 	;
M_.params(	244	)	=	recruitsharePercent	;
M_.params(	245	)	=	SppScaled          	;
M_.params(	246	)	=	muzPercent         	;
M_.params(	247	)	=	mupsiPercent       	;
M_.params(	248	)	=	work_cap_para       ;
M_.params(	249	)	=	wesharePercent     	;
M_.params(	250	)	=	hours_target       	;
M_.params(	251	)	=	iy_target          	;
M_.params(	252	)	=	xy_target          	;
M_.params(	253	)	=	nk_target          	;
M_.params(	254	)	=	avg_Nash           	;
M_.params(	255	)	=	ylessd             	;
M_.params(	256	)	=	sigma_a            	;
M_.params(	257	)	=	epscall0           	;
M_.params(	258	)	=	epscall1           	;
M_.params(	259	)	=	epscall2           	;
M_.params(	260	)	=	epscall3           	;
M_.params(	261	)	=	fcall0             	;
M_.params(	262	)	=	fcall1             	;
M_.params(	263	)	=	fcall2             	;
M_.params(	264	)	=	fcall3             	;
M_.params(	265	)	=	intensity0         	;
M_.params(	266	)	=	intensity1         	;
M_.params(	267	)	=	intensity2         	;
M_.params(	268	)	=	intensity3         	;
M_.params(	269	)	=	mcall0             	;
M_.params(	270	)	=	mcall1             	;
M_.params(	271	)	=	mcall2             	;
M_.params(	272	)	=	mcall3             	;
M_.params(	273	)	=	a0                 	;
M_.params(	274	)	=	a1                 	;
M_.params(	275	)	=	a2                 	;
M_.params(	276	)	=	a3                 	;
M_.params(	277	)	=	SW                 	;
M_.params(	278	)	=	SE                 	;
M_.params(	279	)	=	gcallprime0        	;
M_.params(	280	)	=	gcallprime1        	;
M_.params(	281	)	=	gcallprime2        	;
M_.params(	282	)	=	gcallprime3        	;
M_.params(	283	)	=	fcallprime0        	;
M_.params(	284	)	=	fcallprime1        	;
M_.params(	285	)	=	fcallprime2        	;
M_.params(	286	)	=	fcallprime3        	;
M_.params(	287	)	=	Jtildezplus0       	;
M_.params(	288	)	=	Jtildezplus1       	;
M_.params(	289	)	=	Jtildezplus2       	;
M_.params(	290	)	=	Jtildezplus3       	;
M_.params(	291	)	=	rhosigma_a         	;
M_.params(	292	)	=	sig_sigma_a        	;
M_.params(	293	)	=	BigF               	;
M_.params(	294	)	=	BigFPercent        	;
M_.params(	295	)	=	ab  	;
M_.params(	296	)	=	sigmaLScaled  	;
M_.params(	297	)	=	rRstar  	;



do_steady_print=0;
if do_steady_print ==1 && check==0,
    fprintf(' \n\nsteady state properties of open economy model \n\n')
    fprintf(' Rf = %6.4f, Rx = %6.4f, R = %6.4f, mc = %6.4f, mcx = %6.4f, mcmc = %6.4f, pibar = %6.4f\n',Rf,Rx,R,1/lambdad,mcx,mcmc,pibar)
    fprintf(' pix = %6.4f, pic = %6.4f, pii = %6.4f, pimc = %6.4f, pimi = %6.4f, pc = %6.4f, px = %6.4f\n',pix,pic,pii,pimc,pimi,pc,px)
    fprintf(' pinvest = %6.4f, pmx = %6.4f, pmc = %6.4f,pkprime = %6.4f, kbar = %6.4f, sigmab = %6.4f\n',pinvest,pmx,pmc,pkprime,kbar,sigmab)
    fprintf(' q*pc = %6.4f, H = %6.4f, q = %6.4f, i = %6.4f, c = %6.4f, x = %6.4f \n',q*pc,H,q,i,c,x)
    fprintf(' a = %6.4f, psizplus = %6.4f, ybar = %6.4f, ystar = %6.4f, rkbar = %6.4f, Rk = %6.4f\n',a,psizplus,ybar,ystar,rkbar,Rk)
    fprintf(' ky = %6.4f, cy = %6.4f, iy = %6.4f, xy = %6.4f, yystar = %6.4f, real R = %6.6f\n\n',kbar/ybar,c/ybar,i/ybar,x/ybar,ybar/ystar,R/pibar)
    
    fprintf(' ****************************\n');
    fprintf(' k*Pk/y = %6.4f, c*Pc/y = %6.4f, i*Pinvest/y = %6.4f, x*Px/y = %6.4f, n/(pkprime*k) = %6.4f \n\n',kbar*pkprime/ybar,c*pc/ybar,i*pinvest/ybar,x*px/ybar,n/(pkprime*kbar))
    fprintf(' x*capPx/(y*P) = %6.4f',px*q*pc*x/(ybar*s));
    
    fprintf(' \n model parameter values\n\n')
    fprintf(' mupsi = %6.4f, alphaa = %6.4f, pibar = %6.4f, nustar = %6.4f, nuf = %6.4f, nux = %6.4f\n',mupsi,alphaa,pibar,nustar,nuf,nux)
    fprintf(' lambdamx = %6.4f, lambdax = %6.4f, lambdad = %6.4f, omegax = %6.4f, etax = %6.4f, etag = %6.4f\n',lambdamx,lambdax,lambdad,omegax,etax,etag)
    fprintf(' AL = %6.4f, sigmaL = %6.4f, Upsilon = %6.4f, epsilon = %6.4f\n',AL,sigmaL,Upsilon,epsilon)
    fprintf(' zetac = %6.4f, zetah = %6.4f, sigmaa = %6.4f\n',zetac,zetah,sigmaa)
    fprintf(' kappaw = %6.4f, varkappaw = %6.4f, thetaw = %6.4f\n\n',kappaw,varkappaw,thetaw)
    
    fprintf('\n*** Financial frictions block of s.s. values below *** \n')
    fprintf('sigma = %6.4f, omegabar = %6.4f, n = %6.4f, spread = %6.4f \n\n', sigma,omegabar,n,spread );
    
    fprintf('*** Financial frictions block parameter values below *** \n')
    fprintf('Fomegabar = %6.4f, mu = %6.4f, we = %6.4f, gammaa = %6.4f\n\n', Fomegabar,mu,we,gammaa)
    
    
    fprintf('\n*** Labor block of s.s. values below *** \n\n')
    fprintf('m = %6.4f, f = %6.4f, unemp = %6.4f\n', m,f,unemp  )
    fprintf('H = %6.4f, Jw = %6.4f, Vw = %6.4f \n',H,Jw,Vw )
    fprintf('BigQ = %6.4f, v = %6.4f, eta = %6.4f, w = %6.4f\n',BigQ,v,eta,w)
    fprintf('Jzplus = %6.4f, Vxzplus = %6.4f, Uzplus = %6.4f\n',Jzplus,Vxzplus,Uzplus)
    fprintf('kappa = %6.4f, wbar = %6.4f, bu = %6.4f\n\n', kappa ,wbar, bu)
        
    fprintf('epscall0 = %6.4f, epscall1 = %6.4f, epscall2 = %6.4f, epscall3 = %6.4f\n', epscall0,epscall1,epscall2,epscall3 )
    fprintf('fcall0 = %6.4f, fcall1 = %6.4f, fcall2 = %6.4f, fcall3 = %6.4f\n', fcall0,fcall1,fcall2,fcall3 )
    fprintf('a0 = %6.4f, a1 = %6.4f, a2 = %6.4f, a3 = %6.4f\n', a0,a1,a2,a3 )
    fprintf('l0 = %6.4f, l1 = %6.4f, l2 = %6.4f, l3 = %6.4f\n', l0,l1,l2,l3 )
    fprintf('v0 = %6.4f, v1 = %6.4f, v2 = %6.4f, v3 = %6.4f\n\n', v0,v1,v2,v3 )
    
    fprintf('intensity1 = %6.4f, intensity1 = %6.4f, intensity2 = %6.4f, intensity3 = %6.4f\n\n',intensity0,intensity1, intensity2 ,intensity3 )
    
    fprintf('\n*** Labor block of parameter values below *** \n\n')
    fprintf('sigmam = %6.4f, unemp = %6.4f, BigN = %6.4f, rho = %6.4f\n',sigmam,unemp,BigN,rho)
    fprintf('iota = %6.4f, sigma_match = %6.4f, sigma_a = %6.4f\n\n',iota,sigma_match,sigma_a )
    
    disp('*** key parts of endog. sep FOC for cut-off below *** ');
    intensity0_prime=1/sigmaL*intensity0^(1-sigmaL)*w*wbar*psizplus*(1-tauy)/(1+tauw)*gcallprime0/AL;
    fprintf('intensity0_prime = %6.4f, gcallprime0 = %6.4f\n\n',intensity0_prime, gcallprime0 )
     fprintf('LHS SE-term: wbar*((epscall0/(1-fcall0)-w)*intensity0_prime + gcallprime0*intensity0)= %6.4f\n',wbar*((epscall0/(1-fcall0)-w)*intensity0_prime + gcallprime0*intensity0) );
     fprintf('LHS SW-term: (w*wbar*(1-tauy)/(1+tauw)-AL*intensity0^sigmaL/psizplus)*intensity0_prime= %6.4f\n', (w*wbar*(1-tauy)/(1+tauw)-AL*intensity0^sigmaL/psizplus)*intensity0_prime );   
     fprintf('Both these are multiplied by the factor (1-fcall0) / fcallprime0 = %6.4f\n', (1-fcall0) / fcallprime0 );  
     
     recruiting_costs_to_y=recruiting/ybar
    
    total_separation_rate=1-rho*(1-fcall0)
    job_lasting_in_quarters=1/(1-rho*(1-fcall0))
    job_finding_duration_in_quarters=1/f
    
    theoretical_replacement_ratio=bu*(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)
    
    
    varnames=[  ' ydiffU            '
        ' pidU              '
        ' RU                '
        ' wdiffU            '
        ' cdiffU            '
        ' idiffU            '
        ' qdiffU            '
        ' HdiffU            '
        ' xdiffU            '
        ' impdiffU          '
        ' picU              '
        ' piiU              '
        ' ystardiffU        '
        ' pistarU           '
        ' RstarU            '];
    
    means=[  100*(muzplus-1)
        400*log(pibar)
        400*(R-1)
        100*(muzplus-1)
        100*(muzplus-1)
        100*(muzplus*mupsi-1)
        0
        0
        100*(muzplus-1)
        100*(muzplus-1)
        400*log(pic)
        400*log(pii)
        100*(muzplus-1)
        400*log(pistar)
        400*(Rstar-1)];
    
    disp(' ');
    disp('Inspect steady state...');
    keyboard
end
