function [ys,check] = No_exp_steadystate(ys,exo)
% function [ys,check] = NK_baseline_steadystate(ys,exo)
% computes the steady state for the NK_baseline.mod and uses a numerical
% solver to do so
% Inputs: 
%   - ys        [vector] vector of initial values for the steady state of
%                   the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%
% Output: 
%   - ys        [vector] vector of steady state values fpr the the endogenous variables
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impos restriction on parameters)

global M_ 

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = deblank(M_.param_names(ii,:));
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;


%% Enter model equations here

options=optimset(); % set options for numerical solver

% the steady state computation follows FV (2006), section 4.1
pi=piss;
% u=1;
% q=1;
% d=1;
% phi=1;
% m=0;
% zeta=1;
% LambdaYd= (LambdaA+alppha*Lambdamu)/(1-alppha);
% mu_z=exp(LambdaYd);
% mu_I=exp(Lambdamu);
% mu_A=exp(LambdaA);

%set the parameter Lambdax
% Lambdax=mu_z;
% 
%set the parameter gammma1
% L = 1 - u_s - u_u - n_s - n_u;
% if L>0.5 % parameter violates restriction; Preventing this cannot be implemented via prior restriction as it is a composite of different parameters and the valid prior region has unknown form
%     check=1; %set failure indicator
%     return; %return without updating steady states
% end

% 
% r=1*gammma1;
% R=1+(PI*mu_z/betta-1);
% 
% set Rbar
% Rbar=R;
% 
% PIstar=((1-thetap*PI^(-(1-epsilon)*(1-chi)))/(1-thetap))^(1/(1-epsilon));
% PIstarw=((1-thetaw*PI^(-(1-chiw)*(1-eta))*mu_z^(-(1-eta)))/(1-thetaw))^(1/(1-eta));

% model;
% ////*********** PATIENT HHs ********************************************************
lam_p = (eps_c*beta_p*(c_p - zeta*c_p)^(-sigma_c) )/P; % 1. 

