var c n I k x u eta lambda mu delta_u y g z ng nz;
varexo eg ez e_ng e_nz ;

parameters beta theta alpha gamma phi uParam psi omega delta_0;

%
% Parameters:
%

beta = 0.985; % subjective discount factor
theta = 1.4; % labor elasticity
alpha = 0.64; % labor share in production
gamma = 0.001; % The KEY Parameter
phi = 1.3; % Adjustment cost param
uParam = 0.15; %Util Param.
psi = 1; %CHECK CHECK
omega = 1 + uParam;
delta_0 = 1;

%
% Steady State Values
%

u_bar = (((1/beta)-1)/(delta_0-(delta_0/omega)))^(1/omega);
delta_u_bar = delta_0*(u_bar^omega)/omega;

AX = (1/beta - (1-delta_u_bar))/(1-alpha);
BX = (theta*(AX-delta_u_bar)/(alpha*AX)) - (gamma/((beta*(1-gamma))-1));

k_bar = (1/(psi*BX)^(1/theta))*(((u_bar^(1-alpha))/AX)^(1/alpha));
n_bar = ((AX/u_bar^(1-alpha))^(1/alpha))*k_bar;
c_bar = (AX-delta_u_bar)*k_bar;
I_bar = delta_u_bar*k_bar;
x_bar = c_bar;
mu_bar = (1/c_bar)*(1/(1-(psi*(n_bar^theta))))*(psi*(n_bar^theta))/((beta*(1-gamma))-1);
lambda_bar = ((((beta*(1-gamma))-1)/(psi*(n_bar^theta))) + gamma)*mu_bar;
eta_bar = lambda_bar;
y_bar = c_bar + I_bar;

g_bar = 1;
z_bar = 1;
model;

(1/(c-(psi*(n^theta)*x))) + (mu*gamma*(x(-1)/(g*c))^(1-gamma)) = lambda;
(psi*(n^theta)/(c-(psi*(n^theta)*x))) + mu = beta*(mu(+1)/g(+1))*(1-gamma)*((g(+1)*c(+1)/x)^gamma);
(theta*psi*x*(n^(theta-1))/(c-(psi*(n^theta)*x))) = lambda*alpha*((u*k(-1)/(g*z^(1/(1-alpha))*n))^(1-alpha));
eta = (beta*(lambda(+1)/g(+1)/z(+1)^(1/(1-alpha)))*(1-alpha)*(u^(1-alpha))*((n(+1)*g(+1)/k)^alpha)) + (beta*(eta(+1)/g(+1)/z(+1)^(1/(1-alpha)))*(1-delta_u(+1)));
lambda*(1-alpha)*(u^(-1*alpha))*((k(-1)/g)^(-alpha))*(n^alpha) = eta*delta_0*(u^(omega-1));
lambda/z = eta - (phi*eta*0.5*(1-(I*g*z^(1/(1-alpha))/(I(-1))))^2) + (phi*eta*(1-(I*g*z^(1/(1-alpha))/(I(-1))))*(I*g*z^(1/(1-alpha))/(I(-1)))) - (beta*phi*(eta(+1)/g(+1))*(1-(I(+1)*g(+1)*z(+1)^(1/(1-alpha))/I))*((I(+1)*g(+1)/I)^2));
x = (c^gamma)*((x(-1)/g)^(1-gamma));
y = c + I/z;
c + I/z = ((u*k(-1)/g)^(1-alpha))*(n^alpha);
k = (((1-delta_u)*k(-1))/g) + I - (0.5*phi*((1 - (I*g/(I(-1))))^2)*I);
delta_u = delta_0*(u^omega)/omega;
g=(eg+ng(-4));
ng=e_ng;
z=(ez+nz(-4));
nz=e_nz;
end;

initval;

c = c_bar;
n = n_bar;
I = I_bar;
k = k_bar;
x = x_bar;
u = u_bar;
eta = eta_bar;
lambda = lambda_bar;
mu = mu_bar;
delta_u = delta_u_bar;
y = y_bar;

g = g_bar;
z = z_bar;
end;

//maxit_ = 1000;

steady;
check;

//shocks;
//var g;
//periods 4;
//values 1.01;
//var g;
//periods 1:3;
//values 1;
//var g;
//periods 5:1000;
//values 1;
//end;
schocks;
var eng;stderr 0.01;
var enz;stderr 0.01;
end;

stoch_simul(periods=1000,irf=20);

startFigure = 1;
lineWidthOfPlot = 1;

% Set this Parameters:
periodOfShock = 4;
periodOfSimulations = 1000;
numberOfPeriodToPlot = 1000;

A = ones(periodOfSimulations+2,1);
A(1) = 1;
A(periodOfShock+1) = 1.01;
A(periodOfShock+2:periodOfSimulations+2) = 1.01;

%y_actual = log((y.*A)./y(1))*100;
%I_actual = log((I.*A)./I(1))*100;
%n_actual = log((n)./n(1))*100;
%x_actual = log((x.*A)./x(1))*100;
%k_actual = log((k.*A)./k(1))*100;
%c_actual = log((c.*A)./c(1))*100;

y_actual = (((y.*A)-y(1))./y(1))*100;
I_actual = (((I.*A)-I(1))./I(1))*100;
n_actual = (((n)-n(1))./n(1))*100;
x_actual = (((x.*A)-x(1))./x(1))*100;
k_actual = (((k.*A)-k(1))./k(1))*100;
c_actual = (((c.*A)-c(1))./c(1))*100;


% Figure 1
figure(startFigure)
subplot(2,3,1);
hold on;
plot(y_actual(1:numberOfPeriodToPlot),'r-o','LineWidth',lineWidthOfPlot);
%title('Output w/ trend')
title('Output - Y_{t}')

subplot(2,3,2);
hold on;
plot(n_actual(1:numberOfPeriodToPlot),'r-o','LineWidth',lineWidthOfPlot);
%title('Labor')
title('Labor - N_{t}')

subplot(2,3,3);
hold on;
plot(x_actual(1:numberOfPeriodToPlot),'r-o','LineWidth',lineWidthOfPlot);
%title('Durable Good Consumption w/ trend')
title('X Param - X_{t}')

subplot(2,3,4);
hold on;
plot(I_actual(1:numberOfPeriodToPlot),'r-o','LineWidth',lineWidthOfPlot);
%title('Investment w/ trend')
title('Investment - I_{t}')

subplot(2,3,5);
hold on;
plot(k_actual(1:numberOfPeriodToPlot),'r-o','LineWidth',lineWidthOfPlot);
%title('Capital w/ trend')
title('Capital - K_{t}')

subplot(2,3,6);
hold on;
plot(c_actual(1:numberOfPeriodToPlot),'r-o','LineWidth',lineWidthOfPlot);
%title('Durable Good Produced w/ trend')
title('Consumption - C_{t}')