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

%----------------------------------------------------------------
% 1. Defining variables
%----------------------------------------------------------------
var  c r r_n r_k p pm w y ym k l i q b dep d pstar f f1 lambda uc 
     e_S e_A e_R;

%predetermined_variables b;
varexo eta_S eta_A eta_R;

parameters  alpha beta gamma gammaf gammap delta rho phi_p phi_y h sigmaC sigmaL veps kappa
            Rss LAMBDAss Pss Lss Qss Iss Kss Yss Css PMss Fss F1ss MUss Wss Bss Dss
            IKss KYss IYss KLss RKss
            % shocks
			rho_S rho_A rho_R;

%----------------------------------------------------------------
% 2. Calibration
%----------------------------------------------------------------
alpha       = 0.3;              % Share of capital in production
beta    	= 0.99;			    % Discount factor savers
gammaf      = 0.8;              % Discount factor intermediary firms
gamma       = 0.779;            % Proba of keeping prices fixed (GK)
gammap      = 0.24;             % Price indexation (GK)        
delta       = 0.025;            % Capital depreciation
rho			= 0.80;				% Monetary policy rule, Smoothing, autocorr parameter
phi_p		= 1.5;				% Monetary policy rule, Inflation (=1.5)
phi_y		= 0.5;              % Monetary policy rule, GDP (=0.1)
h           = 0.7;              % Consumption habit (=0.7)
sigmaC	    = 1.5;		        % Risk aversion consumption
sigmaL	    = 2;		        % Labor disutility
veps        = 4.2;              % from final firms prices (GK)
kappa       = 1.8;              % investment adjustement cost

% Implied by steady state
LAMBDAss   = 1;
Pss        = 1;
Rss        = 1/beta;
Lss        = 1/3;
Qss        = 1;
Css        = LAMBDAss/((1-h)^(-sigmaC));
IKss       = 0.025;
KYss       = 0.2;              %(=0.149 in GK)
IYss       = IKss/KYss; 
KLss       = KYss^(1/(alpha-1));
Kss        = KLss * Lss;
Yss        = Kss^alpha*Lss^(1-alpha);
Iss        = Yss - Css;
PMss       = (veps-1)/veps;
RKss       = (alpha*Yss/Kss*PMss)/Qss + (1-delta);
Fss        = Yss*PMss/(1-beta*gamma);
F1ss       = Yss/(1-beta*gamma);
Wss        =(1-alpha)*Yss/Lss*Pss;
Bss        = Qss*Kss;
MUss       = 1;                %(perfect subsitutability => epsilon -> infty)
Dss        = 1;


% shock processes (autoregressive parameter)
rho_S   = 0.3; 
rho_A   = 0.8;
rho_R 	= 0.4;

%----------------------------------------------------------------
% 3. Model
%----------------------------------------------------------------
model(linear); 
    %% Household
    % Euler equation
    c = (h/(1+h))*c(-1) + (1/(1+h))*c(+1) - (1-h)/((1+h)*sigmaC)*(r-p(+1));
    % Labour supply equation 
    w = sigmaC*c - sigmaC*h*c(-1)+ sigmaL*l + p;
    

    %% Intermediary firms
    % Intermediary firm production function 
    ym = alpha*k(-1) + (1-alpha)*l + e_A;
    % Labour demand
    w = ym - l + pm;
	% Capital demand
    %r_k(+1) = (pm + ym - k)*(1-(1-delta)/RKss) + (1-delta)/RKss * q(+1) - q;
    r_k = (pm(-1) + ym(-1) - k(-1))*(1-(1-delta)/RKss) + (1-delta)/RKss * q - q(-1);

    % Asset market equilibrium
    k = b;

    %% Capital producers
    % Capital supply
    q = kappa*delta*(i - k(-1));
    % Capital accumulation
    k = delta*i + (1-delta)*k(-1);

    %% Final firms
    % Final firms prices
    Fss*f = Yss*PMss*(y+pm) + beta*gamma*LAMBDAss*Fss*(lambda(1)-gammap*veps*p+veps*p(1)+f(1));
    lambda(1) + r = 0;
    F1ss*f1 = Yss*y + beta*gamma*LAMBDAss*Pss^(gammap-1)*F1ss*(lambda(1)+gammap*(1-veps)*p-(1-veps)*p(1)+f1(1));
    pstar = f - f1 + p;
    p = gamma*gammap*p(-1) + (1-gamma)* pstar + e_S;
    uc = -1/((1-beta*h)*(1-h))*(c-h*c(-1)-beta*h*(c(1)-h*c));
    lambda = uc - uc(-1);
    Dss*d=gamma*Dss*Pss^(veps*(1-gammap))*(-gammap*veps*p(-1)+veps*p+d(-1))+veps/(1-gamma)*((1-gamma*Pss^((gamma-1)*(1-gammap)))/(1-gamma))^(-veps/(1-gamma))*gamma*Pss^((gamma-1)*(1-gammap))*((gamma-1)*p+gammap*(1-gamma)*p(-1));
    y = ym + d;

    %% Financial intermediary
    % Bank balance sheet evolution 
    %r(-1)*(1+b(-1))+ b(-1) - b = r(-1) * (1 + dep(-1))+ dep(-1) - dep;
    b=dep;

	% Monetary Policy Rule
	r_n = rho*r_n(-1) +  (1-rho)*(phi_p*p + phi_y*y ) + e_R  ;
    % Ressource constraint
    y = Css/Yss * c + Iss/Yss * i;
    % Fisher equation
    r_n = r + p(+1);
       
	% Supply shock
    e_S = rho_S * e_S(-1) + eta_S;
    % Monetary policy shock
    e_A = rho_A * e_A(-1)+ eta_A;
    % Supply shock
    e_R = rho_R* e_R(-1) + eta_R;
	
	end;

resid(1);
steady;
check;

%----------------------------------------------------------------
% 4. Shocks
%----------------------------------------------------------------
shocks;
var eta_S;  stderr 0.1;
var eta_A;  stderr 0.1;
var eta_R;  stderr 0.1;
end;


stoch_simul (order=1, irf=10);