
function [ys, check] = tc(ys,exe)
global M_
check=0;
% Chortares, Katirzidis and Tsoukalas CKT(2008)
bet=0.99; %value ss anual real interst rate 4%
hds = 0.333;
hdu = 0.333;
hms = 0.333;
hmu = 0.333;
sigma = 1.4925;  % Elasticity between skilled & unskilled domestic sigma>0 sigma goes to inf (perfect subst)
sigma1 = 1.4925; % Elasticity between skilled & unskilled immigrants
%ze =1; %1.3;   % Productivity difference between  skilled & unskilled, Canova-Ravn: ze=2.
omegau = 1; % Productivity difference between  natives and immigrants, UNSKILLED Labour efficiency. 
omegas = 1; % Productivity difference between  natives and immigrants, SKILLED Labour efficiency.
            % omegau > 1 productive diff in high-skilled worker
            % Kidland & Prescott (1984) omegai=2, Canova-Ravn omegai = 1.3, 
rho = 0.67; % Elasticity between skilled & unskilled labour rho=0 PS, so here ImpS sigma = 1/rho
            % Ottaviano-Peri(2005) [0.5,0.67], Canova-Ravn(1997) [0, 0.25,
            % 0.5]
            % Ben and Pesarides(1991) rh0 = [0.2,2]
rho1 = 0.1; % Elasticity between skilled natives and immigrants rho=0 PS.  
rho2 = 0.1; % Elasticity between skilled natives and immigrants rho=0 PS.  
            % Ottaviano-Peri(2005) for skilled & unskilled rhok = 5 and
            % rhok = 10. No perfect substitutes.
epsi = 6.1;   % Elasticity of subs Intermediate-goods constant 20% markup.
            % CKT(2008) and me 2nd chapter, some authors epsi = 10.

%Production Function
delt = 0.025;       % Depreciation rate
alf = 0.64;         % The Labour Share

%Shocks
zetabar = 1.0;   % US data 1948 - 2003  
zeta = zetabar;

%Dynamic population 
%phibar = 1.008;
gd = 1;             % % Domestic Population Growth
gT = 1;
gm = 1;
gms = 1;
gmu = 1;
gds = 1;
gdu = 1;
gammt = 0.75; %
phim = 0.5; %0.46;    % EU Skilled Population/Total EU Population PwC(2016)
                % phim = Nsm/Nm
phi_t = 0.5; %0.55;   % EU Skilled Population/Total EU Population PwC(2016)
                % phi_t = Nds/Nd

%Shocks 
%rho_zeta = 0.60;  % Technology 
rho_zeta = 0.60;  % Technology 
rho_ad = 0.9048;  % Domestic Consumption 
rho_am = 0.9048;  % Immigration Consumption 

%Prices
varphi_h = 77; %Rotemberg cost

%Monetary Policy
rho_r = 0.82;
phi_pai = 1.5;
phi_y = 0.1; 

pai = 1; 
paih = 1;
Q = bet/(pai*gd);
rk = (gd/bet)-1 +delt; 
choque=1;
rr = choque;
% Dynamics Population
%Households
hd = ((hds*phi_t)^((sigma+1)/sigma) + (hdu*(1-phi_t))^((sigma+1)/sigma))^(sigma/(sigma+1));
hm = ((hms*phim)^((sigma+1)/sigma) + (hmu*(1-phim))^((sigma+1)/sigma))^(sigma/(sigma+1));

%Firms
hs=(omegas*(hds*phi_t*gammt)^(1-rho1) + (hms*phim*(1-gammt))^(1-rho1))^(1/(1-rho1));
hu=(omegau*(hdu*(1-phi_t)*gammt)^(1-rho2) + (hmu*(1-phim)*(1-gammt))^(1-rho2))^(1/(1-rho2));
he=((hs)^(1-rho)+ (hu)^(1-rho))^(1/(1-rho));
rmc=(epsi-1)/epsi;
k = ((1/rk)*(1-alf)*zeta)^(1/alf)*he;
kd =(1/gammt)*k;
yh=zeta*k^(1-alf)*he^(alf);
ind = gd*kd-(1-delt)*kd; % gd=1 for now 
in = (gammt)*ind;
wds =(alf*omegas)*zeta*k^(1-alf)*(phi_t*gammt*hds)^(-rho1)*he^(alf+rho-1)*hs^((rho1-rho));  
wdu =(alf*omegau)*zeta*k^(1-alf)*((1-phi_t)*gammt*hdu)^(-rho2)*he^(alf+rho-1)*hu^(rho2-rho);  
wmu = alf*zeta*k^(1-alf)*((1-phim)*(1-gammt)*hmu)^(-rho2)*he^(alf+rho-1)*hu^(rho2-rho);  
wms = alf*zeta*k^(1-alf)*(phim*(1-gammt)*hms)^(-rho1)*he^(alf+rho-1)*hs^(rho1-rho);  
bd=0;
b = gammt*bd;
dd =yh- wds*hds*phi_t*gammt - wdu*hdu*(1-phi_t)*gammt - kd*rk*gammt - phim*hms*wms*(1-gammt) - (1-phim)*hmu*wmu*(1-gammt) - gammt*yh*varphi_h/2*(pai-1)^(2);
dd = 0;
%%dd = 1/gammt*(yh-rmc*yh);
d = gammt*dd;
cm=phim*hms*wms+(1-phim)*hmu*wmu;
yh1= wds*hds*phi_t*gammt + wdu*hdu*(1-phi_t)*gammt + phim*hms*wms*(1-gammt)+ (1-phim)*hmu*wmu*(1-gammt) + kd*rk*gammt + (1/gd)*gammt*delt*k+gammt*dd;
cd= wds*hds*phi_t + wdu*hdu*(1-phi_t)+kd*rk-ind + dd +(1-Q)*bd/pai;
c =  cd*gammt + cm*(1-gammt);
ad = (cd*(phi_t*hds)^(1/sigma))/((1-hd)*hd^(1/sigma)*wds);
am = (cm*(phim*hms)^(1/sigma))/((1-hm)*hm^(1/sigma)*wms);
lambd =  ad/cd;
lambd1 = am/cm;
GDP = yh;
DOMESTIC_DEBT = b/GDP;
cd;
hd;
hds;
hdu;
k;
b;
hd;
ind;
dd;
cm;
hm;
hms;
hmu;
c;
kd;
he;
hs;
hu;
yh;
in;
d;
phim;
phi_t;
%gd;
%gds;
%gdu;
%gm;
%gms;
%gmu;
gT;
gammt;
zetabar;
eps;
Q;
rk;
wds;
wdu;
wms;
wmu;
paih;
pai;
rmc;
rr;
zeta;
am;
yh1;
lambd;
lambd1;


