function [ys,check] = baseline70_steadystate(ys,junk);

%dbstop if warning
%dbstop if error
%dbstop at 423

global M_ check;
persistent LastSuccessFsolve;


%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	);
lambdaw     	=	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	);
xiw         	=	M_.params(	117	);
rkbar       	=	M_.params(	118	);
Rf          	=	M_.params(	119	);
Rnustar     	=	M_.params(	120	);
x           	=	M_.params(	121	);
Rx          	=	M_.params(	122	);
varkappad   	=	M_.params(	123	);
varkappax   	=	M_.params(	124	);
varkappamx  	=	M_.params(	125	);
varkappami  	=	M_.params(	126	);
varkappamc  	=	M_.params(	127	);
varkappaw   	=	M_.params(	128	);
thetaw      	=	M_.params(	129	);
pibreve     	=	M_.params(	130	);
pitildedpi  	=	M_.params(	131	);
Kd          	=	M_.params(	132	);
Fd          	=	M_.params(	133	);
phalo       	=	M_.params(	134	);
pitildexpi  	=	M_.params(	135	);
Kx          	=	M_.params(	136	);
Fx          	=	M_.params(	137	);
phalox      	=	M_.params(	138	);
pitildemcpi 	=	M_.params(	139	);
Kmc         	=	M_.params(	140	);
Fmc         	=	M_.params(	141	);
phalomc     	=	M_.params(	142	);
cm          	=	M_.params(	143	);
pitildemipi 	=	M_.params(	144	);
Kmi         	=	M_.params(	145	);
Fmi         	=	M_.params(	146	);
phalomi     	=	M_.params(	147	);
im          	=	M_.params(	148	);
pitildemxpi 	=	M_.params(	149	);
Kmx         	=	M_.params(	150	);
Fmx         	=	M_.params(	151	);
phalomx     	=	M_.params(	152	);
xm          	=	M_.params(	153	);
pitildewpi  	=	M_.params(	154	);
Kw          	=	M_.params(	155	);
Fw          	=	M_.params(	156	);
whalo       	=	M_.params(	157	);
smallh      	=	M_.params(	158	);
aofu        	=	M_.params(	159	);
taud        	=	M_.params(	160	);
taux        	=	M_.params(	161	);
taumc       	=	M_.params(	162	);
taumi       	=	M_.params(	163	);
taumx       	=	M_.params(	164	);
rhotaud     	=	M_.params(	165	);
rhotaux     	=	M_.params(	166	);
rhotaumc    	=	M_.params(	167	);
rhotaumi    	=	M_.params(	168	);
rhotaumx    	=	M_.params(	169	);
piw         	=	M_.params(	170	);
taub        	=	M_.params(	171	);
pitarget    	=	M_.params(	172	);
etaa        	=	M_.params(	173	);
a11         	=	M_.params(	174	);
a12         	=	M_.params(	175	);
a13         	=	M_.params(	176	);
a21         	=	M_.params(	177	);
a22         	=	M_.params(	178	);
a23         	=	M_.params(	179	);
a24         	=	M_.params(	180	);
a31         	=	M_.params(	181	);
a32         	=	M_.params(	182	);
a33         	=	M_.params(	183	);
a34         	=	M_.params(	184	);
sig_ystar   	=	M_.params(	185	);
c21         	=	M_.params(	186	);
sig_pistar  	=	M_.params(	187	);
c24         	=	M_.params(	188	);
c31         	=	M_.params(	189	);
c32         	=	M_.params(	190	);
sig_Rstar   	=	M_.params(	191	);
c34         	=	M_.params(	192	);
SppScaled   	=	M_.params(	193	);
muzPercent  	=	M_.params(	194	);
% mupsiPercent	=	M_.params(	195	);
intensity_target=   M_.params(195);
iy_target       =   M_.params(196);
xy_target       =   M_.params(197);
hours_target    =   M_.params(198);
ylessd          =   M_.params(199);
sigmaLScaled    =   M_.params(200);
work_cap_para   =   M_.params(201);

% 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;
% get mupsi residually, so that muzplus is always the same (and must be given in .mod-file)
mupsi=(muzplus/muz)^((1-alphaa)/alphaa);

sigmaL=sigmaLScaled*10;
Spp=SppScaled*10;   % this obviously makes passing Spp redundant

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Steady State Calculations%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Balanced growth factor

