%%%%% DSGE implementation of:
%     The Liquidity-Augmented Model of Macroeconomic Aggregates
%       Athanasios Geromichalos (UC Davis) and Lucas Herrenbrueck (Simon Fraser University)
%       First version: 2017; this version: 2020

%     This file uses code contributed by Johannes P. Strobel (Bundesbank)
%       and Marieh Azizirad (Simon Fraser University)


@#define interestraterule=2

@#if interestraterule==1
    title_string='Economy with strict money growth rule and endogenous interest rates' 
@#endif
@#if interestraterule==2
    title_string='Economy with Taylor rule and persistent money growth (random walk)'
@#endif
@#if interestraterule==3
    title_string='Economy with fixed interest rates (ZLB) and stationary money growth'
@#endif


%%%%%%%%%%%%
%1. Declare model variables, parameters, ss values
%%%%%%%%%%%%

var 
Y,C,H,K,I,           % output and friends
w,r,                 % factor prices
Q,                   % wedge between PM and CM price of output
pi,                  % inflation
ell,                 % aggregate liquidity premium
j,sB,                % policy rate and its inverse (*s*econdary market price of liquid *b*onds)
b,                   % bond-to-money ratio B(t)/M(t); we assume B(t) grows steadily at rate PIBAR
m,                   % gross money growth from AM to PM := (N+M)/M, standing in for new money N
Z,                   % CM real balances := M/P, standing in for price level P
tfp,                 % total factor productivity shock
nrs,                 % natural rate shock
lambda,              % fraction of shoppers (shock)
eta,                 % capital tradability (shock)
jshock,              % policy rate shock
mshock,              % AM money injection shock (only valid when interest rates are endogenous)
                     % display variables; the first four in percent, the latter two in percentage points:
output,consumption,investment,moneygrowth,inflation,policyrate;

% Exogenous innovations
varexo eps_a, eps_b, eps_l, eps_e, eps_j, eps_m; 


parameters           % See the parameter declaration below for definitions
BETA SIGMA KAPPA KAPPA0 ALPHA DELTA 
JBAR PIBAR 
TAU_A TAU_PI TAU_Y

%St. Dev's and persistence params of exogenous LOM's
STD_A STD_B STD_L STD_E STD_J STD_M
RHO_A RHO_B RHO_L RHO_E RHO_J RHO_M

%Steady State Computation
Y_ss, C_ss, H_ss, K_ss, I_ss, w_ss, r_ss, Q_ss, pi_ss, 
ell_ss, sB_ss, j_ss, Z_ss, m_ss, b_ss, 
A_ss, lambda_ss, eta_ss, jshock_ss, mshock_ss, tfp_ss, nrs_ss;


%%%%%%%%%%%%
%2. Parameterization
%%%%%%%%%%%%

%Quarterly calibration:

%%%%%%%%%%
%Real side
%%%%%%%%%%
BETA = 1/1.0125;     % Patience (3% discounting + 2% GDP growth annually; Herrenbrueck, 2019b)
SIGMA = 1;           % Inverse elasticity of intertemporal substitution (discussed in the paper, Section 4.2)
KAPPA = sqrt(.1);    % Inverse elasticity of labor supply (geom. mean. of two preferred values in Christiano-Trabandt-Walentin)
ALPHA = 0.25;        % Capital elasticity of production function (Christiano-Trabandt-Walentin)
DELTA = 0.025;       % Capital depreciation rate (Christiano-Trabandt-Walentin)
A_ss = 1;

%%%%%%%%%%
%Monetary policy
%%%%%%%%%%
PIBAR = 0.005;       % Inflation target (net, quarterly)
JBAR = 0.0075;       % Nominal interest rate intercept of 3% annually
sB_ss = 1/(1+JBAR);
pi_ss = PIBAR; 
m_ss = 1;            % normalized money growth between AM and PM (money growth across a period equals m*(1+PIBAR) )
TAU_A = 0.75;        % Smets-Wouters prior
TAU_PI = 1.5;        % Smets-Wouters prior
TAU_Y = 0.125;       % Smets-Wouters prior

