%----------------------------------------------------------------
% 0. Housekeeping (close all graphic windows)
%----------------------------------------------------------------

close all;


%----------------------------------------------------------------
% 1. Defining variables
%----------------------------------------------------------------

var
    
ce ch c b inv y k n
omegah omegae
qL qk w
Le Lh LF
dce dch
dc dinv dqL dq db % 观测变量
a fi z q theta; % 15 variables and 5 shocks

varexo
epsa epsf epsz epsq epst; 

parameters
beta_h
beta_e
delta
alpha
n_bar
theta_bar
%Lh_Le_bar
%miu
%by_bar
LV_Lh_LF_bar % 住宅-生产房产面积比值
LV_LV_Lh_bar % 住房空置率

% 估计参数
gamma_h
gamma_e
eta
v
phi
fi_bar
phik % 资本调整成本系数

% 外生参数
rhoa
rhof
rhoz
rhoq
rhot;

%----------------------------------------------------------------
% 2. Set the parameters
%---------------------------------------------------------------

beta_h = 0.99;
beta_e = 0.98;
n_bar = 0.25;
delta = 0.025;
alpha = 0.5;

LV_Lh_LF_bar = 8.05; % 1998Q1_2015Q4
LV_LV_Lh_bar = 0.224; % 空置率

%miu = 0.6433;

%theta_bar = 0.5415; % phi=0.05
%theta_bar = 0.4847; % phi=0.10
theta_bar = 0.4047; % phi=0.20

%theta_bar = 0.5670; % phi=0.0321
%theta_bar = 0.2;  
eta = 1.02251; % 产出的长期稳态增长率1996Q1_2015Q4

%by_bar = 5.12; % 贷款产出比

% 估计参数
gamma_h =  0.5;
gamma_e =  0.5;
phi     = 0.20;  % 土地生产指数
fi_bar  = 0.90; %0.07~0.09
v       =    1; % 劳动供给弹性倒数
phik    =    1; % 资本调整成本系数

%外生参数
rhoa = 0.6;
rhof = 0.6;
rhoz = 0.6;
rhoq = 0.6;
rhot = 0.6;

%----------------------------------------------------------------
% 3. Model
%---------------------------------------------------------------

model(linear);

%-----------------------------------------------------------
% Compute the steady state
%----------------------------------------------------------
# Lh_LV_bar = 1/LV_LV_Lh_bar-1;
# LV_LF_bar = LV_Lh_LF_bar*LV_LV_Lh_bar;
# LF_LV_bar = 1/LV_LF_bar;
# LV_Le_bar = 1/(1+LF_LV_bar);
# miu = LV_Le_bar;

# LF_Lh_bar = 1/LV_LF_bar/Lh_LV_bar;
# LV_Lh_bar = 1/Lh_LV_bar;

# Lh_Le_bar = 1/(LF_Lh_bar+LV_Lh_bar);


# R_bar = eta/beta_h;
# omegab_omegae_bar = (beta_h-beta_e)/eta;
# ik_bar = 1-(1-delta)/eta;
# yk_bar = (eta/beta_e*(1-omegab_omegae_bar*theta_bar)-(1-delta))/alpha/(1-phi)/eta;
# qL_Le_y_bar = alpha*phi/((1-omegab_omegae_bar*theta_bar*eta)/beta_e-1);
# qL_Lh_y_bar = qL_Le_y_bar*Lh_Le_bar;
# omegah_ch_bar = (1-beta_h*gamma_h/eta)/(1-gamma_h/eta);
# by_bar = theta_bar*eta*qL_Le_y_bar+theta_bar/yk_bar;
% theta_bar = by_bar/(eta*qL_Le_y_bar+1/yk_bar);
# omegah_qL_Lh_bar = fi_bar/(1-beta_h);
# omegah_y_bar = omegah_qL_Lh_bar/qL_Lh_y_bar;
# ch_y_bar = omegah_ch_bar/omegah_y_bar;
# cy_bar = 1-ik_bar/yk_bar;
# ce_y_bar = cy_bar-ch_y_bar;
# psi = omegah_y_bar*(1-alpha)/n_bar^(1+v);
# mh = (1-beta_h*gamma_h/eta)*(1-gamma_h/eta);
# me = (1-beta_e*gamma_e/eta)*(1-gamma_e/eta);

%(1)家户关于消费的一阶条件
mh*omegah = (1-beta_h*gamma_h/eta*rhoa)*(1-gamma_h/eta)*a + gamma_h/eta*ch(-1) - (1+beta_h*gamma_h^2/eta^2)*ch + beta_h*gamma_h/eta*ch(+1);

%(2)家户关于劳动的一阶条件
w + omegah = a + v*n;

%(3)家户关于房地产的一阶条件
qL = fi_bar/omegah_qL_Lh_bar*(a + fi - omegah-Lh)+beta_h*(omegah(+1)-omegah+qL(+1));

%(4)企业家关于消费的一阶条件
me*omegae = gamma_e/eta*ce(-1)-(1+beta_e*gamma_e^2/eta^2)*ce+beta_e*gamma_e/eta*ce(+1);

