// est time series are first difference of logs
// addpath C:\dynare\4.3.3\matlab
// C:\Publication and research\Journal Articles\2013 - Carbon taxes and Productivity\analysis\dynare
// dC, dH, dK, dY, dS, dE

var C, H, K, E, S, p, y, A_k, A_e, A_x, A_u, A_s, z_k, z_e, z_x, z_u, z_s;
varexo eps_k, eps_e, eps_u, eps_x, eps_s, tau_marg, tau_c, tau_d, tau_k, tau_x, tau_h;

parameters ga, gs, sig, nu, mu, phi, eta, alpha, beta, chi_h, chi_e, theta_k, theta_e, theta_h, delta, xi, rho_corr, rho_k, rho_e, rho_x, rho_u, rho_s;

// PARAMETERS
alpha=.46; 
mu=.88;
chi_h=.95; // similar to the cobb douglas shares in Golosov ECTA, but calibrated to EUKLEMS data on electricity sector share
chi_e=1-chi_h;
beta=0.98282; 
theta_k=.3;
theta_e=.05;
theta_h=1-theta_k-theta_e;
delta=0.1; // annual depreciation 
nu=-5; // using .05 instead as (.05-1)/.05 p17 of krusell, olovsson, and hassler 2011 
xi=.7; 
ga=.037; // balanced growth
gs=2.26-1.23*ga; // .0098449; // something else p21 krusell paper, gs=2.26-1.23ga
rho_k=0.98; // productivity shock
rho_e=.96; // energy shock
rho_u=.98; // preference
rho_x=.98; // investment
rho_s=.96; // emissions shock
rho_corr=-.2; // correlation between K-H and E shocks
eta=.5; // Hall 2008

// turn off parameter value when doing sensitivity on it in .m file
phi=-1;
sig=1/(1-nu); // (sig-1)/sig=nu


model;

// Environment - this is how energy leads to emissions
S=(E*xi)*A_s;

// FOC 1 intertemporal Euler
// complete specification: A_u*((mu*(C^alpha*(1-H)^(1-alpha))^phi+(1-mu)*(1/S)^phi)^((1-eta-phi)/phi))*((C^alpha*(1-H)^(1-alpha))^(phi-1))*alpha*mu*C^(alpha-1)*(1-H)^(1-alpha)=beta*A_u(+1)*((mu*(C(+1)^alpha*(1-H(+1))^(1-alpha))^phi+(1-mu)*(1/S(+1))^phi)^((1-eta-phi)/phi))*((C(+1)^alpha*(1-H(+1))^(1-alpha))^(phi-1))*alpha*mu*C(+1)^(alpha-1)*(1-H(+1))^(1-alpha)*((1-delta)+((1-tau_k)/(1+tau_x))*(theta_k*A_k(+1)*(A_k(+1)*K)^(theta_k-1)*(A_k(+1)*H(+1))^theta_h*(A_e(+1)*E(+1))^theta_e));
// simple cobb douglas: A_u*C(+1)=beta*C*A_u(+1)*(1-delta+((1-tau_k)/(1+tau_x))*(theta_k*A_k(+1)*(A_k(+1)*K)^(theta_k-1)*(A_k(+1)*H(+1))^theta_h*(A_e(+1)*E(+1))^theta_e));
A_u*((mu*(C^alpha*(1-H)^(1-alpha))^phi+(1-mu)*(1/S)^phi)^((1-eta-phi)/phi))*((C^alpha*(1-H)^(1-alpha))^(phi-1))*alpha*mu*C^(alpha-1)*(1-H)^(1-alpha)=beta*A_u(+1)*((mu*(C(+1)^alpha*(1-H(+1))^(1-alpha))^phi+(1-mu)*(1/S(+1))^phi)^((1-eta-phi)/phi))*((C(+1)^alpha*(1-H(+1))^(1-alpha))^(phi-1))*alpha*mu*C(+1)^(alpha-1)*(1-H(+1))^(1-alpha)*((1-delta)+((1-tau_k)/(1+tau_x))*(theta_k*A_k(+1)*(A_k(+1)*K)^(theta_k-1)*(A_k(+1)*H(+1))^theta_h*(A_e(+1)*E(+1))^theta_e));


// r=theta_k*A_k*(A_k*K(-1))^(theta_k-1)*(A_k*H)^theta_h*(A_e*E)^theta_e;
// w=theta_h*A_k*(A_k*H)^(theta_h-1)*(A_k*K(-1))^theta_k*(A_e*E)^theta_e;
// marginal product p=theta_e*A_e*(A_k*H)^theta_h*(A_k*K(-1))^theta_k*(A_e*E)^(theta_e-1);
// solv for p if just labor input p=(theta_h*A_k*(A_k*H)^(theta_h-1)*(A_k*K(-1))^theta_k*(A_e*E)^theta_e+tau_d*xi*A_k)/A_k

// try 1: if using labor input for energy production
// p=theta_e*A_e*(A_k*H)^theta_h*(A_k*K(-1))^theta_k*(A_e*E)^(theta_e-1);
// p=(theta_h*A_k*(A_k*H)^(theta_h-1)*(A_k*K(-1))^theta_k*(A_e*E)^theta_e+tau_d*xi*A_k)/A_k;


C*(1-alpha)=alpha*(1-H)*(1-tau_marg)*(theta_h*A_k*(A_k*H)^(theta_h-1)*(A_k*K(-1))^theta_k*(A_e*E)^theta_e);

