var EM Q QW QK xW xK cW p cK w n inv r kC kX invC invX k nC nX;
var YC YX A KC KX NC NX;
varexo e;
parameters gamma Qbar h s1 s2 epsilon s delta s3 beta alpha rho;
gamma = 0.01;
Qbar = 1;
//environmental quality without carbon emissions
h = 0.1;
//natural absorption of carbon
s1 = 0.2;
s2 = 0.3;
epsilon = 1;
//the ratio or w*n allocated to workers
s = 0.9;
//ratio of workers in total population 1
delta = 0.025;
s3 = 0.4;
beta = 0.985;
alpha = 0.35;
rho = 0.95;
model;
//Environment sector
//1-carbon emission
EM = (YC)^(1-gamma);
//2-nature environmental quality
Q = h*Qbar + (1-h)*Q(-1) - EM;
//3-Environmental quality for workers
QW = Q + xW;
//4-Environmental quality for capitalists
QK = Q + xK;
//Workers sector//
//5-two goods tradeoff
s1 * cW = p * QW;
//6-labor supply
s2 * cK = epsilon * w * (1-n);
//7-workers' budget constraint
cW + p * xW = epsilon * w * n;
//Capitalists sector//
//8-capitalists' budget constraint
cK + p * xK + inv = r * k(-1) + (1-epsilon)/(1-s) * w * n;
//9-law of motion of capital for C
kC = (1-delta)*kC(-1) + invC;
//10-law of motion of capital for X
kX = (1-delta)*kX(-1) + invX;
//11-total k
k = kC + kX;
//12-capitalists' consumption tradeoff
s3 * cK = p * QK;
//13-euler equation
cK(+1)/cK = beta*(r(+1)+1-delta);
//Production sector//
//14-production of C
YC = A * KC(-1)^alpha * NC^(1-alpha);
//15-production of X
YX = A * KX(-1)^alpha * NX^(1-alpha);
//16-demand for labor for C
w = (1-alpha) * YC / NC;
//17-demand for labor for X
w = p * (1-alpha) * YX / NX;
//18-demand for capital for C
r = alpha * YC / KC(-1);
//19-demand for capital for X
r = p * alpha * YX / KX(-1);
//Labor market equilibrium
//20-labor market equilibrium for C
NC = s * nC;
//21-labor market equilibrium for X
NX = s * nX;
//22-total labor division
n = nC + nX;
//Capital market equilibrium
//23-capital market equilibrium for C
KC = (1-s) * kC;
//24-capital market equilibrium for X
KX = (1-s) * kX;
//Goods market equilibrium
//25-C
YC = s*cW + (1-s)*cK + (1-s)*invC;
//26-X
YX = s*xW + (1-s)*xK + (1-s)*invX;
//Shocks
//27-productivity shocks
log(A) = rho*log(A(-1)) + e;
end;
initval;
A = 1;
p = 1;
r = 1/beta + delta - 1;
w = (1-alpha)*(r/alpha)^(-alpha/(1-alpha));
YC = (epsilon*(1-alpha)/s2*(2*s+s*s1-s1-s*s3-1)*((r/alpha)^(-alpha/(1-alpha))))/(2*epsilon*(1-alpha)/s2*(2-s3+2*s2+s1-s1/s-s2/s-1/s)+(1-s3)*(1-delta*alpha/r));
YX = YC;
NC = (r/alpha)^(alpha/(1-alpha))*YC;
KC = YC*alpha/r;
nC = NC/s;
kC = KC/(1-s);
nX = nC;
n = nX + nC;
kX = kC;
cW = epsilon*w*(1-n)/s2;
cK = (1+2*epsilon*(1-alpha)/s2-delta*alpha/r)/(1-s)*YC-s*epsilon*(1-alpha)/s2/(1-s)*((r/alpha)^(-alpha/(1-alpha)));
xW = epsilon*w*n - cW;
xK = s3*cK - s1*cW + xW;
EM = YC^(1-gamma);
Q = Qbar - EM/h;
NX = NC;
KX = KC;
k = kC + kX;
QW = Q + xW/h;
QK = Q + xK/h;
invC = delta * kC;
invX = delta * kX;
inv = invC + invX;
end;
steady;
check(qz_zero_threshold=1e-20);
model_diagnostics;
model_info;
shocks;
var e;
stderr 0.01;
end;
stoch_simul;