%-----------------------------------
%	NEOCLASSICAL GROWTH MODEL
% 	EXOGENOUS LABOR SUPPLY
%-----------------------------------

%-----------------------------------
%	THIS VERSION: SEPTEMBER 26, 2012
%-----------------------------------

% Load calibration parameters

load china_calib
load china_shocks

%-----------------------------------
%	ENDOGENOUS VARIABLES
%-----------------------------------

var 

c k y x rk w q sav;

%-----------------------------------
%	EXOGENOUS VARIABLES
%-----------------------------------

varexo gamma n tau g;

%-----------------------------------
%	PARAMETERS
%-----------------------------------

parameters

% Structural
beta delta sigma theta psi1k psi2k

% Steady state exogenous
tauss gammass nss gss;

beta  = disc_fact;
delta = depre;
sigma = risk_av;
theta = k_share;
psi1k = k_adjust;

tauss   = tau_ss;
gammass = gamma_ss;
nss     = n_ss;
gss     = g_ss;

psi2k   = gammass * nss;

%------------------------------------
%	MODEL
%------------------------------------

model;

% Euler equation

q = (beta / (gamma(+1) ^ sigma)) * (rk(+1) + (1 - delta) * q(+1)) * ((c(+1) / c) ^ (-sigma));

% Tobin's q
1 = (1 - (psi1k / 2) * (gamma * n * x / x(-1) - psi2k) ^ 2 
    - psi1k * (gamma * n * x / x(-1) - psi2k) * (gamma * n * x / x(-1))) * q
        + (beta / gamma(+1) ^ sigma) * q(+1) * ((c(+1) / c) ^ (-sigma)) 
            * psi1k * (gamma(+1) * n(+1) * x(+1) / x - psi2k) * (gamma(+1) * n(+1) * x(+1) / x)^ 2;

% Resource constraint
c + x = (1 - g) * y;

% Production function
y = k(-1) ^ theta;

% Investment

gamma(+1) * n(+1) * k = (1 - delta) * k(-1) 
    + (1 - (psi1k / 2) * (gamma * n * x / x(-1) - psi2k) ^ 2) * x;

% Rental rate
rk = theta * y / k(-1);

% Wage
w = (1 - theta) * y;

% Saving rate
sav = ((1 - g) * y - c - delta * k(-1)) / (y - delta * k(-1));

end;

%--------------------------------------
%	INITIAL STEADY STATE
%--------------------------------------

initval;

gamma = gammass;

tau = tauss;

n = nss;

g = gss;

k = 1.75;

rk = theta * k ^ (theta - 1);

y = k ^ theta;

x = (gammass * nss + delta - 1) * k;

q = 1;

c = (1 - gss) * y - x;

w = (1 - theta) * y;

sav = ((1 - g) * y - c - delta * k) / (y - delta * k);

end;

shocks(shocks_file=china_shocks);
var gamma;
periods 1;
values 0;
//periods 1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20	21	22	23	24	25	26	27	28	29	30	31	32	33	34	35	36	37	38	39	40	41	42	43	44	45	46	47;
//values 0.902406029	1.113754768	1.054179609	1.047562191	1.07711342	0.787755928	0.901510405	1.105294212	1.234291267	1.000297531	0.949254704	1.054263193	0.957315718	1.058176281	0.940491805	1.02200402	1.111848485	1.067795246	1.044603122	1.02192917	1.103050832	1.078122667	1.116190458	1.094822672	1.068219284	1.080650654	1.095327344	0.979105657	1.014752332	1.077908874	1.135768077	1.185161873	1.113924783	1.088719964	1.081714901	1.058077087	1.04125179	1.041481825	1.037296603	1.067881836	1.09053663	1.099193253	1.08476273	1.134066119	1.13375115	1.090561595	1.069262496;
end;

simul(periods=500);