% Dynare 
% FVGQRRUmodel
% case for Argentina


var Y K C H INV r epsilon_r logsigma_r logY mu w D lambda NX  r   ;  

varexo u_r u_sigma_r  ;

parameters delta alpha ffi nu rho_x phi phi_big_D Dfix beta
r_bar rho_r rho_sigma_r logsigmar_bar eta_r ;

%FVGQRRUmodel_parameters;

    delta = 0.014;    % Depreciation rate 
    alpha = 0.32;     % Capital share
    ffi   = 1000;     % Pameter determining the elasticity of labor to wages
    nu    = 5;        % Inverse of the elasticity of intertemporal substitution
    rho_x = 0.95;     % Autoregressive coefficient of the productivity process

% They estimate those for 4 countries
    phi       = 95;              % The adjustment cost of capital
    Dfix      = 4;               % Parameter that determines average value of debt
    phi_big_D =0.001;            % the holding costs of debt

% Argentina (posteriors)
    r_bar        = 0.0268;


    rho_r        = 0.97;
    rho_sigma_r  = 0.94;
    logsigmar_bar  = -5.71;   %exp(-5.71)=0.0033
    eta_r        = 0.046;



beta=0.980; 

%--------------------------------------------------------------------------
% 3. Model
%--------------------------------------------------------------------------



model;

    lambda=C^(-nu);                          

    lambda/(1+r)=lambda*phi_big_D*(D-Dfix)+beta*lambda(+1);         

    mu=beta*(alpha*lambda(+1)*Y(+1)/K+(1-delta)*mu(+1));

    ffi*H^ffi=w*lambda;

    w=(1-alpha)*Y/H;                                  

    lambda=mu*(1-phi/2*(INV/INV(-1)-1)^2-phi*INV/INV(-1)*(INV/INV(-1)-1))+
            beta*mu(+1)*phi*(INV(+1)/INV)^2*(INV(+1)/INV-1);
                                                            
    Y = C+INV+NX;                                               
    NX = D(-1)-D/(1+r)+phi_big_D/2*(D-Dfix)^2;                  
    K = (1-delta)*K(-1)+(1-phi/2*(INV/INV(-1)-1)^2)*INV;         
    Y = K(-1)^alpha*(H)^(1-alpha);                       

    % Stochastic process for the interest rate
    r=r_bar+epsilon_r;                               

                                                                              
    % Demeaned Country Spread
    epsilon_r=rho_r*epsilon_r(-1)+exp(logsigma_r)*u_r;                                     % u_r=N(0,1)                 
    logsigma_r=(1-rho_sigma_r)*logsigmar_bar+rho_sigma_r*logsigma_r(-1)+eta_r*u_sigma_r;   % u_sigma_r=N(0,1) 
                                                                             
    logY=log(Y);
end;

%--------------------------------------------------------------------------
% 4. Computation
%--------------------------------------------------------------------------



initval;


    Y      =		 1.89192;
    K      =		 31.8219;
    INV    =		 0.445507;
    H      =		 0.50124;
    C      =		 2.20082;
    mu     =		 0.206456;
    w      =		 2.56665;
    D      =		 -77.3417;
    lambda =		 0.206456;
    NX     =		 -0.754401;
    r           =		 r_bar;
    logY=log(Y);
    epsilon_r   =  0;
    logsigma_r  =  logsigmar_bar; 


end;

resid;
check;
steady;

shocks;
    var u_r        = 1; 
    var u_sigma_r  = 1;

end;



stoch_simul(order=3,  irf=40, nograph) ;





%--------------------------------------------------------------------------

ex_ = zeros(40,2);
ex_(1,2) = 1;
irfsetar = simult_(oo_.steady_state,oo_.dr,ex_,3);


irfY_etar         =(irfsetar(1,:) -oo_.steady_state(1))';
irfK_etar         =(irfsetar(2,:) -oo_.steady_state(2))';
irfC_etar         =(irfsetar(3,:) -oo_.steady_state(3))';
irfN_etar         =(irfsetar(4,:) -oo_.steady_state(4))';
irfI_etar         =(irfsetar(5,:) -oo_.steady_state(5))';
irfR_etar         =(irfsetar(6,:) -oo_.steady_state(6))';
irfepsr_etar      =(irfsetar(7,:) -oo_.steady_state(7))';
irflogsigmar_etar  =(irfsetar(8,:) -oo_.steady_state(8))';
irflogY_etar       =(irfsetar(9,:) -oo_.steady_state(9))';

%----------------------------------------------------


NPER = size(irfY_etar,1);
time=(0:1:NPER-1)';



for i=1:9
   subplot(3,3,i)
   set(gca,'FontSize',10)
   title('Impulse Responses')
end


subplot(3,3,1)
AA = plot(time,100*irfY_etar, '-');
title('Y');
grid;

subplot(3,3,2)
AA = plot(time,100*irfK_etar, '-');
title('K')
grid;

subplot(3,3,3)
AA = plot(time,100*irfC_etar, '-');
title('C')
grid;

subplot(3,3,4)
AA = plot(time,100*irfN_etar, '-');
title('N')
grid;


subplot(3,3,5)
AA = plot(time,100*irfI_etar, '-');
title('I')
grid;

subplot(3,3,6)
AA = plot(time,100*irfR_etar, '-');
title('R')
grid;

subplot(3,3,7)
AA = plot(time,(100*irfepsr_etar), '-');
title('epsilon r')
grid;


subplot(3,3,8)
AA = plot(time,(100*irflogsigmar_etar), '-');
title('logsigmar')
grid;

subplot(3,3,9)
AA = plot(time,100*irflogY_etar, '-');
title('logY');
grid;