% eps_h*(h_p^(-sigma_h)) - lam_p * P_h + beta_p * lam_p *(1-delta_h)*P_h = 0; % 2. 
% 
% (1+r_d) = lam_p/(beta_p * lam_p+1) ; //  3. d/dd
% 
% chi*(l_p^(-sigma_l)) = L_p;                 // 4.
% 
% P*c_p + P_h*i_h_p + i_d + ks*n_u = w_s*n_s + w_u*n_u + tau*(u_s+u_u) + J_R/gamma_p - t_p;  // 5. Patient household budget constraint
% 
% ////*********** IMPATIENT HHs ********************************************************
% eps_c*beta_i*(c_i - zeta*c_i(-1))^(-sigma_c) = lam_i*P;  // 6.
% 
% beta_i*eps_h*(gamma_h*h_i)^(-sigma_h) - lam_i*P_h +beta_i*(1-delta_h)*lam_i(+1)*P_h(+1) 
%     - lam_ib*((gamma_h*m_i*(P_h(+1))*(1-delta_h))/(gamma_b*(1+r_ib))) =0;      // 7. Housing demand (d/d_hi)
% beta_i*eps_h*((1-gamma_h)*h_is)^(-sigma_h) - lam_i*p_hs +beta_i*(1-delta_h)*lam_i(+1)*p_hs(+1) 
%     + lam_ib*(((1-gamma_h)*m_f*(P_h(+1))*(1-delta_h))/((1-gamma_b)*(1+r_ibs)))=0; // 8. Housing demand (d/dh_is)
% 
% lam_i = lam_ib - beta_i * lam_i(+1) * (1+r_ib);           // 9. (d/db)
% lam_ib= lam_i - beta_i* ( E(+1)/E )*lam_i(+1) * (1+r_ibs);       // 10. (d/dbs)
% 
% L_i = chi*(l_i^(-sigma_l));   // 11.
% 
% P*c_i + P_h*i_h_i + ks*n_u = w_s*n_s + w_u*n_u + tau*(u_s+u_u) + i_b + E*i_bs - t_i;    // 12. Budget constraint in Gerali manca t_i
% 
% gamma_h*(m_i*(P_h(+1))*(1-delta_h)*h_i)/(1+r_ib) + (1-gamma_h)*(m_f*(p_hs(+1))*(1-delta_h)*h_is)/(1+r_ib) = gamma_b*b + (1-gamma_b)*bs;     // 13. Borrowing constraint impatient household 
% 
% ////************ ENTREPRENEURS *********************************************************
% (c_e - zeta*c_e(-1))^(-sigma_c) = lam_e*P;         // 14. FOC consumption eps_c? beta_e?
% 
% y_e = (y_eu)^(rho_ye)*(y_es)^(1-rho_ye);    // 15. Production technology
% y_es = (n_s)^(alfa)*(k)^(1-alfa);           // 16. Non distinguo tra capitle s ed u
% y_eu = (n_u)^(alfa)*(k)^(1-alfa);           // 17. Non distinguo tra capitle s ed u
% 
% PRu = rho_ye*(y_e/y_eu);           // 18. Prices of intermediate goods (Unskilled)
% PRs = (1-rho_ye)*(y_e/y_es);       // 19. Prices of intermediate goods (Skilled)
% 
% (lam_eb* m_e * p_k * (1-delta_k) )/(1+r_be) + beta_e * lam_e(+1) * ( r_k + p_k(+1)*(1-delta_k) ) = lam_e*p_k;  // 20. FOC capital exp(r_k(+1))*exp(u(+1))  exp(s_e) 
% 
% p_w*r_k = (1-alfa) * k^(-alfa) * ( n_u^rho_ye * n_s^(1-rho_ye) )^alfa ;  // 21. Capital rental rate dy_e/dk
% 
% lam_eb = beta_e * lam_e(+1) * (1+r_be) - lam_e ;  // 22. FOC credit demand
% 
% (1+r_be) * be = m_e * p_k(+1) * k * (1-delta_k);  // 23. Borrowing constraint entrepreneurs
% 
% P*c_e + p_k * k + w_s*n_s + w_u*n_u = y_e*p_w + i_be + (1-delta_k)*k(-1);   // 24. budget constraint entrepreneurs (22) k_e
% 
% /////LABOUR MARKET\\\\\
% L = 1- n_s - n_u - u_s - u_u;                                     // 25. Workers pool
% exp(L) = gamma_p*exp(L_p) + gamma_i*exp(L_i);              //E' CORRETTO SPECIFICARE COSI'? // 26.
% 
% n_s = (1-phi)*(n_s(-1) + p_s(-1)*(u_s(-1) + ks(-1)*n_s(-1)));      // 27. Skilled workers law of motion
% n_u = (1-phi)*(n_u(-1)*(1-p_s(-1)*ks(-1)) + p_u(-1)*u_u(-1));      // 28. Unskilled workers law of motion
% 
% m_s = eta_s*(v_s^(lam_w))*(u_s + ks*n_s)^(1-lam_w);                     // 29. Skilled matching function
% m_u = eta_u*(v_u^(lam_w))*(u_u)^(1-lam_w);                                        // 30. Unskilled matching function
% 
% p_s = m_s/(u_s + ks*n_s);                         // 31. Job finding prob. Skilled
% p_u = m_u/u_u;                                            // 32. Job finding prob. Unskilled
% 
% q_s = m_s/v_s;                                            // 33. Vacancy filling prob. Skilled
% q_u = m_u/v_u;                                            // 34. Vacancy filling prob. Unskilled
% 
% theta_s = v_s/(u_s + ks*n_s);                     // 35. Labour market tightness Skilled
% theta_u = v_u/u_u;                                        // 36. Labour market tightnes Unskilled
% 
% theta_u/theta_s = gamma_s/gamma_u;                               // 37. Equilibrium condition (tightness/costs)
% 
% lam_w = (lam_p+lam_i)/2;                                    // 38.
% 
% ((beta_w*chi*L^(-sigma_l)) - lam_w*tau)/p_s = (1-phi)*lam_w(+1)*w_s(+1);        // 39. (d/du_s) Participating condition skilled (util marg di l - util marg di u)
% (beta_w*chi*(L^(-sigma_l)) - (lam_w*tau))/p_u = (1-phi)*lam_w(+1)*w_u(+1);      // 40. (d/du_u) Participating condition unskilled
% 
% ks = (1-phi)*(lam_w(+1)/lam_w)*p_s*(w_s(+1)-w_u(+1));                      // 41. (d/dks)
% 
% //Nash
% w_u = (1-mu)*(tau + ks) + mu*(alfa*( (y_e/n_u)*PRu) ) + mu*(gamma_u/q_u)-(lam_w(+1)/lam_w)*(1-mu)*(1-phi)*w_u(+1)*((1-p_s*ks)-p_u);           // 42. w_u eq.
% w_s = mu*( (alfa*((y_e/n_s)*PRs) )+(gamma_s/q_s))+(1-mu)*tau -(1-mu)*(1-p_s)*(1-phi)*(lam_w(+1)/lam_w)*w_s(+1);                               // 43. w_s eq.
% 
% //Domanda di lavoro (Entrepreneurs)
% gamma_s/q_s = alfa*((p_w*y_e)/n_s) - w_s + (lam_e(+1)/lam_e)*(1-phi)*(w_s(+1)*(1+m_s(+1)));        // 44. Job creation condition (Skilled) (d/dn_s)
% gamma_u/q_u = alfa*((p_w*y_e)/n_u) - w_u + (lam_e(+1)/lam_e)*(1-phi)*(w_u(+1)*((1-p_s)+m_u));      // 45. Job creation condition (Unskilled) (d/dn_u)
% 
% ////*********** /////PRODUCERS\\\\\ *****************************************************
% k = (1-delta_k) * k(-1) + ( 1 - kappa_i/2 * (I*ee_qk/I(-1) - 1)^2 ) * I ; // 46. Capital pag 39 d/dk
% 
% 1 = p_k * ( 1 -  kappa_i/2 * (I*ee_qk/I(-1) - 1)^2  - kappa_i * (I*ee_qk/I(-1) - 1) * I*ee_qk/I(-1) ) 
%   + beta_e * lam_e(+1) / lam_e * p_k(+1) *   kappa_i * (I(+1)*ee_qk(+1)/I - 1) * ee_qk(+1) * (I(+1)/I)^2 ; // 47. Real price of capital d/di
% 
% i_h_p = h_p - (1-delta_h)*h_p(-1);                                 // 48. Housing
% i_d = d - (1+r_d(-1))*d(-1);                                  //49. Deposit
% i_h_i= h_i - (1-delta_h)*h_i(-1);                                  //50. Impatient Housing
% i_h_is= h_is - (1-delta_h)*h_is(-1);                               //51. Foreign Impatient Housing
% i_bs= bs - (1+r_ibs(-1))*bs(-1);                              //52. Foreign Bonds
% i_b = b - (1+r_ib(-1))*b(-1);                                 //53. Domestic Bonds
% i_be= be - (1+r_be(-1))*be(-1);                               //54. Entrepreneurs Bonds
% 
% ////************* /////BANKS\\\\\ ****************************************************************
% k_b = (1-delta_kb) * k_b(-1) + omega_b*j_B(-1);                         //55. Bank capital accumulation (law of motion)
% 
% //FOCs
% r_ib  = r_cb - ck_i; //56. FOC wholesale branch, where r_cb is policy rate and r_ib the net wholsale loan rate
% r_ibs = r_cb - ck_is; //57. FOC wholesale branch, where r_cb is policy rate and r_ibs the net foreign wholsale loan rate
% r_be  = r_cb - ck_e; //58. FOC wholesale branch, where r_cb is policy rate and r_be the net wholsale loan rate
% //Funzioni di costo del capitale (quadratiche)
% ck_i  = kappa_kb * ( k_b / gamma_i*b          - vi ) * (k_b / gamma_i*b)^2;  // 59.
% ck_is = kappa_kb * ( k_b /(E*gamma_i*bs) - vi)  * (k_b / (E*gamma_i*bs))^2;  // 60.
% ck_e  = kappa_kb * ( k_b / BE           - vi ) * (k_b / BE)^2;  // 61.
% 
% b + E*bs + b_e  = d_b + E*bf + k_b;           //62. Balance sheet constraint (bf=foreign funding)
% 
% //Wholsale Bank's profits
% j_B = (r_ib*b) + (r_ibs*bs) + (r_be*b_e) - (r_d*d) - (r_x*E*bf) - ck_i - ck_is - ck_e;  // 63.
% 
% //Loans = Deposits identities
% gamma_b * d_b = gamma_p * d;              //64.
% gamma_b * b_e = gamma_e * be;             //65.
% gamma_b * bi  = gamma_i * (b + bs);  //66. gamma_b = bank share ?
% 
% ////***********  FINAL GOODS PRODUCERS AND RETAILERS  **************************************************************
% J_d = y_d*(1 -(1/p_w)  - (kappa_p/2) * (( P_d/P_d(-1) ) - ( pi(-1)^ind_p * piss^(1-ind_p) ))^2 ) ; //67. Aggregate retail profits (Patient households)
%  
% J_f = y_f*(1 -(1/p_w)  - (kappa_p/2) * (( P_f/P_f(-1) ) - ( pi(-1)^ind_p * piss^(1-ind_p) ))^2 ) ; //68. Aggregate retail profits (Patient households)
% 
% J_R = J_d + J_f;   //69.
% 
% // New Keynesian Phillips curve domestic:
% 1 - eps_y + eps_y / p_w - kappa_p * ( ( P_d/P_d(-1) )  - ( pi(-1)^ind_p * piss^(1-ind_p) )) * ( P_d/P_d(-1) ) 
%     + beta_p*(lam_p(+1)/lam_p)*kappa_p * (pi(+1) - ( pi^ind_p * piss^(1-ind_p) )) * ( P_d(+1)/P_d ) * (y_d(+1)/y_d) = 0;  //70.
% 
% // New Keynesian Phillips curve import:
% 1 - eps_y + eps_y / p_w - kappa_p * ( P_f/P_f(-1) )  - ( pi(-1)^ind_p * piss^(1-ind_p) ) * ( P_f/P_f(-1) )  
%     + beta_p*(lam_p(+1)/lam_p)*kappa_p * (pi(+1) - ( pi^ind_p * piss^(1-ind_p) )) * ( P_f(+1)/P_f )  * (y_f(+1)/y_f) = 0;  //71.
% 
% /////FOREINGN SECTOR (S-VAR)\\\\\
% //external_function(name = funcname);
% //Copia bvar_dsge code (y_x, r_x, p_x) oppure vedi "codes_dynare_wp_buss_unemployment"
% 
% y_x = rho_y*y_x(-1)+e_y;        //72.
% p_x = rho_p*p_x(-1)+e_p;        //73.
% r_x = rho_r*r_x(-1)+e_r;        //74. Foreign interest rate
% 
% E = P_f/P_d;                    //75. Nominal exchange rate - exp(E)*(exp(p_x)/exp(P_d)); -> E=Nominal exchange rate
% rer = E*(p_x/P);           //76. Real exchange rate
% 
% //UIP
% (1+r_cb)/(1+r_x) = (rer(+1)/rer)*(pi(+1)/(p_x(+1)/p_x));  //77. Secondo Gerali r_BF = r_cb -> Foreign interbank rate
% 
% /////PREZZI\\\\\
% P_f = piss*(1-omega) + omega*P_f(-1);        //79.
% P = piss*(1-omega) + omega*P(-1);            //79.
% P_d = piss*(1-omega) + omega*P_d(-1);        //80.
% p_x = piss*(1-omega) + omega*p_x(-1);        //81.
% 
% P/P_d = ( 1 - psi + psi*(E^(1 - mu_p)) )^(1/1-mu_p);        //82.
% 
% pi = P/P(-1);        //83.
% 
% ////************  AGGREGATION & EQUILIBRIUM  ************************************************
% C         = gamma_p * c_p + gamma_i * c_i + gamma_e * c_e;  //84. Aggregate consumption
% BE        = gamma_b * b_e;  //85. Chiude le identità del settore bancario [exp(BE)  = gamma_e * exp(be) - Entrep.]
% IB        = gamma_b * bi;        //86.
% B         = (IB + BE); //87.
% y_e       = y_h + y_f;        //88.
% J_B       = gamma_b * j_B;  // 89.
% K         = gamma_e * k;    // 90.
% H         = gamma_p * h_p + gamma_i *( h_i + h_is ); // 91.
% Y1        = C +    1     * (K-(1-delta_k)*K(-1)) ; //92. Ausiliaria per CB Taylor rule
% 
% (y_f*P_f) + E*((1+r_x(-1))*bf(-1)) = E*(y_x*p_x) + E*bf;    //93. BOP (in domestic currency)
% 
% Y = C + k - (1-delta_k)*k(-1) + gamma_p*i_h_p + gamma_i*(i_h_i+i_h_is) + G + gamma_s*v_s + gamma_u*v_u + n_u;     //94. Vincolo delle risorse aggregate (Demand)
% 
% G = gamma_p*t_p + gamma_i*t_i + G_e;                                                      //95. Public spending constraint 
% 
% Y = (P_d*y_d)/P + (E*(p_x*y_x))/P - (P_f*y_f)/P;       //96. Vincolo delle risorse aggregate (Supply)
% 
% //Production function
% Y = ((1-psi)^(1/mu_y)*y_h^((mu_y - 1)/mu_y) + psi^(1/mu_y)*y_f^((mu_y - 1)/mu_y))^(mu_y/(mu_y-1));  //97. Final good technology
% 
% ////*********** TAYLOR RULE & PROFITS CB *****************************************************                                                  
% (1+r_ib) = (1+r_ib_ss)^(1 - rho_ib ) * (1+r_ib(-1))^rho_ib * (( pi / piss )^phi_pi * (Y1/Y1(-1))^phi_y  ) ^ ( 1 - rho_ib ) * (1+e_r_cb) ;       //98. Mmonetary policy
% 
% ////*********** EXOGENOUS PROCESSES ****************************************************
% eps_c     = 1 - rho_eps_c   *    1      + rho_eps_c   * eps_c(-1)  + e_z;
% eps_h     = 1 - rho_eps_h   *    1      + rho_eps_h   * eps_h(-1)  + e_j;
% m_i      = (1-rho_mi)     *  m_i_ss     + rho_mi     * m_i(-1)     + e_mi;
% m_e      = (1-rho_me)     *  m_e_ss     + rho_me     * m_e(-1)     + e_me;
% m_f      = (1-rho_mf)     *  1          + rho_mf     * m_f(-1)     + e_mf;  //Foreign housing
% ee_qk    =  1-rho_ee_qk   *    1        + rho_ee_qk  * ee_qk(-1)   + e_qk;
% eps_y    = (1-rho_eps_y)  * eps_y_ss    + rho_eps_y  * eps_y(-1)   + e_y_cb;
% eps_l    = (1-rho_eps_l)  * 1           + rho_eps_l  * eps_l(-1)   + e_l;
% 
% G_e      =  1 - rho_G     *    1        + rho_G   * G_e(-1)        + e_g;
% e_y      =  1 - rho_y     *    1        + rho_y   * e_y(-1)        + ee_y;
% e_p      =  1 - rho_p     *    1        + rho_p   * e_p(-1)        + ee_p;
% e_r      =  1 - rho_r     *    1        + rho_r   * e_r(-1)        + ee_r;
% e_r_cb   =  1 - rho_CB    *    1        + rho_CB  * e_r_cb(-1)     + ee_r_cb;
% 
% ////*********** AUXILIARY VARIABLES *****************************************************
% rr_e         = lam_e - beta_e*lam_e(+1)*(1+r_be)/pi(+1);              //A che serve?!
% aux1         = k_b/B;                                                           //A che serve?!
% bm      = (bi(-1)/(bi(-1)+b_e(-1)) * r_ib(-1) + b_e(-1)/(bi(-1)+b_e(-1)) *r_be(-1)) - r_d(-1);       //A che serve?!
% spr_b   =  0.5*r_ib + 0.5*r_be - r_d;                                 //A che serve?!
%                


