diff --git a/calvo_tr_linear.mod b/calvo_tr_linear.mod new file mode 100644 index 0000000..fe3ee7f --- /dev/null +++ b/calvo_tr_linear.mod @@ -0,0 +1,78 @@ +var c n y r pi yhat ghat g beta; +varexo ub ug; + +parameters alpha gamma chi varphi epsilon betass rhog rhob gbar phipi phig; + +gamma = 2; +chi = 11; +varphi = 1; +epsilon = 7; +betass = (1/1.04)^(1/12); +rhob = 0.95; +rhog = 0.9; +alpha = 0.9; +gbar = 0.07; +phipi = 2; +phig = 0; + +model; +# nu_u = gamma*STEADY_STATE(y)/STEADY_STATE(c); +# nu_v = varphi*STEADY_STATE(y)/STEADY_STATE(n); +# capgamma = nu_u/(nu_u + nu_v); +# kappa = (1-alpha)*(1-alpha*betass)/alpha*(nu_u + nu_v); +y = n; +y = c + g; +yhat-ghat = yhat(+1)-ghat(+1) - 1/nu_u*(r - pi(+1) + log(beta(+1))); +yhat = log(y/STEADY_STATE(y)); +ghat = (g-STEADY_STATE(g))/STEADY_STATE(y); +pi = beta(+1)*pi(+1) + kappa*(yhat-capgamma*ghat); +r = max(-log(betass) + phipi*pi + phig*log(g/gbar),0); +g/gbar = (g(-1)/gbar)^rhog*exp(ug); +beta/betass = (beta(-1)/betass)^rhob*exp(ub); +end; + +shocks; +var ug; stderr .01; +var ub; stderr .01; +end; + +yss = fsolve(@(y) epsilon/(epsilon-1)*chi*y^varphi*(y-gbar)^gamma - 1,0.5); + +initval; +beta = betass; +g = gbar; +r = -log(beta); +y = yss; +c = y - g; +n = y; +end; + +steady; + +shocks; +var ub; +periods 1:1; +values 0.005; +end; + +simul(periods=200, maxit=500, stack_solve_algo=0); + +figure(1) +subplot(311) +plot(1200*r); +title('r') +hold on +xlim([0 50]) +subplot(312) +plot(y/yss-1); +title('y/yss'); +hold on +xlim([0 50]) +subplot(313) +plot(pi); +title('pi'); +hold on +xlim([0 50]) + + + diff --git a/calvo_tr_linear_stoch.mod b/calvo_tr_linear_stoch.mod new file mode 100644 index 0000000..b6fd0f1 --- /dev/null +++ b/calvo_tr_linear_stoch.mod @@ -0,0 +1,88 @@ +var c n y r pi yhat ghat g beta; +varexo ub ug; + +parameters alpha gamma chi varphi epsilon betass rhog rhob gbar phipi phig; + +@# define stochastic = 0 + +gamma = 2; +chi = 11; +varphi = 1; +epsilon = 7; +betass = (1/1.04)^(1/12); +rhob = 0.95; +rhog = 0.9; +alpha = 0.9; +gbar = 0.07; +phipi = 2; +phig = 0; + +model; +# nu_u = gamma*STEADY_STATE(y)/STEADY_STATE(c); +# nu_v = varphi*STEADY_STATE(y)/STEADY_STATE(n); +# capgamma = nu_u/(nu_u + nu_v); +# kappa = (1-alpha)*(1-alpha*betass)/alpha*(nu_u + nu_v); +y = n; +y = c + g; +yhat-ghat = yhat(+1)-ghat(+1) - 1/nu_u*(r - pi(+1) + log(beta(+1))); +yhat = log(y/STEADY_STATE(y)); +ghat = (g-STEADY_STATE(g))/STEADY_STATE(y); +pi = beta(+1)*pi(+1) + kappa*(yhat-capgamma*ghat); +r = max(-log(betass) + phipi*pi + phig*log(g/gbar),0); +%r = -log(betass) + phipi*pi + phig*log(g/gbar); +g/gbar = (g(-1)/gbar)^rhog*exp(ug); +beta/betass = (beta(-1)/betass)^rhob*exp(ub); +end; + +shocks; +var ug; stderr .01; +var ub; stderr .01; +end; + +yss = fsolve(@(y) epsilon/(epsilon-1)*chi*y^varphi*(y-gbar)^gamma - 1, 0.5); + +initval; +beta = betass; +g = gbar; +r = -log(beta); +y = yss; +c = y - g; +n = y; +end; + +steady; + +@# if stochastic==1 + shocks; + var ub; stderr .01; + var ug; stderr .01; + end; + stoch_simul(order=1,nograph,noprint,nomoments,nocorr,nofunctions); + irf_calvo_ub +@# else + shocks; + var ub; + periods 1:1 2:200; + values 0.01 0; + end; + simul(periods=200, maxit=500, stack_solve_algo=0); + figure(1) + subplot(311) + plot(1200*r); + title('r') + hold on + xlim([0 50]) + subplot(312) + plot(y/yss-1); + title('y/yss'); + hold on + xlim([0 50]) + subplot(313) + plot(1200*pi); + title('pi'); + hold on + xlim([0 50]) +@# endif + + + diff --git a/calvo_tr_nl.mod b/calvo_tr_nl.mod new file mode 100644 index 0000000..9ed4009 --- /dev/null +++ b/calvo_tr_nl.mod @@ -0,0 +1,82 @@ +var c n y r pi pstar w Delta K F; +var beta g; +varexo ub ug; + +parameters gamma betass chi varphi epsilon rhob rhog lambda gbar phir phipi phig; + +gamma = 2; +chi = 11; +varphi = 1; +epsilon = 7; +betass = (1/1.04)^(1/12); +rhob = 0.95; +rhog = 0.9; +lambda = 1/10; +gbar = 0.07; +phipi = 2; +phig = 0; +phir = 0; + +model; +c^(-gamma)*w = chi*n^varphi; +1 = r*beta(+1)*(c(+1)/c)^(-gamma)/pi(+1); +n = y*Delta; +y = c + g; +1 = lambda*pstar^(1-epsilon) + (1-lambda)*pi^(epsilon-1); +Delta = lambda*pstar^(-epsilon) + (1-lambda)*pi^epsilon*Delta(-1); +pstar = K/F; +K = c^(-gamma)*epsilon/(epsilon-1)*w*y + beta(+1)*(1-lambda)*pi(+1)^epsilon*K(+1); +F = c^(-gamma)*y + beta(+1)*(1-lambda)*pi(+1)^(epsilon-1)*F(+1); +r = max((1/betass)*(r(-1)/(1/betass))^phir*((pi/1)^phipi*(g/gbar)^phig)^(1-phir),1); +g/gbar = (g(-1)/gbar)^rhog*exp(ug); +beta/betass = (beta(-1)/betass)^rhob*exp(ub); +end; + + +yss = fsolve(@(y) epsilon/(epsilon-1)*chi*y^varphi*(y-gbar)^gamma - 1, .3); + +initval; +y = yss; +g = gbar; +c = y - g; +n = y; +beta = betass; +r = 1/beta; +pi = 1; +pstar = 1; +w = chi*n^varphi*c^gamma; +Delta = 1; +K = c^(-gamma)*epsilon/(epsilon-1)*w*y/(1-beta*lambda); +F = c^(-gamma)*y/(1-beta*lambda); +end; + +steady; + +shocks; +var ub; +periods 1:1 ; +values 0.005; +end; + +simul(periods=200, maxit=500, stack_solve_algo=0); + +figure(1) +subplot(311) +plot(1200*(r-1)); +title('r') +hold on +xlim([0 50]) +subplot(312) +plot(y/yss-1); +title('y/yss'); +hold on +xlim([0 50]) +subplot(313) +plot(pi-1); +title('pi'); +hold on +xlim([0 50]) + + + + diff --git a/calvo_tr_nl_stoch.mod b/calvo_tr_nl_stoch.mod new file mode 100644 index 0000000..75bfc1f --- /dev/null +++ b/calvo_tr_nl_stoch.mod @@ -0,0 +1,90 @@ +var c n y r pi pstar w Delta K F; +var beta g; +varexo ub ug; + +parameters gamma betass chi varphi epsilon rhob rhog lambda gbar phir phipi phig; + +@# define stochastic = 0 + +gamma = 2; +chi = 11; +varphi = 1; +epsilon = 7; +betass = (1/1.04)^(1/12); +rhob = 0.95; +rhog = 0.9; +lambda = 1/10; +gbar = 0.07; +phipi = 2; +phig = 0; +phir = 0; + +model; +c^(-gamma)*w = chi*n^varphi; +1 = r*beta(+1)*(c(+1)/c)^(-gamma)/pi(+1); +n = y*Delta; +y = c + g; +1 = lambda*pstar^(1-epsilon) + (1-lambda)*pi^(epsilon-1); +Delta = lambda*pstar^(-epsilon) + (1-lambda)*pi^epsilon*Delta(-1); +pstar = K/F; +K = c^(-gamma)*epsilon/(epsilon-1)*w*y + beta(+1)*(1-lambda)*pi(+1)^epsilon*K(+1); +F = c^(-gamma)*y + beta(+1)*(1-lambda)*pi(+1)^(epsilon-1)*F(+1); +r = max((1/betass)*(r(-1)/(1/betass))^phir*((pi/1)^phipi*(g/gbar)^phig)^(1-phir),1); +%r = ((1/betass)*(r(-1)/(1/betass))^phir*((pi/1)^phipi*(g/gbar)^phig)^(1-phir)); +g/gbar = (g(-1)/gbar)^rhog*exp(ug); +beta/betass = (beta(-1)/betass)^rhob*exp(ub); +end; + + +yss = fsolve(@(y) epsilon/(epsilon-1)*chi*y^varphi*(y-gbar)^gamma - 1, .3); +initval; +y = yss; +g = gbar; +c = y - g; +n = y; +beta = betass; +r = 1/beta; +pi = 1; +pstar = 1; +w = chi*n^varphi*c^gamma; +Delta = 1; +K = c^(-gamma)*epsilon/(epsilon-1)*w*y/(1-beta*lambda); +F = c^(-gamma)*y/(1-beta*lambda); +end; + +steady; + +@# if stochastic==1 + shocks; + var ub; stderr .01; + var ug; stderr .01; + end; + stoch_simul(order=1,nograph,noprint,nomoments,nocorr,nofunctions); + irf_calvo_ub +@# else + shocks; + var ub; + periods 1:1 2:200; + values 0.01 0; + end; + simul(periods=200, maxit=500, stack_solve_algo=0); + figure(1) + subplot(311) + plot(1200*(r-1)); + title('r') + hold on + xlim([0 50]) + subplot(312) + plot(y/yss-1); + title('y/yss'); + hold on + xlim([0 50]) + subplot(313) + plot(1200*(pi-1)); + title('pi'); + hold on + xlim([0 50]) +@# endif + + + diff --git a/do_cplt.m b/do_cplt.m new file mode 100644 index 0000000..6cabda8 --- /dev/null +++ b/do_cplt.m @@ -0,0 +1,70 @@ +% CONSTANT PRICE LEVEL TARGETING RULE +% (C) Anton Nakov + +% MODEL PARAMETERS (quarterly frequency) + beta = 1/1.005; % quarterly time discount factor + sigma = 2; % relative risk aversion + kappa = 0.024; % slope of the Phillips curve + lx = 0.003; % inflation weight in loss f-n + +% EXOGENOUS SHOCK PROCESS (natural real interest rate) + rnstst = 100*(1/beta-1); % steady-state + stdrn = 0.003; % standard deviation + rho = 0.85; % persistence + vare = stdrn^2*(1-rho^2); % variance of the innovation + +% DECLARE MODEL FUNCTION + model.func = 'func_cplt'; + +% DEFINE APPROXIMATION SPACE + n = [49 49]; % number of grid points + smin = [-2 -0.05]/4; % minimum states (quarterly) + smax = [+4 +0.05]/4; % maximum states (quarterly) + fspace = fundefn('lin',n,smin,smax); % function space + scoord = funnode(fspace); % state collocation grid coordinates + snodes = gridmake(scoord); % state collocation grid points + +% SET OPTIONS + optset('remsolve','nres',1); + optset('arbit','lcpmethod','minmax'); + +% INITIALIZE + nn = length(snodes); + xinit = [zeros(nn,2) snodes(:,1)]; + hinit = zeros(nn,2); + +% GAUSSIAN QUADRATURE + [e,w] = qnwnorm(3,0,vare); + model.e = e; % shocks + model.w = w; % probabilities + +% SOLVE RATIONAL EXPECTATIONS EQULIBRIUM + model.params = {sigma,rho,kappa,lx,beta,rnstst}; + [c,s,xx,p,f,resid] = remsolve(model,fspace,scoord,xinit); + +% HOMOTOPY + for rho = [0.85]; + vare = stdrn^2*(1-rho^2); % variance of the innovation + [e,w] = qnwnorm(3,0,vare); % (number of grid points; mean; variance) + model.e = e; % shocks + model.w = w; + model.params = {sigma,rho,kappa,lx,beta,rnstst}; + xinit = reshape(xx,nn,3); + hinit = reshape(p,nn,2); + [c,s,xx,p,f,resid] = remsolve(model,fspace,scoord,xinit,hinit); + end + +%% +% SIMULATE LIQUIDITY TRAP + + init = [-2/4 0]; % initial state + nper = 40; + [ssim,xxsim] = simultrap(model,init,nper,scoord,xx); + + pi_sim = 4*squeeze(xxsim(:,1,:)); + x_sim = 4*squeeze(xxsim(:,2,:)); + i_sim = 4*squeeze(xxsim(:,3,:)); + rn_sim = 4*squeeze(ssim(:,1,:)); + +% PLOT SIMULATED LIQUIDITY TRAP + plot_paths \ No newline at end of file diff --git a/do_ocp.m b/do_ocp.m new file mode 100644 index 0000000..db871b8 --- /dev/null +++ b/do_ocp.m @@ -0,0 +1,82 @@ +% OPTIMAL COMMITMENT POLICY WITH ZERO FLOOR ON THE NOMINAL INTEREST RATE +% (C) Anton Nakov + +% MODEL PARAMETERS (quarterly frequency) + beta = 1/1.005; % quarterly time discount factor + sigma = 2; % relative risk aversion + kappa = 0.024; % slope of the Phillips curve + lx = 0.003; % inflation weight in loss f-n + +% EXOGENOUS SHOCK PROCESS + rnstst = 100*(1/beta-1); % steady-state (quarterly x 100) + stdrn = 0.5; % standard deviation + rho = 0.65; % persistence + vare = stdrn^2*(1-rho^2); % variance of the innovation + +% DECLARE MODEL FUNCTION + model.func = 'func_ocp'; + +% DEFINE APPROXIMATION SPACE + n = [25 13 5]; % number of grid points + %smin = [-2 0.000 -0.03]/4; % minimum states (quarterly) + %smax = [+4 +0.008 +0.01]/4; % maximum states (quarterly) + smin = [-2 0.000 -0.05]/4; % minimum states (quarterly) + smax = [+4 +0.015 +0.05]/4; % maximum states (quarterly) + + fspace = fundefn('lin',n,smin,smax); % function space + scoord = funnode(fspace); % state collocation grid coordinates + snodes = gridmake(scoord); % state collocation grid points + +% SET OPTIONS + optset('remsolve','nres',1); + optset('arbit','lcpmethod','minmax'); + +% INITIALIZE + nn = length(snodes); + xinit = [zeros(nn,2) max(0,snodes(:,1))]; % [inflation; output gap; nominal interest rate] + hinit = zeros(nn,2); + +% GAUSSIAN QUADRATURE + [e,w] = qnwnorm(3,0,vare); % (number of grid points; mean; variance) + model.e = e; % shocks + model.w = w; % probabilities + +% SOLVE RATIONAL EXPECTATIONS EQULIBRIUM + model.params = {sigma,0.35,kappa,lx,beta,rnstst}; + [c,s,xx,p,f,resid] = remsolve(model,fspace,scoord,xinit,hinit); + +% HOMOTOPY + for rho = [rho] + vare = stdrn^2*(1-rho^2); % variance of the innovation + [e,w] = qnwnorm(3,0,vare); % (number of grid points; mean; variance) + model.e = e; % shocks + model.w = w; + model.params = {sigma,rho,kappa,lx,beta,rnstst}; + xinit = reshape(xx,nn,3); + hinit = reshape(p,nn,2); + [c,s,xx,p,f,resid] = remsolve(model,fspace,scoord,xinit,hinit); + end + +%% +% SIMULATE LIQUIDITY TRAP + +% GET ERGODIC DISTRIBUTION OF ENDOGENOUS STATE + init = [-3/4 0 0]; % initial state (exog. and endogenous) + nper = 40; + + [ssim,xxsim] = simultrap(model,init,nper,scoord,xx); + + pi_sim = 4*(squeeze(xxsim(:,1,:))); + x_sim = 4*(squeeze(xxsim(:,2,:))); + i_sim = 4*(squeeze(xxsim(:,3,:))); + rn_sim = 4*(squeeze(ssim(:,1,:))); + +% PLOT SIMULATED LIQUIDITY TRAP + plot_paths; + +% figure(2) +% surf(4*scoord{2},4*scoord{1},squeeze(4*xx(:,:,1,3))) +% xlabel('\lambda') +% ylabel('r^n') +% zlabel('i') + \ No newline at end of file diff --git a/do_odp.m b/do_odp.m new file mode 100644 index 0000000..bf1d5b8 --- /dev/null +++ b/do_odp.m @@ -0,0 +1,69 @@ +% OPTIMAL DISCRETIONARY POLICY WITH ZERO FLOOR ON THE NOMINAL INTEREST RATE +% (C) Anton Nakov + +% MODEL PARAMETERS (quarterly frequency) + beta = 1/1.005; % quarterly time discount factor + sigma = 2; % relative risk aversion + kappa = 0.024; % slope of the Phillips curve + lx = 0.003; % inflation weight in loss f-n + +% EXOGENOUS SHOCK PROCESS + rnstst = 100*(1/beta-1); % steady-state (quarterly x 100) + stdrn = 0.5; % standard deviation + rho = 0.65; % persistence + vare = stdrn^2*(1-rho^2); % variance of the innovation + +% DECLARE MODEL FUNCTION + model.func = 'func_odp'; + +% DEFINE APPROXIMATION SPACE + n = 101; % number of grid points + smin = -2/4; % minimum natural real interest rate + smax = +4/4; % maximum natural real interest rate + fspace = fundefn('lin',n,smin,smax); % function space + scoord = funnode(fspace); % state collocation grid coordinates + snodes = gridmake(scoord); % state collocation grid points + +% SET OPTIONS + optset('remsolve','nres',1); + optset('arbit','lcpmethod','minmax'); + +% INITIALIZE + nn = length(snodes); + xinit = [zeros(nn,1) zeros(nn,1) zeros(nn,1)+scoord]; % [inflation; output gap; nominal interest rate] + hinit = zeros(nn,2); + +% GAUSSIAN QUADRATURE + [e,w] = qnwnorm(7,0,vare); % (number of grid points; mean; variance) + model.e = e; % shocks + model.w = w; % probabilities + +% SOLVE RATIONAL EXPECTATIONS EQULIBRIUM + model.params = {sigma,rho,kappa,lx,beta,rnstst}; + [c,s,xx,p,f,resid] = remsolve(model,fspace,scoord,xinit,hinit); + +% HOMOTOPY +for rho = [] + vare = stdrn^2*(1-rho^2); % variance of the innovation + [e,w] = qnwnorm(3,0,vare); % (number of grid points; mean; variance) + model.e = e; + model.w = w; + model.params = {sigma,rho,kappa,lx,beta,rnstst}; + xinit = reshape(xx,nn,3); + hinit = reshape(p,nn,2); + [c,s,xx,p,f,resid] = remsolve(model,fspace,scoord,xinit,hinit); +end + +%% +% SIMULATE LIQUIDITY TRAP + init = -3/4; % initial state (exog.) + nper = 40; + [ssim,xxsim] = simultrap(model,init,nper,scoord,xx); + + pi_sim = 4*(squeeze(xxsim(:,1,:))); + x_sim = 4*(squeeze(xxsim(:,2,:))); + i_sim = 4*(squeeze(xxsim(:,3,:))); + rn_sim = 4*(squeeze(ssim(:,1,:))); + +% PLOT SIMULATED LIQUIDITY TRAP + plot_paths \ No newline at end of file diff --git a/do_ttr.m b/do_ttr.m new file mode 100644 index 0000000..49e6480 --- /dev/null +++ b/do_ttr.m @@ -0,0 +1,69 @@ +% TRUNCATED TAYLOR RULE +% (C) Anton Nakov + +% MODEL PARAMETERS + beta = 1/1.005; % quarterly time discount factor + sigma = 3; % relative risk aversion + kappa = 0.024; % slope of the Phillips curve + elb = 0/4; % effective lower bound + +% EXOGENOUS SHOCK PROCESS + rnstst = 100*(1/beta-1); % steady-state (quarterly x 100) + stdrn = 0.5; % standard deviation + rho = 0.65; % persistence + vare = stdrn^2*(1-rho^2); % variance of the innovations + +% DECLARE MODEL FUNCTION + model.func = 'func_ttr'; + +% DEFINE APPROXIMATION SPACE + n = 101; % number of grid points + smin = -2/4; % minimum states (quarterly x 100) + smax = +4/4; % maximum states (quarterly x 100) + fspace = fundefn('lin',n,smin,smax); % function space + scoord = funnode(fspace); % state collocation grid coordinates + snodes = gridmake(scoord); % state collocation grid points + +% SET OPTIONS + optset('remsolve','nres',1); + optset('arbit','lcpmethod','minmax'); + +% INITIALIZE + nn = length(snodes); + xinit = zeros(nn,2); + hinit = zeros(nn,2); + +% GAUSSIAN QUADRATURE + [e,w] = qnwnorm(3,0,vare); + model.e = e; % shocks + model.w = w; % probabilities + +% SOLVE RATIONAL EXPECTATIONS EQULIBRIUM + model.params = {sigma,rho,kappa,beta,rnstst,elb}; + [c,s,xx,p,f,resid] = remsolve(model,fspace,scoord,xinit,hinit); + +% HOMOTOPY +for rho = []; + vare = stdrn^2*(1-rho^2); % variance of the innovation + [e,w] = qnwnorm(3,0,vare); % GAUSSIAN QUADRATURE + model.e = e; % shocks + model.w = w; % probabilities + model.params = {sigma,rho,kappa,beta,rnstst,elb}; + xinit = reshape(xx,nn,2); + hinit = reshape(p,nn,2); + [c,s,xx,p,f,resid] = remsolve(model,fspace,scoord,xinit,hinit); +end + +%% +% SIMULATE LIQUIDITY TRAP + init = -3/4; % initial state (exog.) + nper = 40; + [ssim,xxsim] = simultrap(model,init,nper,scoord,xx); + + pi_sim = 4*(squeeze(xxsim(:,1,:))); + x_sim = 4*(squeeze(xxsim(:,2,:))); + i_sim = taylor(4*rnstst, pi_sim, x_sim,4*elb); + rn_sim = 4*(squeeze(ssim(:,1,:))); + +% PLOT SIMULATED LIQUIDITY TRAP + plot_paths \ No newline at end of file diff --git a/func_cplt.m b/func_cplt.m new file mode 100644 index 0000000..158e863 --- /dev/null +++ b/func_cplt.m @@ -0,0 +1,95 @@ +% Function file for CONSTANT PRICE LEVEL TARGETING RULE +% Program for the article "Optimal and Simple Monetary Policy Rules +% with Zero Floor on the Nominal Interest Rate", +% International Journal of Central Banking, Vol. 4(2), pages 73-127, June +% (C) Anton Nakov +function [out1,out2,out3] = func_cplt(flag,s,xx,ep,e,sigma,rho,kappa,lx,beta,rnstst) + +rn = s(:,1); +p_lag = s(:,2); + +n = length(rn); + +if ~isempty(ep) + E_pi = ep(:,1); + E_x = ep(:,2); +end +if ~isempty(xx) + pi = xx(:,1); + x = xx(:,2); + i = xx(:,3); + p = p_lag + pi/100; +end + +switch flag + +case 'b'; % BOUND FUNCTION + + out1(:,1) = -inf*ones(n,1); % pi low + out1(:,2) = -inf*ones(n,1); % x low + out1(:,3) = zeros(n,1); % i low + + out2(:,1) = inf*ones(n,1); % pi high + out2(:,2) = inf*ones(n,1); % x high + out2(:,3) = inf*ones(n,1); % i high + +case 'f'; % EQUILIBRIUM FUNCTION + + out1(:,1) = pi - beta*E_pi - kappa*x; % f1 (NKPC) + out1(:,2) = x - E_x + 1/sigma*(i - E_pi - rn); % f2 (NIS) + out1(:,3) = i.*(p + lx/kappa*x); % f3 (ZB) + + out2(:,1,1) = ones(n,1); % f1pi + out2(:,1,2) = -kappa*ones(n,1); % f1x + out2(:,1,3) = zeros(n,1); % f1i + + out2(:,2,1) = zeros(n,1); % f2pi + out2(:,2,2) = ones(n,1); % f2x + out2(:,2,3) = 1/sigma*ones(n,1); % f2i + + out2(:,3,1) = i/100; % f3pi p = p_lag + pi/100; + out2(:,3,2) = lx/kappa*i; % f3x + out2(:,3,3) = p + lx/kappa*x; % f3i + + out3(:,1,1) = -beta*ones(n,1); % f1E_pi + out3(:,1,2) = zeros(n,1); % f1E_x + + out3(:,2,1) = -ones(n,1)/sigma; % f2E_pi + out3(:,2,2) = -ones(n,1); % f2E_x + + out3(:,3,1) = zeros(n,1); % f3E_pi + out3(:,3,2) = zeros(n,1); % f3E_x + +case 'g'; % STATE TRANSITION FUNCTION + + out1(:,1) = rho*rn + (1-rho)*rnstst + e; % g1 + out1(:,2) = p; % g2 + + out2(:,1,1) = zeros(n,1); % g1pi + out2(:,1,2) = zeros(n,1); % g1x + out2(:,1,3) = zeros(n,1); % g1i + + out2(:,2,1) = ones(n,1)/100; % g2pi + out2(:,2,2) = zeros(n,1); % g2x + out2(:,2,3) = zeros(n,1); % g2i + +case 'h'; % EXPECTATION FUNCTION + + out1(:,1) = pi; % E_pi + out1(:,2) = x; % E_x + + out2(:,1,1) = ones(n,1); % E_pi_pi + out2(:,1,2) = zeros(n,1); % E_pi_x + out2(:,1,3) = zeros(n,1); % E_pi_i + + out2(:,2,1) = zeros(n,1); % E_x_pi + out2(:,2,2) = ones(n,1); % E_x_x + out2(:,2,3) = zeros(n,1); % E_x_i + + out3(:,1,1) = zeros(n,1); % E_pi_rn + out3(:,2,1) = zeros(n,1); % E_x_rn + + out3(:,1,2) = zeros(n,1); % E_pi_p + out3(:,2,2) = zeros(n,1); % E_x_p + +end \ No newline at end of file diff --git a/func_ocp.m b/func_ocp.m new file mode 100644 index 0000000..58b589b --- /dev/null +++ b/func_ocp.m @@ -0,0 +1,103 @@ +function [out1,out2,out3] = func_ocp(flag,s,xx,ep,e,sigma,rho,kappa,lx,beta,rnstst) +% Function file for OPTIMAL COMMITMENT POLICY WITH ZERO FLOOR ON THE NOMINAL INTEREST RATE +% Program for the article "Optimal and Simple Monetary Policy Rules +% with Zero Floor on the Nominal Interest Rate", +% International Journal of Central Banking, Vol. 4(2), pages 73-127, June +% (C) Anton Nakov + +rn = s(:,1); +phi1_lag = s(:,2); +phi2_lag = s(:,3); + +n = length(rn); + +if ~isempty(ep) + E_pi = ep(:,1); + E_x = ep(:,2); +end + +if ~isempty(xx) + pi = xx(:,1); + x = xx(:,2); + i = xx(:,3); +end + +switch flag + +case 'b'; % BOUND FUNCTION + + out1(:,1) = -inf*ones(n,1); % pi low + out1(:,2) = -inf*ones(n,1); % x low + out1(:,3) = zeros(n,1); % i low + + out2(:,1) = inf*ones(n,1); % pi high + out2(:,2) = inf*ones(n,1); % x high + out2(:,3) = inf*ones(n,1); % i high + +case 'f'; % EQUILIBRIUM FUNCTION + out1(:,1) = pi - beta*E_pi - kappa*x; % f1 (NKPC) + out1(:,2) = x - E_x + 1/sigma*(i - E_pi - rn); % f2 (NIS) + out1(:,3) = i.*((kappa/sigma +1)/beta*phi1_lag + kappa*phi2_lag - kappa*pi - lx*x); % f3 (ZB) + + out2(:,1,1) = ones(n,1); % f1pi + out2(:,1,2) = -kappa*ones(n,1); % f1x + out2(:,1,3) = zeros(n,1); % f1i + + out2(:,2,1) = zeros(n,1); % f2pi + out2(:,2,2) = ones(n,1); % f2x + out2(:,2,3) = 1/sigma*ones(n,1); % f2i + + out2(:,3,1) = -kappa*i; % f3pi + out2(:,3,2) = -lx*i; % f3x + out2(:,3,3) = ((kappa/sigma +1)/beta*phi1_lag + kappa*phi2_lag - kappa*pi - lx*x); % f3i + + + out3(:,1,1) = -beta*ones(n,1); % f1E_pi + out3(:,1,2) = zeros(n,1); % f1E_x + + out3(:,2,1) = -ones(n,1)/sigma; % f2E_pi + out3(:,2,2) = -ones(n,1); % f2E_x + + out3(:,3,1) = zeros(n,1); % f3E_pi + out3(:,3,2) = zeros(n,1); % f3E_x + + +case 'g'; % STATE TRANSITION FUNCTION + out1(:,1) = rho*rn + (1-rho)*rnstst + e; % g1 + out1(:,2) = (kappa/sigma +1)/beta*phi1_lag + kappa*phi2_lag - kappa*pi - lx*x; % g2 + out1(:,3) = phi2_lag + phi1_lag/(beta*sigma) - pi; % g3 + + out2(:,1,1) = zeros(n,1); % g1pi + out2(:,1,2) = zeros(n,1); % g1x + out2(:,1,3) = zeros(n,1); % g1i + + out2(:,2,1) = -kappa*ones(n,1); % g2pi + out2(:,2,2) = -lx*ones(n,1); % g2x + out2(:,2,3) = zeros(n,1); % g2i + + out2(:,3,1) = -ones(n,1); % g3pi + out2(:,3,2) = zeros(n,1); % g3x + out2(:,3,3) = zeros(n,1); % g3i + +case 'h'; % EXPECTATION FUNCTION + out1(:,1) = pi; % E_pi + out1(:,2) = x; % E_x + + out2(:,1,1) = ones(n,1); % E_pi_pi + out2(:,1,2) = zeros(n,1); % E_pi_x + out2(:,1,3) = zeros(n,1); % E_pi_i + + out2(:,2,1) = zeros(n,1); % E_x_pi + out2(:,2,2) = ones(n,1); % E_x_x + out2(:,2,3) = zeros(n,1); % E_x_i + + out3(:,1,1) = zeros(n,1); % E_pi_rn + out3(:,2,1) = zeros(n,1); % E_x_rn + + out3(:,1,2) = zeros(n,1); % E_pi_phi1_lag + out3(:,2,2) = zeros(n,1); % E_x_phi1_lag + + out3(:,1,3) = zeros(n,1); % E_pi_phi2_lag + out3(:,2,3) = zeros(n,1); % E_x_phi2_lag + +end \ No newline at end of file diff --git a/func_odp.m b/func_odp.m new file mode 100644 index 0000000..88fa3a3 --- /dev/null +++ b/func_odp.m @@ -0,0 +1,67 @@ +% Function file for OPTIMAL DISCRETIONARY POLICY WITH ZERO FLOOR ON THE NOMINAL INTEREST RATE +% Program for the article "Optimal and Simple Monetary Policy Rules +% with Zero Floor on the Nominal Interest Rate", +% International Journal of Central Banking, Vol. 4(2), pages 73-127, June +% (C) Anton Nakov + +function [out1,out2,out3] = func_odp(flag,rn,xx,ep,e,sigma,rho,kappa,lx,beta,rnstst) +n = length(rn); +pi = xx(:,1); +x = xx(:,2); +i = xx(:,3); + +switch flag +case 'b'; % BOUND FUNCTION + out1(:,1) = -inf*ones(n,1); % pi low + out1(:,2) = -inf*ones(n,1); % x low + out1(:,3) = zeros(n,1); % i low + + out2(:,1) = inf*ones(n,1); % pi high + out2(:,2) = inf*ones(n,1); % x high + out2(:,3) = inf*ones(n,1); % i high + +case 'f'; % EQUILIBRIUM FUNCTION + out1(:,1) = pi - beta*ep(:,1) - kappa*x; % f1 + out1(:,2) = x - ep(:,2) + 1/sigma*(i - ep(:,1) - rn); % f2 + out1(:,3) = -i.*(lx*x + kappa*pi); % f3 + + out2(:,1,1) = ones(n,1); % f1pi + out2(:,1,2) = -kappa*ones(n,1); % f1x + out2(:,1,3) = zeros(n,1); % f1i + + out2(:,2,1) = zeros(n,1); % f2pi + out2(:,2,2) = ones(n,1); % f2x + out2(:,2,3) = 1/sigma*ones(n,1); % f2i + + out2(:,3,1) = -kappa*i; % f3pi + out2(:,3,2) = -lx*i; % f3x + out2(:,3,3) = -(lx*x + kappa*pi); % f3i + + out3(:,1,1) = -beta*ones(n,1); % f1e1 + out3(:,1,2) = zeros(n,1); % f1e2 + out3(:,2,1) = -ones(n,1)/sigma; % f2e1 + out3(:,2,2) = -ones(n,1); % f2e2 + out3(:,3,1) = zeros(n,1); % f3e1 + out3(:,3,2) = zeros(n,1); % f3e2 + +case 'g'; % STATE TRANSITION FUNCTION + out1 = rho*rn + (1-rho)*rnstst + e; % g + out2(:,1,1) = zeros(n,1); % gpi + out2(:,1,2) = zeros(n,1); % gx + out2(:,1,3) = zeros(n,1); % gi + +case 'h'; % EXPECTATION FUNCTION + out1(:,1) = pi; % e1 + out1(:,2) = x; % e2 + + out2(:,1,1) = ones(n,1); % e1pi + out2(:,1,2) = zeros(n,1); % e1x + out2(:,1,3) = zeros(n,1); % e1i + + out2(:,2,1) = zeros(n,1); % e2pi + out2(:,2,2) = ones(n,1); % e2x + out2(:,2,3) = zeros(n,1); % e2i + + out3(:,1,1) = zeros(n,1); % e1rn + out3(:,2,1) = zeros(n,1); % e2rn +end \ No newline at end of file diff --git a/func_ttr.m b/func_ttr.m new file mode 100644 index 0000000..0290671 --- /dev/null +++ b/func_ttr.m @@ -0,0 +1,64 @@ +function [out1,out2,out3] = func_ttr2(flag,s,xx,ep,e,sigma,rho,kappa,beta,rnstst,elb) +% FUNCTION FILE FOR TRUNCATED TAYLOR RULE +% Program for the article "Optimal and Simple Monetary Policy Rules +% with Zero Floor on the Nominal Interest Rate", +% International Journal of Central Banking, Vol. 4(2), pages 73-127, June +% (C) Anton Nakov + +rn = s(:,1); +n = length(rn); + +if ~isempty(ep) + E_pi = ep(:,1); + E_x = ep(:,2); +end +if ~isempty(xx) + pi = xx(:,1); + x = xx(:,2); +end + +switch flag + +case 'b'; % BOUND FUNCTION + + out1(:,1) = -inf*ones(n,1); % pi low + out1(:,2) = -inf*ones(n,1); % x low + + out2(:,1) = inf*ones(n,1); % pi high + out2(:,2) = inf*ones(n,1); % x high + +case 'f'; % EQUILIBRIUM FUNCTION + i = taylor(rnstst, pi, x, elb); + + out1(:,1) = pi - (beta*(E_pi) + kappa*x); % f1 (NKPC) + out1(:,2) = x - E_x + 1/sigma*(i - E_pi - rn); % f2 (NKIS) + + out2(:,1,1) = (ones(n,1)); % f1_pi + out2(:,1,2) = -kappa/ones(n,1); % f1_x + + out2(:,2,1) = zeros(n,1); % f2_pi + out2(:,2,2) = ones(n,1); % f2_x + + out3(:,1,1) = -beta/ones(n,1); % f1E_pi + out3(:,1,2) = zeros(n,1); % f1E_x + out3(:,2,1) = -ones(n,1)/sigma; % f2E_pi + out3(:,2,2) = -ones(n,1); % f2E_x + +case 'g'; % STATE TRANSITION FUNCTION + out1(:,1) = rho*rn + (1-rho)*rnstst + e; % g1 + out2(:,1,1) = zeros(n,1); % g1_pi + out2(:,1,2) = zeros(n,1); % g1_x + +case 'h'; % EXPECTATION FUNCTION + out1(:,1) = pi; % E(pi) + out1(:,2) = x; % E(x) + + out2(:,1,1) = ones(n,1); % E(pi)_pi + out2(:,1,2) = zeros(n,1); % E(pi)_x + + out2(:,2,1) = zeros(n,1); % E(x)_pi + out2(:,2,2) = ones(n,1); % E(x)_x + + out3(:,1,1) = zeros(n,1); % E(pi)_rn + out3(:,2,1) = zeros(n,1); % E(x)_rn +end \ No newline at end of file diff --git a/plot_ocp.m b/plot_ocp.m new file mode 100644 index 0000000..1176e86 --- /dev/null +++ b/plot_ocp.m @@ -0,0 +1,44 @@ +% PLOT OPTIMAL POLICY FUNCTIONS + + PI = 4*xx(:,:,:,1); % optimal inflation + X = 4*xx(:,:,:,2); % optimal output gap + I = 4*xx(:,:,:,3); % optimal nominal interest rate + + rn = 4*s{1}; + phi1 = 4*s{2}; + phi2 = 4*s{3}; + + pi = reshape(PI(:,:,2),length(rn),length(phi1)); + x = reshape(X(:,:,2),length(rn),length(phi1)); + i = reshape(I(:,:,2),length(rn),length(phi1)); + +% PLOT OPTIMAL INFLATION + figure + surf(phi1,rn,pi); + title('Inflation Policy'); + xlabel('\phi_{t-1}'); + ylabel('r^{n}_{t}'); + zlabel('\pi_{t}'); + view(60,20); + axis tight + +% PLOT OPTIMAL OUTPUT GAP + figure + surf(phi1,rn,x); + title('Output Gap'); + xlabel('\phi_{t-1}'); + ylabel('r^{n}_{t}'); + zlabel('x_{t}'); + view(60,20); + axis tight + +% PLOT OPTIMAL NOMINAL INTEREST RATE + figure + surf(phi1,rn,i); + title('Nominal interest rate','FontSize',14,'interpreter','latex'); + xlabel('$\phi_{t-1}$','FontSize',14,'interpreter','latex'); + ylabel('$r^{n}_{t}$','FontSize',14,'interpreter','latex'); + zlabel('$i_{t}$','FontSize',14,'interpreter','latex'); + %view(60,20); + axis tight + diff --git a/plot_odp.m b/plot_odp.m new file mode 100644 index 0000000..644af37 --- /dev/null +++ b/plot_odp.m @@ -0,0 +1,29 @@ +% PLOT OPTIMAL POLICY FUNCTIONS + pi = 4*xx(:,1); + x = 4*xx(:,2); + i = 4*xx(:,3); + s = 4*s; + xinit = 4*xinit; % solution without zero lower bound + + figure +% inflation + subplot(3,1,1) + plot(s,pi) + hold on + plot(s,xinit(:,1),'k:'); + ylabel('\pi'); + +% output gap + subplot(3,1,2) + plot(s,x) + hold on + plot(s,xinit(:,2),'k:'); + ylabel('x'); + +% nominal interest rate + subplot(3,1,3) + plot(s,i) + hold on + plot(s,xinit(:,3),'k:'); + xlabel('r_t^n'); + ylabel('i'); diff --git a/plot_paths.m b/plot_paths.m new file mode 100644 index 0000000..46236ec --- /dev/null +++ b/plot_paths.m @@ -0,0 +1,21 @@ +nperirf = 40; +figure(2) + subplot(1,3,1) + plot(pi_sim); + hold on + title('Inflation') + xlim([1 nperirf]) + subplot(1,3,2) + plot(x_sim); + hold on + title('Output gap') + xlim([1 nperirf]) + subplot(1,3,3) + plot(i_sim); + hold on + title('Interest rate') + hold on + plot(rn_sim,'k:'); + legend('Nominal','Natural') + xlim([1 nperirf]) + xlabel('Quarters') \ No newline at end of file diff --git a/plot_paths_do.m b/plot_paths_do.m new file mode 100644 index 0000000..d088a4e --- /dev/null +++ b/plot_paths_do.m @@ -0,0 +1,30 @@ + +figure(1) + +cut = [1 40]; + +subplot(231) +plot(rn_sim) +title('$r^n_t$','FontSize',14,'interpreter','latex') +set(gcf,'Name','Responses to natural real rate shock') +hold on +xlim(cut) + +subplot(232) +plot(i_sim) +title('$i_t$','FontSize',14,'interpreter','latex') +hold on +xlim(cut) + +subplot(233) +plot(pi_sim) +title('$\pi_t$','FontSize',14,'interpreter','latex') +hold on +xlim(cut) + +subplot(235) +plot(x_sim) +title('$y_t$','FontSize',14,'interpreter','latex') +hold on +xlim(cut) + diff --git a/plot_ttr.m b/plot_ttr.m new file mode 100644 index 0000000..634c7e8 --- /dev/null +++ b/plot_ttr.m @@ -0,0 +1,18 @@ +% PLOT TRUNCATED TAYLOR RULE + + figure + subplot(3,1,1) + plot(pi_sim); + title('Inflation') + xlim([1 nper]) + subplot(3,1,2) + plot(x_sim); + title('Output gap') + xlim([1 nper]) + subplot(3,1,3) + plot(i_sim); + title('Interest rate') + hold on + plot(rn_sim,'k:'); + legend('Nominal','Natural') + xlim([1 nper]) \ No newline at end of file diff --git a/ramsey_lin.mod b/ramsey_lin.mod new file mode 100644 index 0000000..8474c4a --- /dev/null +++ b/ramsey_lin.mod @@ -0,0 +1,83 @@ +var c n y r pi yhat ghat g beta; +varexo ub ug; + +parameters alpha gamma chi varphi epsilon b rhog rhob gbar lambda; + +gamma = 2; +chi = 11; +varphi = 1; +epsilon = 7; +b = (1/1.04)^(1/12); +rhob = 0.9; +rhog = 0.9; +alpha = 0.9; +gbar = 0.07; +phipi = 2; +phig = 0; +lambda = 0.003; + +model; +# nu_u = gamma*STEADY_STATE(y)/STEADY_STATE(c); +# nu_v = varphi*STEADY_STATE(y)/STEADY_STATE(n); +# capgamma = nu_u/(nu_u + nu_v); +# kappa = (1-alpha)*(1-alpha*b)/alpha*(nu_u + nu_v); +y = n; +y = c + g; +yhat-ghat = yhat(+1)-ghat(+1) - 1/nu_u*(r - pi(+1) + log(beta(+1))); +yhat = log(y/STEADY_STATE(y)); +ghat = (g-STEADY_STATE(g))/STEADY_STATE(y); +pi = beta(+1)*pi(+1) + kappa*(yhat-capgamma*ghat); +g/gbar = (g(-1)/gbar)^rhog*exp(ug); +beta/b = (beta(-1)/b)^rhob*exp(ub); +end; + +yss = fsolve(@(y) epsilon/(epsilon-1)*chi*y^varphi*(y-gbar)^gamma - 1,0.5); + +initval; +beta = b; +g = gbar; +r = -log(beta); +y = yss; +c = y - g; +n = y; +end; + +planner_objective pi^2 + lambda*yhat^2; + +ramsey_model(planner_discount=(1/1.04)^(1/12)); + +shocks; +var ub; +periods 1:1; +values 0.02; +end; + +ramsey_constraints; +r > 0; +%r < 0.005; +end; + +perfect_foresight_setup(periods=200); + +options_.stack_solve_algo = 7; +options_.solve_algo = 10; + +perfect_foresight_solver; + +figure(1) +cut = [0 200] +subplot(311) + plot(100*(exp(r).^12-1)); + title('r') + xlim(cut) + hold on +subplot(312) + plot(100*(yhat)); + title('yhat'); + xlim(cut) + hold on +subplot(313) + plot(100*(exp(pi).^12-1)); + title('pi'); + xlim(cut) + hold on \ No newline at end of file diff --git a/ramsey_nk.mod b/ramsey_nk.mod index 3556f25..25a4288 100644 --- a/ramsey_nk.mod +++ b/ramsey_nk.mod @@ -2,16 +2,14 @@ var y, pi, i, p, rr, rn; varexo e; parameters beta, gamma, sigma, kappa, r, rho, lambda; -% Program execution parameters - -% MODEL PARAMETERS +% PARAMETERS beta = 1/1.005; % quarterly time discount factor - gamma = 0.9; % indexation to past inflation + gamma = 0; % 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 +% NATURAL REAL RATE PROCESS r = 100*(1/beta-1); % steady-state (quarterly x 100) rho = 0.95; % persistence @@ -23,7 +21,7 @@ pi = p - p(-1); rn = r + rho*(rn(-1)-r) + e; end; -planner_objective pi^2 + lambda*y^2; +planner_objective pi^2 + lambda*y^2 + lambda*(i-rn)^2; ramsey_model(planner_discount=1/1.005); diff --git a/ramsey_nk_foc.mod b/ramsey_nk_foc.mod new file mode 100644 index 0000000..aba9d74 --- /dev/null +++ b/ramsey_nk_foc.mod @@ -0,0 +1,41 @@ +var y, pi, i, rr, rn, 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; % effective lower bound + +% 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); +[mcp = 'pi < 0.01'] +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); +rn = r + rho*(rn(-1)-r) + e; +end; + +initval; +rn = -0.5; +phi1 = 0; +end; + +endval; +rn = r; +end; + +perfect_foresight_setup(periods=200); +perfect_foresight_solver(lmmcp); + +do_irf; \ No newline at end of file diff --git a/ramsey_nk_michel.mod b/ramsey_nk_michel.mod new file mode 100644 index 0000000..2e755e1 --- /dev/null +++ b/ramsey_nk_michel.mod @@ -0,0 +1,45 @@ +var y, pi, i, rr, rn; +varexo e; +parameters beta, sigma, kappa, elb, r, rho, lambda; + +% 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: 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); +[mcp = 'pi < 0.01'] +pi = beta*pi(+1) + kappa*y; +rr = i - pi(+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 = -0.5; +end; + +ramsey_constraints; +i > 0; +end; + +perfect_foresight_setup(periods=200); + +perfect_foresight_setup(periods=200); +// useful to check what lmmcp is doing +options_.lmmcp.Display = 'iter'; +perfect_foresight_solver(lmmcp); + +do_irf; diff --git a/ramsey_nl.mod b/ramsey_nl.mod new file mode 100644 index 0000000..8278a95 --- /dev/null +++ b/ramsey_nl.mod @@ -0,0 +1,85 @@ +var c n y yhat r pi pstar w D K F; +var beta; +varexo e; + +parameters b gamma chi phi eps rho alpha g lambda; + +b = (1/1.04)^(1/12); +gamma = 2; +chi = 11; +phi = 1; +eps = 7; +rho = 0.9; +alpha = 0.10; +g = 0.07; +lambda = 0.003; + +model; +y = c + g; +y = n/D; +yhat = log(y/STEADY_STATE(y)); +w = chi*c^gamma*n^phi; +1 = r*beta(+1)*(c(+1)/c)^(-gamma)/pi(+1); +1 = alpha*pstar^(1-eps) + (1-alpha)*pi^(eps-1); +D = alpha*pstar^(-eps) + (1-alpha)*pi^(eps)*D(-1); +K = c^(-gamma)*eps/(eps-1)*w*y + beta(+1)*(1-alpha)*pi(+1)^(eps)*K(+1); +F = c^(-gamma)*y + beta(+1)*(1-alpha)*pi(+1)^(eps-1)*F(+1); +pstar = K/F; +log(beta)= (1-rho)*log(b) + rho*log(beta(-1)) + e; +end; + +planner_objective (pi-1)^2 + lambda*yhat^2; +% planner_objective -log(c) + chi*n^(1+phi)/(1+phi); + +ramsey_model(planner_discount=(1/1.04)^(1/12)); + +yss = fsolve(@(y) eps/(eps-1)*chi*y^phi*(y-g)^gamma - 1,0.475018591498115); + +initval; +y = yss; +beta = b; +r = 1/beta; +pi = 1; +pstar = 1; +D = 1; +c = y - g; +n = y; +w = chi*n^phi*c^gamma; +K = c^(-gamma)*eps/(eps-1)*w*y/(1-beta*alpha); +F = c^(-gamma)*y/(1-beta*alpha); +end; + +shocks; +var e; +periods 1:1; +values 0.02; +end; + +ramsey_constraints; +r > 1; +end; + +perfect_foresight_setup(periods=200); + +options_.stack_solve_algo = 7; +options_.solve_algo = 10; + +perfect_foresight_solver; + +figure(1) +cut = [0 200]; +subplot(311) + plot(100*(r.^12-1)); + title('r') + xlim(cut) + hold on +subplot(312) + plot(100*(yhat)); + title('yhat'); + xlim(cut) + hold on +subplot(313) + plot(100*(pi.^12-1)); + title('pi'); + xlim(cut) + hold on \ No newline at end of file diff --git a/ramsey_villaverde.mod b/ramsey_villaverde.mod new file mode 100644 index 0000000..e33b68e --- /dev/null +++ b/ramsey_villaverde.mod @@ -0,0 +1,88 @@ +var Y, mcHAT, yHAT, C, L, D, N, w, PI, PIstar, R, DP, mc, beta; + +varexo e; + +parameters alpha, b, theta, epsilon, psi, lambda, PI_bar, R_bar, rho; + +PI_bar = 1^(1/4); +R_bar = 1.04^(1/4); +b = PI_bar/R_bar; +alpha = 0.667; +theta = 0.75; +epsilon = 6; +psi = 1; +rho = 0.9; +lambda = 0.003; + +model; +C = Y; +w = C*L^psi; +1 = beta(+1)*R*C/C(+1)/PI(+1); +N = mc*Y/C + theta*beta(+1)*PI(+1)^epsilon*N(+1); +D = PIstar*(Y/C + theta*beta(+1)*PI(+1)^(epsilon-1)/PIstar(+1)*D(+1)); +N/D = (epsilon-1)/epsilon; +1 = theta*(PI)^(epsilon-1) + (1-theta)*(PIstar)^(1-epsilon); +DP = theta*(PI)^epsilon*DP(-1) + (1-theta)*(PIstar)^(-epsilon); +w*L = alpha*mc*Y*DP; +Y = L^alpha/DP; +mcHAT = log(mc/STEADY_STATE(mc)); +yHAT = log(Y/STEADY_STATE(Y)); +log(beta)= (1-rho)*log(b) + rho*log(beta(-1)) + e; +end; + +initval; +beta = b; +PI = PI_bar; +R = PI/beta; +PIstar = ((1 - theta*(PI)^(epsilon-1))/(1-theta))^(1/(1-epsilon)); +DP = (1-theta)*PIstar^(-epsilon)/(1-theta*PI^(epsilon)); +mc = (epsilon-1)/epsilon; % approx.(assuming flex. markup) +L = (alpha*mc*DP)^(1/(1+psi)); % approx.(assuming Y = C) +Y = L^alpha/DP; +C = Y; +w = C*L^psi; +D = PIstar*Y/C/(1-beta*theta*PI^(epsilon-1)/PIstar); +N = (epsilon-1)*D/epsilon; +end; + + planner_objective - log(C) + L^(1+psi)/(1+psi); +% planner_objective (PI-1)^2 + lambda*(mcHAT)^2; + +ramsey_model(planner_discount=1.0^(1/4)/1.04^(1/4)); + +shocks; +var e; +periods 1:1; +values 0.02; +end; + +ramsey_constraints; +R > 1; +end; + +perfect_foresight_setup(periods=200); + +options_.stack_solve_algo = 7; +options_.solve_algo = 10; + +perfect_foresight_solver; + +figure(1) +cut = [0 200] +subplot(311) + plot(100*(R.^4-1)); + title('R') + xlim(cut) + hold on +subplot(312) + plot(100*(PI.^4-1)); + title('PI') + xlim(cut) + hold on +subplot(313) + plot(100*(yHAT)); + title('yHAT') + xlim(cut) + hold on + + diff --git a/ramsey_villaverde_crra.mod b/ramsey_villaverde_crra.mod new file mode 100644 index 0000000..a8338e4 --- /dev/null +++ b/ramsey_villaverde_crra.mod @@ -0,0 +1,89 @@ +var Y, mcHAT, yHAT, C, L, D, N, w, PI, PIstar, R, DP, mc, beta; + +varexo e; + +parameters alpha, b, gamma, theta, epsilon, psi, lambda, PI_bar, R_bar, rho; + +PI_bar = 1.0^(1/4); +R_bar = 1.04^(1/4); +b = PI_bar/R_bar; +alpha = 1; +theta = 0.75; +epsilon = 6; +psi = 1; +rho = 0.9; +gamma = 1; +lambda = 0.003; + +model; +C = Y; +w = C^gamma*L^psi; +1 = beta(+1)*R*(C/C(+1))^gamma/PI(+1); +N = mc*Y/C^gamma +theta*beta(+1)*PI(+1)^epsilon*N(+1); +D = PIstar*(Y/C^gamma + theta*beta(+1)*PI(+1)^(epsilon-1)/PIstar(+1)*D(+1)); +N/D = (epsilon-1)/epsilon; +1 = theta*(PI)^(epsilon-1) + (1-theta)*(PIstar)^(1-epsilon); +DP = theta*(PI)^epsilon*DP(-1) + (1-theta)*(PIstar)^(-epsilon); +w*L = alpha*mc*Y*DP; +Y = L^alpha/DP; +mcHAT = log(mc/STEADY_STATE(mc)); +yHAT = log(Y/STEADY_STATE(Y)); +log(beta)= (1-rho)*log(b) + rho*log(beta(-1)) + e; +end; + +initval; +beta = b; +PI = PI_bar; +R = PI/beta; +PIstar = ((1 - theta*(PI)^(epsilon-1))/(1-theta))^(1/(1-epsilon)); +DP = (1-theta)*PIstar^(-epsilon)/(1-theta*PI^(epsilon)); +mc = (epsilon-1)/epsilon; % approx.(assuming flex. markup) +L = (alpha*mc*DP)^(1/(1+psi)); % approx.(assuming Y = C) +Y = L^alpha/DP; +C = Y; +w = C^gamma*L^psi; +D = PIstar*Y/C^gamma/(1-beta*theta*PI^(epsilon-1)/PIstar); +N = (epsilon-1)*D/epsilon; +end; + +% planner_objective - C^(1-gamma)/(1-gamma) + L^(1+psi)/(1+psi); + planner_objective - log(C) + L^(1+psi)/(1+psi); +% planner_objective (PI-1)^2 + lambda*(mcHAT)^2; + +ramsey_model(planner_discount=1.0^(1/4)/1.04^(1/4)); + +shocks; +var e; +periods 1:1; +values 0.02; +end; + +ramsey_constraints; +R > 1; +end; + +perfect_foresight_setup(periods=200); + +options_.stack_solve_algo = 7; +options_.solve_algo = 10; + +perfect_foresight_solver; + +figure(1) +subplot(311) + plot(100*(R.^4-1)); + title('R') + xlim([0 50]) + hold on +subplot(312) + plot(100*(PI.^4-1)); + title('PI') + xlim([0 50]) + hold on +subplot(313) + plot(100*(yHAT)); + title('yHAT') + xlim([0 50]) + hold on + + diff --git a/remsimul.m b/remsimul.m new file mode 100644 index 0000000..17c0f07 --- /dev/null +++ b/remsimul.m @@ -0,0 +1,50 @@ +% REMSIMUL Monte Carlo simulation of discrete time controlled Markov process +% USAGE +% st = remsimul(model,s0,nper,sres,xres,) +% INPUTS +% model : name of model structure +% fspace : name of projection space structure +% s0 : k by 1 vector of initial states +% nper : number of simulated time periods +% sres : coordinates of the evaluation grid (from remsolve) +% xres : optimal control function values +% at grid defined by sres (from remsolve) +% OUTPUT +% st : k by nyrs vector of simulated states +% +% USES: DISCRAND + +% Copyright (c) 1997-2000, Paul L. Fackler & Mario J. Miranda +% paul_fackler@ncsu.edu, miranda.4@osu.edu + +function [ssim,xsim] = remsimul(model,s0,nper,sres,xres) + +% DUMMY SHOCKS/WEIGHTS IF MODEL DETERMINISTIC + if isfield(model,'e'), e=model.e; else, e=0; end; + if isfield(model,'w'), w=model.w; else, w=1; end; + + nrep = size(s0,1); + ds = size(s0,2); + if ds>1 + st = gridmake(sres); + dx = ds*length(xres(:))/length(st(:)); + else + dx = length(xres(:))/length(sres(:)); + end + ssim = zeros(nrep,ds,nper+1); + xsim = zeros(nrep,dx,nper+1); + + nx = prod(size(xres))/dx; + xres = reshape(xres,nx,dx); + for t=1:nper+1 + xx = minterp(sres,xres,s0); + ssim(:,:,t) = s0; + xsim(:,:,t) = xx; + ee = e(discrand(nrep,w),:); + s0 = feval(model.func,'g',s0,xx,[],ee,model.params{:}); + end + + ssim = ssim(:,:,1:nper+1); + xsim = xsim(:,:,1:nper+1); + if ds==1; ssim=squeeze(ssim(:,1,:)); end + if dx==1; xsim=squeeze(xsim(:,1,:)); end \ No newline at end of file diff --git a/remsolve.m b/remsolve.m new file mode 100644 index 0000000..ff4a716 --- /dev/null +++ b/remsolve.m @@ -0,0 +1,176 @@ +% REMSOLVE Solves rational expectations models +% USAGE +% [c,scoord,x,h,f,resid] = remsolve(model,fspace,scoord,x,h); +% INPUTS +% model : rational expectations model (structure array) +% fspace : projection space structure (structure array) +% scoord : n-vector or grid (cell array) of state nodes +% x : initial guess for equilibrium responses at nodes +% h : initial guess for expectation variable at nodes +% OUTPUTS +% c : expectation function approximation basis coefficients +% scoord : residual evaluation coordinates (cell array for n>1) +% x : equilibrium responses at evaluation points +% h : expectation variables at evaluation points +% f : marginal arbitrage benefits at evaluation points +% resid : rational expectation residuals at evaluation points +% MODEL STRUCTURE FIELDS: +% func : model function file (see below) +% e : shocks +% w : probabilities +% params : additional parameters to function file +% MODEL FUNCTION FILE FORMAT: +% [out1,out2,out3] = func(flag,s,x,e,additional parameters) +% if flag = 'b' returns bound function +% xl:n.dx, xu:n.dx +% if flag = 'f' returns reward function and derivatives +% f:n.dx, fx:n.dx.dx, feh:n.dx.dh +% if flag = 'g' returns transition function and derivatives +% g:n.ds, gx:n.ds.dx +% if flag = 'h' returns expectation function and derivatives +% h:n.dh, hx:n.dh.dx, hs:n.dh.ds +% where n = number of collocation nodes +% ds = state space dimension +% dx = response space dimension +% dh = expectation space dimension +% USER OPTIONS FOR REMSOLVE (SET WITH optset('remsolve',option,value)): +% tol : convergence tolerance +% maxit : maximum number of iterations +% nres : nres*fspace.n uniform nodes to evaluate residual +% showiters : 0/1, 1 to display iteration results +% USER OPTIONS FOR ARBIT (SET WITH optset('arbit',option,value)): +% maxit : maximum number of iterations for CP +% tol : convergence tolerance for CP +% lcpmethod : 'minmax' or 'smooth' + +% Copyright (c) 1997-2002, Paul L. Fackler & Mario J. Miranda +% paul_fackler@ncsu.edu, miranda.4@osu.edu + +function [c,scoord,x,h,f,resid] = remsolve(model,fspace,scoord,x,h) + +% SET PARAMETER DEFAULTS + tol = optget('remsolve','tol',sqrt(eps)); + maxit = optget('remsolve','maxit',500); + nres = optget('remsolve','nres',10); + showiters = optget('remsolve','showiters',1); + arbmaxit = optget('arbit','maxit',300); + arbtol = optget('arbit','tol',sqrt(eps)); + lcpmethod = optget('arbit','lcpmethod','minmax'); + +% UNPACK MODEL STRUCTURE + func=model.func; + params=model.params; + if ~isfield(model,'e'); e=0; else, e=model.e; end; + if ~isfield(model,'w'); w=1; else, w=model.w; end; + +% DETERMINE NUMBER OF DIMENSIONS & COORDINATES + ni = fspace.n; % number of collocation coordinates by state dimension + N = prod(ni); % number of collocation states + n = length(ni); % dimension of state space + m = size(x,2); % dimension of response space + +% COMPUTE COLLOCATION NODES & INTERPOLATION MATRIX +% scoord = funnode(fspace); % state collocation coordinates + s = gridmake(scoord); % state collocation nodes + if nargin<5 | isempty(h) + h = feval(func,'h',s,x,[],[],params{:}); + end + p = size(h,2); % dimension of expectation space + Phi = funbasx(fspace); % collocation matrix + c = funfitxy(fspace,Phi,h); % initial basis coefficients + +% SOLVE EULER EQUATION + tic + for it=1:maxit % perform function iterations + cold = c; % store old basis coefficients + [f,x,h] = arbit(s,x,c,fspace,func,params,... + e,w,n,m,p,arbmaxit,arbtol,lcpmethod); % update expectation and response variables + c = funfitxy(fspace,Phi,h); % update basis coefficient + change = norm(c-cold,inf); % compute change + if showiters, fprintf ('%4i %10.1e\n',it,change), end % print progress + if changefspace.b+eps), disp('Warning: extrapolating beyond smax'), end; + +% COMPUTE RESIDUAL + if nargout>5 + ind = ones(1,n); + if isfield(model,'discretestates') + ind(model.discretestates) = 0; + end + if n==1, + if ind + ni = nres*ni; + scoord = linspace(fspace.a,fspace.b,ni)'; + end + else + for i=1:n, + if ind(i) + ni(i) = nres*ni(i); + scoord{i} = linspace(fspace.a(i),fspace.b(i),ni(i))'; + end + end + end + s = gridmake(scoord); + x = funeval(funfitxy(fspace,Phi,x),fspace,scoord); % rough guess for responses at evaluation points + [f,x,h] = arbit(s,x,c,fspace,func,params,... + e,w,n,m,p,arbmaxit,arbtol,lcpmethod); % shadow prices and actions at evaluation points + resid = h-funeval(c,fspace,s); % residual at evaluation points + resid = reshape(resid,[ni p]); % reshape residual for plotting + end + +% RESHAPE OUTPUT + x = reshape(x,[ni m]); + h = reshape(h,[ni p]); + f = reshape(f,[ni m]); +return + + + +% ARBIT - Solves arbitrage equation at specified nodes + function [f,x,h] = arbit(s,x,c,fspace,func,params,e,w,n,m,p,maxit,tol,lcpmethod) + +% COMPUTE BOUNDS + [xl,xu] = feval(func,'b',s,x,[],[],params{:}); + + N=size(s,1); +% SOLVE FIRST ORDER CONDITIONS + for it=1:maxit + xold = x; + [f,fx]=equilibrium(s,x,c,fspace,func,params,e,w,n,m,p); + [lcpf,deltax] = lcpstep(lcpmethod,x,xl,xu,reshape(f,N,m),fx); + x = x + deltax; + if norm(deltax(:))< tol, break, end; + if any(isnan(x)), error('NaNs encountered'); end + end + h = feval(func,'h',s,x,[],[],params{:}); +return + + +% EQUILIBRIUM Computes the equilibrium condition at arbitrary s and x +function [f,fx]=equilibrium(s,x,c,fspace,func,params,e,w,n,m,p) + K = length(w); + N=size(s,1); + eh = zeros(N,p); + ehder = zeros(N,p,m); + for k=1:K + kk = k+zeros(N,1); + [g,gx] = feval(func,'g',s,x,[],e(kk,:),params{:}); + [hnext,hnextder] = fund(c,fspace,g,1); + eh = eh + w(k)*hnext; + ehder = ehder + w(k)*arraymult(hnextder,gx,N,p,n,m); + clear hnext hnextder gx g + end + [f,fx,feh] = feval(func,'f',s,x,eh,[],params{:}); + fx = fx + arraymult(feh,ehder,N,m,p,m); +return \ No newline at end of file diff --git a/simultrap.m b/simultrap.m new file mode 100644 index 0000000..0038db3 --- /dev/null +++ b/simultrap.m @@ -0,0 +1,49 @@ +% Monte Carlo simulation of discrete time controlled Markov process +% Modified by Anton Nakov (2005) + +% INPUTS +% model : model structure variable +% ss : k by d vector of initial states +% nper : number of simulated time periods +% sres : coordinates of the evaluation grid (from dpsolve) +% xres : optimal control function values +% at grid defined by sres +% OUTPUTS +% ssim : k by d by nper+1 vector of simulated states +% xsim : k by d by nper+1 vector of simulated actions +% For finite horizon problems, xsim is k by d by nper +% If d=1, ssim and xsim are k by nper+1 +% +% USES: DISCRAND + +% Copyright (c) 1997-2002, Paul L. Fackler & Mario J. Miranda +% paul_fackler@ncsu.edu, miranda.4@osu.edu + + +function [ssim,xsim] = simultrap(model,ss,nper,sres,xres) + + func = model.func; + params = model.params; + + nrep = size(ss,1); + ds = size(ss,2); + st = gridmake(sres); + + dx = ds*length(xres(:))/length(st(:)); + + ssim = zeros(nrep,ds,nper+1); + xsim = zeros(nrep,dx,nper+1); + + nx = numel(xres)/dx; + xres = reshape(xres,nx,dx); + + for t=1:nper+1 + xx = minterp(sres,xres,ss); + ssim(:,:,t) = ss; + xsim(:,:,t) = xx; + ee = 0; + ss = feval(func,'g',ss,xx,[],ee,params{:}); + end + + if ds==1; ssim=squeeze(ssim); end + if dx==1; xsim=squeeze(xsim); end \ No newline at end of file diff --git a/taylor.m b/taylor.m new file mode 100644 index 0000000..23aa8de --- /dev/null +++ b/taylor.m @@ -0,0 +1,3 @@ +function i = taylor(iss, pi, x, elb) + +i = max(elb, iss + 1.5*pi + 0.5*x); \ No newline at end of file diff --git a/ttr_2_constraints.mod b/ttr_2_constraints.mod new file mode 100644 index 0000000..e60b42e --- /dev/null +++ b/ttr_2_constraints.mod @@ -0,0 +1,53 @@ +var y, pi, i, rr, rn, u; +varexo e, eps; +parameters beta, sigma, kappa, elb, ub, 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 + ub = 3/4; % upper 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 = min(max(elb, r + phi_pi*pi + phi_y*y + u),ub); +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:21 22:200; +values 0 0.1 -0.1 0; +end; + +simul(periods=200, maxit=500, stack_solve_algo=0); +do_irf; \ No newline at end of file diff --git a/zlb.mod b/zlb.mod new file mode 100644 index 0000000..4a72581 --- /dev/null +++ b/zlb.mod @@ -0,0 +1,64 @@ +var Y, C, L, D, N, w, PI, PIstar, R, DP, mc, beta; + +varexo e; + +parameters alpha, beta_ss, theta, epsilon, psi, phi_pi, PI_bar, R_bar, rho; + +PI_bar = 1.0^(1/4); +R_bar = 1.04^(1/4); +beta_ss = PI_bar/R_bar; +alpha = 1; +theta = 0.75; +epsilon = 6; +psi = 1; +phi_pi = 1.5; +rho = 0.9; + +model; +C = Y; +w = C*L^psi; +1 = beta(+1)*R*C/C(+1)/PI(+1); +N = mc*Y/C +theta*beta(+1)*PI(+1)^epsilon*N(+1); +D = PIstar*(Y/C + theta*beta(+1)*PI(+1)^(epsilon-1)/PIstar(+1)*D(+1)); +N/D = (epsilon-1)/epsilon; +1 = theta*(PI)^(epsilon-1) + (1-theta)*(PIstar)^(1-epsilon); +DP = theta*(PI)^epsilon*DP(-1) + (1-theta)*(PIstar)^(-epsilon); +w*L = alpha*mc*Y*DP; +Y = L^alpha/DP; +R = max(1, PI_bar/beta_ss + phi_pi*(PI-PI_bar)); +log(beta)= (1-rho)*log(beta_ss) + rho*log(beta(-1)) + e; +end; + +initval; +beta = beta_ss; +PI = PI_bar; +R = PI/beta; +PIstar = ((1 - theta*(PI)^(epsilon-1))/(1-theta))^(1/(1-epsilon)); +DP = (1-theta)*PIstar^(-epsilon)/(1-theta*PI^(epsilon)); +mc = (epsilon-1)/epsilon; % approximation (assuming flex. markup) +L = (alpha*mc*DP)^(1/(1+psi)); % approximation (assuming Y = C) +Y = L^alpha/DP; +C = Y; +w = C*L^psi; +D = PIstar*Y/C/(1-beta*theta*PI^(epsilon-1)/PIstar); +N = (epsilon-1)*D/epsilon; +end; + +steady; check; + +shocks; +var e; +periods 1:1 2:200 ; +values 0.01 0; +end; + +simul(periods=200, maxit=500, stack_solve_algo=0); + +figure(1) +subplot(211) +plot(100*(R.^4-1)); +xlim([0 50]) +subplot(212) +plot(100*(PI.^4-1)); +xlim([0 50]) + diff --git a/zlb_lin.mod b/zlb_lin.mod new file mode 100644 index 0000000..77aad90 --- /dev/null +++ b/zlb_lin.mod @@ -0,0 +1,59 @@ +var y, pi, i, rr, rn; +varexo e; +parameters beta, sigma, kappa, r, rho, phi_pi, pibar, R_bar; + +% Program execution parameters + +% MODEL PARAMETERS +PI_bar = 1.0^(1/4); +R_bar = 1.04^(1/4); +pibar = 100*(PI_bar - 1); +beta = PI_bar/R_bar; +theta = 0.75; +epsilon = 6; +psi = 1; +phi_pi = 1.5; +r = 100*(1/beta-1); % steady-state natural real rate (quarterly x100) +rho = 0.9; +sigma = 1; % relative risk aversion +phi_pi = 1.5; % reaction to inflation + +%kappa = (1-theta)*(1-theta*beta)/theta*(psi+sigma)/(1+psi*epsilon); +kappa = (1-theta)*(1-theta*beta)/theta*(psi+sigma); +% slope of the Phillips curve + +model; +y = y(+1) - 1/sigma*(rr - rn); +pi-pibar = beta*(pi(+1)-pibar) + kappa*y; +i = max(0, r + pibar + phi_pi*(pi-pibar)); +rr = i - pi(+1); +rn = r + rho*(rn(-1)-r) + e; +end; + +initval; +rn = r; +pi = pibar; +rr = r; +i = 100*(R_bar-1); +y = 0; +end; + +steady; check; + +shocks; +var e; +periods 1:1 2:200 ; +values -0.9 0; +end; + +simul(periods=200, maxit=500, stack_solve_algo=0); + +figure(1) +subplot(211) +hold on +plot(((1+i/100).^4-1)*100) +xlim([0 50]) +subplot(212) +hold on +plot(((1+pi/100).^4-1)*100) +xlim([0 50])