%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;
Rk=pibar*muzplus/betaa; %not valid in financial frictions model

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


whalo=((1-xiw*(pitildew/piw)^(1/(1-lambdaw)))/(1-xiw))^(1-lambdaw)...
    /((1-xiw*(pitildew/piw)^(lambdaw/(1-lambdaw)))/(1-xiw))^((1-lambdaw)/lambdaw);

%find parameters AL, delta and varphi to match targets for intensity
%(smallh), pinvest*i/y and px*x/y

%it might be possible to get these parameters in closed form. However, this
%would require rewriting the most of the steady state program....thats why
%we go for fsolve instead


%set options
opts=optimset('Display','off','MaxIter',2000,'TolFun',1e-6,'TolX',1e-6);

%initial guess
xxx=[4      %AL
    0.025   %delta
    0.1];     %varphi

max_fsolve_tries=100;
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 == 2,
        if ~isempty(LastSuccessFsolve)
            xxx=LastSuccessFsolve;
        end

    elseif fsolve_tries>2
        %disp('Trouble finding a valid steady state.');
        %disp(['I am trying a random guess for fsolve: ',num2str(fsolve_tries),' out of ',num2str(max_fsolve_tries)]);
        %disp('Panic: LastSuccessFsolve did not work, I try random guesses now. Keep your fingers crossed');
        %Generate random guesses from the interval [a, b].
        ten=10;
        lowerbound=0.01;
        % yields implicit upperbound of 10/2+lowerbound, roughly 5 in
        xxx(1)=4*(ten*rand*rand/2+lowerbound);%AL
        xxx(2)=0.025*(ten*rand*rand/2+lowerbound);%delta
        xxx(3)=0.2*(ten*rand*rand/2+lowerbound);%varphi
    end

    try

         [yyy,fval,exitflag]=fsolve(@baseline_steadystate_find_AL_delta_varphi,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,piw,phalomi,phalomc,...
            phalomx,Rstar,lambdami,lambdamx,phalox,lambdax,taumx,taumc,taumi,mcmi,mcmx,mcmc,taux,mcx,...
            intensity_target,iy_target,xy_target,Rk,xiw,lambdaw);

        fsolve_tries=fsolve_tries+1;
 
        if  exitflag~=1,
            %disp('Cannot find a steady state');
            ys=ones(size(ys))*1e+2;
            check=1;

        elseif exitflag==1
            AL=yyy(1)*10;
            delta=yyy(2);
            varphi=yyy(3);
            LastSuccessFsolve=yyy;
        
            
            %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
            rkbar=(Rk*mupsi*pkprime/pibar-(1-delta)*pkprime-tauk*delta*mupsi*pkprime/pibar)/(1-tauk);

            %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
            kh=muzplus*((taud*mupsi^alphaa*wbar*Rf)/((1-alphaa)*mc))^(1/alphaa);

            %Abbreviation ThetacH
            ThetacH=1/(zetah*AL)*((1-xiw*(pitildew/piw)^(1/(1-lambdaw)))/(1-xiw))^(1-lambdaw*(1+sigmaL))...
                *(wbar*(1-tauy)/(lambdaw*(1+tauw)))*(1-betaa*xiw*(pitildew/piw)^(lambdaw*(1+sigmaL)/(1-lambdaw)))...
                *(zetac/(muzplus-b)*(muzplus-betaa*b)/(pc*(1+tauc)))/(1-betaa*xiw*(pitildew/piw)^(1/(1-lambdaw)));

            %Abbreviation Thetay
            Thetay=phalo^(lambdad/(lambdad-1))*epsilon*(kh/(mupsi*muzplus))^alphaa...
                /(1+phalo^(lambdad/(lambdad-1))*(1/mc-phalo^(lambdad/(1-lambdad))));

            %alphaa abbreviations
            alpha1=omegac*(pc/pmc)^etac*ThetacH*phalomc^(lambdamc/(1-lambdamc));

            alpha2=(1-Rstar*s/(pibar*muzplus))*etaa*Thetay/varphi...
                +omegai*(pinvest/pmi)^etai*kh*(1-(1-delta)/(muzplus*mupsi))/Upsilon*phalomi^(lambdami/(1-lambdami));

            alpha3=px^(1-etaf)/Rnustar-omegax*((omegax*pmx^(1-etax)+1-omegax)^(1/(1-etax))/pmx)^etax...
                *phalox^(lambdax/(1-lambdax))*px^(-etaf)*phalomx^(lambdamx/(1-lambdamx));

            alpha4=(1-etag)*Thetay-pinvest^etai*kh*(1-(1-delta)/(muzplus*mupsi))/Upsilon*(1-omegai);

            alpha5=(1-omegac)*pc^etac*ThetacH;

            alpha6=(omegax*pmx^(1-etax)+1-omegax)^(etax/(1-etax))*(1-omegax)...
                *phalox^(lambdax/(1-lambdax))*px^(-etaf);

            %Hours
            H=((alpha5+alpha6*alpha1/alpha3)/(alpha4-alpha6*alpha2/alpha3))^(1/(1+sigmaL));


            %Further steady states
            kbar=kh*H;
            ybar=Thetay*H;
            g=etag*ybar;
            q=varphi/pc;
            i=kbar*(1-(1-delta)/(muzplus*mupsi))/Upsilon;
            ystar=1/alpha3*(alpha1*H^(-sigmaL)+alpha2*H);
            x=px^(-etaf)*ystar;
            c=ThetacH*H^(-sigmaL);
            phi=ybar*(1/mc-phalo^(lambdad/(1-lambdad)));
            a=etaa*ybar;
            k=kbar*u;
            psizplus=(zetac*(muzplus-betaa*b)/(muzplus-b))/(c*pc*(1+tauc));

            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)));

            Kw=(zetah*H^(1+sigmaL))/(1-betaa*xiw*(pitildew/piw)^(lambdaw*(1+sigmaL)/(1-lambdaw)));
            Fw=(psizplus/lambdaw*H*(1-tauy)/(1+tauw))/(1-betaa*xiw*(pitildew/piw)^(1+lambdaw/(1-lambdaw)));


            smallh=whalo^(lambdaw/(1-lambdaw))*H;

            %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;
            pitildewpi=pitildew/piw;


            %Data series steady states
            data_pid=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=400*log(pic);
            data_pii=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_gdiff=0;
            
            data_Hhat=0;
            
            ylessd=ybar;
            
            Slev=1;
            nxdivGDP = 0;
            GDP12q_diff=12*100*log(muzplus);
            GDP4q_diff=4*100*log(muzplus);
            pic4q=data_pic;

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %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);
                %old Calvo wage setting variables
                log(pitildewpi); log(Kw); log(Fw); log(whalo);
                %log(whalohalo);
                log(smallh);
                aofu; log(piw);
                %log(intensity); 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); log(gamma);
                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_spreaddiff;
                data_gdiff;
                log(ylessd); nxdivGDP; log(Slev); data_Hhat; GDP12q_diff;
                GDP4q_diff; pic4q];
            
           %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,smallh,q,i,c,x, ...
                    psizplus,ybar,ystar,cm,xm,im,delta,varphi,AL];

                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
            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.');
    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, dude');
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	)=	lambdaw     	;
M_.params(	59	)=	wbar        	;
M_.params(	60	)=	pic         	;
M_.params(	61	)=	psizplus    	;
M_.params(	62	)=	H           	;
M_.params(	63	)=	q           	;
M_.params(	64	)=	pimi        	;
M_.params(	65	)=	mcmi        	;
M_.params(	66	)=	mc          	;
M_.params(	67	)=	pix         	;
M_.params(	68	)=	mcmx        	;
M_.params(	69	)=	mcx         	;
M_.params(	70	)=	pimx        	;
M_.params(	71	)=	muz         	;
M_.params(	72	)=	ybar        	;
M_.params(	73	)=	u           	;
M_.params(	74	)=	kappaw      	;
M_.params(	75	)=	AL          	;
M_.params(	76	)=	sigmaL      	;
M_.params(	77	)=	Rk          	;
M_.params(	78	)=	varphi      	;
M_.params(	79	)=	ystar       	;
M_.params(	80	)=	mcmc        	;
M_.params(	81	)=	px          	;
M_.params(	82	)=	pimc        	;
M_.params(	83	)=	pmc         	;
M_.params(	84	)=	pmi         	;
M_.params(	85	)=	kbar        	;
M_.params(	86	)=	i           	;
M_.params(	87	)=	c           	;
M_.params(	88	)=	pc          	;
M_.params(	89	)=	pkprime     	;
M_.params(	90	)=	pinvest     	;
M_.params(	91	)=	pmx         	;
M_.params(	92	)=	pii         	;
M_.params(	93	)=	pistar      	;
M_.params(	94	)=	Rstar       	;
M_.params(	95	)=	rhomuz      	;
M_.params(	96	)=	rhoepsilon  	;
M_.params(	97	)=	rhoUpsilon  	;
M_.params(	98	)=	rhozetac    	;
M_.params(	99	)=	rhozetah    	;
M_.params(	100	)=	rhog        	;
M_.params(	101	)=	rhotauy     	;
M_.params(	102	)=	rhophitilde 	;
M_.params(	103	)=	phitilde    	;
M_.params(	104	)=	sig_muz     	;
M_.params(	105	)=	sig_epsilon 	;
M_.params(	106	)=	sig_upsilon 	;
M_.params(	107	)=	sig_zetac   	;
M_.params(	108	)=	sig_zetah   	;
M_.params(	109	)=	sig_phitilde	;
M_.params(	110	)=	sig_epsR    	;
M_.params(	111	)=	sig_pitarget	;
M_.params(	112	)=	sig_g       	;
M_.params(	113	)=	sig_tauy    	;
M_.params(	114	)=	s           	;
M_.params(	115	)=	Stilde      	;
M_.params(	116	)=	Sprimetilde 	;
M_.params(	117	)=	xiw         	;
M_.params(	118	)=	rkbar       	;
M_.params(	119	)=	Rf          	;
M_.params(	120	)=	Rnustar     	;
M_.params(	121	)=	x           	;
M_.params(	122	)=	Rx          	;
M_.params(	123	)=	varkappad   	;
M_.params(	124	)=	varkappax   	;
M_.params(	125	)=	varkappamx  	;
M_.params(	126	)=	varkappami  	;
M_.params(	127	)=	varkappamc  	;
M_.params(	128	)=	varkappaw   	;
M_.params(	129	)=	thetaw      	;
M_.params(	130	)=	pibreve     	;
M_.params(	131	)=	pitildedpi  	;
M_.params(	132	)=	Kd          	;
M_.params(	133	)=	Fd          	;
M_.params(	134	)=	phalo       	;
M_.params(	135	)=	pitildexpi  	;
M_.params(	136	)=	Kx          	;
M_.params(	137	)=	Fx          	;
M_.params(	138	)=	phalox      	;
M_.params(	139	)=	pitildemcpi 	;
M_.params(	140	)=	Kmc         	;
M_.params(	141	)=	Fmc         	;
M_.params(	142	)=	phalomc     	;
M_.params(	143	)=	cm          	;
M_.params(	144	)=	pitildemipi 	;
M_.params(	145	)=	Kmi         	;
M_.params(	146	)=	Fmi         	;
M_.params(	147	)=	phalomi     	;
M_.params(	148	)=	im          	;
M_.params(	149	)=	pitildemxpi 	;
M_.params(	150	)=	Kmx         	;
M_.params(	151	)=	Fmx         	;
M_.params(	152	)=	phalomx     	;
M_.params(	153	)=	xm          	;
M_.params(	154	)=	pitildewpi  	;
M_.params(	155	)=	Kw          	;
M_.params(	156	)=	Fw          	;
M_.params(	157	)=	whalo       	;
M_.params(	158	)=	smallh      	;
M_.params(	159	)=	aofu        	;
M_.params(	160	)=	taud        	;
M_.params(	161	)=	taux        	;
M_.params(	162	)=	taumc       	;
M_.params(	163	)=	taumi       	;
M_.params(	164	)=	taumx       	;
M_.params(	165	)=	rhotaud     	;
M_.params(	166	)=	rhotaux     	;
M_.params(	167	)=	rhotaumc    	;
M_.params(	168	)=	rhotaumi    	;
M_.params(	169	)=	rhotaumx    	;
M_.params(	170	)=	piw         	;
M_.params(	171	)=	taub        	;
M_.params(	172	)=	pitarget    	;
M_.params(	173	)=	etaa        	;
M_.params(	174	)=	a11         	;
M_.params(	175	)=	a12         	;
M_.params(	176	)=	a13         	;
M_.params(	177	)=	a21         	;
M_.params(	178	)=	a22         	;
M_.params(	179	)=	a23         	;
M_.params(	180	)=	a24         	;
M_.params(	181	)=	a31         	;
M_.params(	182	)=	a32         	;
M_.params(	183	)=	a33         	;
M_.params(	184	)=	a34         	;
M_.params(	185	)=	sig_ystar   	;
M_.params(	186	)=	c21         	;
M_.params(	187	)=	sig_pistar  	;
M_.params(	188	)=	c24         	;
M_.params(	189	)=	c31         	;
M_.params(	190	)=	c32         	;
M_.params(	191	)=	sig_Rstar   	;
M_.params(	192	)=	c34         	;
M_.params(	193	)=	SppScaled   	;
M_.params(	194	)=	muzPercent  	;
% M_.params(	195	)=	mupsiPercent	;
M_.params(	195	)=  intensity_target ;
M_.params(	196	)=  iy_target ;
M_.params(	197	)=  xy_target ;
M_.params(198)   =  hours_target ;
M_.params(199)   =  ylessd ;
M_.params(200)   =  sigmaLScaled ;
M_.params(201)   =  work_cap_para ;


