%Here we solve the core equations of the nonlinear steady state. By
%substitution we have reduced the system of equation to the set of
%equations listed below


function ff=mythird_steadystate_solvecore(xxx,tauk,muzplus_p,alpha,lambdad,...
            etag,etax,omegax,etaf,omegai,etai,beta,b, sigmaL,omegac,...
            etac, etaa, lambdamc,taub,pitarget_p, delta,tau,...
          lambdami,lambdamx,lambdax,mu,Fomegabar,xiw,lambdaw,xid,xix,ximc,ximi,ximx,muz_p,...
            xy_target,nk_target,AL) 

        

%backing out the guesses for the variables to be solved for
sigma=xxx(1);
omegabar=xxx(2);
Rk=xxx(3);
kbar=xxx(4);
wbar=xxx(5);
H=xxx(6);
ystar=xxx(7);
rkbar=xxx(8);
varphi=xxx(9);

%Relative import prices
Rstar=(muzplus_p/beta-taub)/(1-taub)*pitarget_p/(1+tau);

mcmx=1/lambdamx*((1-ximx*(1)^(1/(1-lambdamx)))/(1-ximx))^((1-lambdamx))...
    *(1-beta*ximx*(1)^(lambdamx/(1-lambdamx)))...
    /(1-beta*ximx*(1)^(1/(1-lambdamx)));
mcmc=1/lambdamc*((1-ximc*(1)^(1/(1-lambdamc)))/(1-ximc))^((1-lambdamc))...
    *(1-beta*ximc*(1)^(lambdamc/(1-lambdamc)))...
    /(1-beta*ximc*(1)^(1/(1-lambdamc)));
mcmi=1/lambdami*((1-ximi*(1)^(1/(1-lambdami)))/(1-ximi))^((1-lambdami))...
    *(1-beta*ximi*(1)^(lambdami/(1-lambdami)))/(1-beta*ximi*(1)^(1/(1-lambdami)));

piw=muzplus_p*pitarget_p;
R=(muzplus_p/beta-taub)/(1-taub)*pitarget_p;
mcx=1/lambdax*((1-xix*(1)^(1/(1-lambdax)))/(1-xix))^((1-lambdax))...
    *(1-beta*xix*(1)^(lambdax/(1-lambdax)))...
    /(1-beta*xix*(1)^(1/(1-lambdax)));
%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));

mupsi=(muzplus_p/muz_p)^((1-alpha)/alpha);
mc=1/lambdad*(1-beta*xid*(1)^(lambdad/(1-lambdad)))...
    /(1-beta*xid*(1)^(1/(1-lambdad)))...
    *((1-xid*(1)^(1/(1-lambdad)))/(1-xid))^(1-lambdad);

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

phalo=((1-xid*(1)^(1/(1-lambdad)))/(1-xid))^(1-lambdad)...
    /((1-xid*(1)^(lambdad/(1-lambdad)))/(1-xid))^((1-lambdad)/lambdad);
phalomi=((1-ximi*(1)^(1/(1-lambdami)))/(1-ximi))^(1-lambdami)...
    /((1-ximi*(1)^(lambdami/(1-lambdami)))/(1-ximi))^((1-lambdami)/lambdami);
s=1;
phalomc=((1-ximc*(1)^(1/(1-lambdamc)))/(1-ximc))^(1-lambdamc)...
    /((1-ximc*(1)^(lambdamc/(1-lambdamc)))/(1-ximc))^((1-lambdamc)/lambdamc);

phalomx=((1-ximx*(1)^(1/(1-lambdamx)))/(1-ximx))^(1-lambdamx)...
    /((1-ximx*(1)^(lambdamx/(1-lambdamx)))/(1-ximx))^((1-lambdamx)/lambdamx);
phalox=((1-xix*(1)^(1/(1-lambdax)))/(1-xix))^(1-lambdax)...
    /((1-xix*(1)^(lambdax/(1-lambdax)))/(1-xix))^((1-lambdax)/lambdax);


%Rental rate on capital
equ1=rkbar-(Rk*mupsi*((1-omegai+omegai*( Rstar*varphi/mcmi)^(1-etai))^(1/(1-etai)))/pitarget_p-(1-delta)*( (1-omegai+omegai*( Rstar*varphi/mcmi)^(1-etai))^(1/(1-etai)))-tauk*delta*mupsi*((1-omegai+omegai*( Rstar*varphi/mcmi)^(1-etai))^(1/(1-etai)))/pitarget_p)/(1-tauk);

%Real Wage
equ2=wbar-(mc/( (1/(1-alpha))^(1-alpha)*(1/alpha)^alpha*rkbar^alpha*R^(1-alpha)/1))^(1/(1-alpha));

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

% Capital labor ratio
equ4=kbar/H-muzplus_p*(( mupsi^alpha*wbar*R)/((1-alpha)*mc))^(1/alpha);