NumberOfEndogenousVariables = M_.endo_nbr;                    % Number of endogenous variables.
ys = zeros(NumberOfEndogenousVariables,1);                    % Initialization of ys (steady state).
for ii = 1:NumberOfEndogenousVariables                         % Loop...
  varname = deblank(M_.endo_names(ii,:));                      %    Get the name of endogenous variable i.                     
  eval(['ys(' int2str(ii) ') = ' varname ';']);                %    Get the steady state value of this variable.
end        

id1 = strmatch('bet',M_.param_names,'exact');          
M_.params(id1) = bet;
id2 = strmatch('delt',M_.param_names,'exact');          
M_.params(id2) = delt;
%id3 = strmatch('eta0',M_.param_names,'exact');          
%M_.params(id3) = eta0;
id4 = strmatch('rho_ad',M_.param_names,'exact');          
M_.params(id4) = rho_ad;
id5 = strmatch('sigma',M_.param_names,'exact');          
M_.params(id5) = sigma;
id6 = strmatch('sigma1',M_.param_names,'exact');          
M_.params(id6) = sigma1;
id7 = strmatch('rho_am',M_.param_names,'exact');          
M_.params(id7) = rho_am;
id8 = strmatch('omegau',M_.param_names,'exact');          
M_.params(id8) = omegau;
id9 = strmatch('omegas',M_.param_names,'exact');          
M_.params(id9) = omegas;
id10 = strmatch('rho',M_.param_names,'exact');          
M_.params(id10) = rho;
id11 = strmatch('rho1',M_.param_names,'exact');          
M_.params(id11) = rho1;
id12 = strmatch('rho2',M_.param_names,'exact');          
M_.params(id12) = rho2;
id13 = strmatch('alf',M_.param_names,'exact');          
M_.params(id13) = alf;
id14 = strmatch('rho_zeta',M_.param_names,'exact');          
M_.params(id14) = rho_zeta;
%id15 = strmatch('rho_ystar',M_.param_names,'exact');          
%M_.params(id15) = rho_ystar;
id16 = strmatch('varphi_h',M_.param_names,'exact');          
M_.params(id16) = varphi_h;
id17 = strmatch('epsi',M_.param_names,'exact');          
M_.params(id17) = epsi;
%id18 = strmatch('Gam',M_.param_names,'exact');          
%M_.params(id18) = Gam;
%id19 = strmatch('nuh',M_.param_names,'exact');          
%M_.params(id19) = nuh;
id20 = strmatch('rho_r',M_.param_names,'exact');          
M_.params(id20) = rho_r;
id21 = strmatch('phi_pai',M_.param_names,'exact');          
M_.params(id21) = phi_pai;
id22 = strmatch('phi_y',M_.param_names,'exact');          
M_.params(id22) = phi_y;
id21 = strmatch('lambd',M_.param_names,'exact');          
M_.params(id21) = lambd;
id22 = strmatch('lambd1',M_.param_names,'exact');          
M_.params(id22) = lambd1;
id23 = strmatch('choque',M_.param_names,'exact');          
M_.params(id23) = choque;
id24 = strmatch('phi_t',M_.param_names,'exact');          
M_.params(id24) = phi_t;
id25 = strmatch('gd',M_.param_names,'exact');          
M_.params(id25) = gd;
id26 = strmatch('gammt',M_.param_names,'exact');          
M_.params(id26) = gammt;
id27 = strmatch('phim',M_.param_names,'exact');          
M_.params(id27) = phim;
id28 = strmatch('gT',M_.param_names,'exact');          
M_.params(id28) = gT;
id27 = strmatch('am',M_.param_names,'exact');          
M_.params(id27) = am;
id28 = strmatch('ad',M_.param_names,'exact');          
M_.params(id28) = ad;
phi_t, gd
%id23 = strmatch('bias',M_.param_names,'exact');          
%M_.params(id22) = bias;
%id23 = strmatch('thetd',M_.param_names,'exact');          
%M_.params(id22) = thetd;
%id24 = strmatch('thetm',M_.param_names,'exact');          
%M_.params(id22) = thetm;
%id25 = strmatch('thet3',M_.param_names,'exact');          
%M_.params(id22) = thet3;
%id26 = strmatch('thet4',M_.param_names,'exact');          
%M_.params(id22) = thet4;



