var C, H, K, E, S;

// varexo S;

parameters beta, iota, sigma_c, sigma_s, chi_s, eta_h, sigma_h, theta, delta, gamma, xi, rho, nsigma_h, nsigma_c, nsigma_s;

// PARAMETERS
chi_s=0.95;
iota=0.64;
beta=0.992; 
sigma_c=5.32; 
sigma_s=1.33; 
eta_h=0.97;
sigma_h=0.02; 
theta=0.33;
delta=0.0253; 
gamma=0.0; 
xi=0.8; 
rho=0.8; 

// abbreviations
nsigma_h=(sigma_h-1)/sigma_h;
nsigma_c=(sigma_c-1)/sigma_c; 
nsigma_s=(sigma_s-1)/sigma_s; 

model;

// FOC 1 - Intertemporal Euler equation
(C^(-1/sigma_c))*((iota*(C^nsigma_c)+(1-iota)*((chi_s*(S^nsigma_s)+(1-chi_s)*((1-H)^nsigma_s))^(nsigma_c/nsigma_s)))^(1/(sigma_c-1)))=beta*(C(+1)^(-1/sigma_c))*((iota*(C(+1)^nsigma_c)+(1-iota)*((chi_s*(S(+1)^nsigma_s)+(1-chi_s)*((1-H(+1))^nsigma_s))^(nsigma_c/nsigma_s)))^(1/(sigma_c-1)))*(1+(eta_h*theta*((H/K)^(1-theta))*(((K^theta)*(H^(1-theta)))^(-1/sigma_h))*((eta_h*((K^theta)*(H^(1-theta)))^nsigma_h+(1-eta_h)*(E^nsigma_h))^(1/sigma_h-1))));

// FOC 2 - Intratemporal Euler equation
(1/(iota*(C^(-1/sigma_c))))*((1-iota)*(1-chi_s)*(H^(-1/sigma_s))*((chi_s*(S^nsigma_s)+(1-chi_s)*((1-H)^nsigma_h))^((sigma_c-sigma_s)/(sigma_c*sigma_s-sigma_c))))=(eta_h*(1-theta)*((K/H)^theta)*(((K^theta)*(H^(1-theta)))^(-1/sigma_h))*((eta_h*((K^theta)*(H^(1-theta)))^nsigma_h+(1-eta_h)*(E^nsigma_h))^(1/sigma_h-1)));

// FOC 3 - Intratemporal on environment and consumption
(1/((1+gamma)*(C^(-1/sigma_c))))*(1-iota)*chi_s*(S^(-1/sigma_s))*((1-iota)*(1-chi_s)*(H^(-1/sigma_s))*((chi_s*(S^nsigma_s)+(1-chi_s)*((1-H)^nsigma_h))^((sigma_c-sigma_s)/(sigma_c*sigma_s-sigma_c))))=(-1)*(1/(xi+rho))*(1-eta_h)*(E^(-1/sigma_h))*((eta_h*(((K^theta)*(H^(1-theta)))^nsigma_h)+(1-eta_h)*(E^nsigma_h))^(1/sigma_h-1));
// note the -1 multiplied bc lambda equals -1 multiplied by...

// Aggregate resource constraint
/* k' = (1-delta)*k + y - c */ 
K=(1-delta)*K(-1)+((eta_h*(((K^theta)*(H^(1-theta)))^nsigma_h)+(1-eta_h)*(E^nsigma_h))^(1/nsigma_h))-C; 

// environment
S=S(-1)*(1+gamma)-E*(xi+rho);

end;

initval;
K=100.0;
C=100.0;
H=0.3;
E=100.0;
S=100.0;

end;

homotopy_setup;
chi_s, 1.0, 0.95; 
iota, 0.64; 
end;

steady(homotopy_mode = 1, homotopy_steps = 50);