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


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


% var ty tc tg p_k tj ti tzk tzy tvk tvy n_y tak tay thk thy x tk l tok toy tpik tpiy tpks tpke n_k tjy tjk Lambda tie tis lambda_y lambda_k p_k_bar chi u b_k ttfp;
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_k tfp ty tc tg tj ti tis tie toy tok thk thy tzk tzy tak tay tpke tvk tvy tpik tpiy tjk tjy tpks ttfp tk;

varexo eps sigma varrho eps_st;

parameters kappa gchi 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 rho eta nu rho_st b_y gy gay gak gzy gzk gvk gjk gvy gjy gpke lambda_y_bar lambda_k_bar;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Give your parameters values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
beta = 0.98;
delta = 0.015;
alpha = 0.35;
gamma = 1 - 0.17/alpha;
zeta = 1;
theta = 1.4; 
vartheta = 1.25;
theta_bar = 0.8;
mu = 1.1;
mu_w = 1.2;
mu_k = 1.15;
% mu_k_bar = (mu_k*theta)/(theta*gamma + 1 - gamma);
phi = 0.99;
%chi_k_bar = 0.0304/4;
rho_lambda = 0.95;
ksi_k = 1;
ksi_y = 0.6;
% lambda_y_bar = 1.041388535907322;
% lambda_k_bar = 1.041388535907322;
lambda_k_bar = 1.067688464348893;
lambda_y_bar = 1.067688464348893;
rho = 1;
eta = 0.9;
nu = 0.9;
rho_st = 1;
gy = 1.006;
delta_der_u = (gy/beta + delta - 1)/0.8;
b_y = 1;
% 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;
gpke = 0.99125;
gak = (gpke)^(1/(1-theta));
gay = (gy)^((1-alpha)/(vartheta-1));
gzk = gak;
gzy = gay;
ghk = gy/gak;
ghy = gy/gay;
gvk = ghk;
gvy = ghy;
gjk = ghk;
gjy = ghy;
%chi_k_bar = (gzk-phi)*(((gzy-phi)/(chi_y_bar))^(-ksi_k/ksi_y));
gchi = 1;
%chi_y_bar = 0.0202/4;
%chi_k_bar = ((gzk-phi)*(((gzy-phi)/chi_y_bar)^(-ksi_k/ksi_y)))/4;
chi_y_bar = gzy - phi;
chi_k_bar = gzk - phi;
kappa = (1-alpha)/(vartheta-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define your model equations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model;
% Resource constraint
ty = tc + tg + p_k*tj/mu_k + (mu-1)*ty/(mu*vartheta) + (mu_k - 1)*ti + (tzk(-1) - tak(-1)*(theta^(1/(1-theta))))*thk + (tzy(-1) - tay(-1))*thy; %1
% aggregate production function
ty = x*((tay(-1))^(vartheta-1))*((n_y)^(mu-1))*((u*tk(-1))^alpha)*(l^(1-alpha));%2
% capital accumulation
gy*tk = (1 - delta)*tk(-1) + tj; %3
% product amount
gak*tak = lambda_k*(tzk(-1)*(theta^(1/(theta-1))) - tak(-1)) + phi*tak(-1); %4
(gy^kappa)*tay = lambda_y*(tzy(-1) - tay(-1)) + phi*tay(-1); %5
% stock of new technologies
gak*tzk = (chi_k_bar*(chi^ksi_k) + phi)*tzk(-1); %6
gay*tzy = (chi_y_bar*(chi^ksi_y) + phi)*tzy(-1); %7
% factor market equilibria
(1-alpha)*ty = vartheta*mu*mu_w*tc*(l^(zeta+1)); %8
alpha*ty = vartheta*mu*delta_der_u*u*p_k*tk(-1); %9
% new capital 
(p_k*tj)/mu_k = theta*tie + tis; %10
(theta*tie)/tis = (1 - gamma)/gamma; %11
% euler equation
(beta*Lambda(+1)/p_k)*((alpha*ty(+1))/(mu*vartheta*tk) + (1 - delta)*p_k(+1)) = 1; %12
Lambda = tc/(gy*tc(+1)); %13
% optimal adoption
1 = phi*beta*rho_lambda*Lambda(+1)*(lambda_k*gy*(tvk(+1) - tjk(+1))/(thk*gak)); %14
1 = phi*beta*rho_lambda*Lambda(+1)*(lambda_y*(gy^(1-kappa))/thy)*(tvy(+1) - tjy(+1)); %15
tjk = -thk + gy*phi*beta*Lambda(+1)*(lambda_k*tvk(+1) + (1 - lambda_k)*tjk(+1))/gak; %16
tjy = -thy + (gy^(1-kappa))*phi*beta*Lambda(+1)*(lambda_y*tvy(+1) + (1 - lambda_y)*tjy(+1)); %17
lambda_k = lambda_k_bar*((tak(-1)*thk*(theta^(1/(1-theta)))/tok)^rho_lambda); %18
lambda_y = lambda_y_bar*((tay(-1)*thy/toy)^rho_lambda); %19
tvk = tpik + gy*phi*beta*Lambda(+1)*tvk(+1)/gak; %20
tvy = tpiy + (gy^(1-kappa))*phi*beta*Lambda(+1)*tvy(+1); %21
tpik = ((theta-1)/theta)*(1-gamma)*ti*(theta^(1/(1-theta)))/tak(-1); %22
tpiy = ((vartheta-1)/(mu*vartheta))*ty/tay(-1); %23
(mu-1)*ty/(mu*vartheta*n_y) = toy; %24
(mu_k - 1)*ti/n_k = tok; %25
% operation costs
tok = b_k*p_k_bar*tk(-1); %26
toy = b_y*p_k_bar*tk(-1); %27
p_k = mu_k*((n_k)^(1-mu_k))*(tpks^gamma)*(tpke^(1-gamma)); %28
tpke = theta*((tak(-1))^(1-theta)); %29
p_k_bar = (tpke^(1-gamma))*(tpks^gamma); %30
% stochastic processes
log(tg) =  log(tg(-1)) + sigma; %31
log(tpks) = rho_st*log(tpks(-1)) + eps_st; %32
log(chi) = rho*log(chi(-1)) + eps; %33
log(x) = log(x(-1)) + varrho; %34
% investments vs new capital
p_k*tj/mu_k = ti; %35
% equation for b_k
% b_k = ((b_y*thk*(theta^(1/(1-theta))))/(thy))*((((gy^kappa) - phi)*(tzk*(theta^(1/(theta-1)))-tak(-1))*(tay(-1)^(1-rho_lambda)))/((tzy(-1)-tay(-1))*(gak-phi)*(tak(-1)^(1-rho_lambda))))^rho_lambda;
b_k = (b_y*tak(-1)*thk*(theta^(1/(1-theta)))/(tay(-1)*thy));
% tfp
ttfp = x*(n_y^(mu-1))*(tay(-1)^(vartheta-1)); %37
% equations for all stationarized tilde variables
ty = 1.006;
%ty = y/y(-1); %38
tc = c/y(-1); %39
tg = 0.2012;
%tg = g/y(-1); %40
tj = j/y(-1); %41
ti = i/y(-1); %42
tis = i_s/y(-1); %43
tie = i_e/y(-1); %44
toy = o_y/y(-1); %45
tok = o_k/y(-1); %46
thk = h_k*a_k(-2)/y(-1); %47
thy = h_y/(y(-1)^(1 - kappa)); %48
tzk(-1) = z_k(-1)/a_k(-2); %49
tzy(-1) = z_y(-1)/a_y(-2); %50
tak(-1) = (theta/gpke)^(1/(theta-1)); %51
tay(-1) = a_y(-1)/(y(-1)^kappa); %52
tpke = 0.99125; %53
tvk = v_k*a_k(-2)/y(-1); %54
tvy = v_y/(y(-1)^(1-kappa)); %55
tpik = pi_k*a_k(-2)/(y(-1)); %56
tpiy = pi_y/(y(-1)^(1-kappa)); %57
tjy = j_y/(y(-1)^(1-kappa)); %58
tjk = j_k*a_k(-2)/(y(-1)); %59
tpks = p_k_st/(p_k_et^((gamma-1)/gamma)); %60
ttfp = tfp/(y(-1)^(1-alpha)); %61
tk = k/y; % 62

end;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set initial values for some variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
initval;
ty = 1.006; % 1
tg = 0.2012; % 2
tpke = 0.99125; % 3
tak = (theta/gpke)^(1/(theta-1)); % 4
Lambda = 1/gy; % 5
u = 0.8; % 6
lambda_k = 0.05; % 7
lambda_y = 0.05; % 8
% p_k_st = 1;
%chi = ((chi_k_bar)/(chi_y_bar))^(1/(ksi_y - ksi_k));
chi = 1; % 9
y = 5; % 10 
g = y/5; % 11
ti = ((delta + gy - 1)*alpha*ty)/(vartheta*mu_k*mu*u*delta_der_u); % 12
tzk = (tak*(gak - phi + lambda_k))/(lambda_k*(theta^(1/(theta-1)))); % 13
%z_k = 1;
%a_k = lambda_k*z_k/(gak - phi + lambda_k);
tay = ty*tak*theta*(vartheta - 1)/(vartheta*mu*ti*(theta - 1)*(1 - gamma)*(theta^(1/(1-theta)))); % 14
tzy = tay*(gy^kappa - phi + lambda_y)/lambda_y; % 15
tpik = (theta - 1)*(1 - gamma)*ti*(theta^(1/(1-theta)))/(theta*tak); % 16
tpiy = (vartheta - 1)*ty/(vartheta*tay*mu); % 17
tvk = gak*tpik/(gak - phi*beta*Lambda*gy); % 18
tvy = tpiy/(1 - phi*beta*Lambda*(gy^(1-kappa))); % 19
thk = (phi*beta*Lambda*lambda_k*tvk*rho_lambda*gy*(phi*beta*Lambda*gy-gak))/(gak*(-gak + phi*beta*gy*Lambda*(1-lambda_k + rho_lambda*lambda_k))); % 20
thy = (phi*beta*Lambda*lambda_y*tvy*rho_lambda*(gy^(1-kappa))*(phi*beta*Lambda*(gy^(1-kappa))-1))/(phi*beta*Lambda*(gy^(1-kappa))*(1 + lambda_y*(rho_lambda-1))-1); % 21
% h_k = (phi*beta*Lambda*lambda_k*gvk*rho_lambda*v_k*(phi*beta*Lambda*gjk - 1))/(phi*beta*Lambda*gjk*(1+lambda_k*(rho_lambda - 1))-1);
% h_y = (phi*beta*Lambda*lambda_y*gvy*rho_lambda*v_y*(phi*beta*Lambda*gjy - 1))/(phi*beta*Lambda*gjy*(1+lambda_y*(rho_lambda - 1))-1);
% b_k = (((alpha*y/(u*mu*vartheta*delta_der_u*mu_k*(mu_k-1)*i))^(-1/mu_k))*((lambda_k_bar/lambda_y_bar)^(-1/rho_lambda))*a_k*h_k)^(mu_k/(1-mu_k));
% b_k = (b_y*tak*thk*(theta^(1/(1-theta)))/(tay*thy))*((lambda_k_bar/lambda_y_bar)^(1/rho_lambda)); % 22
% b_k = ((b_y*thk*(theta^(1/(1-theta))))/(thy))*((((gy^kappa) - phi)*(tzk*(theta^(1/(theta-1)))-tak)*(tay^(1-rho_lambda)))/((tzy-tay)*(gak-phi)*(tak^(1-rho_lambda))))^rho_lambda;
b_k = (b_y*tak*thk*(theta^(1/(1-theta)))/(tay*thy));
% p_k_et = theta*(a_k^(1-theta));
tpks = (gpke)^((gamma-1)/gamma); % 23
p_k_bar = (tpke^(1 - gamma))*(tpks^gamma); % 24
tk = ((alpha*ty/(u*mu*vartheta*delta_der_u*mu_k))^(1/mu_k))*(((mu_k-1)*ti/b_k)^((mu_k-1)/mu_k))/p_k_bar; % 25
%k = (lambda_k_bar^(1/rho_lambda))*a_k*h_k/(p_k_bar*b_k*(lambda_k^(1/rho_lambda)));
n_k = (mu_k-1)*ti/(b_k*p_k_bar*tk); % 26
toy = p_k_bar*tk*b_y; % 27
n_y = (mu - 1)*ty/(mu*toy*vartheta); % 28
p_k = alpha*ty/(u*mu*vartheta*delta_der_u*tk); % 29
tok = p_k_bar*tk*b_k; % 30
tjk = (-thk*gak + phi*beta*Lambda*lambda_k*tvk*gy)/(gak - phi*beta*Lambda*(1-lambda_k)*gy); % 31
tjy = (-thy + phi*beta*Lambda*lambda_y*tvy*(gy^(1-kappa)))/(1 - phi*beta*Lambda*(1-lambda_y)*(gy^(1-kappa))); % 32
tis = gamma*ti; % 33
tie = (ti - tis)/theta; % 34
tj = (gy - 1 + delta)*tk; % 35
l = (((1-alpha)*ty)/(mu*mu_w*vartheta*(ty - tg - p_k*tj/mu_k - (mu-1)*ty/(mu*vartheta) - (mu_k - 1)*ti - thk*(tzk - tak*(theta^(1/(1-theta)))) - thy*(tzy - tay))))^(1/(zeta+1)); % 36
tc = (1 - alpha)*ty/(mu_w*mu*vartheta*(l^(zeta + 1))); % 37
x = ty*(tay^(1-vartheta))*(n_y^(1-mu))*((u*tk)^(-alpha))*(l^(alpha-1)); % 38
ttfp = x*(n_y^(mu - 1))*(tay^(vartheta - 1)); % 39
c = tc*y; % 40
j = tj*y; % 41
i = ti*y; % 42
i_s = tis*y; % 43
i_e = tie*y; % 44
o_y = toy*y; % 45
o_k = tok*y; % 46
v_y = tvy*(y^(1-kappa)); % 47
pi_y = tpiy*(y^(1-kappa)); % 48
j_y = tjy*(y^(1-kappa)); % 49
h_y = thy*(y^(1-kappa)); % 50
a_y = tay*(y^kappa); % 51
z_y = tzy*a_y; % 52
tfp = ttfp*(y^(1-alpha)); % 53
a_k = tak; % 54
v_k = tvk*y/a_k; % 55
h_k = thk*y/a_k; % 56
j_k = tjk*y/a_k; % 57
pi_k = tpik*y/a_k; % 58
p_k_et = theta*(a_k^(1-theta)); % 59
p_k_st = tpks*(p_k_et^((gamma-1)/gamma)); % 60
k = tk*y; % 61
z_k = tzk*a_k; % 62





end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the steady state, check for existence
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options_.Schur_vec_tol = 1e-6;
resid;
% model_diagnostics;
steady(maxit = 2000);
check(qz_zero_threshold=1e-15);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Add shocks
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
shocks;
var eps = 0.049;
% var sigma = 0.07;
% var varrho = 0.07;
% var eps_st = 0.07;
end;



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