// simpler (1-alpha)*(mu*(1-H)^phi+(1-mu)*(1/((E*xi)*A_s))^phi)^((1-phi)/phi)*(1-H)^(phi-1)*mu*C/((mu*(1-H)^phi+(1-mu)*(1/((E*xi)*A_s))^phi)^(1/phi))=alpha*(1-tau_marg)*(theta_h*A_k*(A_k*H)^(theta_h-1)*(A_k*K(-1))^theta_k*(A_e*E)^theta_e);

// FOC 2 intratemporal euler; (dU/d1-h)/(dU/c)=w*tax wedge
// C*(1-alpha)=alpha*(1-H)*(1-tau_marg)*(theta_h*A_k*(A_k*H)^(theta_h-1)*(A_k*K(-1))^theta_k*(A_e*E)^theta_e);


// Output for final goods
y=(A_k*K(-1))^theta_k*(A_k*H)^theta_h*(A_e*E)^theta_e;


// just have labor input (p-tau_d*xi)*A_k=theta_h*A_k*(A_k*H)^(theta_h-1)*(A_k*K(-1))^theta_k*(A_e*E)^theta_e;

// energy demand
// E=((p^(-1)*theta_e*A_e*(A_k*H)^theta_h*(A_k*K(-1))^theta_k)^(1/(1-theta_e)))*(1/A_e);

// try 2: solve for p and substitute into final good's energy demand
// (theta_e-1)*(1/A_e)*(p^(theta_e-2))*(theta_e*A_e*(A_k*H)^theta_h*(A_k*K(-1))^theta_k)^(1/(1-theta_e))*(p-tau_d*xi-1)=tau_d*xi*((p^(-1)*theta_e*A_e*(A_k*H)^theta_h*(A_k*K(-1))^theta_k)^(1/(1-theta_e)))*(1/A_e);

// with return on capital included as cost: 1=(1/(1-theta_e))*(1/p)^theta_e*(theta_k*A_k*(A_k*K(-1))^(theta_k-1)*(A_k*H)^theta_h*(A_e*E)^theta_e+tau_d*xi-p);
(1/(1-theta_e))*(1/p)^theta_e*(p-tau_d*xi-1)+1=0;
E=A_e^(-1)*(1/p)^(1/(1-theta_e))*((theta_e*A_e*(A_k*H)^theta_h*(A_k*K(-1))^theta_k)^(1/(1-theta_e)));


// Resource constraint: final output = consumption + investment
y=C+K-(1-delta)*K(-1);

// note: correlated shocks as per krusell, hassler, olovsson
z_k=rho_k*z_k(-1)+rho_e*z_e(-1)+eps_k; A_k=exp(z_k);
z_e=rho_e*z_e(-1)+rho_k*z_k(-1)+eps_e; A_e=exp(z_e);
z_u=rho_u*z_u(-1)+eps_u; A_u=exp(z_u);
z_x=rho_x*z_x(-1)+eps_x; A_x=exp(z_x);
z_s=rho_s*z_s(-1)+eps_s; A_s=exp(z_s);


end;



// no shocks in steady state
//steady state initial values - use fsolve to do this
initval;
p=.2;
// r=.03; w=.15;
C=0.5;
H=0.3;
K=0.8158;
E=0.4;
S=0.2;
y=0.3853; 
A_k=1; A_e=1; A_u=1; A_x=1; A_s=1;
eps_k=0; eps_e=0; eps_u=0; eps_x=0; eps_s=0;
tau_c=.076; tau_x=.0276; tau_h=.32; tau_k=.293;
tau_marg=.32; // .32;
tau_d=0.05; // .0449; 

end;


steady;

shocks;
var eps_k; stderr 0.03;
var eps_e; stderr .03;
var eps_u; stderr .03;
var eps_x; stderr .03;
var eps_s; stderr .05;
var eps_k, eps_e = rho_corr*.03*.03; // correlated shock
// var tau_d;
// periods 1:10;
// values 0;
end;



resid(1);

steady;
check;

// observation_trends;
// y (ga);
// S (gs);
// end;

estimated_params;
phi, normal_pdf, -1, 2, -2, 1;
xi, beta_pdf, .7, .2;
mu, beta_pdf, .9, .1;
alpha, beta_pdf, .5, .1;
rho_k, beta_pdf, .95, .03;
rho_e, beta_pdf, .95, .03;
rho_u, beta_pdf, .95, .03;
rho_x, beta_pdf, .95, .03;
rho_s, beta_pdf, .95, .05;
rho_corr, normal_pdf, -.1, .2;
stderr eps_k, inv_gamma_pdf, .01, inf;
stderr eps_e, inv_gamma_pdf, .01, inf;
stderr eps_u, inv_gamma_pdf, .01, inf;
stderr eps_x, inv_gamma_pdf, .01, inf;
stderr eps_s, inv_gamma_pdf, .01, inf;
end;

varobs H, C, y, K, E;
write_latex_dynamic_model;

// adjust mh_jscale so that acceptance rate is between 25-50%
// estimation(datafile=kdata2,mh_replic=6000,mh_nblocks=4,mh_drop=.5,mh_jscale=.6,conf_sig=.95);
// estimation(datafile=kdata2,mode_compute=4,mh_replic=0);
stoch_simul(periods=1000, drop=200, solve_algo=0, order=2, hp_filter=6.25,nograph, noprint); // 
// extended_path(periods=42);