%(5)企业家关于劳动的一阶条件
w = y-n;

%(6)企业家关于投资的一阶条件
qk + q = (1+beta_e)*phik*eta^2*inv - phik*eta^2*inv(-1) - beta_e*phik*eta^2*inv(+1);

%(7)企业家关于资本存量的一阶条件
qk=(beta_e/eta*(alpha*(1-phi)*eta*yk_bar+1-delta)-omegab_omegae_bar*theta_bar*beta_e/eta/(1/R_bar-beta_e/eta))*(omegae(1)-omegae)+beta_e/eta*alpha*(1-phi)*eta*yk_bar*(y(1)-k)+(beta_e/eta*(1-delta)+theta_bar*omegab_omegae_bar)*qk(1)-omegab_omegae_bar*theta_bar/R_bar/(1/R_bar-beta_e/eta)*(omegah-omegah(1))+omegab_omegae_bar*theta_bar*theta;

%(8)企业家关于房地产的一阶条件
qL=(beta_e*(alpha*phi/qL_Le_y_bar+1)-omegab_omegae_bar*theta_bar*beta_e/(1/R_bar-beta_e/eta))*(omegae(1)-omegae)+beta_e*alpha*phi/qL_Le_y_bar*(y(1)-Le)+(beta_e+theta_bar*omegab_omegae_bar*eta)*qL(1)+eta*omegab_omegae_bar*theta_bar/R_bar/(1/R_bar-beta_e/eta)*(omegah(1)-omegah)+eta*omegab_omegae_bar*theta_bar*theta;


%(9)生产函数
y = z + alpha*phi*LF(-1)+alpha*(1-phi)*k(-1)+(1-alpha)*n;

%(10)资本运动方程    
k-(1-delta)/eta*k(-1)=ik_bar*inv;

%(11)信贷约束条件
b=theta+eta*qL_Le_y_bar*yk_bar/(eta*qL_Le_y_bar*yk_bar+1)*(qL(1)+Le)+1/(eta*qL_Le_y_bar*yk_bar+1)*(qk(1)+k);

%(12)企业月算约束条件
ce_y_bar*ce+qL_Le_y_bar*(Le-Le(-1))+ik_bar/yk_bar*(inv-q)+(1-alpha)*(w+n)+1/eta*by_bar*b(-1)=y+1/R_bar*by_bar*(b-omegah+omegah(1));

%(13)总体资源约束
cy_bar*c+ik_bar/yk_bar*(inv-q)=y;

%(14)土地约束
Lh_Le_bar*Lh+Le=0;

%(15)总消费
cy_bar*c = ch_y_bar*ch+ce_y_bar*ce;

%(16)生产性房产
LF=Le-eta*miu/(1-miu)*(qL-qL(-1));

% shock
a=rhoa*a(-1)+epsa;
fi = rhof*fi(-1)+epsf;
z = rhoz*z(-1)+epsz;
q = rhoq*q(-1)+epsq;
theta = rhot*theta(-1)+epst;

dce = ce-ce(-1);
dch = ch-ch(-1);
% 观测方程
dc = c-c(-1);
dinv = inv-inv(-1);
dqL = qL-qL(-1);
dq = q-q(-1);
db = b-b(-1);
 
end;

steady;
check;

shocks;
var epsa;
stderr .01;
var epsf;
stderr .01;
var epsz;
stderr .01;
var epsq;
stderr .01;
var epst;
stderr .01;
end;

estimated_params;

stderr epsa,   3.3767,inv_gamma_pdf, 0.5,1;
stderr epsf,   6.8809,inv_gamma_pdf, 0.5,1;
stderr epsz,   0.9462,inv_gamma_pdf, 0.5,1;
stderr epsq,   0.9508,inv_gamma_pdf, 0.5,1;
stderr epst,   0.8211,inv_gamma_pdf, 0.5,1;

gamma_h,     beta_pdf,    0.50, 0.15;
gamma_e,     beta_pdf,    0.50, 0.15;        
v,           gamma_pdf,   4.00, 1.50;
phi,         beta_pdf,    0.20, 0.15;
fi_bar,      gamma_pdf,   0.90, 0.50;
phik,        gamma_pdf,   0.20, 0.10;

rhoa,        0.8263,beta_pdf,   0.60, 0.20;
rhof,        0.9968,beta_pdf,   0.60, 0.20;
rhoz,        0.7139,beta_pdf,   0.60, 0.20;
rhoq,        0.9405,beta_pdf,   0.60, 0.20;
rhot,        0.9447,beta_pdf,   0.60, 0.20;

end;

varobs
%dq 
dc dinv dqL db; 

estimation(
nodisplay,
nograph,
datafile=china_data,
%conf_sig = 0.68,
mh_replic=0,
mode_compute=5)
%dq
dc dinv dqL db;


stoch_simul(irf=20,order=2,conditional_variance_decomposition=[5 10 20]) dc dce dch dqL dinv db;

%shock_decomposition dqL;