%------------------------------------------------------------
close all;
%------------------------------------------------------------
% 1. Variables (15) 
%------------------------------------------------------------ 
var C_t0,
lambda1_t0,
lambda2_t0,
lambda3_t0,
K_t0,
U_t0,
mu_t0,
X_t0,
A1_t0,
A2_t0,
A1_z,
A2_z,
I_t0;



predetermined_variables K_t0, X_t0; 
%------------------------------------------------------------
% 2. Exogenous variables (2)
%------------------------------------------------------------
varexo eps_A1,
       eps_A2;   

%------------------------------------------------------------
% 3. Parameters (16)
%------------------------------------------------------------
parameters alpha 
           beta  
           delta 
           rhoa1 
           rhoa2
           omega 
           psai 
           eta  
           theta1 
           theta2 
           d2 
           d1 
           d0 
           gamma
           phic
           Mrow;

% set parameter values
alpha   =  0.36;          
beta    = 0.98267;       
delta   = 0.025;         
rhoa1   = 0.95;         
rhoa2   = 0.95;        
omega   = 4;              
psai    = 5;            
eta     = 0.9979;         
theta1  = 0.05607;        
theta2  = 2.8;
d2      =  1.4647E-08;      
d1      = -6.6722E-06;
d0      =  1.3950E-03;
gamma   = 1-0.696;        
phic    = 2;           
Mrow    = 4.7727;         
%------------------------------------------------------------
% 4. Model (13 equations)
%------------------------------------------------------------

model;

 
//1. foc_C =
 
1/C_t0^phic = lambda1_t0;
 
 
//2. foc_U =
 
(A1_t0*alpha*lambda1_t0*(K_t0*U_t0)^alpha*(mu_t0^theta2*theta1 - 1)*(d2*X_t0^2 + d1*X_t0 + d0 - 1))/U_t0 = K_t0*U_t0^(omega - 1)*lambda2_t0 + (A1_t0*K_t0*alpha*lambda3_t0*(K_t0*U_t0)^(alpha - 1)*(gamma - 1)*(mu_t0 - 1)*(d2*X_t0^2 + d1*X_t0 + d0 - 1))/(-A1_t0*(K_t0*U_t0)^alpha*(d2*X_t0^2 + d1*X_t0 + d0 - 1))^gamma;
 
 
//3. foc_mu =
 
 A1_t0*lambda1_t0*mu_t0^(theta2 - 1)*theta1*theta2*(K_t0*U_t0)^alpha*(d2*X_t0^2 + d1*X_t0 + d0 - 1) = lambda3_t0*(-A1_t0*(K_t0*U_t0)^alpha*(d2*X_t0^2 + d1*X_t0 + d0 - 1))^(1 - gamma);
 
 
//4. foc_I =
 
lambda1_t0=(A2_t0(+1)*I_t0(+1)^2*beta*lambda2_t0(+1)*psai*(I_t0(+1)/I_t0 - 1))/I_t0^2 - lambda2_t0*(A2_t0*((psai*(I_t0/I_t0(-1) - 1)^2)/2 - 1) + (A2_t0*I_t0*psai*(I_t0/I_t0(-1) - 1))/I_t0(-1));
 
 
//5. foc_K =
 
lambda2_t0=beta*((lambda2_t0(+1)*(omega - U_t0(+1)^omega))/omega + (A1_t0(+1)*alpha*lambda1_t0(+1)*(K_t0(+1)*U_t0(+1))^alpha*(mu_t0(+1)^theta2*theta1 - 1)*(d2*X_t0(+1)^2 + d1*X_t0(+1) + d0 - 1))/K_t0(+1) - (A1_t0(+1)*U_t0(+1)*alpha*lambda3_t0(+1)*(K_t0(+1)*U_t0(+1))^(alpha - 1)*(gamma - 1)*(mu_t0(+1) - 1)*(d2*X_t0(+1)^2 + d1*X_t0(+1) + d0 - 1))/(-A1_t0(+1)*(K_t0(+1)*U_t0(+1))^alpha*(d2*X_t0(+1)^2 + d1*X_t0(+1) + d0 - 1))^gamma);
 
 
 
//6. foc_X =
 