do_steady_print=0;
if check==0,
    if do_steady_print==1,
        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.4f\n\n',kbar/ybar,c/ybar,i/ybar,x/ybar,ybar/ystar,R/pibar)
        fprintf(' smallh = %6.4f \n\n',smallh)
        fprintf(' ****************************\n');
        fprintf(' k*Pk/y = %6.4f, c*Pc/y = %6.4f, i*Pinvest/y = %6.4f, x*Px/y = %6.4f \n\n',kbar*pkprime/ybar,c*pc/ybar,i*pinvest/ybar,x*px/ybar);%,n/(pkprime*kbar))

        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('sigmaa = %6.4f, Fomegabar = %6.4f, mu = %6.4f, we = %6.4f, gamma = %6.4f\n\n', sigmaa,Fomegabar,mu,we,gamma)


        %fprintf('\n*** Labor block of s.s. values below *** \n\n')
        %fprintf('m = %6.4f, f = %6.4f, unemp = %6.4f, intensity = %6.4f\n', m,f,unemp,intensity )
        %fprintf('H = %6.4f, Jw = %6.4f, Vw = %6.4f, vtilde = %6.4f\n',H,Jw,Vw,vtilde )
        %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('rho = %6.4f, iota = %6.4f, kappa = %6.4f, BigN = %6.4f\n', rho,iota,kappa,BigN )
        %fprintf('wbar = %6.4f\n\n',wbar)

        %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, kappa = %6.4f, eta = %6.4f, bu = %6.4f\n\n',iota,sigma_match,kappa,eta,bu )

        %recruiting_costs_to_y=recruiting/ybar
        %job_lasting_in_quarters=1/(1-rho)
        %job_finding_duration_in_quarters=1/f
        %if varkappaw~=0 | thetaw~=1
        %replacement_ratio=(1-tauy)*bu/(((1-tauy)/(1+tauw)*w*wbar*intensity*1/BigN*(G0+G1+G2+G3)...
        %            -(zetah*AL*intensity^(1+sigmaL)/((1+sigmaL)*psizplus))))
        %else
        %replacement_ratio=(1-tauy)*bu/((((1-tauy)/(1+tauw))*w*wbar*intensity-(zetah*AL*intensity^(1+sigmaL)/((1+sigmaL)*psizplus))))
        %end

        varnames=[  ' ydiffU            '
            ' pidU              '
            ' RU                '
            ' wdiffU            '
            ' cdiffU            '
            ' idiffU            '
            ' qdiffU            '
            ' HdiffU            '
            ' xdiffU            '
            ' impdiffU          '
            ' picU              '
            ' piiU              '
            ' ystardiffU        '
            ' pistarU           '
            ' RstarU            '];

        means=[  100*(muzplus-1)
            400*log(pibar)
            400*log(R)
            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*log(Rstar)];

        %disp('   ');
        %disp('Model means');
        %disp([varnames, num2str(means)]);
        %pause
        keyboard
    end
end