%%%%% TRIAL JP %%%%%

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

var c y i q s pi piH piF psiF a mc phi pi_star y_star i_star eps_a eps_g eps_cp delta_e y_obs i_obs pi_obs s_obs q_obs y_star_obs i_star_obs pi_star_obs r delta_s; 
varexo gshock mshock phishock pi_starshock y_starshock i_starshock ashock cpshock sshock;

parameters h sigma alpha eta delta thetaH beta thetaF chi rhoi psi_pi psi_y psi_dy psi_de varphi rhog rhoa rhocp rhophi rhopi_star rhoy_star rhoi_star;


%----------------------------------------------------------------
% 2. Calibration (of parameters)
%----------------------------------------------------------------

h = 0.5; 
sigma = 1.2; 
alpha = 0.185; 
eta = 1.5; 
delta = 0.5; 
thetaH = 0.5;
beta = 0.99;
thetaF = 0.5;
chi = 0.01;
rhoi = 0.5;
psi_pi = 2.5;     
psi_y = 0.25;
psi_dy = 0.25;
psi_de = 0.25;
varphi = 1.5;   
rhog = 0.8;
rhoa = 0.8;
rhocp = 0.5;
rhophi = 0.8;
rhopi_star = 0.8;
rhoy_star = 0.8;  
rhoi_star = 0.8; 

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

model(linear);

c-h*c(-1) = (c(+1)-h*c) - sigma^(-1)*(1-h)*(i-pi(+1)) + sigma^(-1)*(1-h)*(eps_g - eps_g(+1));            % eq (14)
(1-alpha)*c = y - alpha*eta*(2-alpha)*s - alpha*eta*psiF - alpha*y_star;                                 % eq (15)
s-s(-1) = piF - piH;                                                                                     % eq (16)
q = psiF + (1-alpha)*s;                                                                                  % eq (17)
piH-delta*piH(-1) = thetaH^(-1)*(1-thetaH)*(1-thetaH*beta)*mc + beta*(piH(+1)-delta*piH);                % eq (18)
mc = varphi*y - (1+varphi)*eps_a + alpha*s + sigma*(1-h)^(-1)*(c-h*c(-1));                               % equation for marginal cost
piF-delta*piF(-1) = thetaF^(-1)*(1-thetaF)*(1-thetaF*beta)*psiF + beta*(piF(+1)-delta*piF) + cpshock;    % eq (19)
pi = piH +alpha*(s-s(-1));                                                                               % eq (20)
(i-pi(+1)) - (i_star-pi_star(+1)) = (q(+1)-q) - chi*a - phi;                                             % eq (21)
c+a = beta^(-1)*a(-1) - alpha*(s+psiF) + y;                                                              % eq (22)
i = rhoi*i(-1) + psi_pi*pi + psi_y*y + psi_dy*(y-y(-1))+ +psi_de*delta_e + mshock;                       % eq (23)

delta_e + pi_star - piF =(psiF - psiF(-1));                                                              % equation for including delta_e

eps_g = rhog*eps_g(-1) + gshock;                                                                         % equation specifying gshock
eps_a = rhoa*eps_a(-1) + ashock;                                                                         % equation specifying eps_a
eps_cp = rhocp*eps_cp(-1) + cpshock;                                                                     % equation specifying eps_cp
phi = rhophi*phi(-1) + phishock;                                                                         % equation specifying phi


pi_star = rhopi_star*pi_star(-1)+ pi_starshock;                                                          % AR(1) process for pi_star
y_star = rhoy_star*y_star(-1) + y_starshock;                                                             % AR(1) process for y_star
i_star = rhoi_star*i_star(-1) + i_starshock;                                                             % AR(1) process for i_star


y_obs = y - steady_state(y);
pi_obs = pi - steady_state(pi);
i_obs = i - steady_state(i);
q_obs = q - steady_state(q);
s_obs = s - steady_state(s);
y_star_obs = y_star - steady_state(y_star);
i_star_obs = i_star -steady_state(i_star);
pi_star_obs = pi_star - steady_state(pi_star);                                                                                                         

r = i - pi;
delta_s = s-s(-1);

end;


%----------------------------------------------------------------
% 4. Computation (Steady state and the 1st-order approx. solution
%----------------------------------------------------------------

% Exact steady state
initval; 
    c = 0; 
    i = 0; 
    pi = 0; 
    y = 0; 
    s = 0; 
    psiF = 0; 
    piH = 0; 
    q = 0; 
    a = 0; 
    mc = 0; 
    eps_g = 0; 
    phi = 0; 
    pi_star = 0; 
    y_star = 0; 
    i_star = 0;
    eps_a = 0;
    delta_e = 0;
    y_obs = 0;
    i_obs = 0;
    pi_obs = 0;
    s_obs = 0;
    q_obs= 0;
    y_star_obs = 0;
    i_star_obs = 0;
    pi_star_obs = 0;
end;


steady;

check;

shocks;
var gshock; stderr 1;
var ashock; stderr 1;
var phishock; stderr 1;
var mshock; stderr 1;
var cpshock; stderr 1;
var sshock; stderr 1;
var pi_starshock; stderr 1;
var y_starshock; stderr 1;
var i_starshock; stderr 1;
end;

varobs y_obs pi_obs s_obs q_obs i_obs y_star_obs pi_star_obs i_star_obs;

estimated_params;
h, beta_pdf, 0.5, 0.25;
sigma, gamma_pdf, 1.2, 0.4;
eta, gamma_pdf, 1.50, 0.75;
delta, beta_pdf, 0.5, 0.25;
thetaH, beta_pdf, 0.5, 0.1;
varphi,gamma_pdf, 1.5, 0.75;
thetaF,  beta_pdf, 0.5, 0.1;
rhoi,  beta_pdf, 0.5, 0.25;
psi_pi,  gamma_pdf, 1.50, 0.30;     
psi_y, gamma_pdf, 0.25, 0.13;     
psi_dy, gamma_pdf, 0.25, 0.13;
psi_de, gamma_pdf, 0.25, 0.13;
end;


%BAYNESIAN ESTIMATION

%Step 1
//estimation(datafile=JPdata2,mode_compute=4,mh_replic=0);

%Step 2
//estimation(datafile=JPdata2,mode_compute=0,mode_file=TRIAL_JP4_MATT_mode,mh_replic=5000,mh_nblocks=2,mh_drop=0.4,mh_jscale=0.7);

%Step 3
estimation(datafile=JPdata2,mode_compute=6,mh_replic=10000,mh_nblocks=2,mh_drop=0.1,conf_sig=0.95,mh_conf_sig=0.95,moments_varendo,bayesian_irf,irf=6);

% SOLVE THE MODEL - USING CALIBRATIONS
//stoch_simul(irf=6, order=1);