clc; clear all; close all


%% Structural Parameters (quarterly calibration)
beta=0.99;                    % discount factor
alpha=0.33;                   % elasticity of production wrt capital
epsilon=6;                    % elasticity of substitution btw differentiated goods
delta=0.025;                  % depreciation rate
sigma=2;                      % relative risk aversion
phi=1;                        % inverse of Frisch elasticity
g=0.2;                        % share of public spending in ss
eta=1.5;                      % elasticity of intratemporal substitution
gamma=0.3;                    % trade openness
D=0.25;                       % external debt/GDP ratio
yz=1;
tot=1;
%% Chang params
parameter=calibration1;
DELTA=parameter(1);
G=1;
OMEGAK=parameter(3);
PSI=parameter(4);
PSIY=parameter(5);
ETA=parameter(6);
BETA=parameter(7);
EPS=parameter(8);
OMEGAP=parameter(9);
PISS=parameter(10);
ALPHA=parameter(11);
THETA=parameter(12);
XIS=parameter(13);
XIP=parameter(14);
K=parameter(15);
OMEGAM=parameter(16);
AS=parameter(17);
AP=parameter(18);
LS=parameter(19);
LP=parameter(20);
TAU=parameter(21);
GSS=parameter(22);
SIGMAM=parameter(23);
MS=parameter(24);
MP=parameter(25);
SIGS=parameter(26);
SIGG=parameter(27);
SIGAP=parameter(28);
SIGAS=parameter(29);
SIGDA=parameter(30);

GDPSS = gdpsss(DELTA, G, OMEGAK, PSI, ETA, BETA, EPS, OMEGAP, PISS, ALPHA, THETA, MS,MP, XIS,XIP, K, OMEGAM, AS,AP, LS, LP, TAU,GSS,SIGMAM,PSIY);
CSS   = csss(DELTA, G, OMEGAK, PSI, ETA, BETA, EPS, OMEGAP, PISS, ALPHA, THETA, MS,MP, XIS,XIP, K, OMEGAM, AS,AP, LS, LP, TAU,GSS,SIGMAM,PSIY);
DSS   = dsss(DELTA, G, OMEGAK, PSI, ETA, BETA, EPS, OMEGAP, PISS, ALPHA, THETA, MS,MP, XIS,XIP, K, OMEGAM, AS,AP, LS, LP, TAU,GSS,SIGMAM,PSIY);

PSIM  = (1 - BETA/PISS/G)/(CSS)*(DSS/0.95*0.05);%currency stands for 5% of M2

RHOS = 0.95;
RHOG = 0.95;
RHOAP = 0.95;
RHOAS = 0.95;
RHODA = 0.95;

RR = 0;
RP = -0.397/(1-0.391);
RY = 0.183/(1-0.391);
TP = 0; 
TY = 0; 
TT = 0;
%% Chang part steady state
hes = 1;
hep = 1;
I_K = 1-(1-DELTA)/G;% Eq. 4
rk = G/BETA-(1-DELTA); % Eq. 9; Eq. 8 gives uc = uk

r = G*PISS/BETA; % Eq. 7
rr=1/beta; 
x = EPS/(EPS-1); % Eq. 15
rs = (r-1)/(1-TAU)+1; % Eq. 31
rp = r; % Eq. 34
pi = PISS; % steady state inflation 
tau = TAU; % steady state RRR
q=1;

% solve for bankruptcy cutoffs omegas omegap 
% using Eq. 16-28.
options = optimset('MaxFunEvals',100000,'MaxIter',10000,'TolFun',1e-16);
FF =fsolve(@(x)solveff(AS,AP,K,OMEGAM,LS,LP,MS,MP,XIS,XIP,ALPHA,THETA,G,PISS,BETA,TAU,OMEGAM*(1+exp(x(1))),OMEGAM*(1+exp(x(2)))),[log(1.02),log(1.01)],options);

omegas = OMEGAM*(1+exp(FF(1)));
omegap = OMEGAM*(1+exp(FF(2)));

