// addpath C:\dynare\4.3.3\matlab

var C, H, K, z;
varexo eps, tau_c, tau_k, tau_d, tau_h, tau_x;

parameters alpha, beta, sigma_s, chi_s, theta, delta, xi, nsigma_s, varsig, se;

// PARAMETERS
alpha=1.5;
chi_s=0.08;
beta=0.992; 
sigma_s=1.33; 
theta=0.33;
delta=0.0253; 
xi=0.8; 
varsig=0.985;
se=0.007;


// abbreviations
nsigma_s=(sigma_s-1)/sigma_s; 

model;

// FOC 1 Euler
1=beta*(C/C(+1))*(1+(exp(z(+1))*theta*((H(+1)/K)^(1-theta)))*(1-tau_k)/(1+tau_x)-delta-xi*tau_d)+(alpha*chi_s*xi*(((chi_s*(1/(xi*K(-1)))^nsigma_s)+((1-chi_s)*((1-H)^nsigma_s)))^(1/(sigma_s-1)))*(1/(xi*K))^((2*sigma_s-1)/sigma_s)/(((chi_s*(1/(xi*K(-1)))^nsigma_s)+((1-chi_s)*((1-H)^nsigma_s)))^nsigma_s));

// FOC 2
C*alpha*(1-chi_s)*((1-H)^(-1/sigma_s))*((chi_s*((1/(xi*K(-1)))^nsigma_s)+(1-chi_s)*((1-H)^nsigma_s))^(1/(sigma_s-1)))=((1-tau_h)*(1+tau_x)/(1+tau_c))*(exp(z)*(1-theta)*((K(-1)/H)^theta));

// Resource constraint
K=(1-delta)*K(-1)+(K(-1)^theta)*(H^theta)-C; 

z=varsig*z(-1)+se*eps;

end;

// no shocks in steady state


initval;
K=11.16;
C=0.99;
H=0.3;
tau_c=.076; 
tau_x=.0276;
tau_d=.0449;
tau_h=.32; 
tau_k=.293;
end;
steady;

endval;
K=11.16;
C=.99;
H=.3;
tau_c=.076; 
tau_x=.0276;
tau_d=.1;
tau_h=.32; 
tau_k=.293;
end;
steady;

shocks;
var eps; 
stderr 1;
end;

resid(1);

steady;
check;

stoch_simul(periods=500, irf=40, order=1);
dynasave('simudata.mat');