%Kung(2014)
var U_t u_t expUsigma gammaN_t C_t L_t SDF_t K_t LambdaK_t LambdaN_t
LambdaPrimeK_t LambdaPrimeN_t qK_t qN_t rK_t rN_t I_t S_t PP Pi_t w_t
Y_t A_t mc_t R_t;


varexo epsA; 

parameters beta psi gamma nu alpha deltaK deltaN kappaP zetaK eta zetaN rhoA
seA rho_r phi_p phi_y tau utilityconstant alphaK1 alphaN1 alphaK2 alphaN2 

A_ss Pi_ss L_ss mc_ss qK_ss qN_ss Y_ss w_ss rN_ss gammaN_ss SDF_ss R_ss rK_ss K_ss 
I_ss S_ss C_ss  U_ss u_ss LambdaK_ss LambdaN_ss LambdaPrimeK_ss LambdaPrimeN_ss; 


beta=0.997;       %discount factor
psi=2;            %IES
gamma=10;         %Risk Aversion
nu=6;             %Demand Elasticity
alpha=0.33;       %Capital Share
deltaK=0.02;      %Capital Depreciation
deltaN=0.0375;    %R&D Depreciation
kappaP=30;        %Price Adjustment Costs

zetaK=4.8;        %Capital Adjustment Costs
eta=0.1;          %Degree of technological spillovers
zetaN=3.3;        %R&D Adjustment Costs

rhoA=0.983;       %Persistence of TFP shock  
seA=0.012;        %TFP volatility

rho_r=0.7;     
phi_p=1.5;
phi_y=0.1;

%-----------------------
%  Steady State Values
%-----------------------
A_ss=1;
Pi_ss=1;
L_ss=0.3;
mc_ss=(nu-1)/nu;
qK_ss=1;
qN_ss=1;


ksssolution=fsolve(@(kssguess) alpha*kssguess^(alpha-1)-eta*(1-alpha)*kssguess^alpha-(deltaK-deltaN)/(mc_ss*L_ss^(1-alpha)),1,optimoptions('fsolve','Display','off','OptimalityTolerance', 1e-16));
K_ss=ksssolution;
Y_ss=K_ss^alpha*L_ss^(1-alpha);

w_ss=(1-alpha)*mc_ss*Y_ss/L_ss;

rN_ss=eta*(1-alpha)*mc_ss*Y_ss;
gammaN_ss=(beta*(rN_ss+1-deltaN))^(psi);
SDF_ss=beta*gammaN_ss^(-1/psi);
R_ss=Pi_ss/SDF_ss;
rK_ss=1/SDF_ss-(1-deltaK);
I_ss=K_ss*(gammaN_ss-(1-deltaK));
S_ss=gammaN_ss-(1-deltaN);
C_ss=Y_ss-I_ss-S_ss;

tau=w_ss*(1-L_ss)/C_ss;
U_ss=1;
u_ss=C_ss*(1-L_ss)^tau;
utilityconstant=((1-beta)*u_ss^(1-1/psi)+beta*gammaN_ss^(1-1/psi))^(1/(1/psi-1));

LambdaK_ss=gammaN_ss-(1-deltaK);
LambdaN_ss=gammaN_ss-(1-deltaN);
LambdaPrimeK_ss=1;
LambdaPrimeN_ss=1;

alphaK1=(I_ss/K_ss)^(1/zetaK);
alphaN1=(S_ss)^(1/zetaN);

alphaK2=(I_ss/K_ss)/(1-zetaK);
alphaN2=(S_ss)/(1-zetaN);




model;

%eq1
U_t=utilityconstant*((1-beta)*u_t^(1-1/psi)+beta*expUsigma^((1-1/psi)/(1-gamma))*gammaN_t^(1-1/psi))^(1/(1-1/psi));

%eq2
expUsigma=U_t(+1)^(1-gamma);

%eq3
u_t=exp(C_t)*(1-exp(L_t))^tau;

