clear all
close all
clc

% set parameters
beta = 0.993;
kappa = 1-40^(-1);
spF = 0.03;
psi = 0.6;
epsp = 11;
epsw = 11;
alpha = 1/3;
delta0 = 0.1;
delta2 = 0.2; % squared term capital adjustment cost
g = 0.2;
bb = 0.7; 
sigma = 0.9;

upse = 0.96;
upsc = 0.96;
v = 0.96;
Delta = 0.5;
etap = 1.5;
m=0.15;
deltaT = 0.04;

% parameters unrelated to steady state
phip = 0;
phiw = 0;
psik = 3; % capital adjustment cost
eta = 1; % Frisch elasticity

% Taylor rule
rhor = 0.7;
phipi = 1.5;
phiy = 0.25;
sr = 0.025;

% Shock processes
rhof = 0.97;
sf = 0.01;
rhoA = 0.9;
sA = 0.01;
rhot = 0.95;
st = 0.01;
indp = 0.3;
indw = 0.3;



% solve for steady state
Rds = beta^(-1);
Rfs = (1+spF)^(1/4)*Rds;
Rvs = Rds + Delta*(Rfs-Rds);
Rtrs = Rds;
Qs = (Rfs - kappa)^(-1);
M1s = beta/(Qs*(1-beta*kappa));
M2s = 1 + (M1s-1)*psi;
pws = (epsp-1)/epsp;
Ks = (alpha*pws/(M2s*(1/beta - (1-delta0))))^(1/(1-alpha));
Ys = Ks^(alpha);
Is = delta0*Ks;
Gs = g*Ys;
ws = (1-alpha)*pws*Ks^(alpha);
delta1 = (alpha*pws*Ks^(alpha-1)/M2s);

fcbs = 0; 
fws = psi*Is/(Qs*(1-kappa));
mrss = ((epsw-1)/epsw)*ws;
fs = fws - fcbs;
Pes=1;

%syms nes phies lambdaes W1 thetaes Omegaes rers Pcs Vs ncs dcs phics lambdacs W2 thetacs Omegacs Ccs Ces deltaT
%eqns1 = [nes == (sigma*(Rfs-Rvs)+deltaT)*Qs*fs/(1-sigma*Rvs),...
%phies== Qs*fs/nes, lambdaes == (Rfs-Rvs)*phies/Rvs,W1 == (1+lambdaes)*beta*(Rfs-Rvs),...
%thetaes == (1-sigma)*W1/(lambdaes-sigma*phies*W1),Omegaes == 1-sigma+sigma*thetaes*phies,...
%rers == ((2*v-1)/(1-v)*(v/(2*v-1)-Pes^(1-etap)))^(1/(1-etap)),...
%Pcs == (v/(2*v-1)-((1-v)/(2*v-1))*rers^(etap-1))^(1/(1-etap)),...
%Vs == Qs*fs-nes,ncs == (sigma*((Rfs-Rds)*Qs*fs+(m/(1-m))*(Rvs-Rds)*Vs)+deltaT*Qs*fs)/(1-Rds*sigma),...
%dcs == (m/(1-m))*Vs+Qs*fs-ncs, phics == (Qs*fs+Delta*Vs*(m/1-m))/ncs,lambdacs == (phics*(Rfs-Rds))/Rds,...
%W2 == (1+lambdacs)*(Rfs-Rds)*beta, thetacs == (1-sigma)/(lambdacs/W2-sigma*phics),...
%Omegacs == 1-sigma+sigma*phics*thetacs,mrss+(Rvs-1)*dcs-deltaT*Qs*fs==Ccs,...
%Ces==(Ys-(1-v)*(Ccs+Is+Gs))/v-Is-Gs,Ccs==Ces];
%S1 = solve(eqns1,[nes phies lambdaes W1 thetaes Omegaes rers Pcs Vs ...
 %   ncs dcs phics lambdacs W2 thetacs Omegacs Ccs Ces deltaT]);
%deltaT = double(S1.deltaT);

nes = (sigma*(Rfs-Rvs)+deltaT)*Qs*fs/(1-sigma*Rvs);
phies= Qs*fs/nes;
lambdaes = (Rfs-Rvs)*phies/Rvs;
W1 = (1+lambdaes)*beta*(Rfs-Rvs);
thetaes = (1-sigma)*W1/(lambdaes-sigma*phies*W1);
Omegaes = 1-sigma+sigma*thetaes*phies;

rers = ((2*v-1)/(1-v)*(v/(2*v-1)-Pes^(1-etap)))^(1/(1-etap));
Pcs = (v/(2*v-1)-((1-v)/(2*v-1))*rers^(etap-1))^(1/(1-etap));

Vs = Qs*fs-nes;
ncs = (sigma*((Rfs-Rds)*Qs*fs+(m/(1-m))*(Rvs-Rds)*Vs)+deltaT*Qs*fs)/(1-Rds*sigma);
dcs = (m/(1-m))*Vs+Qs*fs-ncs;
phics = (Qs*fs+Delta*Vs*(m/(1-m)))/ncs;
lambdacs = (phics*(Rfs-Rds))/Rds;
W2 = (1+lambdacs)*(Rfs-Rds)*beta;
thetacs = (1-sigma)/(lambdacs/W2-sigma*phics);
Omegacs = 1-sigma+sigma*phics*thetacs;

%Ccs = mrss+(Rvs-1)*dcs-deltaT*Qs*fs;
%Ces = (Ys-(1-v)*(Ccs+Is+Gs))/v-Is-Gs;
%bes = (Ces-mrss+deltaT*Qs*fs)/(rers*(Rvs-1));

%Ces = Ys - Is - Gs;
%Ccs = Ys - Is - Gs;
%O1 = (1-v)/v*(1-2*m)/(m*(1-m));
%Ces = (Ys-Is-Gs-(1-m)*O1*(Is+Gs))/(m+(1-m)*(1+O1));
%Ccs = (1+O1)*Ces+O1*(Is+Gs);

syms Ccs Ces;
eqns2 = [Ys==v*(Ces+Is+Gs)+(1-v)*(1-m)/m*(Ccs+Is+Gs),...
         Ys==(1-v)*m/(1-m)*(Ces+Is+Gs)+v*(Ccs+Is+Gs)];
S2 = solve(eqns2,[Ccs Ces]);
Ccs = double(S2.Ccs);
Ces = double(S2.Ces);
     

bes = (Ces - mrss + deltaT*Qs*fs-(1-sigma)*nes)/((Rds-1)*rers);

mues = (1/Ces)*(1-beta*bb)/(1-bb);
mucs = (1/Ccs)*(1-beta*bb)/(1-bb);

chie = mues*mrss;
chic = mucs*mrss;


save param_sw beta kappa spF psi epsp epsw alpha delta0 delta1 delta2 g bb sigma...
     upse upsc v Delta deltaT etap m phip phiw psik eta rhor phipi phiy sr...
     rhof sf rhoA sA rhot st indp indw ...
     Rds Rfs Rvs Rtrs Qs M1s M2s pws Ks Ys Is Gs ws fws fcbs mrss fs...
     Pes nes phies lambdaes thetaes Omegaes rers Pcs Vs ncs dcs phics...
     lambdacs thetacs Omegacs Ccs Ces bes mues mucs chie chic

%%
dynare tomeme4