%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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;
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 u delta_der_u rho eta nu rho_st chi_k_bar chi_y_bar phi lambda_y_bar lambda_k_bar;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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_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);
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)/u;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
j = (mu_k_bar*i)/p_k; %12
% euler equation
p_k = beta*Lambda(+1)*((alpha*y(+1))/(mu*k) + (1 - delta)*p_k(+1)); %13
Lambda(+1) = c/c(+1); %14
% optimal adoption
1 = phi*beta*Lambda(+1)*(v_y(+1) - j_y(+1))*(lambda_y/h_y)*rho_lambda; %15
1 = phi*beta*Lambda(+1)*(v_k(+1) - j_k(+1))*(lambda_k/h_k)*rho_lambda; %16
v_y = pi_y + phi*beta*(Lambda(+1)*v_y(+1)); %17
v_k = pi_k + phi*beta*(Lambda(+1)*v_k(+1)); %18
pi_k = ((theta - 1)/theta)*(1 - gamma)*(i/(a_k*mu_k)); %19
pi_y = ((vartheta - 1)/vartheta)*y/(a_y*mu); %20
j_y = -h_y + phi*beta*Lambda(+1)*(lambda_y*v_y(+1) + (1 - lambda_y)*j_y(+1)); %21
j_k = -h_k + phi*beta*Lambda(+1)*(lambda_k*v_k(+1) + (1 - lambda_k)*j_k(+1)); %22
lambda_y = lambda_y_bar*((a_y*h_y/o_y)^(rho_lambda)); %23
lambda_k = lambda_k_bar*((a_k*h_k/o_k)^(rho_lambda)); %24
((mu - 1)/mu)*(y/n_y) = o_y; %25
((mu_k - 1)/mu_k)*(i/n_k) = o_k; %26
p_k = mu_k*(n_k^(1 - mu_k))*(p_k_st^gamma)*(p_k_et^(1 - gamma)); %27
p_k_et = theta*(a_k^(1-theta)); %28
p_k_bar = (theta^(1 - gamma))*(a_k^((gamma - 1)*(theta - 1)))*(p_k_st^gamma); %29
i = theta*i_e + i_s; %30
% stochastic processes
log(chi) = rho*log(chi(-1)) + eps; %31
log(x) = eta*log(x(-1)) + sigma; %32
log(g) = nu*log(g(-1)) + varrho; %33
log(p_k_st) = rho_st*log(p_k_st(-1)) + eps_st; %34




end;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set initial values for some variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
initval;
Lambda = 1;
k = 0.5;
a_y = (u*delta_der_u*(vartheta - 1)*k)/(alpha*vartheta);
a_k = (delta*(theta - 1)*(1 - gamma)*k)/(mu_k_bar*theta*mu_k);
z_k = a_k*(1 - phi + 0.05)/0.05;
z_y = a_y*(1 - phi + 0.05)/0.05;
n_k = (delta*(mu_k - 1)*k)/(mu_k_bar*mu_k);
n_y = (u*delta_der_u*(mu - 1)*k)/(alpha);
lambda_k = (1 - phi)*a_k/(z_k - a_k);
lambda_y = (1 - phi)*a_y/(z_y - a_y);
chi = (chi_k_bar/chi_y_bar)^(1/(ksi_y - ksi_k));
p_k_et = theta*((a_k)^(1 - theta));
j = delta*k;
i = (phi*beta*Lambda*rho_lambda*lambda_k*mu_k*a_k*n_k)*((lambda_k_bar/lambda_k)^(1/rho_lambda))/((1 - phi*beta*Lambda*(1 - lambda_k))*(mu_k - 1));
o_k = (mu_k - 1)*i/(mu_k*n_k);
h_k = ((lambda_k/lambda_k_bar)^(1/rho_lambda))*(o_k/a_k);
%pi_k = h_k*(1 - phi*beta*Lambda)*(1/(phi*beta*Lambda) + (1 - rho_lambda)/(rho_lambda*(1 - phi*beta*Lambda)));
pi_k = (theta - 1)*(1 - gamma)*i/(a_k*mu_k*theta);
v_k = pi_k/(1 - phi*beta*Lambda);
%j_k = (h_k*(1/rho_lambda - 1))/(1 - phi*beta*Lambda);
j_k = (-h_k + phi*beta*Lambda*lambda_k*v_k)/(1 - phi*beta*Lambda*(1 - lambda_k));
y = (mu_k - 1)*i*mu*n_y/(mu_k*n_k*(mu - 1));
p_k = (mu_k_bar*i*(1 - beta*Lambda*(1 - delta)))/(j*beta*Lambda*u*delta_der_u);
p_k_st = (p_k*(n_k^(mu_k - 1))*(p_k_et^(gamma - 1))/mu_k)^(1/gamma);
pi_y = (vartheta - 1)*y/(vartheta*a_y*mu);
o_y = (mu - 1)*y/(mu*n_y);
v_y = pi_y/(1 - phi*beta*Lambda);
%h_y = (phi*beta*Lambda*lambda_y*v_y)/((1 - rho_lambda)*(1 - lambda_y) + 1);
h_y = ((lambda_y/lambda_y_bar)^(1/rho_lambda))*(o_y/a_y);
j_y = (-h_y + phi*beta*Lambda*lambda_y*v_y)/(1 - phi*beta*Lambda*(1 - lambda_y));
%j_y = (h_y*(1 - rho_lambda))/(rho_lambda*(1 - phi*beta*Lambda));
i_s = gamma*i;
i_e = (i - i_s)/theta;
p_k_bar = (p_k_et^(1 - gamma))*(p_k_st^gamma);
g = y/5;
c = y/mu - g - (2*mu_k - 1)*i/mu_k - (z_y - a_y)*h_y - (z_k - a_k)*h_k;
l = ((1 - alpha)*y/(mu*mu_w*c))^(1/(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;
steady(maxit = 2000);
check;


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



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