clear all;
clc;

global sigma1 beta xi theta gamma sstar; 

% PARAMETERS

% Banks

theta = 0.475;      % divertable proportion of assets
gamma = 6.4;        % home bias in funding
sigma1 = 0.94;       % survival probability
xi = 0.000588;      % fraction of total assets brought by new banks

% Households 


beta = 0.985;       % discount rate
zeta = 0.2;         % inverse of Frisch elasticity of labor supply
zeta0 = 5.89;       % inverse of labor supply capacity
varkappa = 0.000985;    % cost parameter of direct finance

% Producers

alphaK = 0.3;       % cost share of capital
alphaM = 0.15;      % cost share of imported intermediate goods
lambda = 0.98;      % one minus depreciation rate
eta = 9;            % elasticity of demand
omega = 0.66;       % fraction of non-adjusters
kappa = (eta-1)*omega/((1-omega)*(1-beta*omega));       % slope of Phillips curve
kappaI = 1;         % cost of adjusting investment goods production
varphi = 1;         % price elasticity of export demand

Q = 1; %1
R = 1/beta; % 2 
Rstar = 1.005;
A = 1;
eps = 1;
Lambda = beta;
mC = (eta-1)/eta;
pi = 1;
i = 1/beta-1;

% find s
sstar = 1-beta*Rstar; %3


global s;

myfun = @(s,sigma1,beta,xi,theta) (1-sigma1)*(beta*s+xi*(1+s))*(sigma1*s+xi*(1+s))-theta*(beta-sigma1)*(sigma1*(1-beta)*s+(1-sigma1)*xi*(1+s));  % parameterized function
%c = 2;                    % parameter
fun = @(s) myfun(s,sigma1,beta,xi,theta);    % function of x alone
x0 = [0.0001 1];
s = fzero(fun,x0)
%fprintf('\nx = %10.6f\n',s)

%All below is OK^

%Here is missed code START

%options = optimset('Display','iter');
%s = fzero('H',[0.000001 1],options);   

%function f = fun(n)
 %   f = (1-sigma1)*(beta*s+xi*(1+s))*(sigma1*s+xi(1+s))-theta*(beta-sigma1)*(sigma1*(1-beta)*s+(1-sigma1)*xi*(1+s));
%end

%fun = (1-sigma1)*(beta*s+xi*(1+s))*(sigma1*s+xi(1+s))-theta*(beta-sigma1)*(sigma1*(1-beta)*s+(1-sigma1)*xi*(1+s));
%x0 = [0.000001 1];
%swork = fzero('fun',x0,options);
%s = swork;

%myfun = @(s,sigma1,beta,xi,theta) (1-sigma1)*(beta*s+xi*(1+s))*(sigma1*s+xi(1+s))-theta*(beta-sigma1)*(sigma1*(1-beta)*s+(1-sigma1)*xi*(1+s));  % parameterized function
%c = 2;                    % parameter
%fun1 = @(s) myfun(s,sigma1,beta,xi,theta);    % function of x alone
%swork = fzero(myfun,[0.000001 1])

%s = swork

%Here is missed code FINISH

x = s/sstar*(-1+sqrt(1+2/gamma*(sstar/s)^2)); %4 ???
Theta = theta*(1+gamma/2*x^2);  %5

Z = (s+1)/beta-lambda; %6
Kh = s/varkappa; %7

K2Y = (1-1/eta)*alphaK/Z; %8

Y = ((1-1/eta)^(1+zeta*(alphaK+alphaM))*(A/(eps^alphaM*Z^alphaK))^(1+zeta))^(1/(zeta*(1-alphaK-alphaM)))/((1-alphaK-alphaM)*zeta0^(1/zeta)); %9

K = K2Y*Y;

I = (1-lambda)*K; %10

Kb = K-Kh; %15

mustar = sstar/s;

phi = (beta-sigma1)/(sigma1*(s+sstar*x)+xi*(1+s)); %13

psi = ((1-sigma1)*((s+sstar*x)*phi+1))/(1-sigma1*((s+sstar*x)*phi+1)); %14

Omega = Lambda*(1-sigma1+sigma1*psi); %16

mu = Omega*(Z+lambda-R); %17 ???

mudstar = Omega*(R-Rstar);  %21

nu = Omega*R;

Ystar2Y = (alphaM*(1-1/eta)+(Rstar-1)*x*Kb/Y)/(eps^varphi); %12 ???

Ystar = Ystar2Y*Y;

Ex = eps^varphi*Ystar; %18

Chi = varkappa/2*Kh^2;

C2Y = 1-(1-lambda)*K2Y-eps^varphi*Ystar2Y-Chi/Y;  %11 ???

C = C2Y*Y;

M = alphaM/eps*(1-1/eta)*Y;

Dstar = x*Kb/eps;

N = Kb/phi; %19

D = Kb-N-Dstar;

L = (1-alphaK-alphaM)*(Y/A*(alphaK/K)^alphaK*(alphaM/M)^alphaM)^(1/(1-alphaK-alphaM));

w = (1-alphaK-alphaM)/alphaK*Z*K/L; %20

fprintf('\nSteady state values:\n');

fprintf('\nPrices:\n');

fprintf('\nmC = %10.6f\n',mC);
fprintf('\npi = %10.6f\n',pi);
fprintf('\nZ = %10.6f\n',Z);
fprintf('\nw = %10.6f\n',w);
fprintf('\ni = %10.6f\n',i);
fprintf('\neps = %10.6f\n',eps);
fprintf('\nQ = %10.6f\n',Q);

fprintf('\nQuantities:\n');

fprintf('\nY = %10.6f\n',Y);
fprintf('\nM = %10.6f\n',M);
fprintf('\nL = %10.6f\n',L);
fprintf('\nC = %10.6f\n',C);
fprintf('\nI = %10.6f\n',I);
fprintf('\nK = %10.6f\n',K);
fprintf('\nEx = %10.6f\n',Ex);
fprintf('\nN = %10.6f\n',N);
fprintf('\nKb = %10.6f\n',Kb);
fprintf('\nKh = %10.6f\n',Kh);
fprintf('\nD = %10.6f\n',D);
fprintf('\nDstar = %10.6f\n',Dstar);

fprintf('\nBank variables:\n');

fprintf('\nx = %10.6f\n',x);
fprintf('\npsi = %10.6f\n',psi);
fprintf('\nphi = %10.6f\n',phi);
fprintf('\nnu = %10.6f\n',nu);
fprintf('\nmu = %10.6f\n',mu);
fprintf('\nmudstar = %10.6f\n',mudstar);

fprintf('\nBank variables:\n');

fprintf('\nA = %10.6f\n',A);
fprintf('\nYstar = %10.6f\n',Ystar);
fprintf('\nRstar = %10.6f\n',Rstar);
fprintf('\nR = %10.6f\n',R);
fprintf('\nOmega = %10.6f\n',Omega);
fprintf('\nLambda = %10.6f\n',Lambda);
fprintf('\nTheta = %10.6f\n',Theta);
fprintf('\nmustar = %10.6f\n',mustar);
fprintf('\nChi = %10.6f\n',Chi);