FS = 1/(K-1)*OMEGAM^K*omegas^(-K+1); % fraction of revenue income going to the SOE borrower f(omegas)
GS = LS*omegas+(1-LS)*(1-(1-MS)*K/(K-1))*OMEGAM^K*omegas^(-K+1)+(1-MS)*(1-LS)*K/(K-1)*OMEGAM; % fraction of revenue income going to the SOE lender g(omegas)
G1S = LS+(1-LS)*(1-MS*K)*OMEGAM^K*omegas^(-K); % first order derivative of g(omegas)
F1S = -OMEGAM^K*omegas^(-K);% first order derivative of f(omegas)
BS_YS = GS/rs; % Eq. 26
NS_YS = -G1S/F1S*FS/rs*PISS*G; % Eq. 27
KS_YS = (1-ALPHA)*(NS_YS/PISS/G+BS_YS)/rk; % Eq. 22

FP = 1/(K-1)*OMEGAM^K*omegap^(-K+1); % fraction of revenue income going to the POE borrower f(omegap)
GP = LP*omegap+(1-LP)*(1-(1-MP)*K/(K-1))*OMEGAM^K*omegap^(-K+1)+(1-MP)*(1-LP)*K/(K-1)*OMEGAM;% fraction of revenue income going to the SOE lender g(omegap)
G1P = LP+(1-LP)*(1-MP*K)*OMEGAM^K*omegap^(-K);% first order derivative of g(omegap)
F1P = -OMEGAM^K*omegap^(-K);% first order derivative of g(omegap)
BP_YP = GP/rp; % Eq. 26
NP_YP = -G1P/F1P*FP/rp*PISS*G;% Eq. 27
KP_YP = (1-ALPHA)*(NP_YP/PISS/G+BP_YP)/rk; % Eq. 22

HSP1 =(NS_YS/PISS/G+BS_YS)/(NP_YP/PISS/G+BP_YP);
YSP = ((PSIY/(1-PSIY))^(SIGMAM/(SIGMAM-1))*AS/AP*(KS_YS/KP_YP)^(1-ALPHA)*HSP1^(ALPHA*THETA))^(1/(ALPHA-ALPHA*THETA+1/(SIGMAM-1)));%YSPP = YS/YP*PS/PP

%WSP = ((AS/AP)*(MU/(1-MU))^((1-THETA)*ALPHA/XIL)*(NS_YS/PISS/G+BS_YS)/(NP_YP/PISS/G+BP_YP))^(1/(ALPHA+(1-THETA)*ALPHA/XIL));
%YSP = WSP^((1-ALPHA+ALPHA*THETA)/XIL+1-ALPHA)/(MU/(1-MU))^((1-ALPHA+ALPHA*THETA)/XIL)*AS/AP;

YP_Y = 1/(1+YSP)/x;% YP_Y = YP*PP/(Y*P) = YP*PP/(YS*PS+YP*PP)*Pw/P
YS_Y = YSP/(1+YSP)/x;% YP_Y = YP*PP/(Y*P) = YP*PP/(YS*PS+YP*PP)*Pw/P
ps = (YS_Y/PSIY^SIGMAM/x^(-SIGMAM))^(1/(1-SIGMAM));
pp = (YP_Y/(1-PSIY)^SIGMAM/x^(-SIGMAM))^(1/(1-SIGMAM));

K_Y = (KS_YS*YS_Y+KP_YP*YP_Y)*G; % Eq. 46 
GDP_Y = 1-YS_Y*MS*K/(K-1)*(OMEGAM - OMEGAM^K*omegas^(-K+1))-YP_Y*MP*K/(K-1)*(OMEGAM - OMEGAM^K*omegap^(-K+1)); % Eq. 37
G_Y = GSS*GDP_Y; % Eq. 39
C_Y = GDP_Y-I_K*K_Y-G_Y; % Eq. 38

LIQ=YS_Y*MS*K/(K-1)*(OMEGAM - OMEGAM^K*omegas^(-K+1))+YP_Y*MP*K/(K-1)*(OMEGAM - OMEGAM^K*omegap^(-K+1));

WAGE_Y = ALPHA*THETA*((NP_YP/PISS/G+BP_YP)*YP_Y+(NS_YS/PISS/G+BS_YS)*YS_Y); % Eq. 20
h = (WAGE_Y/C_Y/PSI)^(1/(ETA+1)); % Eq. 5,6
HS_HP = (NS_YS/PISS/G+BS_YS)/(NP_YP/PISS/G+BP_YP)*(YSP);
hhp = h/(1+(HS_HP));
hhs = h/(1+1/(HS_HP)); % Eq. 20 ???

