% New Keynesian Model with lower bound
% Michael Paetz, November 15, 2011
% You find the complete algorithm in the attachment of Tom Holden's paper: http://www.dynare.org/wp-repo/dynarewp004.pdf
%In line 6 and 7 you have to define the number of periods you would like to simulate the IRFs (TIRF), and the number of periods, after which you believe the constraint will no longer bind (TShadowI)

@#define TshadowI = 4
@#define TIRF = 4

var pi, y, i, shock_pi;

@#for index in 0:( TshadowI - 1 )
	varexo epsilon_shadowI_@{index};
	var LAG0_shadowI_@{index};
	@#for offset in 1:index
		var LAG@{offset}_shadowI_@{index};
	@#endfor
	@#define offset = 10000000000
@#endfor
@#define index = 10000000000

varexo epsilon_pi;

parameters alpha, beta, kappa, lambda, phi, sigma, epsilon, omega, bigtheta, phi_pi, rho_pi, i_bar;

beta = 0.95;
omega = 0.8598;
epsilon = 6;
sigma = 5;
phi_pi = 1.5;
i_bar = 1/beta;
alpha = 0.33;
phi = 0.33;
bigtheta = ((1 - alpha)/(1 - alpha + alpha * epsilon));
lambda = (((1-omega)*(1-beta*omega))/omega)*bigtheta;
kappa = lambda * (sigma + ((phi + alpha)/(1 - alpha)));
rho_pi = 0.7;

//Model

model(linear);

	@#for index in 0:( TshadowI - 1 )
		LAG0_shadowI_@{index} = epsilon_shadowI_@{index};
		@#for offset in 1:index
			LAG@{offset}_shadowI_@{index}(0) = LAG@{offset-1}_shadowI_@{index}(-1);
		@#endfor
		@#define offset = 10000000000
	@#endfor
	@#define index = 10000000000
	
	#shadowI_price = 0
		@#for index in 0:( TshadowI - 1 )
			+ LAG@{index}_shadowI_@{index}(0)
		@#endfor
		@#define index = 10000000000
	;

y=y(+1)-sigma^(-1)*(i - log(i_bar) - pi(+1));

i = log(i_bar) + phi_pi * pi + shadowI_price;

pi = 0.8 * beta * pi(+1) + 0.2 * pi(-1)+  kappa* y - shock_pi;

shock_pi = rho_pi * shock_pi(-1) + epsilon_pi;

end;

initval;
y = 0;
pi = 0;
end;

check;
steady;

shocks;
	@#for index in 0:( TshadowI - 1 )
		var epsilon_shadowI_@{index} = 1;
	@#endfor
	@#define index = 10000000000

var epsilon_pi; stderr 0.2;
end;

stoch_simul(order = 1, irf = @{TIRF}, noprint, nograph);

solve_shadowfun(@{TshadowI},@{TIRF});