% Climate DSGE Habit Model SMM Estimation
%
% =========================================================================

var a c ih il io uc w em cm ltau sdf anew yh yl y ph pl ra x qr qx L 
kh kl kd e qh ql Dh Dl DPh DPl Vh Vl n nh nl rh rl rf p pnew premh preml premhl theta 
cy dc ivy div iv dy ca logca dca hml rm;

varexo ea etau el;

parameters mu rhotau ltauss alpha beta gamma delta rho Khss Klss Yhss Ylss Yss Ihss 
Ilss Ioss Css Ws Anewss cms ems Ns Nhss Nlss ass phi1 varphi xsi rhocm 
 nu pi1 pi2 Kdss Rhss Rlss Vhss Vlss Dhss Dlss DPhss DPlss Pss Pnewss rhoa alpha1 alpha2 eta premhss premlss premhlss
sigma lambda xi eis gammah gammal rhoL Lss Ess omega;

alpha = 0.40;                                    %input share for capital
beta = 0.995;                                    %time discount
delta = 0.06;                                    %capital depreciation
rho = 0.9;                                       %logrithmic of growth autocorrelation 
mu=0.00000005;                                      %growth sensitivity to carbon
nu=0.6;                                             %share of dirty capital cumulation
rhoa=0.1;                                            %robert ready 2019 JME

omega=10;                                         %consumption risk aversion
varphi=6.6;                                         %preference for emission in consumption changed varphi 4-5 is sensible
phi1=0.033;                                       %emission generation 
xsi=7;                                            %aversion to leasure

rhotau=0.66;                                      %autocorrelation for climate policy
rhocm=0.93;                                        %autocorrelation for carbon accumulation
eta=1/3;                                          %convex adjustment cost power

pi1=0.15;                                            
pi2=0.01;                                            
kappa=0.95;                                       %solution to the relative shares
rhoL=0.9;                                             
ass=0.0018;  
Lss=0.6;                                          %Krussell 2022
Lss2=rhoa*ass/(1-rhoL)+Lss;
gammah=0.34;                                      %Fried 2019
gammal=1-gammah;                                  %Fried 2019
sigma=0.05;
lambda=0.5;
xi=0.15;
eis=0.05;
    
Ns=1/3;                                               

Plss=(gammal/(kappa^((eis-1)/eis)*gammah+gammal)^(1/(1-eis)));
Phss=(Plss*gammah/gammal*(kappa)^(-1/eis));
Nhss=Ns*Phss/Plss*kappa/(Phss*kappa/Plss*gammah+gammal);
Nlss=(Ns-Nhss*gammah)/gammal;
sdfss=beta;

%Kdss Solution
options = optimset('TolX',1e-15, 'TolFun',1e-15,'MaxIter',10000,'MaxFunEvals',100000,'Display','off');
FUN1=@(x) find_ss_kd_habit(x, log(Phss), log(Plss), alpha, nu, rhocm, mu, phi1, beta, delta, log(Nhss), ass, sigma);
Kdssolution=fsolve(FUN1,0,options);

if imag(Kdssolution)<10^(-10)
    Kdss=exp(real(Kdssolution));
    disp('no complex');
else
    Kdss=Kdssolution;
    disp('complex dirty capital');
end;

% Extraction Rate in Asset Pricing Model
options = optimset('TolX',1e-15, 'TolFun',1e-15,'MaxIter',10000,'MaxFunEvals',100000,'Display','off');
FUN2=@(x) find_ss_theta_habit1(x, log(Kdss), pi1, beta, xi, ass, lambda, Lss, rhoa, rhoL, gammah);
ThetaSolution=fsolve(FUN2,0.01,options);
thetass=ThetaSolution;

Ioss=(pi1/2*lambda*xi*thetass^2/((1/sdfss-1)*(1/sdfss-1+lambda)))^(1/(1-xi))*exp(Lss2);
Xss=Ioss^xi*exp(Lss2)^(1-xi)/lambda;
Qxss=(Ioss/exp(Lss2))^(1-xi)/xi;
Qrss=(1/sdfss-1+lambda)/lambda*Qxss;
Rass=Xss*lambda/thetass;
Pss=pi1*thetass+pi2+Qrss;

ems=Kdss*phi1*gammah;
cms=ems/(1-rhocm);
Anewss=exp(ass-mu*cms);   %interpret the cms accumulatio as long run effect 
Ws=Plss*(1-alpha)*Anewss*((1/sdfss-1+delta)/Plss/alpha)^(alpha/(alpha-1));
Klss=(Nlss*Anewss*((1/sdfss-1+delta)/(Plss*alpha))^(1/(alpha-1)));
Ylss=(1/sdfss-1+delta)*Klss/alpha/Plss;