%eq4
#eUlag=expUsigma(-1)^(1/(1-gamma));
SDF_t=beta*(u_t/u_t(-1))^(1-1/psi)*(exp(C_t(-1))/exp(C_t))*gammaN_t(-1)^(-1/psi)*(U_t/eUlag)^(1/psi-gamma);

%eq5
exp(K_t)*gammaN_t=(1-deltaK+LambdaK_t)*exp(K_t(-1));

%eq6
gammaN_t=(1-deltaN)+LambdaN_t;

%eq7
exp(w_t)=tau*exp(C_t)/(1-exp(L_t));

%eq8
1=R_t*SDF_t(+1)/Pi_t(+1);

%eq9
1=exp(qK_t)*LambdaPrimeK_t;

%eq10
1=exp(qN_t)*LambdaPrimeN_t;

%eq11
exp(qK_t)=SDF_t(+1)*(rK_t(+1)+exp(qK_t(+1))*
(1-deltaK-LambdaPrimeK_t*exp(I_t(+1))/exp(K_t)+LambdaK_t(+1)));

%eq12
exp(qN_t)=SDF_t(+1)*(rN_t(+1)+exp(qN_t(+1))*
(1-deltaN-LambdaPrimeN_t*exp(S_t(+1))+LambdaN_t(+1)));

%eq13
PP=Pi_t/Pi_ss;

%eq14
(1-nu)+nu*exp(mc_t)=kappaP*(PP-1)*PP-SDF_t(+1)*(PP(+1)-1)*PP(+1)*(exp(Y_t(+1))/exp(Y_t))*gammaN_t;

%eq15
exp(w_t)=(1-alpha)*exp(mc_t)*exp(Y_t)/exp(L_t);

%eq16
rK_t=alpha*exp(mc_t)*exp(Y_t)/exp(K_t(-1));

%eq17
rN_t=eta*(1-alpha)*exp(mc_t)*exp(Y_t);

%eq18
R_t/R_ss=(R_t(-1)/R_ss)^rho_r*((PP^phi_p)*((exp(Y_t)/exp(Y_t(-1)))^phi_y))^(1-rho_r);

%eq19
exp(Y_t)=exp(C_t)+exp(I_t)+exp(S_t)+0.5*kappaP*(PP-1)^2*exp(Y_t);

%eq20
exp(Y_t)=exp(K_t(-1))^alpha*(exp(A_t)*exp(L_t))^(1-alpha);

%eq21
LambdaK_t=(alphaK1/(1-1/zetaK))*(exp(I_t)/exp(K_t(-1)))^(1-1/zetaK)+alphaK2;

%eq22
LambdaN_t=(alphaN1/(1-1/zetaN))*exp(S_t)^(1-1/zetaN)+alphaN2;

%eq23
LambdaPrimeK_t=alphaK1*(exp(I_t)/exp(K_t(-1)))^(-1/zetaK);

%eq24
LambdaPrimeN_t=alphaN1*exp(S_t)^(-1/zetaN);

%eq25
A_t=rhoA*A_t(-1)+seA*epsA;


end;

initval;
U_t=U_ss; 
u_t=u_ss; 
expUsigma=U_ss^(1-gamma); 
gammaN_t=gammaN_ss; 
C_t=log(C_ss); 
L_t=log(L_ss); 
SDF_t=SDF_ss; 
K_t=log(K_ss); 
LambdaK_t=LambdaK_ss; 
LambdaN_t=LambdaN_ss;
LambdaPrimeK_t=LambdaPrimeK_ss; 
LambdaPrimeN_t=LambdaPrimeN_ss; 
qK_t=log(qK_ss); 
qN_t=log(qN_ss); 
rK_t=rK_ss; 
rN_t=rN_ss; 
I_t=log(I_ss); 
S_t=log(S_ss); 
PP=1; 
Pi_t=Pi_ss; 
R_t=R_ss;
Y_t=log(Y_ss); 
A_t=log(A_ss); 
mc_t=log(mc_ss);
w_t=log(w_ss);
end;

steady;
resid;
shocks;
var epsA=1;

end;

check;
stoch_simul(order = 1, noprint, nomoments, irf = 20);