% Open-economy New Keynesian model by Alpanda et al (2010) with Rotemberg (1982) quaratic costs of price adjustments
%
% Only 10a definition is proper
%q,e-real (nominal) exchange rates refers to inverse definition of exchange rate (the higher is the exchange rate, the lower is the value of domestic currency)
pkg load io;
pkg load struct;
pkg load optim;


%----------------------------------------------------------------
% 0. Housekeeping (close all graphic windows)
%----------------------------------------------------------------
options_.qz_criterium=1;%Tolerance for unit eigenvalue 
close all;

%----------------------------------------------------------------
% 1. Defining variables
%----------------------------------------------------------------
var c,i,y,r,ne,s_D,cpi,c_h,m,s_MI,pi,pi_f,mc,s_FCP,q,d,pi_for,y_for,i_for,Psi_F,s,s_I,s_RP,a,s_Z,y_obs,pi_obs,cpi_obs,ir_obs,pi_f_obs,de_obs,y_for_obs,pi_for_obs,ir_for_obs;
varobs y_obs,pi_obs,cpi_obs,ir_obs,pi_f_obs,de_obs,y_for_obs,pi_for_obs,ir_for_obs;
varexo e_d,e_I,e_z,e_MI,e_FCP,e_rp,e_pi_for,e_i_for,e_y_for;
%s_D-preference shock (domestic)
%cpi- (domestic) CPI inflation
%c_h-exports
%m-imports
%s_MI- home produced goods cost-push shock (mark-up shock)
%pi- domestic goods inflation (GDP deflator)
%pi_f-foreign goods inflation
%mc-marginal cost of domestic firms
%s_FCP-imported goods cost-push shock
%q-Real exchange rate
%d-nominal depreciation rate
%pi_for-foreign inflation (inflation abroad, GDP deflator)
%y_for-foreign output
%i_for-foreign interest rate
%Psi_F-deviation from the law-of-one-price
%s-terms of trade
%s_I-monetary policy shock
%s_RP-risk premium shock
%a-foreign bonds held by domestic households (ratio to trend GDP)
%s_Z-productivity shock
%-----------------------------------------
% 2. Calibration
%----------------------------------------------------------------
parameters beta, zeta, sigma, gamma,ni, alpha, chi, theta, delta, phi, phi_f, kappa, kappa_f, alpha_pi, phi_i, alpha_y, alpha_d, rho_z, rho_d, rho_rp, rho_i_for_1, rho_i_for_2, rho_y_for_1, rho_y_for_2, rho_pi_for_1, rho_pi_for_2;
beta    	= 0.994535;				% discount factor calibrated to match realized real interest rate
%zeta=0.831; %habit persistence parameter in consumption
%sigma=0.78;% Inverse of intertemporal elasticity of substitution
%gamma=1.478; % Inverse of labour supply elasticity
%ni=0.572; %Elasticity of subtitution between home-produced and imported goods
alpha=0.4216; % Share of imported goods in consumption (calibrated to match average share of imports in GDP)
chi=0.01; %Sensitivity of risk premium to foreign bond holdings (calibrated)
theta=6; % Elasticity of substitution between intermediate goods (values 6 implies 20% mark-up), calibrated
%delta=0.928; % Persistence of export demand
%phi=0.474;%Price indexation of home-produced goods
%phi_f=0.608; % Price indexation of imported goods
%kappa=96.56; % Price adjustment cost of home-produced goods
%kappa_f=112.796; % Price adjustment cost of foreign produced goods
%alpha_pi=1.422; %Taylor rule coefficient for inflation
%phi_i=0.916;%Taylor rule smoothing parameter
%alpha_y=0.294; % Taylor rule coefficient for output 
%alpha_d=0.249; %Taylor rule coefficient by depreciation (adjusted to quarterly interest rate)
%rho_z=0.823; %Productitivity shock persistence
%rho_d=0.812;%Demand (preference) shock persistence
%rho_rp=0.912;%Risk premium shock persistence
%rho_i_for_1=1.2;
rho_i_for_2=0;
%rho_y_for_1=0.983;
rho_y_for_2=0;
%rho_pi_for_1=0.332;
rho_pi_for_2=0;


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