Ess=(Ws*Anewss^(alpha-1)*Nhss^alpha/Phss/(1-alpha))^(1/alpha);
Khss=(Ess^((sigma-1)/sigma)/(1-nu)-nu*Kdss^((sigma-1)/sigma)/(1-nu))^(sigma/(sigma-1));

Yhss=Ylss*kappa;
Yss=(gammah*Yhss^((eis-1)/eis)+gammal*Ylss^((eis-1)/eis))^(eis/(eis-1));
ltauss=log((Phss*Yhss*nu*alpha/Kdss^(1/sigma)/Ess^(1-1/sigma)-Pss)/phi1);

Ihss=delta*Khss;
Ilss=delta*Klss;
Css=Yss-gammah*Ihss-gammal*Ilss-Ioss-(pi1/2*thetass^2+pi2*thetass)*Rass;
ucss=(Css-varphi*ems)^(-omega);
gamma=ucss*Ws*(1-Ns)^xsi;

alpha1=(delta)^(eta);
alpha2=delta-(delta)^(1-eta)/(1-eta)*alpha1;
Qhss=1;
Qlss=1;

Pss=pi1*thetass+pi2+Qrss;
Pnewss=Pss+phi1*exp(ltauss);
Dhss=Phss*Yhss-Ihss-Ws*Nhss-Pss*Kdss-exp(ltauss)*Kdss*phi1;
Dlss=Plss*Ylss-Ilss-Ws*Nlss;
Vhss=Dhss/(1-beta);
Vlss=Dlss/(1-beta);
Rhss=Vhss/(Vhss-Dhss);
Rlss=Vlss/(Vlss-Dlss);
Rfss=1/beta;
DPhss=Dhss/Vhss;
DPlss=Dlss/Vlss;
premhss=Rhss-Rfss;
premlss=Rlss-Rfss;
premhlss=Rhss-Rlss;

model;    
 anew=a-mu*cm; %1
 exp(uc) = (exp(c)-varphi*em)^(-omega); %2
 sdf=exp(uc)/exp(uc(-1))*beta; %3
 gamma-exp(uc)*exp(w)*((1-exp(n))^xsi)=0; %4

 //H Sector
 exp(e)=((1-nu)*exp(kh*((sigma-1)/sigma))+nu*exp(kd*((sigma-1)/sigma)))^(sigma/(sigma-1)); %5
 exp(w)=(1-alpha)*exp(yh+ph)/exp(nh); %6
 nu*alpha*exp(yh+ph)/(exp(kd*(1/sigma)+e*(1-1/sigma)))=exp(p)+phi1*exp(ltau); //7
 nu*alpha*exp(yh+ph)/(exp(kd*(1/sigma)+e*(1-1/sigma)))=exp(pnew); //8
 exp(qh)=sdf(1)*(alpha*(1-nu)*exp(yh(1)+ph(1))/exp(kh(1)*1/sigma+e(1)*(1-1/sigma))+exp(qh(1))*(1-delta+alpha1*eta*(exp((ih(1)-kh(1))*(1-eta)))/(1-eta)+alpha2)); //9 
 exp(qh)=exp((ih-kh)*eta)/alpha1; //10
 exp(kh)=((1-delta)*exp(kh(-1))+(alpha1/(1-eta)*exp((ih(-1)-kh(-1))*(1-eta))+alpha2)*exp(kh(-1))); //11
 em=exp(kd)*phi1*gammah; //12

 //L Sector
 exp(w)=(1-alpha)*exp(yl+pl)/exp(nl); //13
 exp(ql)=sdf(1)*(alpha*exp(yl(1)+pl(1))/exp(kl(1))+exp(ql(1))*(1-delta+alpha1*eta*(exp((il(1)-kl(1))*(1-eta)))/(1-eta)+alpha2));//14
 exp(kl)=((1-delta)*exp(kl(-1))+(alpha1/(1-eta)*exp((il(-1)-kl(-1))*(1-eta))+alpha2)*exp(kl(-1))); //15
 exp(ql)=exp((il-kl)*eta)/alpha1; //16

