diff --git a/do_irf.m b/do_irf.m
new file mode 100644
index 0000000..142eef1
--- /dev/null
+++ b/do_irf.m
@@ -0,0 +1,32 @@
+
+figure(1)
+set(gcf,'Position',[500 300 900 600])
+
+cut = 2:60;
+
+subplot(231)
+plot(4*rn(cut))
+title('$r^n_t$','FontSize',14,'interpreter','latex')
+set(gcf,'Name','Responses to natural real rate shock')
+hold on
+
+subplot(232)
+plot(4*i(cut))
+title('$i_t$','FontSize',14,'interpreter','latex')
+hold on
+
+subplot(233)
+plot(4*pi(cut))
+title('$\pi_t$','FontSize',14,'interpreter','latex')
+hold on
+
+subplot(234)
+plot(4*rr(cut))
+title('$r_t$','FontSize',14,'interpreter','latex')
+hold on
+
+subplot(235)
+plot(y(cut))
+title('$y_t$','FontSize',14,'interpreter','latex')
+hold on
+
diff --git a/ocp.mod b/ocp.mod
new file mode 100644
index 0000000..6db6ba8
--- /dev/null
+++ b/ocp.mod
@@ -0,0 +1,90 @@
+var y, pi, i, rr, rn, p, phi1, phi2; 
+varexo e;
+parameters beta, sigma, kappa, lambda, elb, r, rho;
+
+% Program execution parameters
+
+% MODEL PARAMETERS 
+  beta    = 1/1.005;                % quarterly time discount factor 
+  sigma   = 2;                      % relative risk aversion 
+  kappa   = 0.024;                  % slope of the Phillips curve 
+  lambda  = 0.003;                  % weight on output gap in loss function
+  elb     = -0.0/4;                % effective lower bound
+  
+% EXOGENOUS SHOCK PROCESS 
+  r   = 100*(1/beta-1);             % steady-state (quarterly x 100)  
+  rho = 0.95;                       % persistence
+
+model;      
+y  = y(+1) - 1/sigma*(rr - rn);
+pi = beta*pi(+1) + kappa*y;
+phi1 = (kappa/sigma +1)/beta*phi1(-1) + kappa*phi2(-1) - kappa*pi - lambda*y; 
+phi2 = phi2(-1) + phi1(-1)/(beta*sigma) - pi;        
+0 = min(i-elb, phi1);           
+rr = i - pi(+1);
+pi = p - p(-1);
+rn = r + rho*(rn(-1)-r) + e;
+end;
+
+initval;
+rr = r;
+rn = r;
+i = r;
+phi1 = 0;
+phi2 = 0;
+end;
+
+steady; check;
+
+initval;
+rn = - r;
+phi1 = 0;
+end;
+
+endval;
+rn = r;
+end;
+
+perfect_foresight_setup(periods=200);
+perfect_foresight_solver(lmmcp);
+
+cut = [2 60];
+
+figure(1)
+set(gcf,'Name','Responses to natural real rate shock')
+
+subplot(231)
+plot(4*rn)
+title('$r^n_t$','FontSize',14,'interpreter','latex')
+hold on
+xlim([1 cut(2)])
+
+subplot(232)
+plot(4*i)
+title('$i_t$','FontSize',14,'interpreter','latex')
+hold on
+xlim(cut)
+
+subplot(233)
+plot(4*pi)
+title('$\pi_t$','FontSize',14,'interpreter','latex')
+hold on
+xlim(cut)
+
+subplot(234)
+plot(4*rr)
+title('$r_t$','FontSize',14,'interpreter','latex')
+hold on
+xlim(cut)
+
+subplot(235)
+plot(y)
+title('$y_t$','FontSize',14,'interpreter','latex')
+hold on
+xlim(cut)
+
+subplot(236)
+plot(p)
+title('$p_t$','FontSize',14,'interpreter','latex')
+hold on
+xlim(cut)
diff --git a/odp.mod b/odp.mod
new file mode 100644
index 0000000..0ffbff1
--- /dev/null
+++ b/odp.mod
@@ -0,0 +1,51 @@
+var y, pi, i, rr, rn; 
+varexo e;
+parameters beta, sigma, kappa, lambda, elb, r, rho;
+
+% Program execution parameters
+
+% MODEL PARAMETERS 
+  beta    = 1/1.005;                % quarterly time discount factor 
+  sigma   = 2;                      % relative risk aversion 
+  kappa   = 0.024;                  % slope of the Phillips curve 
+  lambda  = 0.003;                  % weight on output gap in loss function
+  elb     = 0/4;                    % effective lower bound  
+  
+% EXOGENOUS SHOCK PROCESS 
+  r   = 100*(1/beta-1);             % steady-state (quarterly x 100)  
+  rho = 0.85;                       % persistence
+
+model;      
+y  = y(+1) - 1/sigma*(rr - rn);
+pi = beta*pi(+1) + kappa*y;
+0 = min(i-elb, -(pi + lambda/kappa*y));           
+rr = i - pi(+1);
+rn = r + rho*(rn(-1)-r) + e;
+end;
+
+initval;
+rr = r;
+rn = r;
+i = r;
+end;
+
+steady; check;
+
+initval;
+rn = (-2/4-r)/rho + r;
+%rn = r;
+end;
+
+endval;
+rn = r;
+end;
+
+%shocks; 
+%var e;
+%periods  1:5  6:24  26:200 ; 
+%values   0   -0.05     0;
+%end;
+
+
+simul(periods=200, maxit=500, stack_solve_algo=0);
+do_irf;
\ No newline at end of file
diff --git a/ramsey_nk.mod b/ramsey_nk.mod
new file mode 100644
index 0000000..3556f25
--- /dev/null
+++ b/ramsey_nk.mod
@@ -0,0 +1,80 @@
+var y, pi, i, p, rr, rn; 
+varexo e;
+parameters beta, gamma, sigma, kappa, r, rho, lambda;
+
+% Program execution parameters
+
+% MODEL PARAMETERS 
+  beta    = 1/1.005;                % quarterly time discount factor 
+  gamma   = 0.9;                      % indexation to past inflation
+  sigma   = 2;                      % relative risk aversion 
+  kappa   = 0.024;                  % slope of the Phillips curve 
+  lambda  = 0.003;                  % weight on output gap in loss function
+  
+% EXOGENOUS SHOCK PROCESS: NATURAL REAL RATE
+  r   = 100*(1/beta-1);             % steady-state (quarterly x 100)  
+  rho = 0.95;                       % persistence
+
+model;      
+y  = y(+1) - 1/sigma*(rr - rn);
+pi-gamma*pi(-1) = beta*(pi(+1)-gamma*pi) + kappa*y;
+rr = i - pi(+1);
+pi = p - p(-1);
+rn = r + rho*(rn(-1)-r) + e;
+end;
+
+planner_objective pi^2 + lambda*y^2; 
+
+ramsey_model(planner_discount=1/1.005);
+
+initval;
+rn = -r;
+end;
+
+ramsey_constraints;
+i > 0;
+end;
+
+perfect_foresight_setup(periods=200);
+perfect_foresight_solver(lmmcp);
+
+cut = [2 60];
+
+figure(1)
+set(gcf,'Name','Responses to natural real rate shock')
+
+subplot(231)
+plot(4*rn)
+title('$r^n_t$','FontSize',14,'interpreter','latex')
+hold on
+xlim([1 cut(2)])
+
+subplot(232)
+plot(4*i)
+title('$i_t$','FontSize',14,'interpreter','latex')
+hold on
+xlim(cut)
+
+subplot(233)
+plot(4*pi)
+title('$\pi_t$','FontSize',14,'interpreter','latex')
+hold on
+xlim(cut)
+
+subplot(234)
+plot(4*rr)
+title('$r_t$','FontSize',14,'interpreter','latex')
+hold on
+xlim(cut)
+
+subplot(235)
+plot(y)
+title('$y_t$','FontSize',14,'interpreter','latex')
+hold on
+xlim(cut)
+
+subplot(236)
+plot(p)
+title('$p_t$','FontSize',14,'interpreter','latex')
+hold on
+xlim(cut)
diff --git a/tplt.mod b/tplt.mod
new file mode 100644
index 0000000..9e38e3e
--- /dev/null
+++ b/tplt.mod
@@ -0,0 +1,103 @@
+var        y, pi, i, rr, rn, p, dy, u; 
+varexo     eps, ups;
+parameters beta, gamma, sigma, lambda, kappa, epsilon, rho, r;
+
+beta    = 1/1.005;        % quarterly time discount factor 
+gamma   = 0;              % indexation to past inflation
+sigma   = 2;              % coefficient of relative risk aversion 
+kappa   = 0.024;          % slope of the Phillips curve 
+epsilon = 8;              % elasticity of substitution among varieties
+rho     = 0.95;           % shock persistence parameter
+r       = 100*(1/beta-1); % steady-state (net, quarterly, x100)  
+lambda  = kappa/epsilon;  % weight on output gap in CB loss function 
+
+model;      
+
+rn-r = rho*(rn(-1)-r) + eps;
+u = rho*u(-1) + ups;
+
+pi = p - p(-1);
+dy = y - y(-1);
+rr = i - pi(+1);
+
+y  = y(+1) - (rr-rn)/sigma;
+pi-gamma*pi(-1) = beta*(pi(+1)-gamma*pi) + kappa*y + u;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+[mcp = 'i > 0'] %% ALTERNATIVE RULES:
+%0 = (pi + y/epsilon);                                       % IT
+ 0 = (p + y/epsilon);                                        % PLT
+%  0 = (p + y/epsilon)*(p<0)     + (pi + y/epsilon)*(p>0);      % TPLT
+% 0 = (p + y/epsilon)*(rn<0.25) + (pi + y/epsilon)*(rn>0.25); % XTPLT
+% i = r + 800*p + 100*y;                                      % IRRRPL
+% i = r + 2*p + 0.25*y;                                       % IRRRPL
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+end;
+
+% planner_objective pi^2 + lambda*y^2; 
+% ramsey_model(planner_discount=1/1.005);
+
+initval;
+rn = -r;
+end;
+
+shocks; 
+var eps;
+periods   1;
+values    0;
+end;
+
+% ramsey_constraints;
+% i > 0;
+% end;
+
+perfect_foresight_setup(periods=200);
+perfect_foresight_solver(lmmcp);
+
+cut = [2 40];
+
+color = 'r--';
+
+figure(1)
+set(gcf,'Name','Responses to natural real rate shock')
+
+subplot(231)
+plot(4*rn-2,color)
+title('Natural real rate shock','FontSize',12)
+hold on
+xlim([1 cut(2)])
+
+subplot(232)
+plot(4*i,color)
+title('Nominal interest rate','FontSize',12)
+hold on
+xlim(cut)
+
+subplot(233)
+plot(4*pi+2,color)
+title('Inflation','FontSize',12)
+hold on
+xlim(cut)
+
+subplot(234)
+plot(4*rr-2,color)
+title('Real interest rate','FontSize',12)
+hold on
+xlim(cut)
+xlabel('Quarters')
+
+subplot(235)
+plot(y,color)
+title('Output gap','FontSize',12)
+hold on
+xlim(cut)
+xlabel('Quarters')
+
+subplot(236)
+trend = [1:(0.02/4):2.005]';
+plot(p+trend,color)
+title('Price level','FontSize',12)
+hold on
+xlim([1 cut(2)])
+xlabel('Quarters')
\ No newline at end of file
diff --git a/ttr.mod b/ttr.mod
new file mode 100644
index 0000000..473c23f
--- /dev/null
+++ b/ttr.mod
@@ -0,0 +1,52 @@
+var y, pi, i, rr, rn, u; 
+varexo e, eps;
+parameters beta, sigma, kappa, elb, r, rho, phi_pi, phi_y;
+
+% Program execution parameters
+
+% MODEL PARAMETERS 
+  beta    = 1/1.005;                % quarterly time discount factor 
+  sigma   = 2;                      % relative risk aversion 
+  kappa   = 0.024;                  % slope of the Phillips curve 
+  elb     = 0/4;                    % effective lower bound  
+  phi_pi  = 1.5;                    % reaction to inflation
+  phi_y   = 0.5;                    % reaction to output gap
+  
+% EXOGENOUS SHOCK PROCESS: NATURAL REAL RATE
+  r   = 100*(1/beta-1);             % steady-state (quarterly x 100)  
+  rho = 0.95;                       % persistence
+
+model;      
+y  = y(+1) - 1/sigma*(rr - rn);
+pi = beta*pi(+1) + kappa*y;
+i  = max(elb, r + phi_pi*pi + phi_y*y + u);                         
+rr = i - pi(+1);
+rn = r + rho*(rn(-1)-r) + e;
+u  = rho*u(-1) + eps;
+end;
+
+initval;
+rr = r;
+rn = r;
+i = r;
+end;
+
+steady; check;
+
+initval;
+rn = (-2/4-r)/rho + r;
+%rn = r;
+end;
+
+endval;
+rn = r;
+end;
+
+%shocks; 
+%var e;
+%periods  1:1  2:10  11:200 ; 
+%values   0   -0.06     0;
+%end;
+
+simul(periods=200, maxit=500, stack_solve_algo=0);
+do_irf;
\ No newline at end of file
diff --git a/ttrp.mod b/ttrp.mod
new file mode 100644
index 0000000..70d61da
--- /dev/null
+++ b/ttrp.mod
@@ -0,0 +1,97 @@
+var y, pi, i, rr, rn, tr, p, phi_pi, phi_p; 
+varexo e;
+parameters beta, sigma, kappa, elb, r, rho, phi_y;
+
+% Model Parameters
+  beta    = 1/1.005;             % quarterly time discount factor 
+  sigma   = 2;                   % relative risk aversion 
+  kappa   = 0.024;               % slope of the Phillips curve 
+  phi_y   = 0.5;                 % reaction to output gap
+  elb     = 0/4;                 % effective lower bound    
+
+% Natural Real Interest Rate Process 
+  r   = 100*(1/beta-1);          % steady-state (quarterly x 100)  
+  rho = 0.95;                    % persistence
+
+model;      
+
+y  = y(+1) - 1/sigma*(rr - rn);
+pi = beta*pi(+1) + kappa*y;
+i  = max(elb, tr);                         
+tr = r + phi_pi*pi + phi_p*p + phi_y*y;
+
+pi = p - p(-1);
+rr = i - pi(+1);
+
+rn = r + rho*(rn(-1)-r) + e;
+
+% phi_pi = 1.5;
+% phi_p  = 0;
+phi_pi = 1.5*(rn>0.25);
+phi_p  = 1.5*(rn<0.25);
+
+end;
+
+initval;
+phi_p  = 1.5;
+phi_pi = 1.5;
+rr = r;
+rn = r;
+tr = r;
+i  = r;
+p  = 0;
+pi = 0;
+end;
+
+steady; check;
+
+initval;
+rn = -r;
+end;
+
+endval;
+rn = r;
+end;
+
+simul(periods=200, maxit=500, stack_solve_algo=0);
+
+figure(1)
+cut = [2 60];
+
+subplot(231)
+plot(4*rn)
+title('$r^n_t$','FontSize',14,'interpreter','latex')
+set(gcf,'Name','Responses to natural real rate shock')
+hold on
+xlim(cut)
+
+subplot(232)
+plot(4*i)
+title('$i_t$','FontSize',14,'interpreter','latex')
+hold on
+xlim(cut)
+
+subplot(233)
+plot(4*pi)
+title('$\pi_t$','FontSize',14,'interpreter','latex')
+hold on
+xlim(cut)
+
+subplot(234)
+plot(4*rr)
+title('$r_t$','FontSize',14,'interpreter','latex')
+hold on
+xlim(cut)
+
+subplot(235)
+plot(y)
+title('$y_t$','FontSize',14,'interpreter','latex')
+hold on
+xlim(cut)
+
+subplot(236)
+plot(p)
+title('$p_t$','FontSize',14,'interpreter','latex')
+hold on
+xlim(cut)
+