clear all
close all
clc

% set parameters
beta = 0.995;
kappa = 1-40^(-1);
spF = 0.03;
psi = 0.8;
epsp = 11;
epsw = 11;
alpha = 1/3;
delta0 = 0.025;
delta1 = 0.03;
delta2 = 0.01; % squared term capital adjustment cost
g = 0.2;
bb = 0.2; 
sigma = 0.9;
psid = 1.3;
gamma = 0.42;

upse = 0.9;
upsc = 0.9;
v = 0.9;
Delta = 1;
etap = 1.5;
m=0.3;
deltaT = 0.04;
Omegab = 1;
Omega2 = 1000;

% parameters unrelated to steady state
phip = 0.75;        % Calvo price setting
phiw = 0.75;        % Calvo wage setting
psik = 3; % capital adjustment cost
eta = 1; % Frisch elasticity

% Taylor rule
rhor = 0.6;
phipi = 1.1;
phiy = 0.1;
sr = 0.0025;

% Shock processes
rhof = 0.9;
sf = 0.01;
rhoA = 0.9;
sA = 0.01;
rhot = 0.95;
st = 0.01;
indp = 0;       % Price indexation
indw = 0;       % Wage indexation



% solve for steady state
Rds = beta^(-1);
Rfs = Rds;
Rvs = Rds;
Rtrs = Rds;
Qs = (Rfs - kappa)^(-1);

pws = (epsp-1)/epsp;
pks = 1;
v1s = pks;

Ks = (pws*alpha/(v1s*delta1))^(1/(1-alpha));
Ys = Ks^(alpha);
Is = delta0*Ks;
Gs = g*Ys;
ws = (1-alpha)*pws*Ks^(alpha);

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

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));

Rves = Rfs*psid;
Omegaes = (1-sigma*beta*(Rves/psid))^(-1)*(1-sigma);
alphaes = beta*Omegaes*(Rves/psid);
nes = (1-sigma*Rves/psid)^(-1)*(sigma*(Rfs-Rves/psid)+deltaT)*Qs*fs;
Vs = (Qs*fs-nes)/(psid*rers);

Omegacs = (1-sigma*beta*(Rds))^(-1)*(1-sigma);
alphacs = beta*Omegaes*Rds;
ncs = (1-sigma*Rds)^(-1)*((sigma*(Rfs-Rds)+deltaT)*Qs*fs+sigma*(Rvs-Rds)*(m/(1-m)*Vs));
dcs = Qs*fs+m/(1-m)*Vs-ncs;

phies = Qs*fs/nes;
phics = ((m/(1-m))*Vs + Qs*fs)/ncs;

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);

divkes = pks*Is*(1 - (psik/2)*(Is/Is - 1)^2) - Is;

des = (psid-1)*rers*Vs;
bes = (Ces - mrss - (Rds-1)*des + deltaT*Qs*fs-gamma*(1-sigma)*nes-divkes)/((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;

Ues = log(Ces-bb*Ces) - chie/(1+eta);
Ucs = log(Ccs-bb*Ccs) - chic/(1+eta);

RLes = 1/Qs + kappa;
RLcs = 1/Qs + kappa;

rles = log(RLes);
rlcs = log(RLcs);

xis = des/(rers*bes+des);


save param_sw beta kappa spF psi epsp epsw alpha delta0 delta1 delta2 g bb sigma psid gamma...
     upse upsc v Delta deltaT etap m phip phiw psik eta rhor phipi phiy sr...
     rhof sf rhoA sA rhot st indp indw Omegab Omega2 v1s...
     Rds Rfs Rvs Rves Rtrs Qs pws Ks Ys Is Gs ws fws fcbs mrss fs xis...
     Pes nes Omegaes alphaes rers Pcs Vs ncs dcs des phies phics...
     Omegacs alphacs Ccs Ces mues mucs chie chic divkes Ues Ucs...
     RLes RLcs rles rlcs
 
 dynare frictionless.mod
  