%----------------------------------------------------------------
% 0. Housekeeping (close all graphic windows)
%----------------------------------------------------------------
close all;

%----------------------------------------------------------------
% 1. Defining variables
%----------------------------------------------------------------
var y pi pistar k i q d b l c w r a pm f f1 lambda lambdaC uc
    e_D e_A e_R;

varexo eta_D eta_A eta_R;

parameters  alpha beta gamma delta sigmaC sigmaL veps varphi epsilon 
            rho phi_pi phi_y eta h 
            Yss Kss Lss Iss Css Rss PMss LAMBDAss PMNss PIss Fss F1ss
            % shocks
			rho_D rho_A rho_R;

%----------------------------------------------------------------
% 2. Calibration
%----------------------------------------------------------------
alpha       = 0.3;              % Share of capital in production
beta    	= 0.99;			    % Discount factor savers
gamma       = 0.779;            % Proba of keeping prices fixed (Gertler Karadi)        
delta       = 0.025;            % Capital depreciation
sigma		= 1;				% Risk aversion consumption
veps        = 4.2;              % from Gertler Karadi
varphi		= 1;				% Labor disutility
epsilon     = 4.2;              % Retail firms elasticity of substitution (Gertler Karadi)
rho			= 0.80;				% Monetary policy rule, Smoothing, autocorr parameter
phi_pi		= 1.5;				% Monetary policy rule, Inflation (=1.5)
phi_y		= 0.5;              % Monetary policy rule, GDP (=0.1)
eta         = 0;                % Inflation indexation (=0.241 in Gertler Karadi)
h           = 0;                % Consumption habit (=0.7)
sigmaC	    = 1.5;		        % risk aversion consumption
sigmaL	    = 2;		        % labor disutility

% Implied by steady state
Kss = alpha/Rss * Yss;
Iss = delta * Kss;
Lss = (Yss/Kss ^alpha)^(1/(1-alpha));
Rss = 1/beta;
Yss = Kss^alpha*Lss^(1-alpha);
Css = Yss - Iss;
PMss = (veps-1)/veps;
LAMBDAss = 1;
Fss = Yss*PMss/(1-beta*gamma);
F1ss = Yss/(1-beta*gamma);
PIss = 1;

% shock processes (autoregressive parameter)
rho_D   = 0.80; 
rho_A   = 0.95;
rho_R 	= 0.40;

%----------------------------------------------------------------
% 3. Model
%----------------------------------------------------------------
model(linear); 
    %% Household
    % Euler equation savers
    sigmaC*c(+1)+ h*c = r-pi(+1);
    % Labour supply equation savers
    w = (1/sigmaL)*l + sigmaC*c;
    % Stochastic discount factor
	lambdaC = -sigmaC*c;
   	% IS curve
	y = (h/(1+h))*y(-1) + (1/(1+h))* y(+1)- 1/(sigmaC*(1+h))*(r-pi(+1)) + e_D;

    %%Intermediary firms
    % Capital accumulation
    k = i + (1-delta)*k(-1);
    % Intermediary firm production function
    y = a + alpha*k(-1) + (1-alpha)*l; 
    % Labour demand 
    w = (1-alpha)* Yss/Lss *(1 + y - l);
	% Capital demand
    r = alpha * Yss/Kss * (1 + y - k(-1));

    %% Final firms
    % Final firms prices
    Fss*f = Yss*PMss*(y+pm) + beta*gamma*LAMBDAss*Fss*(lambda(1)-eta*veps*pi+veps*pi(1)+f(1));
    lambda(1) + r = 0;
    F1ss*f1 = Yss*y + beta*gamma*LAMBDAss*PIss^(eta-1)*F1ss*(lambda(1)+eta*(1-veps)*pi-(1-veps)*pi(1)+f1(1));
    pistar = f - f1 + pi;
    pi = gamma*eta*pi(-1) + (1-gamma)* pistar;
    uc = -1/((1-beta*h)*(1-h))*(c-h*c(-1)-beta*h*(c(1)-h*c));
    lambda = uc - uc(-1);


	% Monetary Policy Rule
	r = rho*r(-1) +  (1-rho)*( phi_pi*pi + phi_y*y ) + e_R  ;
    % Ressource constraint
    y = c + i ;
    % Bank funding constraint
    r(-1)* b(-1) + d = r(-1)*d(-1) + b;
    % Investment equation
    q*i = b - r(-1)*b(-1);

    % Monetary policy shock
    e_A = rho_A * e_A(-1)+ eta_A;
    % Demand shock
    e_D = rho_D * e_D(-1) + eta_D;
    % Supply shock
    e_R = rho_R* e_R(-1) + eta_R;
	
	end;

steady;

check;

%----------------------------------------------------------------
% 4. Shocks
%----------------------------------------------------------------
shocks;
var eta_R;  stderr 1;
var eta_D;  stderr 1;
var eta_A;  stderr 1;
end;


stoch_simul (order=1, irf=10);