%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Dmitri Bogonos
% dmitri.bogonos@rwi-essen.de
% 05.08.2015
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define variables, shocks and parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


var y c g p_k j i z_k z_y v_k v_y n_y a_k a_y h_k h_y x k l o_k o_y pi_k pi_y p_k_st p_k_et n_k j_y j_k Lambda i_s i_e lambda_y lambda_k p_k_bar chi u b_y;
varexo eps sigma varrho eps_st;
parameters mu_k_bar mu mu_k beta rho_lambda vartheta alpha delta theta gamma ksi_k ksi_y mu_w zeta delta_der_u chi_k_bar chi_y_bar phi lambda_y_bar lambda_k_bar rho eta nu rho_st b_k;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Give your parameters values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
beta = 0.98;
delta = 0.015;
gy = 0.2;
alpha = 0.35;
gamma = 1 - 0.17/alpha;
zeta = 1;
theta = 1.4; 
vartheta = 1.25;
theta_bar = 0.8;
%u = 0.8;
mu = 1.1;
mu_w = 1.2;
mu_k = 1.15;
mu_k_bar = (mu_k*theta)/(theta*gamma + 1 - gamma);
chi_k_bar = 0.02;
chi_y_bar = 0.02;
phi = 0.98;
%chi_y_bar = 0.0202;
%chi_k_bar = 0.0304;
%phi = 0.99;
rho_lambda = 0.9;
ksi_k = 1;
%ksi_y = ksi_k*log((1 - phi)/chi_y_bar)/log((1 - phi)/chi_k_bar);
ksi_y = 0.6;
lambda_y_bar = 0.1;
lambda_k_bar = 0.1;
rho = 0.9;
eta = 0.9;
nu = 0.9;
rho_st = 0.85;
delta_der_u = (1/beta + delta - 1)/0.8;
b_k = (((mu_k-1)/mu_k)^(1-mu_k))*((lambda_k_bar/0.05)^(mu_k/rho_lambda))*((phi*beta*1*rho_lambda*0.05*(theta-1)*(1-gamma)/(theta*mu_k*(1-phi*beta*1*(1-0.05*(1-rho_lambda)))))^mu_k)*delta*mu_k/mu_k_bar;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define your model equations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model;
% Resource constraint
y/mu = c + g + i*(2*mu_k - 1)/mu_k + (z_k - a_k)*h_k + (z_y - a_y)*h_y; %1
% aggregate production function
y = x*((a_y)^(vartheta - 1))*((n_y)^(mu - 1))*((u*k(-1))^alpha)*(l^(1 - alpha)); %2
% capital accumulation
k = (1 - delta)*k(-1) + (mu_k_bar*i)/p_k; %3
% product amount
a_y(+1) = lambda_y*(z_y - a_y) + phi*a_y; %4
a_k(+1) = lambda_k*(z_k - a_k) + phi*a_k; %5
% stock of new technologies
z_y(+1) = (chi_y_bar*(chi^ksi_y) + phi)*z_y; %6
z_k(+1) = (chi_k_bar*(chi^ksi_k) + phi)*z_k; %7
% factor market equilibria
(1 - alpha)*y/l = mu*mu_w*(l^zeta)*c; %8
alpha*y/u = mu*delta_der_u*p_k*k(-1); %9
% new capital 
(p_k*j)/mu_k_bar = i; %10
(theta*i_e)/i_s = (1 - gamma)/gamma; %11
% euler equation
p_k = beta*Lambda(+1)*((alpha*y(+1))/(mu*k) + (1 - delta)*p_k(+1)); %12
Lambda(+1) = c/c(+1); %13
% optimal adoption
1 = phi*beta*(Lambda(+1)*(v_y(+1) - j_y(+1))*(lambda_y/h_y)*rho_lambda); %14
1 = phi*beta*(Lambda(+1)*(v_k(+1) - j_k(+1))*(lambda_k/h_k)*rho_lambda); %15
v_y = pi_y + phi*beta*(Lambda(+1)*v_y(+1)); %16
v_k = pi_k + phi*beta*(Lambda(+1)*v_k(+1)); %17
pi_k = ((theta - 1)/theta)*(1 - gamma)*(i/(a_k*mu_k)); %18
pi_y = ((vartheta - 1)/vartheta)*y/(a_y*mu); %19
j_y = -h_y + phi*beta*(Lambda(+1)*(lambda_y*v_y(+1) + (1 - lambda_y)*j_y(+1))); %20
j_k = -h_k + phi*beta*(Lambda(+1)*(lambda_k*v_k(+1) + (1 - lambda_k)*j_k(+1))); %21
lambda_y = lambda_y_bar*((a_y*h_y/o_y)^(rho_lambda)); %22
lambda_k = lambda_k_bar*((a_k*h_k/o_k)^(rho_lambda)); %23
((mu - 1)/mu)*(y/n_y) = o_y; %24
((mu_k - 1)/mu_k)*(i/n_k) = o_k; %25
p_k = mu_k*(n_k^(1 - mu_k))*(p_k_st^gamma)*(p_k_et^(1 - gamma)); %26
p_k_et = theta*(a_k^(1-theta)); %27
p_k_bar = (theta^(1 - gamma))*(a_k^((gamma - 1)*(theta - 1)))*(p_k_st^gamma); %28
i = theta*i_e + i_s; %29
% stochastic processes
log(chi) = rho*log(chi(-1)) + eps; %30
log(x) = log(x(-1)) + sigma; %31
log(g) = nu*log(g(-1)) + varrho; %32
log(p_k_st) = rho_st*log(p_k_st(-1)) + eps_st; %33
% operation costs
o_y = p_k_bar*k*b_y; %34
o_k = p_k_bar*k*b_k; %35
a_k*h_k/(p_k_bar*b_k) = a_y*h_y/(p_k_bar*b_y); % 36


