%% delay_models_jan24.m
% 1/7/2023: written by Carl walsh
% cleaned up version of delay_models_dec23.m
% loops over 'RE' and "CD" models, values of rhoi, phi_pi and k

%% housekeeping
close all
clear all
%clc
%% options
savefigs   = 'no'; % or 'no'
cdprograms = cd;
saveoutput   = 'C:\Users\cewjl\Dropbox\ResearchTopics\MonetaryPolicy\Hakamada\programs_jan24\output_jan24';
savefigures  = 'C:\Users\cewjl\Dropbox\ResearchTopics\MonetaryPolicy\Hakamada\programs_jan24\figures_jan24';
TK         = 'yes'; % if 'yes', loss as T x K, if 'no', loss as function of k for given T
%% model name
models  = {'RE' 'CD'}; % not yet set up to run model = 'CD', see delay_model_cd.m
T       = 20; % length of IRFs, baseline 20
%% parameters that vary: k, phi_pi, and rhoi
KK      = 7;             % maximum dely -- k runs from 0 to KK
phij    = [ 1.5 3];      % policy reaction coefficient on pi(t-k)
rho_i   = [ 0.85 0.5];   %[ coefficient on lagged policy rate: see run_all_models_dr_dyn.m to loop ove this coefficient
                             % Dennis 2009 estimates rho_i as 0.87,
                             % Coibion, 0.85 ordered first as it is our
                             % baseline
%% parameters to investigate
shk     = 0.16;    % Generates a *% (ann rate) rise in inflation under baseline policy with k = 4
mu      = 0;       %0.3; -- 1 (full indeation) or 0 (no indexation)
lambda  = 0.25;    % weight on putput gap for dual mandate with annualized inflation
%% fixed parameters
h       = 0.761; % Afsar, et al 2023, HNK model, p. 18	
beta    = 0.995;
sigma   = 1;  
disc    = 1.0; 
varphi  = 5;
alpha   = 0.25;
epsp    = 9;
epsw    = 4.5;
thetap  = 0.75; % Calvo parameter
eta     = 1.0; 
rho     = -log(beta)*100;
lamp    = (1-thetap)*(1-beta*thetap)*(1-alpha)/(thetap*(1-alpha+alpha*epsp));
varkap  = alpha*lamp/(1-alpha);
kappa   = lamp*(sigma/(1-h) + (varphi+alpha)/(1-alpha));
psiya   = (1+varphi)/(sigma*(1-alpha)+varphi+alpha);
psiwa   = (sigma+varphi)/(sigma*(1-alpha)+varphi+alpha);
rnss    = (1/beta) - 1;
elbss   = -rnss;
lambda_x= (sigma+(varphi+alpha)/(1-alpha)) ; % weight on inflation objective; can be set here to differ from theory based value of (sigma+(varphi+alpha)/(1-alpha)) = 7.56 
lambda_l= (sigma+(varphi+alpha)/(1-alpha))/(epsp/lamp);% or kappa/theta;
k0      = lamp;
k1      = kappa;
k2      = -lamp*sigma*h/(1-h);
phi_x   = 0.5/4;
% cognitive discounting
mbar     = 0.85; 
theta_cd = thetap;
Mf       = mbar*(theta_cd+(1-theta_cd)*(1-beta*theta_cd)/(1-beta*theta_cd*mbar));
Phi      =      (theta_cd+(1-theta_cd)*(1-beta*theta_cd)/(1-beta*theta_cd*mbar)); % M^f in paper is mbar*Phi
%% %%%%%%%%%%% Stochastic processes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rho_d   = 0.85;     %0.7;% demand shock
rho_z   = 0.96;     %0.7;% inflation trend shock persistency
rho_v   = 0.85;     % 0.92;%0.92;%0.7; % cost-push shock persistency % 0.92:Ilabaca, F., G. Meggiorini, and F.Milani (2020)
rho_zz  = 0.85;     % TFP shock persistency 
sigd    = 0.027;    % Billi and Walsh;  0.0357; % from Dobrew et al 2023 %0.001;
sigz    = 0.0005;
sigv    = shk;      % Size of inflation shock, set in line 27
sigzz   = 0.000199; % from Dobrew et al 2023
siged   = sqrt(1-rho_d^2)*sigd;
sigez   = sqrt(1-rho_z^2)*sigz;
sigezz  = sqrt(1-rho_zz^2)*sigzz;
sigev   = sqrt(1-rho_v^2)*sigv;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%%%%% Begin loops over models, k, phi_pi and rhoi %%%%%%%%%%%%%%%
for md  = 1:length(models)
    model = char(models(md));
    for k = 1:KK
        for jj = 1:length(phij)
            for m = 1:length(rho_i)
%k       = 5;        
               phi_pi  = phij(jj);
               rhoi    = rho_i(m);
               %% %%%%%%%%%%% Key matrices for solution once policy acts %%%%%%%%%%%%%%%%%
               % See appendix subsection A.2
               A0 = [rho_v 0 0 0 0
                   0 rhoi 0 0 0
                   0  0   0 0 0
                   0  0  0 mu k2
                   0 0 0 0 h];
               if strcmp(model,'RE')
                   A1 = [1 0 0 0 0
                       0 1 0 0 0
                       0 -1 1 0 0
                       -k0 0 0 1+beta*mu -k1
                       0 (1-h)/sigma 0 0 1+h];
                   A2 = [0 0 0 0 0
                       0 0 0 0 0
                       0 0 0 -1 0
                       0 0 0 beta 0
                       0 0 0 (1-h)/sigma 1];
               elseif strcmp(model,'CD')
                   A1 = [1 0 0 0 0
                       0 1 0 0 0
                       0 -1 1 0 0
                       -k0 0 0 1+beta*mu*Mf -k1
                       0 (1-h)/sigma 0 0 1+h];
                   A2 = [0 0 0 0 0
                       0 0 0 0 0
                       0 0 0 -mbar 0
                       0 0 0 beta*Mf 0
                       0 0 0 mbar*(1-h)/sigma mbar];
               end
               A3 = [0; 1-rhoi; 0; 0; 0];
               A4 = [1; 0; 0; 0; 0];
               F  = [0 0 0 phi_pi phi_x];
               %% %%%%%%%%%%%%% Solve for t+k, t+k+1, ...  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
               %  solves the model when policy follows the policy rule
               %  solution takes the form Z(t) = H0*Z(t-1) + H1*z(t)

               [H0, H1] = solvemodel(A0,A1,A2,A3,A4,F);

               %% %%%%%%%%%%%%% IRFs for no delay case (k = 0) %%%%%%%%%%%%%%%%%%%%%%%%%%%
               kk    = 0;
               N     = size(A1,1);
               IRF_0 = zeros(T,N);

               %% %%%%%%%%%%%%%% Solving for IRFs when policy follows the policy rule %%%%
               %   for time t when shock occurs, lagged values = 0 so Z(t) = H1*shk 
               %   for t+k, t+k+1, ...
               IRF_0(1,:) = sigv*H1';
               for i = 2:T
                   IRF_0(i,:) = IRF_0(i-1,:)*H0';
               end
               % convert i, r, and pi to annual rates
               IRF_0(:,[2 3 4]) = 4*IRF_0(:,[2 3 4]);

               %% Loss function choices: TK = 'yes' if function of T AND K, 'no' if function of k for given T
               % volatility as function of k for given T
               SS         = 1:T; SS = SS';
               if strcmp(TK,'no')
                   % V1_0       = [ave squared, scaled loss, contributions from inflation and output gap] 
                   % V2_0       = [ave cumulative scaled deviation, contributions from inflation and output gap] 
                   Vall_0       = cumsum((IRF_0.^2).*Disc);
                   V1_0         = Vall_0(T,:)/T;
                   Vloss1_0     =  1/2*epsp/lamp*V1_0*[0 0 0 1 lambda]';
                   V1_0         = [ Vloss1_0 1/2*epsp/lamp*V1_0(1,4) 1/2*epsp/lamp*lambda*V1_0(1,5)];
                   VLevel_0     = cumsum(IRF_0.*Disc); 
                   V2_0         = VLevel_0(T,:)./T;
                   VLevelloss_0 = 1/2*epsp/lamp*V2_0*[0 0 0 1 -lambda]';
                   V2_0         = [VLevelloss_0  1/2*epsp/lamp*V2_0(1,4) -1/2*epsp/lamp*lambda*V2_0(1,5)]; 
               end
               %% %%%%%%%%%%%%%%% volatility as function of T AND k  %%%%%%%%%%%%%%%%%%%%%
               % %%%%%%%%% Matric to discount future values from IRFs %%%%%%%%%%%%%%%%%%%
               N    = size(A1,1);
               Disc = zeros(T,N);
               for mm = 1:T
                   Disc(mm,:) = ones(1,N)*beta^(mm-1);
               end
               if strcmp(TK,'yes')
                   VallTK_0    = cumsum((IRF_0.^2).*Disc)./SS;
                   VlossTK_0   = 1/2*epsp/lamp*VallTK_0*[0 0 0 1 lambda]';
                   VLevelTK_0  = cumsum(IRF_0.*Disc)./SS;
                   VLevellossTK_0 =  1/2*epsp/lamp*VLevelTK_0*[0 0 0 1 -lambda]';
                   V1_0     = [VallTK_0(:,2:5) VlossTK_0 1/2*epsp/lamp*VallTK_0(:,4) 1/2*epsp/lamp*lambda*VallTK_0(:,5)];
                   V2_0    = [VLevelTK_0,VLevellossTK_0  1/2*epsp/lamp*VLevelTK_0(:,4) -1/2*epsp/lamp*lambda*VLevelTK_0(:,5)];
               end
               %% %%%%% creates numerica for phi_pi to use in file names %%%%%%%%%%%%%%%%%
               if phi_pi == phij(1)
                   phipi = 15;
               else
                   phipi = 30;
               end
               if rhoi == rho_i(1)
                   ri = 85;
               else
                   ri = 5;
               end

               %% %%%%%% change directory to location for saving output %%%%%%%%%%%%%%%%%%
               cd(saveoutput)
               eval(['IRF_',model,'_',num2str(kk),'_',num2str(phipi),'_',num2str(ri),'= IRF_0;']);
               saveIRFname     = ['IRF_',model,'_',num2str(kk),'_',num2str(phipi),'_',num2str(ri)];
               savefilename    = ['IRF_',model,'_',num2str(kk),'_',num2str(phipi),'_',num2str(ri)];
               save(savefilename,saveIRFname);
               
               eval(['V_',model,'_',num2str(kk),'_',num2str(phipi),'_',num2str(ri),'= V1_0;']);
               saveVname      = ['V_',model,'_',num2str(kk),'_',num2str(phipi),'_',num2str(ri)];
               savefilename    = ['V_',model,'_',num2str(kk),num2str(phipi),'_',num2str(ri)];
               save(savefilename,saveVname);
               
               eval(['VL_',model,'_',num2str(kk),'_',num2str(phipi),'_',num2str(ri),'= V2_0;']);
               saveVLname      = ['VL_',model,'_',num2str(kk),'_',num2str(phipi),'_',num2str(ri)];
               savefilename    = ['VL_',model,'_',num2str(kk),'_',num2str(phipi),'_',num2str(ri)];
               save(savefilename,saveVLname);

               % change directory back to location of the program code
               cd(cdprograms)

               %% %%%%%%%%%%% solve for pre-policy reaction equilibrium for t to t_k-1 %%%
               [R, ZNP] = nopolicy(k,A0,A1,A2,A4,H0,H1,sigv); 
               eval(['R_',num2str(k),'= R;']);
               disp('Unstacking ZNP')
               ZNPU = [ZNP(1:5,1)]; 
               eval(['ZNPU_',num2str(k),'= ZNPU;']);
               
               %% %%%%%%%%%%%%%%% IRFs for k > 0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
               IRF  = zeros(T,N);
               for s = 1:k
                   IRF(s,:)   = ZNP(1+(s-1)*N:s*N,1);
               end
               for j = k+1:T
                   IRF(j,:) = IRF(j-1,:)*H0';
               end
               
               %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
               % convert interest rates and inflation to annual rates

               IRF(:,[2 3 4]) = 4*IRF(:,[2 3 4]);
               
               %% %%%%%%% change directory to location for saved results
               cd(saveoutput)
               eval(['IRF_',model,'_',num2str(k),'_',num2str(phipi),'_',num2str(ri),'= IRF;']);
               saveIRFname = ['IRF_',model,'_',num2str(k),'_',num2str(phipi),'_',num2str(ri)];
               savefilename = ['IRF_',model,'_',num2str(k),'_',num2str(phipi),'_',num2str(ri)];
               save(savefilename,saveIRFname);
               
               %% volatility
               if strcmp(TK,'no')
                   Vall       = cumsum((IRF.^2).*Disc);
                   V1         = Vall(T,:)/T;
                   Vloss1     =  1/2*epsp/lamp*V1*[0 0 0 1 lambda]'; 
                   V1        = [ Vloss1 1/2*epsp/lamp*V1(1,4) 1/2*epsp/lamp*lambda*V1(1,5)]; 
                   V2         = VLevel(T,:)./T;
                   VLevelloss = 1/2*epsp/lamp*V2*[0 0 0 1 -lambda]'; 
                   V2         = [VLevelloss  1/2*epsp/lamp*V2(1,4) -1/2*epsp/lamp*lambda*V2(1,5)]; 
               end
               
               if strcmp(TK,'yes')
                   Vall  = cumsum((IRF.^2).*Disc)./SS;
                   Vloss = 1/2*epsp/lamp*Vall*[0 0 0 1 lambda]';
                   V1     = [Vall(:,2:5) Vloss 1/2*epsp/lamp*Vall(:,4) 1/2*epsp/lamp*lambda*Vall(:,5)];
                   VLevel  = cumsum(IRF.*Disc)./SS;
                   VLevelloss =  1/2*epsp/lamp*VLevel*[0 0 0 1 -lambda]';
                   V2     = [VLevel VLevelloss 1/2*epsp/lamp*VLevel(:,4) -1/2*epsp/lamp*lambda*VLevel(:,5)];
               end
               
               eval(['V_',model,'_',num2str(k),'_',num2str(phipi),'_',num2str(ri),'= V1;']);
               saveVname    = ['V_',model,'_',num2str(k),'_',num2str(phipi),'_',num2str(ri)];
               savefilename = ['V_',model,'_',num2str(k),'_',num2str(phipi),'_',num2str(ri)];
               save(savefilename,saveVname);
               
               eval(['VL_',model,'_',num2str(k),'_',num2str(phipi),'_',num2str(ri),'= V2;']);
               saveVLname   = ['VL_',model,'_',num2str(k),'_',num2str(phipi),'_',num2str(ri)];
               savefilename = ['VL_',model,'_',num2str(k),'_',num2str(phipi),'_',num2str(ri)];
               save(savefilename,saveVLname);
               
               %% %%% change directory back to location of program code %%%%%%%%%%%%%%%%%%
               cd(cdprograms)
               
               %% %%%%% end of loops %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            end % end loop over rhoi
        end % end of loop over phi_pi
    end % end of loop over k
    
    %% %%%%%%%%%%%%% call code for figures    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if strcmp(model,'RE')
        refigures;
    end
    if strcmp(model,'CD')
        cdfigures;
    end
    
    %% %%%%%%%%%%%%%%%%%%%

end % end loop over models

%% %%%%%%%%%%%%%%%%    extra figures       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if strcmp(model,'CD')
    figure(21)
    subplot(2,2,1)
    plot(t,IRF_RE_0_30_5(:,2),'k-o',t,IRF_RE_4_30_5(:,2),'o--k',t,IRF_CD_0_30_5(:,2),'d-b',t,IRF_CD_4_30_5(:,2),'d--b',t,zt,'k','LineWidth',1.5,'MarkerSize',markersize)
    title('Nominal rate')
    legend('RE: k=0','CD: k=0','Location','NorthEast')
    
    subplot(2,2,2)
    plot(t,IRF_RE_0_30_5(:,3),'k-o',t,IRF_RE_4_30_5(:,3),'o--k',t,IRF_CD_0_30_5(:,3),'d-b',t,IRF_CD_4_30_5(:,3),'d--b',t,zt,'k','LineWidth',1.5,'MarkerSize',markersize)
    title('Real rate')
    
    subplot(2,2,3)
    plot(t,IRF_RE_0_30_5(:,5),'k-o',t,IRF_RE_4_30_5(:,5),'o--k',t,IRF_CD_0_30_5(:,5),'d-b',t,IRF_CD_4_30_5(:,5),'d--b',t,zt,'k','LineWidth',1.5,'MarkerSize',markersize)
    title('Output gap')
    
    subplot(2,2,4)
    plot(t,IRF_RE_0_30_5(:,4),'k-o',t,IRF_RE_4_30_5(:,4),'o--k',t,IRF_CD_0_30_5(:,4),'d-b',t,IRF_CD_4_30_5(:,4),'d--b',t,zt,'k','LineWidth',1.5,'MarkerSize',markersize)
    title('Inflation rate')
    
    if strcmp(savefigs,'yes')
        cd(savefigures)
        print -dpng 'IRF_RE_CD_compare';
        savefig('IRF_RE_CD_compare'); 
        cd(saveoutput)
    end
end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(model,'CD')
    figure(22)
    subplot(2,2,1)
    plot(t,IRF_RE_0_15_85(:,2),'k-o',t,IRF_CD_0_15_85(:,2),'d--b','LineWidth',1.5,'MarkerSize',markersize)
    title('Nominal rate')
    legend('RE: k=0','CD: k=0','Location','NorthEast')
    
    subplot(2,2,2)
    plot(t,IRF_RE_0_15_85(:,3),'k-o',t,IRF_CD_0_15_85(:,3),'d--b',t,zt,'k','LineWidth',1.5,'MarkerSize',markersize)
    title('Real rate')
    
    subplot(2,2,3)
    plot(t,IRF_RE_0_15_85(:,5),'k-o',t,IRF_CD_0_15_85(:,5),'d--b',t,zt,'k','LineWidth',1.5,'MarkerSize',markersize)
    title('Output gap')
    
    subplot(2,2,4)
    plot(t,IRF_RE_0_15_85(:,4),'k-o',t,IRF_CD_0_15_85(:,4),'d--b','LineWidth',1.5,'MarkerSize',markersize)
    title('Inflation rate')
       
    if strcmp(savefigs,'yes')
        cd(savefigures)
        print -dpng 'IRF_RE_CD_0_15_85';
        savefig('IRF_RE_CD_0_15_85'); 
        cd(saveoutput)
    end
end
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(18)
subplot(2,1,1)
plot(t,V_RE_0_15_85(:,6),'or','LineWidth',1.5,'MarkerSize',markersize+2)
hold on
plot(t,V_RE_1_30_5(:,6),'d-b',t,V_RE_2_30_5(:,6),'o--b',t,V_RE_3_30_5(:,6),'s-.b',t,V_RE_4_30_5(:,6),'d:b',t,zt,'k','LineWidth',1.5,'MarkerSize',markersize)
axis([-Inf Inf 0 0.4]);
legend('k=0, baseline policy', 'k=1', 'k=2', 'k=3','k=4','Location','NorthEast')
title('Rational expectations: Delay and aggression')

subplot(2,1,2)

plot(t,V_CD_0_15_85(:,6),'or','LineWidth',1.5,'MarkerSize',markersize+2)
hold on
plot(t,V_CD_1_30_5(:,6),'d-b',t,V_CD_2_30_5(:,6),'o--b',t,V_CD_3_30_5(:,6),'s-.b',t,V_CD_4_30_5(:,6),'d:b',t,zt,'k','LineWidth',1.5,'MarkerSize',markersize)
axis([-Inf Inf 0 0.4]);
title('Cognitive discounting: Delay and aggression')

if strcmp(savefigs,'yes')
    cd(savefigures)
    print -dpng 'loss_RE_CD_V1_comp';
    savefig('loss_RE_CD_V1_comp'); 
    cd(saveoutput)
end
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(19)
subplot(2,1,1)
plot(t,VL_RE_0_15_85(:,6),'or','LineWidth',1.5,'MarkerSize',markersize+2)
hold on
plot(t,VL_RE_1_30_5(:,6),'d-b',t,VL_RE_2_30_5(:,6),'d--b',t,VL_RE_3_30_5(:,6),'d-.b',t,zt,'k','LineWidth',1.5,'MarkerSize',markersize)
axis([-Inf Inf 0 8])
legend('k=0, baseline policy', 'k=1', 'k=2', 'k=3','Location','NorthEast')
title('Rational expectations: Delay and aggression')

subplot(2,1,2)
plot(t,VL_CD_0_15_85(:,6),'or','LineWidth',1.5,'MarkerSize',markersize+2)
hold on
plot(t,VL_CD_1_30_5(:,6),'o-b',t,VL_CD_2_30_5(:,6),'d--b',t,VL_CD_3_30_5(:,6),'d-.b',t,zt,'k','LineWidth',1.5,'MarkerSize',markersize)
axis([-Inf Inf 0 8])
title('Cognitive discounting: Delay and aggression')

if strcmp(savefigs,'yes')
    cd(savefigures)
    print -dpng 'loss_RE_CD_V2_comp';
    savefig('loss_RE_CD_V2_comp'); 
    cd(saveoutput)
end

%% %%%%%%% change location back to location of program code %%%%%%%%%%%%%%%
cd(cdprograms)