%%%%%%%%%%
%Shock means, standard deviations, autoregressive parameters
%%%%%%%%%%
tfp_ss = log(A_ss);  % total factor productivity
nrs_ss = 0;          % natural rate shock
lambda_ss = 0.22;    % fraction of shoppers; mean calibrated to empirical short-run money demand function (discussed in the paper)
eta_ss = 0.5;        % fraction of tradable capital; mean calibrated to long-run return on capital (discussed in the paper)
jshock_ss = 0;       % policy rate shock
mshock_ss = 0;       % AM money growth shock (only valid when interest rates are endogenous)

STD_A = .01; STD_B = .01; STD_L = .01;STD_E = .01; STD_J = .01; STD_M = .01;

RHO_A = .5; RHO_B = .5; RHO_L = .5; RHO_E = .5; 
RHO_J = 0;           % zero since the Taylor rule is already autoregressive
RHO_M = 0;           % much better for interpretation; recall the m-shock affects AM-to-PM money growth, not CM-to-CM money growth!

%Steady state value computation
ell_ss = (1+PIBAR)/(BETA*(1+JBAR)) - 1;
j_ss = 1/sB_ss-1;
r_ss = 1 / (BETA*(1+eta_ss*ell_ss)) - 1 + DELTA;
Q_ss = 1 / ((1 + ell_ss/lambda_ss)*(1+j_ss));
w_ss = (1 - ALPHA) * 
    (A_ss * Q_ss)^(1 /(1 - ALPHA)) * 
    (ALPHA * BETA * (1 + eta_ss * ell_ss) / ( 1 - (1 - DELTA)* BETA * (1 + eta_ss * ell_ss) ) )^(ALPHA /(1 - ALPHA));
H_ss = 1/3;
K_ss = (Q_ss * ALPHA / r_ss )^(1 /(1 - ALPHA)) * H_ss; 
Y_ss = A_ss * H_ss^(1 - ALPHA) * K_ss^ALPHA;
b_ss = (1 - lambda_ss)/lambda_ss  * (1+j_ss) - eta_ss / (BETA * (1 + eta_ss*ell_ss)) * K_ss/(Q_ss*Y_ss);
Z_ss = Q_ss * Y_ss * (1+PIBAR) / m_ss;
C_ss = Y_ss - DELTA*K_ss;
KAPPA0 = w_ss / (C_ss^SIGMA * H_ss^KAPPA);
I_ss = DELTA*K_ss;

%%%%%%%%%%%%
%2. Declare equilibrium relations
%%%%%%%%%%%%
model;

%%%%%%%%%DEFINITIONS%%%%%%%

ell = lambda / ((1+j) * Q) - lambda;

j = 1/sB-1;

%%%%%%%%%HOUSEHOLDS%%%%%%%
1 / C^SIGMA = BETA / (1+nrs) / C(+1)^SIGMA / (1+pi(+1)) * ((1+ell(+1)) * (1+j(+1)));

1 / C^SIGMA = BETA / (1+nrs) / C(+1)^SIGMA * (r(+1) + 1 - DELTA) * (1+ell(+1) * eta(+1));

w / C^SIGMA = KAPPA0 * H^KAPPA;


%%%%%%%%%Production%%%%%%%

Y = exp(tfp) * K(-1)^ALPHA * H^(1-ALPHA);

r * K(-1) = ALPHA * Q * Y;

w * H = (1 - ALPHA) * Q * Y;

%%%%%%%%%Market Clearing%%%%%%%

m = lambda * (1 + sB * b(-1) + sB * (r + 1 - DELTA) *  (eta*K(-1) * m / (Q*Y)) ); 

(1+pi) * Q * Y = m * Z(-1); % LHS = spending (in last period's prices), RHS = PM real balances (in last period's prices)

C + K - (1-DELTA) * K(-1) = Y;

I = K - (1-DELTA) * K(-1);

@#if interestraterule==2
    Z = m * (1+pi_ss) / (1+pi) * Z(-1);    % with bygones-be-bygones persistent money growth rule (M.2)
@#else
    Z = (1+pi_ss) / (1+pi) * Z(-1);        % with stationary money growth rule (M.1)
@#endif