end;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set initial values for some variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
initval;
Lambda = 1;
u = 0.8;
lambda_k = 0.05;
lambda_y = 0.05;
p_k_st = 1;
chi = 1;
y = 5;
g = y/5;
%x = 1;
i = (delta*alpha*y)/(mu_k_bar*mu*u*delta_der_u);
z_k = 1.56;
a_k = lambda_k*z_k/(1 - phi + lambda_k);
%z_k = a_k*(1 - phi + lambda_k)/lambda_k;
a_y = y*a_k*mu_k*theta*(vartheta - 1)/(vartheta*mu*(theta - 1)*(1 - gamma));
z_y = a_y*(1 - phi + lambda_y)/lambda_y;
p_k_bar = (theta^(1 - gamma))*(a_k^((1 - theta)*(1 - gamma)))*(p_k_st^gamma);
pi_k = (theta - 1)*(1 - gamma)*i/(theta*a_k*mu_k);
pi_y = (vartheta - 1)*y/(vartheta*a_y*mu);
v_k = pi_k/(1 - phi*beta*Lambda);
v_y = pi_y/(1 - phi*beta*Lambda);
h_k = ((1 - phi*beta*Lambda*(1 - lambda_k))*phi*beta*Lambda*rho_lambda*lambda_k*v_k - phi*beta*Lambda*rho_lambda*lambda_k*phi*beta*Lambda*lambda_k*v_k)/(1 - phi*beta*Lambda*(1 - lambda_k) - phi*beta*Lambda*rho_lambda*lambda_k);
h_y = ((1 - phi*beta*Lambda*(1 - lambda_y))*phi*beta*Lambda*rho_lambda*lambda_y*v_y - phi*beta*Lambda*rho_lambda*lambda_y*phi*beta*Lambda*lambda_y*v_y)/(1 - phi*beta*Lambda*(1 - lambda_y) - phi*beta*Lambda*rho_lambda*lambda_y);
b_y = b_k*a_y*h_y/(a_k*h_k);
k = (lambda_k_bar^(1/rho_lambda))*a_k*h_k/(p_k_bar*b_k*(lambda_k^(1/rho_lambda)));
o_k = p_k_bar*k*b_k;
o_y = p_k_bar*k*b_y;
n_y = (mu - 1)*y/(mu*o_y);
n_k = (mu_k - 1)*i/(mu_k*o_k);
p_k_et = theta*(a_k^(1-theta));
p_k = alpha*y/(u*mu*delta_der_u*k);
%p_k = mu_k*(n_k^(1 - mu_k))*(p_k_st^gamma)*(p_k_et^(1 - gamma)); 
j_k = (phi*beta*Lambda*rho_lambda*lambda_k*v_k/h_k - 1)/(phi*beta*Lambda*rho_lambda*lambda_k/h_k);
j_y = (phi*beta*Lambda*rho_lambda*lambda_y*v_y/h_y - 1)/(phi*beta*Lambda*rho_lambda*lambda_y/h_y);
i_s = gamma*i;
i_e = (i - i_s)/theta;
j = delta*k;
l = ((1 - alpha)*y/(mu*mu_w*(y/mu - g - i*(2*mu_k - 1)/mu_k - h_k*(z_k - a_k) - h_y*(z_y - a_y))))^(1/(zeta + 1));
c = (1 - alpha)*y/(mu_w*mu*(l^(zeta + 1)));
x = y*(a_y^(1-vartheta))*(n_y^(1-mu))*((u*k)^(-alpha))*(l^(alpha-1));


end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the steady state, check for existence
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
resid;
model_diagnostics;
steady(maxit = 5000);
check;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Add shocks
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
shocks;
var eps = 0.007;
var sigma = 0;
var varrho = 0;
var eps_st = 0;
end;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start simulation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
options_.debug=1;
stoch_simul(order=1,irf=20,periods = 200,nograph);