ys = ((AS)*KS_YS^(1-ALPHA)*(hes^(1-THETA)*hhs^THETA)^ALPHA*ps^(1-ALPHA))^(1/ALPHA); % Eq. 16
yp = ((AP)*KP_YP^(1-ALPHA)*(hep^(1-THETA)*hhp^THETA)^ALPHA*pp^(1-ALPHA))^(1/ALPHA);% Eq. 16


y = ys/PSIY^SIGMAM/(ps*x)^(-SIGMAM);% Eq. 18
b=4*D;
s=1;
pH=1; % domestic price
gdp=GDP_Y *y;
k = K_Y*y;
i = I_K*k;
g = G_Y*y;
iz=gamma*s^(-eta)*i;
ih=(1-gamma)*pH^(-eta)*i;
ch=y - LIQ*y - ih - g- b*(1/BETA-1);
c=ch*((1-gamma)*pH^(-eta))^(-1);
cz=gamma*s^(-eta)*c;
uc = 1/c; % Eq. 5
w =PSI*h^(ETA)*c;
gammaz=1-(1-gamma)*(c+i)-g;   % foreign parameter
tb=pH*gammaz*yz*(tot)^eta-s*(cz+iz);
gdp=c+i+g +tb;
uk=uc;
WE_YP = ALPHA*(1-THETA)*(NP_YP/PISS/G+BP_YP);
WE_YS = ALPHA*(1-THETA)*(NS_YS/PISS/G+BS_YS);

wes =  WE_YS*ys*ps;
wep =  WE_YP*yp*pp;

ks = KS_YS*ys*ps;
ns = NS_YS*ys*ps;
bs = BS_YS*ys*ps;
as = ys*ps/(ns/PISS/G+bs);

kp = KP_YP*yp*pp;
np = NP_YP*yp*pp;
bp = BP_YP*yp*pp;
ap = yp*pp/(np/PISS/G+bp);

%ytfp = y/((ks+kp)^(1-ALPHA)*((h)^(THETA))^ALPHA);


%gtfp = gdp/((ks+kp)^(1-ALPHA)*(h^(THETA))^ALPHA);% Eq. 43
%utility = (log(c)-PSI*h^(1+ETA)/(1+ETA));
rzs  = omegas*as*(ns/PISS/G+bs)/bs;% Eq. 23
rzp  = omegap*ap*(np/PISS/G+bp)/bp;% Eq. 23
bcp = 1-(OMEGAM/omegap)^K;% Transform bankruptcy cutoff to bankruptcy ratio
bcs = 1-(OMEGAM/omegas)^K;
x=log(x);
ih=log(ih);
iz=log(iz);
ch=log(ch);
cz=log(cz);
s=log(s);
k=log(k);
i=log(i);
uc=log(uc);
c=log(c);
h=log(h);
w=log(w);
r=log(r);
pi=log(pi);
uk=log(uk);
rk=log(rk);
rr=log(rr);
pH=log(pH);
y=log(y);
wes=log(wes);
wep=log(wep);
ys=log(ys);
ks=log(ks);
hes=log(hes);
hhs=log(hhs);
ns=log(ns);
bs=log(bs);
as=log(as);
rs=log(rs);
omegas=log(omegas);
yp=log(yp);
kp=log(kp);
hep=log(hep);
hhp=log(hhp);
np=log(np);
bp=log(bp);
ap=log(ap);
rp=log(rp);
omegap=log(omegap);
gdp=log(gdp);
%ytfp=log(ytfp);
%gtfp=log(gtfp);
rzs=log(rzs);
rzp=log(rzp);
bcs=log(bcs);
bcp=log(bcp);
g=log(g);
pp = log(pp);
ps = log(ps);
tb= log(tb);

b=log(b);
yz=log(yz);

save par beta alpha delta sigma phi eta gamma gammaz k i uc c h w r pi uk rk y wes wep    x rr...
ys ks hes hhs ns bs as rs omegas ...
yp kp hep hhp np bp ap rp omegap ...
 tau gdp  rzs rzp  g bcs bcp ps pp tb s b yz  % save parameters to use them in Dynare

