function [varargout] = ss3(pi,firstcall)

%initiate params
global params Zhat wedge

pi
rho    = 0.112;  %  rate of time preference of entrepreneurs
rho_hh = 0.025;    %  rate of time preference of HH
delta  = 0.065;    %  Capital depreciation rate
gamma  = 1.4;     %  borrowing constraint parameter
alpha  = 0.3;     %  capital share in production
thetaP = 100;      % price adjustment cost
epsilon= 10;        % markup
phiTR  = 1.25;     % Taylor rule coefficient 
rhoTR  = 0.8 ;     % Quarterly DT = 0.8, CT (1-0.8)/(1/4)=0.8
pibar  = 0;         % SS inflation
eta_q   =.2;        % TFP shock persistence, yearly DT = 0.8, CT (1-0.8)/1=0.2
Deltat  = 1/4;       % Time period
share_d  = rho;      % 0 is model where entrepreneurs consume but have 0 welfare weight. rho is model where they pay dividends.
alphaZ   = alpha;   %exponent on exogenous and endogenous productivity in production function.

psi      = 1;%0.00001;   %  Inverse Frisch elasticity
phiK     = 10;%10;%0;     % capital adj costs

timeless = 1;      % timeless or time 0

N_period  = 300/Deltat; % Number of periods to have 200 years


Upsilon=0.801719796182658;
Zhat     = 1.016451505882546;
wedge    = 0.585586117016045;

%  Zhat     = 1;
% wedge    = 1;
%Initialize
K        = 0;
w        = 0;
Y        = 0;
T        = 0;
C        = 0;
I        = 0;


 

 m        = (epsilon-1)/epsilon + rho_hh*pi*thetaP/epsilon ; % CT 
 mbar     = (m/((epsilon-1)/epsilon)); % = m if no subsidy; =(m/((epsilon-1)/epsilon)) if subsidy
 R        = (rho_hh+delta);
 r        = (R-delta);
 i        = r + pi;
 iota     = delta;
 q        = 1;
 rho_hh_t = rho_hh;
 Q        = 1;
 Welf     = 0;
 Aux_Va   = 0;


L        = 1;
SS_given_L;


if firstcall==0
    opt2 = optimoptions('fsolve','FunctionTolerance',1e-14,'OptimalityTolerance',1e-14,'Display','off');
    sol = fsolve(@residual,1,opt2);   
    residual(sol)
end


function [out] =SS_given_L
    %K        = ((r+delta)/(alpha*(m/((epsilon-1)/epsilon))*L^(1-alpha)))^(1/(alpha-1));
    K        = (R/(wedge*alpha*mbar*Zhat*L^(1-alpha)))^(1/(alpha-1));
    %w        = (1-alpha)*K^alpha*L^(-alpha)*(m/((epsilon-1)/epsilon));
    w        = (1-alpha)*Zhat*K^alpha*L^(-alpha)*(mbar);
    Y        = Zhat*K^alpha*L^(1-alpha);
    T        = (1-mbar)*Y-thetaP/2*pi^2*Y+(R/wedge-R)*K;
    %T        = (1-m)*Y-thetaP/2*pi^2*Y;
    
    C         = ((R-delta*q)*K+w*L+T);
    %C        = Y-(q*iota*K+thetaP*1/2*pi^2*Y);
    %C        = Y-(q*iota*K);

    I        = delta*K;  
end

function out = residual(x_in) % to find L; calls the SS_give_L function
        L        = x_in;
        SS_given_L;
        out      = L - (Upsilon/(w)*C)^(1/-psi);
end

varargout{1}=L;
varargout{2}=r;
varargout{3}=m;
varargout{4}=K;
varargout{5}=w;
varargout{6}=i;
varargout{7}=Y;
varargout{8}=T;
varargout{9}=C;
varargout{10}=rho_hh_t;
varargout{11}=Q;
varargout{12}=Welf;
varargout{13}=I;
varargout{14}=q;
varargout{15}=R;
varargout{16}=iota;
varargout{17}=Aux_Va;
varargout{18}=mbar;

varargout{end+1}=w/C;;


if firstcall==1
    varargout{1}=cell2mat(varargout);
end


        





end