//Energy Sector
 exp(p)=pi1*theta+pi2+exp(qr); //17
 exp(qr)=sdf(1)*(exp(p(1))*theta(1)-pi1/2*theta(1)^2-pi2*theta(1)+exp(qr(1))*(1-theta(1))); //18
 exp(qx)=sdf(1)*(lambda*exp(qr(1))+(1-lambda)*exp(qx(1))); //19
 exp(ra(-1))*(1-theta)=exp(ra)-lambda*exp(x(-1)); //20
 exp(x(-1))*(1-lambda)+exp(io(-1))^xi*exp(L(-1))^(1-xi)-exp(x)=0; //21
 exp(qx)=exp((io-L)*(1-xi))/xi; //22
 exp(ra)*theta=gammah*exp(kd); //23

 //Final Good Sector
 exp(ph-pl)=exp((yh-yl)*(-1/eis))*gammah/gammal; //24
 exp(ph*(1-eis))*gammah^(eis)+exp(pl*(1-eis))*gammal^(eis)=1; //25
 exp(yh) = exp(anew*(1-alpha))*(exp(nh*(1-alpha))*exp(e*alpha)); //26
 exp(yl) = exp(anew*(1-alpha))*(exp(nl*(1-alpha))*exp(kl*alpha)); //27
 exp(y)=(gammah*exp(yh*((eis-1)/eis))+gammal*exp(yl*((eis-1)/eis)))^(eis/(eis-1)); //28

//dynamic process for state variable
 a = rho*a(-1)+(1-rho)*ass + ea; //29  
 ltau=rhotau*ltau(-1)+ltauss*(1-rhotau)+etau; //30
 L=(1-rhoL)*Lss+rhoL*L(-1)+rhoa*a+el;  //31

//Equity Performance of Aggregate Market
 Dh=exp(yh+ph)-exp(ih)-exp(w+nh)-exp(p+kd)-exp(ltau+kd)*phi1; //32
 Dl=exp(yl+pl)-exp(il)-exp(w+nl); //33
 exp(uc(-1)) = beta*exp(uc) * (exp(rf)); //34
 Vh = Dh+sdf(1) * Vh(1); //35
 Vl = Dl+sdf(1) * Vl(1); //36
 DPh = Dh/Vh; //37
 DPl = Dl/Vl; //38

 exp(rh)=Vh/(Vh(-1)-Dh(-1)); //39
 exp(rl)=Vl/(Vl(-1)-Dl(-1)); //40

//addition of the energy sector into the economy exogenously
 exp(n)=gammal*exp(nl)+gammah*exp(nh); //41
 exp(c)=exp(y)-gammah*exp(ih)-gammal*exp(il)-exp(io)-(pi1/2*theta^2+pi2*theta)*exp(ra); //42

//carbon dynamics
 cm=rhocm*cm(-1)+em(-1); //43

//premium expression
 premh=rh-rf(-1); //44
 preml=rl-rf(-1); //45
 premhl=rh-rl;    //46

//boring variables for moment matching
cy=exp(c-y);
dc=c-c(-1);
iv=log(exp(ih)*gammah+exp(il)*gammal+exp(io));
div=iv-iv(-1);
ivy=exp(iv-y);
dy=y-y(-1);
ca=em;
logca=log(ca);
dca=logca-logca(-1);
hml=premhl;
rm=gammah*premh+gammal*preml;

end;

initval;
c=log(Css);
em=ems;
uc=log(ucss);
kh=log(Khss);
kl=log(Klss);
kd=log(Kdss);
yh=log(Yhss);
yl=log(Ylss);
y=log(Yss);
nh=log(Nhss);
nl=log(Nlss);
n=log(gammah*Nhss+gammal*Nlss);
ih=log(Ihss);
il=log(Ilss);
io=log(Ioss);
a=ass;
anew=log(Anewss);
w=log(Ws);
ltau=ltauss;
Dh=(Dhss);
Dl=(Dlss);
sdf=beta;
cm=cms;
rh=log(Rhss);
rl=log(Rlss);
rf=log(Rfss);
Vh=(Vhss);
Vl=(Vlss);
qh=log(Qhss);
ql=log(Qlss);
p=log(Pss);
x=log(Xss);
pnew=log(Pnewss);
DPh=DPhss;
DPl=DPlss;
premhl=premhlss;
premh=premhss;
preml=premlss;
ph=log(Phss);
pl=log(Plss);
ra=log(Rass);
theta=thetass;
qr=log(Qrss);
qx=log(Qxss);
L=Lss2;
e=log(Ess);

cy=Css/Yss;
dc=0;
iv=log(Ihss*gammah+Ilss*gammal+Ioss);
div=0;
ivy=(Ihss*gammah+Ilss*gammal+Ioss)/Yss;
dy=0;
ca=ems;
logca=log(ems);
dca=0;
hml=0;
rm=0;
end;