b = b(-1);           % here, we assume fixed bond/money ratio. Alternative fiscal assumptions could be coded here


%%%%%% Monetary policy rule variants %%%%%

@#if interestraterule==1
    m = 1+mshock;    % Strictly exogenous money growth (I.1)
@#endif

@#if interestraterule==2
    %% Taylor rule (I.2a); delayed 1 period
    log(1+j) = (1-TAU_A) * log(1 + j_ss) 
        +  TAU_A * log(1 + j(-1)) + (1-TAU_A) * TAU_Y * log(Y(-1) / Y_ss) 
        + (1-TAU_A) * TAU_PI * log( (1+pi(-1)) / (1+pi_ss) ) + jshock;
    
    %% Taylor rule (I.2b); delayed and responding to the year
    % log(1 / sB) = (1-TAU_A) * log(1 / sB_ss) 
    %     +  TAU_A * log(1 / sB(-1)) + (1-TAU_A) * TAU_Y * log(Y(-1) / Y_ss) 
    %     + (1-TAU_A) * TAU_PI * log( (pi(-1) * pi(-2) * pi(-3) * pi(-4)) / pi_ss^4) + jshock;
    
    %% Taylor rule (I.2c); contemporaneous
    % log(1+j) = (1-TAU_A) * log(1 + j_ss) 
    %     +  TAU_A * log(1 + j(-1)) + (1-TAU_A) * TAU_Y * log(Y / Y_ss) 
    %     + (1-TAU_A) * TAU_PI * log( (1+pi) / (1+pi_ss) ) + jshock;
@#endif

@#if interestraterule==3
    j = j_ss;        % Fixed interest rates, like at the ZLB (I.3)
@#endif


%%%%%%%%%SHOCKS%%%%%%%

tfp = (1 - RHO_A) * tfp_ss + RHO_A * tfp(-1) + eps_a; 

nrs = RHO_B * nrs(-1) + eps_b;

lambda = (1 - RHO_L) * lambda_ss + RHO_L * lambda(-1) + eps_l; 

eta = (1 - RHO_E) * eta_ss + RHO_E * eta(-1) + eps_e; 

jshock = RHO_J * jshock(-1) + eps_j;

mshock = RHO_M * mshock(-1) + eps_m;

%%%%%%%%%DISPLAY%%%%%%%

output = 100*log(Y);
consumption = 100*log(C);
investment = 100*log(I);
moneygrowth = 100*log( Z/Z(-1) * (1+pi) ); %NOTE: this is CM-to-CM money growth. CM-to-AM-to-PM different things can happen
inflation = 100*pi;
policyrate = 100*j;

end;

steady_state_model;
Y = Y_ss;
C = C_ss;
H = H_ss;
K = K_ss;
w = w_ss;
r = r_ss;
Q = Q_ss;
pi = pi_ss;
sB = sB_ss;
j = j_ss;
ell = ell_ss;
Z = Z_ss;
m = m_ss;
b = b_ss;
lambda = lambda_ss;
eta = eta_ss;
jshock = jshock_ss; %epsilon_j in the paper
mshock = mshock_ss;
I = I_ss;
nrs = nrs_ss;
tfp = tfp_ss;
output = 100*log(Y_ss);
consumption = 100*log(C_ss);
investment = 100*log(I_ss);
inflation = 100*pi_ss;
policyrate = 100*j_ss;
moneygrowth = 100*log(1+PIBAR);
end;

shocks;
    var eps_a; stderr 0.01;
    var eps_b; stderr 0.01; 
    var eps_l; stderr 0.01; 
    var eps_e; stderr 0.01;
    var eps_j; stderr 0.01; 
    var eps_m; stderr 0.01;
end;

steady;
check;

stoch_simul(order=1,irf=20) output,inflation,consumption,moneygrowth,investment,policyrate;
%stoch_simul(order=1,irf=20) output,inflation,consumption,policyrate,investment;

%% for welfare: 
% stoch_simul(order=1,irf=0,simul_replic=100,periods=1000,nodisplay) C,H;
% U = log(C) + KAPPA0/(1+KAPPA)*H.^(1+KAPPA)
% mean(U)