% [ld,fval,exitflag]=fzero(@(ld)(1-betta*thetaw*mu_z^(eta-1)*PI^(-(1-chiw)*(1-eta)))/(1-betta*thetaw*mu_z^(eta*(1+gammma))*PI^(eta*(1-chiw)*(1+gammma)))...
% -(eta-1)/eta*wstar/(varpsi*PIstarw^(-eta*gammma)*ld^gammma)*((1-h*mu_z^(-1))^(-1)-betta*h*(mu_z-h)^(-1))*...
% ((mu_A*mu_z^(-1)*vp^(-1)*tempvaromega^alppha-tempvaromega*(1-(1-delta)*(mu_z*mu_I)^(-1)))*ld-vp^(-1)*Phi)^(-1),0.25,options);
% if exitflag <1
%     %indicate the SS computation was not sucessful; this would also be detected by Dynare
%     %setting the indicator here shows how to use this functionality to
%     %filter out parameter draws
%     check=1; %set failure indicator
%     return; %return without updating steady states
% end


% l=vw*ld;
% k=tempvaromega*ld;
% x=(1-(1-delta)*(mu_z*mu_I)^(-1))*k;
% yd=(mu_A/mu_z*k^alppha*ld^(1-alppha)-Phi)/vp;
% c=(mu_A*mu_z^(-1)*vp^(-1)*tempvaromega^alppha-tempvaromega*(1-(1-delta)*(mu_z*mu_I)^(-1)))*ld-vp^(-1)*Phi;
% lambda=(1-h*betta*mu_z^(-1))*(1-h/mu_z)^(-1)/c;
% F=yd-1/(1-alppha)*w*ld;
% f=(eta-1)/eta*wstar*PIstarw^(-eta)*lambda*ld/(1-betta*thetaw*mu_z^(eta-1)*PI^(-(1-chiw)*(1-eta)));
% f2=varpsi*d*phi*PIstarw^(-eta*(1+gammma))*ld^(1+gammma)/(1-betta*thetaw*(PI^chiw/PI)^(-eta*(1+gammma))*(wstar/wstar*mu_z)^(eta*(1+gammma)));
% 
% g1=lambda*mc*yd/(1-betta*thetap*PI^((1-chi)*epsilon));
% g2=epsilon/(epsilon-1)*g1;

%% end own model equations

for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'M_.params(' num2str(iter) ') = ' M_.param_names(iter,:) ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
  varname = deblank(M_.endo_names(ii,:));
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end