%Trade balance
equ5=(1-Rstar*s/pitarget_p/muzplus_p)*etaa*phalo^(lambdad/(lambdad-1))*(kbar/H/(mupsi*muzplus_p))^alpha...%assets
    /(1+phalo^(lambdad/(lambdad-1))*(1/mc-phalo^(lambdad/(1-lambdad))))*H...
    +varphi*Rstar*(... %imports
    omegac*(((1-omegac+omegac*( Rstar*varphi/mcmc)^(1-etac))^(1/(1-etac)))/( Rstar*varphi/mcmc))^(etac)*phalomc^(lambdamc/(1-lambdamc))...%import consumption
    *1/(AL)*((1-xiw*(pitildew/piw)^(1/(1-lambdaw)))/(1-xiw))^(1-lambdaw*(1+sigmaL))...
    *(wbar*(1)/(lambdaw*(1)))*(1-beta*xiw*(pitildew/piw)^(lambdaw*(1+sigmaL)/(1-lambdaw)))...
    *(1/(muzplus_p-b)*(muzplus_p-beta*b)/(( (1-omegac+omegac*( Rstar*varphi/mcmc)^(1-etac))^(1/(1-etac)))*(1)))/(1-beta*xiw*(pitildew/piw)^(1/(1-lambdaw)))*H^(-sigmaL)...
    +phalomi^(lambdami/(1-lambdami))*omegai*(((1-omegai+omegai*( Rstar*varphi/mcmi)^(1-etai))^(1/(1-etai)))/( Rstar*varphi/mcmi))^(etai)*kbar*(1-(1-delta)/(muzplus_p*mupsi))...%import investment
    +phalomx^(lambdamx/(1-lambdamx))*omegax*(1/(Rstar*varphi/mcmx)*(omegax*(Rstar*varphi/mcmx)^(1-etax)+1-omegax)^(1/(1-etax)))^etax...%import export
    *phalox^(-lambdax/(lambdax-1))*( ( R/(mcx*varphi))*(omegax*(Rstar*varphi/mcmx)^(1-etax)+1-omegax)^(1/(1-etax)))^(-etaf)*ystar )...
    -varphi*(( R/(mcx*varphi))*(omegax*(Rstar*varphi/mcmx)^(1-etax)+1-omegax)^(1/(1-etax)))^(1-etaf)*ystar; %exports


%Feasibility
equ6=(1-etag)*(phalo^(lambdad/(lambdad-1))*(kbar/H/(mupsi*muzplus_p))^alpha)...
    /(1+phalo^(lambdad/(lambdad-1))*(1/mc-phalo^(lambdad/(1-lambdad))))*H...    %(1-etag)*( (phalo^(lambdad/(lambdad-1))*(kbar/H/(mupsi*muzplus_p))^alpha)/(1+phalo^(lambdad/(lambdad-1))*(1/mc-phalo^(lambdad/(1-lambdad))))*H)
    -(omegax*(Rstar*varphi/mcmx)^(1-etax)+1-omegax)^(etax/(1-etax))*(1-omegax)*phalox^(-lambdax/(lambdax-1))*( ( R/(mcx*varphi))*(omegax*(Rstar*varphi/mcmx)^(1-etax)+1-omegax)^(1/(1-etax)))^(-etaf)*ystar...  % exports
    -(Gamma-Gam_muG)*Rk*((1-omegai+omegai*( Rstar*varphi/mcmi)^(1-etai))^(1/(1-etai)))*kbar/(pitarget_p*muzplus_p)...%financial frictions monitoring costs
    -(1-omegai)*( (1-omegai+omegai*( Rstar*varphi/mcmi)^(1-etai))^(1/(1-etai)))^(etai)*kbar*(1-(1-delta)/(muzplus_p*mupsi))...%investment
    -(1-omegac)*(( (1-omegac+omegac*( Rstar*varphi/mcmc)^(1-etac))^(1/(1-etac)))^etac)*1/(AL)*((1-xiw*(pitildew/piw)^(1/(1-lambdaw)))/(1-xiw))^(1-lambdaw*(1+sigmaL))... %consumption
    *(wbar*(1)/(lambdaw*(1)))*(1-beta*xiw*(pitildew/piw)^(lambdaw*(1+sigmaL)/(1-lambdaw)))...
    *(1/(muzplus_p-b)*(muzplus_p-beta*b)/(( (1-omegac+omegac*( Rstar*varphi/mcmc)^(1-etac))^(1/(1-etac)))*(1)))/(1-beta*xiw*(pitildew/piw)^(1/(1-lambdaw)))*H^(-sigmaL);


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

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


% old incorrect version: equ11=x*(( R/(mcx*varphi))*(omegax*(Rstar*varphi/mcmx)^(1-etax)+1-omegax)^(1/(1-etax)))/ybar-xy_target, where
% x=(( R/(mcx*varphi))*(omegax*(Rstar*varphi/mcmx)^(1-etax)+1-omegax)^(1/(1-etax)))^(-etaf)*ystar (latter statement true); 
equ9=(( R/(mcx*varphi))*(omegax*(Rstar*varphi/mcmx)^(1-etax)+1-omegax)^(1/(1-etax)))^(-etaf)*ystar*(( R/(mcx*varphi))*(omegax*(Rstar*varphi/mcmx)^(1-etax)+1-omegax)^(1/(1-etax)))*varphi/((phalo^(lambdad/(lambdad-1))*(kbar/H/(mupsi*muzplus_p))^alpha)/(1+phalo^(lambdad/(lambdad-1))*(1/mc-phalo^(lambdad/(1-lambdad))))*H)-xy_target; 

%constructing the matrix to be solved with e.g. fsolve
ff=[equ1
    equ2
    equ3
    equ4
    equ5
    equ6
    equ7
    equ8
    equ9];