lambda3_t0=beta*(lambda3_t0(+1)*(eta - (A1_t0(+1)*(K_t0(+1)*U_t0(+1))^alpha*(gamma - 1)*(mu_t0(+1) - 1)*(d1 + 2*X_t0(+1)*d2))/(-A1_t0(+1)*(K_t0(+1)*U_t0(+1))^alpha*(d2*X_t0(+1)^2 + d1*X_t0(+1) + d0 - 1))^gamma) + A1_t0(+1)*lambda1_t0(+1)*(K_t0(+1)*U_t0(+1))^alpha*(mu_t0(+1)^theta2*theta1 - 1)*(d1 + 2*X_t0(+1)*d2));
 
 
//7. foc_lambda1 =
 
A1_t0*mu_t0^theta2*theta1*(K_t0*U_t0)^alpha*(d2*X_t0^2 + d1*X_t0 + d0 - 1)  - A1_t0*(K_t0*U_t0)^alpha*(d2*X_t0^2 + d1*X_t0 + d0 - 1) = C_t0 + I_t0;
 
 
//8. foc_lambda2 =
 
K_t0(+1)=- K_t0*(U_t0^omega/omega - 1) - A2_t0*I_t0*((psai*(I_t0/I_t0(-1) - 1)^2)/2 - 1);
 
 
//9. foc_lambda3 =
 
 X_t0(+1)= Mrow  + X_t0*eta - (mu_t0 - 1)*(-A1_t0*(K_t0*U_t0)^alpha*(d2*X_t0^2 + d1*X_t0 + d0 - 1))^(1 - gamma);
 

// 10. lognormal A1 shock
A1_t0 = exp(A1_z);

//11. Shock process
A1_z = rhoa1*A1_z(-1)+eps_A1;

// 12. lognormal A2 shock
A2_t0 = exp(A2_z);

//13. Shock process
A2_z = rhoa2*A2_z(-1)+eps_A2;



end;
 
initval;
C_t0=2.155941;
lambda1_t0=0.215142;
lambda2_t0=0.215142;
lambda3_t0=-0.002748;
K_t0=22.923560;
U_t0=1.000000;
I_t0=0.573089;
mu_t0=0.209381;
X_t0=3030.278626;
A1_t0=1;
A2_t0=1;
end; 




shocks;
var eps_A1 = 0.007;
var eps_A2 = 0.007;
end;

 

steady;
check;

stoch_simul(order=1,irf=20,nograph);

// Need to obtain the irfs for  Y_t and M_t for each shock:

// OUTPUT: Y_t= A1_t0*(K_t0*U_t0)^alpha*(d2*X_t0^2 + d1*X_t0 + d0 - 1);
//  M_t = (1-mu_t0)*(A1_t0*(K_t0*U_t0)^alpha*(d2*X_t0^2 + d1*X_t0 + d0 - 1))^(1-gamma);


//Constructing the irf for Output and M for the A1 (technology) shock:
Y1_irf_A1 = (K_t0_eps_A1 + oo_.steady_state(5).*ones(size(K_t0_eps_A1))).^(1-alpha);
Y2_irf_A1 = (U_t0_eps_A1 + oo_.steady_state(6).*ones(size(U_t0_eps_A1))).^(1-alpha);
Y3_irf_A1 = (A1_t0_eps_A1 + oo_.steady_state(9).*ones(size(A1_t0_eps_A1)));
Y4_irf_A1 = (ones(size(X_t0_eps_A1))-ones(size(X_t0_eps_A1)).*d0)-(d2*(X_t0_eps_A1 +  oo_.steady_state(8).*ones(size(X_t0_eps_A1)))+d1*(X_t0_eps_A1 +  oo_.steady_state(8).*ones(size(X_t0_eps_A1)))); 
Y_irf_A1 = Y1_irf_A1.*Y2_irf_A1.*Y3_irf_A1.*Y4_irf_A1;
M_irf_A1  = ones(size(mu_t0_eps_A1)) -(mu_t0_eps_A1 + oo_.steady_state(7).*ones(size(mu_t0_eps_A1))).*Y_irf_A1.^(1-gamma);

t_end = size(M_irf_A1,1);
T = [1:t_end];
X0 = zeros(t_end);
X1 = ones(t_end);
figure(1);
subplot(2,1,1);
plot(T,Y_irf_A1,'-k',T,X0,'-r','Linewidth',1.25);
axis auto;
title('Y');
subplot(2,1,2);
plot(T,M_irf_A1,'-k',T,X0,'-r','Linewidth',1.25); 
axis auto;
title('M');
 
//EOF