%First order approximation of the RBC model as described in the slides

clear
clc
syms k2 k1 k z s e alf bet rho
H=sym('1/(exp(z)*k^alf-k1)-bet*(alf*exp(rho*z+s*e)*k1^(alf-1))/(exp(rho*z+s*e)*k1^alf-k2)')

H1=diff(H,k2);
H2=diff(H,k1);
H3=diff(H,k);
H4=diff(H,z);
H5=diff(H,s);

digits(5);

alf=1/2; bet=.9; rho=.9;
k2=(alf*bet)^(1/(1-alf));
k1=(alf*bet)^(1/(1-alf));
k=(alf*bet)^(1/(1-alf));
z=0;
s=0;
e=0; %cheating to cancel the term H_5

H1v=vpa(eval(H1))
H2v=vpa(eval(H2))
H3v=vpa(eval(H3))
H4v=vpa(eval(H4))
H5v=vpa(eval(H5))

syms kk kz

S1=H1v*kk^2+H2v*kk+H3v

Sol1=solve(S1);

%This is our stable kk
kk=min(eval(Sol1(1)),eval(Sol1(2)))

S2=H1v*kk*kz+rho*H1v*kz+H2v*kz+H4v

Sol2=solve(S2);

%This is our kz
kz=Sol2(1)