model(linear); 
c=1/(1+zeta)*c(+1)+(zeta/(1+zeta))*c(-1)-((1-zeta)/(sigma*(1+zeta)))*(i-cpi(+1))+s_D;% New Keynesian IS curve (1)
y=alpha*s+c+alpha*(c_h-m); % National accounting identity (2)
c_h=delta*c_h(-1)+(1-delta)*(ni*q+y_for); %Exports (3)
m=c+(1-ni*(1-alpha))*s;%Imports (4)
pi=(beta/(1+beta*phi))*pi(+1)+(phi/(1+beta*phi))*pi(-1)+((theta-1)/(kappa*(1+beta*phi)))*mc+s_MI;%Domestic-good NKPC (5)
mc=gamma*y-(1+gamma)*s_Z+alpha*s+(sigma/(1-zeta))*(c-zeta*c(-1)); % Marginal cost of domestic firms equation (6)
pi_f=(beta/(1+beta*phi_f))*pi_f(+1)+(phi_f/(1+beta*phi_f))*pi_f(-1)+((theta-1)/(kappa_f*(1+beta*phi_f)))*Psi_F+s_FCP ;%Imported good NKPC (7)
%cpi=(1-alpha)*pi(-1)+alpha*pi_f(-1); %CPI inflation-orignal version (8)
cpi=(1-alpha)*pi+alpha*pi_f;%CPI inflation-version with current domestic and foreign inflation rates (9)
%q-q(-1)=d+pi_for-pi;%Real exchange rate definition (10) wrong definition	
q=(1-alpha)*s+Psi_F; %Real exchange rate changed definition (10a)
s-s(-1)=pi_f-pi; % Terms of trade definition (11)
Psi_F-Psi_F(-1)=d+pi_for-pi_f; % Deviation from the law of one price definition
i=phi_i*i(-1)+(1-phi_i)*(alpha_pi*cpi(-1)+alpha_y*y(-1)+alpha_d*d)+s_I;%Taylor rule (backward-looking version as in original paper)
i-i_for=d(+1)+s_RP-chi*a; % Uncovered interest rate parity condition
a-(1/beta)*a(-1)=alpha*(c_h-m); %Balance of payments
%Additional definitions
r=i-cpi(+1); % Real interest rate (with respect to CPI)
ne=c_h-m; %Net exports
%Domestic shocks processes
s_Z=rho_z*s_Z(-1)+e_z;
s_D=rho_d*s_D(-1)+e_d;
s_MI=e_MI;
s_RP=rho_rp*s_RP(-1)+e_rp;
s_FCP=e_FCP;
s_I=e_I;
%Foreign variables AR(2) processes
y_for=rho_y_for_1*y_for(-1)+rho_y_for_2*y_for(-2)+e_y_for;
pi_for=rho_pi_for_1*pi_for(-1)+rho_pi_for_2*pi_for(-2)+e_pi_for;
i_for=rho_i_for_1*i_for(-1)+rho_i_for_2*i_for(-2)+e_i_for;
%Measurement equations
y_obs=y;
pi_obs=pi;
pi_f_obs=pi_f;
ir_obs=i;
de_obs=d;
y_for_obs=y_for;
ir_for_obs=i_for;
pi_for_obs=pi_for;
cpi_obs=cpi;
end;

%----------------------------------------------------------------
% 4. Computation
%----------------------------------------------------------------
%shocks;
%var e_I; stderr 0.0024;
%var e_z; stderr 0.0048;
%var e_d;stderr 0.0014;
%var e_MI; stderr 0.0084;
%var e_rp;stderr 0.0089;
%var e_FCP; stderr 0.0359;
%var e_i_for; stderr 0.0012;
%var e_y_for; stderr 0.0058;
%var e_pi_for; stderr 0.0021;
%end;
%steady;
%check;

%stoch_simul(irf=20,periods=0) y i r cpi d ne s q pi;
%model_diagnostics;
%stoch_simul(irf=1000,periods=0) y;
%----------------------------------------------------------------
% 5. Estimation
%----------------------------------------------------------------
estimated_params;
stderr e_I,inv_gamma_pdf, 0.0025, Inf;
stderr e_z, inv_gamma_pdf, 0.005, Inf;
stderr e_d,inv_gamma_pdf, 0.005, Inf;
stderr e_MI,inv_gamma_pdf, 0.005, Inf;
stderr e_rp,inv_gamma_pdf, 0.005, Inf;
stderr e_FCP,inv_gamma_pdf, 0.005, Inf;
stderr e_i_for,inv_gamma_pdf, 0.005, Inf;
stderr e_y_for,inv_gamma_pdf, 0.005, Inf;
stderr e_pi_for,inv_gamma_pdf, 0.005, Inf;
zeta, beta_pdf, 0.5,0.25;
sigma, gamma_pdf, 1.2,0.4;
ni, gamma_pdf, 1.5, 0.75;
gamma, gamma_pdf, 1.5, 0.75;
phi, beta_pdf, 0.7, 0.2;
phi_f, beta_pdf, 0.7, 0.2;
kappa, gamma_pdf, 100, 20;
kappa_f, gamma_pdf, 100, 20;
phi_i, beta_pdf, 0.5, 0.25;
alpha_pi, gamma_pdf, 1.5, 0.5;
alpha_y, gamma_pdf, 0.25, 0.125;
alpha_d, gamma_pdf, 0.25, 0.125;
delta, beta_pdf, 0.8, 0.1;
rho_z, beta_pdf, 0.8, 0.1;
rho_d, beta_pdf, 0.8, 0.1;
rho_rp, beta_pdf, 0.8, 0.1;
rho_i_for_1, beta_pdf, 0.6, 0.3, -1,1;
rho_y_for_1, beta_pdf, 0.6, 0.3, -1,1;
rho_pi_for_1, beta_pdf, 0.8, 0.1, -1,1;
end;

estimation(datafile='estimation_data.mat',mode_compute=3,mh_jscale=0.68,prefilter=0,bayesian_irf,irf=20,mh_replic=0,mh_drop=0.65,moments_varendo,conditional_variance_decomposition=[4,8,12,24]) y;
shock_groups(name=group1);
foreign_shocks=e_i_for, e_pi_for, e_y_for;
end;
shock_decomposition(use_shock_groups=group1) y;