% Dekle-Jeong-Kiyotaki Small Open Economy Steady State Calculation for new Shrunk Model (taking real exchange rate as exogenous and dropping foreign bond equations)
% July 12, 2012
% coded by Karrar Hussain


function [ss,check] = djk_shrunk_steadystate(ss,exe)
 
%Parameters of the model (Exogenous)

beta    =0.98;
%phi_k   =10;
mu      =0.75;
k       =0.999999;
alpha   =0.36;
delta   =0.025;
theta   =8;
eta     =1.46;
gamma   =0.00000895;
psi     =0.25;
g       =0.13;
tau_c   =0.000000001;
tau_y   =0.000000001;
a       =1;


% (1) Constants already defined in the paper 

e       =1/(1+tau_c);
f       =1/(1-tau_y);
M       =(theta/(theta-1));
A       =((1-beta)/(beta*f))+delta;
B       =((alpha*a)/A)^(1/(1-alpha));
D       =a*(1-alpha)*(B^alpha);
E       =D-A*B;

N1      =(-beta*tau_c*e)/((1-beta)*eta*f);
N2      =(k-beta)/((1-beta)*(k-1));
N3      =(-beta*tau_y)/(1-beta);
N4      =-beta/(1+beta);
N5      =a*(B^alpha);
N6      =(e*D)/(eta*f);
N7      =delta*B;
N8      =N1*D;
N9      =N1*D+N3*D+N3*B*(A-delta)-N3*E;
N10     =N5/(1-mu);
N11     =N6-N7;
N12     =(1/(1-mu))^(1/(1-theta));
N13     =N12/M;
%N10     =(beta*tau_y)/(1+beta);
  check = 0;
  
%%% (6) Solving for steady state of N, after substituting the steady state values of C and
%%% (wx/N) in equation (32), from part (4) of the code file

x0 = 1;
options=optimset('Display','iter','TolX',1e-20,'TolFun',1e-39) ;
options.MaxFunEvals=100000;
options.MaxIter = 100000;
JN = @(x)paramJN(x,N9,N8,N5,N4,N6,N12,N11,g,N3,N2,psi,beta,mu,theta,N13,alpha,gamma,D,eta,f);
[x,fval,exitflag] = fsolve(JN,x0,options); 

pi = x;

% functions of N in sequential order derived and mentioned above in the
% part (4) of the code file

p_star  = (N12*((1-mu*(x^(theta-1)))^(1/(1-theta))));

phi     = ((N13*((1-mu*(x^(theta-1)))^(1/(1-theta)))*(1-mu*beta*(x^theta)))/(1-mu*beta*(x^(theta-1))));

L       =((1-mu)*(((N12*((1-mu*(x^(theta-1)))^(1/(1-theta)))))^-theta)/(1-mu*(x^theta)));

h       = ((g+N6*((((N13*((1-mu*(x^(theta-1)))^(1/(1-theta)))*(1-mu*beta*(x^theta)))/(1-mu*beta*(x^(theta-1)))))^(1/(1-alpha))))/(((((N13*((1-mu*(x^(theta-1)))^(1/(1-theta)))*(1-mu*beta*(x^theta)))/(1-mu*beta*(x^(theta-1)))))^(alpha/(1-alpha)))*(N5/(((1-mu)*(((N12*((1-mu*(x^(theta-1)))^(1/(1-theta)))))^-theta)/(1-mu*(x^theta)))))+N11*((((N13*((1-mu*(x^(theta-1)))^(1/(1-theta)))*(1-mu*beta*(x^theta)))/(1-mu*beta*(x^(theta-1)))))^(1/(1-alpha)))));

y       = (N5*((((N13*((1-mu*(x^(theta-1)))^(1/(1-theta)))*(1-mu*beta*(x^theta)))/(1-mu*beta*(x^(theta-1)))))^(alpha/(1-alpha)))/(((1-mu)*(((N12*((1-mu*(x^(theta-1)))^(1/(1-theta)))))^-theta)/(1-mu*(x^theta)))))*((g+N6*((((N13*((1-mu*(x^(theta-1)))^(1/(1-theta)))*(1-mu*beta*(x^theta)))/(1-mu*beta*(x^(theta-1)))))^(1/(1-alpha))))/(((((N13*((1-mu*(x^(theta-1)))^(1/(1-theta)))*(1-mu*beta*(x^theta)))/(1-mu*beta*(x^(theta-1)))))^(alpha/(1-alpha)))*(N5/(((1-mu)*(((N12*((1-mu*(x^(theta-1)))^(1/(1-theta)))))^-theta)/(1-mu*(x^theta)))))+N11*((((N13*((1-mu*(x^(theta-1)))^(1/(1-theta)))*(1-mu*beta*(x^theta)))/(1-mu*beta*(x^(theta-1)))))^(1/(1-alpha)))));

d       = y-E*h*(phi^(1/(1-alpha)));

w       =D*(phi^(1/(1-alpha)));

K       =B*(phi^(1/(1-alpha)))*h;

L1      =(a*(K^alpha)*(h^(1-alpha)))/y;

i       =pi/beta;

r       =A;

m       = ( (w*(1-h)/eta*f)*(pi/(pi-beta))*(gamma))^psi;

b       = (beta/(1-beta))*(((tau_c*w*(1-h))/(eta*f))+tau_y*(w*h+(A-delta)*K+d)-g+m*(1-(pi^-1)));

X       =(phi*y*eta*f)/(w*(1-h)*(1-mu*beta*(pi^theta)));

Z       =(y*eta*f)/(w*(1-h)*(1-mu*beta*(pi^(theta-1))));

p_star1 =M*(X/Z);

x1     =delta*K;

c      =((pi-beta)/pi)*(e/gamma)*(m^(1/psi));

c1     =(((1-h)*w*e)/(eta*e*f));

c2     =y-x1-g;

Y1     =w*h+r*K+d;

Y2     =c+delta*K+g;

  ss =[a,tau_c,g,tau_y,c,m,h,b,K,pi,y,x1,d,L,X,Z,p_star,r,w,i,phi];
  
  Equation1=(e/c)-(gamma*m^(-1/psi)+(beta*e/(c*pi)));
  Equation2=(e/c)-(eta*f/(w*(1-h)));
  Equation3=1-(beta*i/pi);
  Equation4=delta*K-x1;
  Equation5=r-A;
  Equation6=L1*y-a*(K^alpha)*(h^(1-alpha));
  Equation7=A-a*alpha*phi*(K^(alpha-1))*(h^(1-alpha));
  Equation8=w-a*(1-alpha)*phi*(K^alpha)*(h^(-alpha));
  Equation9=y-c-delta*K-g;
  Equation10=y-w*h-(r)*K-d+0.6613;
  Equation11=L-((1-mu)*(p_star^-theta))/(1-mu*x^theta);
  Equation12=(p_star^(1-theta))-((1-mu*x^(theta-1))/(1-mu));
  Equation13=X-(e*phi*y/(c*(1-mu*beta*(x^theta))));
  Equation14=Z-(e*y/(c*(1-mu*beta*(x^(theta-1)))));
  Equation15=tau_c*c+tau_y*(w*h+(r-delta)*K+d)-g-(i*b/x)+m*(1-x^-1)+b;
  Equation16=((k-1)*(tau_c*c+tau_y*(w*h+(r-delta)*K+d)-g-(i*b/x)))-m*(1-x^-1);
  Equation17=p_star-M*(X/Z);