clc; clear all; close all

%% This code computes the steady state of an environmental NK model 

% The calibration follows Annicchiarico and Di Dio (2015)
%% Structural Parameters (quarterly calibration)
beta=0.99;                    % discount factor
alpha=0.33;                   % elasticity of production wrt capital
epsilon=6;                    % elasticity of substitution btw differentiated goods
delta=0.025;                  % depreciation rate
sigma=1;                      % relative risk aversion
phi=1;                        % inverse of Frisch elasticity
g=0.0795;                     % public spending
a=1; %1.2480;                     % exogenous TFP in ss
eta=0.9979;                   % pollution decay
gamma1=0.45;                  % shifter in the emission function
gamma2=0;                     % concavity of emission function
theta1=0.1850;                % shifter in the abatement function
theta2=2.8;                   % convexity in the abatement function
d0=1.3950e-3;                 % constant in damage function
d1=-6.6722e-6;                % 1st order coefficient term in damage function
d2=1.4647e-8;                 % 2nd order coefficient term in damage function
d3=1;                         % damage function shifter
kappaL=19.8413;               % labor disutility shifter
erow=1.3331;                  % row emissions
tau= 200.0512;                  % emission tax

%% Steady State
pi=1;                        % inflation targeting 
rr=1/beta;                   % real interest rate
r=pi/beta;                   % nominal interest rate
q=1;                         % marginal value of investment (in terms of lambda)  
rk=1/beta-(1-delta);         % rental rate of capital

options = optimoptions('fsolve','MaxFunEvals',300000,'MaxIter',30000,'TolFun',1e-15);

Z0=[0.7780,800,0.2];                       % Initial guess
Z = fsolve(@(ZZ) find_steady_env(ZZ,eta,erow,gamma1,gamma2,tau,theta1,theta2,d0,d1,d2,d3,a,epsilon,alpha,rk,delta,kappaL,sigma,g,phi),Z0,options);

y=Z(1);                                                                               % output
x=Z(2);                                                                               % atmospheric carbon
h=Z(3);                                                                               % hours

e=x*(1-eta)-erow;                                                                     % emissions
mu=(gamma1*tau/(theta1*theta2)*y^(-gamma2))^(1/(theta2-1));                           % abatement
D=d3*(d0+d1*x+d2*x^2);                                                                % pollution damage
A=(1-D)*a;                                                                            % TFP 
mc=(epsilon-1)/(epsilon)-(1-gamma2)*gamma1*tau*y^(-gamma2)*(1-mu)-theta1*mu^(theta2); % marginal cost net of abatement and tax
k=alpha*mc*y/rk;                                                                      %  capital
i=delta*k;                                                                            % investment
w=(1-alpha)*mc*y/h;                                                                   % wage
c=y-i-g-y*theta1*mu^(theta2);                                                         % consumption
nu=e/y;                                                                               % emission/output ratio

%% Steady state parameters
gss=g;
ass=a;
piss=pi;
rss=r;
yss=y;
xss=x;
hss=h;
tauss=tau;

%% Parameters not affecting the ss
phipi=3;                
phiy=0;               
kappaI=2.48;              
rhoa=0.95;                 
rhog=0.97;                
rhom=0.15;                 
calvo=0.75;               
% adjusment cost coefficient to have the same linear Phillips Curve of the Calvo framework 
kappaP=(epsilon-1)*calvo/(piss^2*(1-calvo)*(1-beta*calvo)); 


% save parameters to use them in Dynare
save par beta alpha delta sigma phi gss ass kappaL piss rss epsilon d0 d1 d2 d3...
    gamma1 gamma2 theta1 theta2 eta erow tauss yss xss hss nu...
     phipi phiy kappaI rhoa rhog rhom kappaP

