%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2021 Alessandro Sardone.
% Permission is granted to copy, distribute and/or modify
% this document under the terms of the GNU Free
% Documentation License, Version 1.3 or any later version
% published by the Free Software Foundation; with no
% Invariant Sections, no Front-Cover Texts, and no
% Back-Cover Texts. A copy of the license is included in
% the section entitled GNU Free Documentation License".
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear, clc

%% Additional Parameters

T=30;             % in T years the tax reaches the new steady state
G=100;             % total periods of simulation

muEINIT=0;                       % abatement effort
muYINIT=0;                       % abatement effort

muEEND=1;                       % abatement effort
muYEND=1;                       % abatement effort

%% The process for tau
% We already have an initial and a final value for mu. We need to specify
% the entire process. We assume that mu move gradually from mu initval to mu endval in T
% period, increasing at decreasing rates

SY=muYEND-muYINIT;             % difference between final and initial abatement
SE=muEEND-muEINIT;             % difference between final and initial abatement

step_stepY=2*(SY)/((T-2)*(T-1)); % in each step, tau grows less compared to the previous step, by step_step
step_stepE=2*(SE)/((T-2)*(T-1));

step=zeros(T-1,2);             % vector of steps
step(T-1,:)=0;
for j=2:1:T-1
    step(T-j,1)=step(T-j+1,1)+step_stepY;
    step(T-j,2)=step(T-j+1,2)+step_stepE;
end
 
mus=zeros(G,2);              % vectors of carbon taxes
mus(1,1)=muYINIT;
mus(1,2)=muEINIT;
for j=2:T                           
    mus(j,1)=mus(j-1,1)+step(j-1,1);
    mus(j,2)=mus(j-1,2)+step(j-1,2);
end
 
mus(T+1:G,1)=ones(G-T,1)*muYEND; % after period T, the carbon tax remains at the final value
mus(T+1:G,2)=ones(G-T,1)*muEEND;
musY=mus(:,1);
musE=mus(:,2);

%% Perfect foresight simulation
dynare tank
