%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Computes foreign SVAR separately from the domestic part
% that is, explicitly assumes the domestic part cannot influence
% the rest of the world.
%
% If these codes are useful, cite: 
% Ginters Buss, 2016. "Financial frictions in Latvia," Empirical Economics, Springer, vol. 51(2), pages 547-575, September.
% 
% Copyright of this version: Latvijas Banka
% These codes are distributed by Latvijas Banka for academic purposes. 
% If in doubt, contact Latvijas Banka.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 

% global structure for parameters and steady states 
global M_;
 
var            

lmuzplusU 
lmupsiU 
 
lRstarU 
lystarU lpistarU 
lmuzU 

data_ystardiffU 
data_pistarU 
data_RstarU 

;


varexo 
mupsi_eps muz_eps 
Rstar_eps 
ystar_eps pistar_eps 
data_ystardiffU_eps 
data_pistarU_eps 

;


parameters  
alphaa 
betaa  
muzplus mupsi 
rhomupsi 

muz
ystar
pistar
Rstar

%further parameters
rhomuz 
sig_muz sig_mupsi 

a11 a12 a13 a21 a22 a23 a24 a31 a32 a33 a34
sig_ystar c21 sig_pistar c24 c31 c32 sig_Rstar c34

muzPercent 

;


% *** Exog. AR1 processes ***	
rhomuz      =	0.5	;
rhomupsi	=	0.5	;

   
% *** Foreign VAR ***	
a11	=	0.5	;
a22	=	0.0	;
a33	=	0.5	;

a12	=	0	;
a13	=	0	;
a21	=	0	;
a23	=	0	;
a24	=	0	;
a31	=	0	;
a32	=	0	;
a34	=	0	;

c21	=	0	;
c31	=	0	;
c32	=	0	;
c24	=	0	;
c34	=	0	;

% ***** not estimated, therefore calib here: ******
%%%%%

% *** calib ***
alphaa=.4; 

betaa=0.995; % GB2013: if real annual interest rate is 3%, then set betaa = 1.03^(-1/4)
              % GB: take EA rate 1.02^(-1/4)
              

muzplus=1.005;      %set to match growth rate of EA GDP as given in datafile

muzPercent=100*(muzplus-1);    % if  = 100*(muzplus-1), then mupsi=1 %
                               % if < 100*(muzplus-1), mupsi>1 
                               % LV: no trend observed 
muz=1+muzPercent/100;
mupsi=(muzplus/muz)^((1-alphaa)/alphaa);

%Steady state inflation rates to be set

pistar=1.005;
Rstar=(muzplus/betaa)*pistar;
ystar=.1531;

% shocks
sig_muz     =	0.15;
sig_mupsi	=	0;

sig_ystar	=	0.5	;
sig_pistar	=	0.5	;
sig_Rstar	=	0.5	;   % different scaling


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%**** begin of model *****%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model;

%105 Composite technology growth
 exp(lmuzplusU)=(exp(lmupsiU))^(alphaa/(1-alphaa))*exp(lmuzU);

%Foreign VAR
lystarU-log(ystar)=a11*(lystarU(-1)-log(ystar))+a12*(exp(lpistarU(-1))-pistar)+a13*(exp(lRstarU(-1))-Rstar)+ystar_eps/100;

exp(lpistarU)-pistar=a21*(lystarU(-1)-log(ystar))+a22*(exp(lpistarU(-1))-pistar)+a23*(exp(lRstarU(-1))-Rstar)
                     +a24*(lmuzU(-1)-log(muz))+a24*alphaa/(1-alphaa)*(lmupsiU(-1)-log(mupsi))
                     +c21*ystar_eps/100+pistar_eps/1000+c24*muz_eps/100+c24*alphaa/(1-alphaa)*mupsi_eps/100;

exp(lRstarU)-Rstar=a31*(lystarU(-1)-log(ystar))+a32*(exp(lpistarU(-1))-pistar)+a33*(exp(lRstarU(-1))-Rstar)
                     +a34*(lmuzU(-1)-log(muz))+a34*alphaa/(1-alphaa)*(lmupsiU(-1)-log(mupsi))
                     +c31*ystar_eps/100+c32*pistar_eps/1000+Rstar_eps/100+c34*muz_eps/100+c34*alphaa/(1-alphaa)*mupsi_eps/100;

%Unit root technology shocks
lmupsiU=(1-rhomupsi)*log(mupsi)+rhomupsi*lmupsiU(-1)+mupsi_eps/100;
lmuzU=(1-rhomuz)*log(muz)+rhomuz*lmuzU(-1)+muz_eps/100;


%Measurement equations 128-142
%linking data and model variables
% all vars are in per capita terms

data_ystardiffU-data_ystardiffU_eps+100*log(muzplus)=100*(lmuzplusU+(lystarU - lystarU(-1))); % foreign output qoq growth %
data_pistarU-data_pistarU_eps=400*log(exp(lpistarU)); % foreign inflation qoq annualized %
data_RstarU=400*(exp(lRstarU)-1); % foreign nominal interest rate annualized %

end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%**** end of model *****%%%%%%%%%%%%%%%%%%%%%%/
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


initval;

lmuzplusU=log(muzplus); 
lmupsiU =log(mupsi);  
   
lRstarU =log(Rstar);  
 
 lystarU =log(ystar);   
lpistarU =log(pistar);  

lmuzU=log(muz);

data_ystardiffU=0;
data_pistarU=400*log(pistar);
data_RstarU=400*log(Rstar);

end;

resid(1);

steady;

check; 

shocks;
var muz_eps=(sig_muz)^2;
var mupsi_eps = 0^2;

var ystar_eps=(sig_ystar)^2;
var pistar_eps=(sig_pistar)^2;
var Rstar_eps=(sig_Rstar)^2;

end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
% *** estim params priors begin ****%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
estimated_params; 

%Initial value,prior shape, mean, std
%shocks
% *** all shock priors re-scaled by some factor; see the above equations
 stderr muz_eps, .316,        inv_gamma_pdf,0.25,inf;  
stderr ystar_eps, .323,      inv_gamma_pdf,0.5,inf;  
stderr pistar_eps, .603,     inv_gamma_pdf,0.5,inf;  % rescaled, dividing by 1000
stderr Rstar_eps,  .061,     inv_gamma_pdf,.075,inf; 


% *** Exog. AR1 processes ***
rhomuz, .591,beta_pdf,0.5, 0.075; 

% *** Foreign VAR ***

% Std setup:
a11,.934, beta_pdf,0.9,0.05; 
a22,.534, beta_pdf,0.5,0.05; 
a33,.964, beta_pdf,0.9,0.05; 

a12,-.18, normal_pdf,-.1,0.1; 
a13,-.059, normal_pdf,-.1,0.1;
a21,.145, normal_pdf,.1,0.1;
a23,-.079, normal_pdf,-.1,0.1;
a24,.09, normal_pdf,.05,0.1;
a31,.0655, normal_pdf,.05,0.1;
a32,-.089,normal_pdf,-.1,0.1; 
a34,.114, normal_pdf,.1,0.1;

c21,.03, normal_pdf,0.05,0.05;
c31,.155, normal_pdf,0.1,0.05;
c32,.377, normal_pdf,0.4,0.05;
c24,.065, normal_pdf,0.05,0.05;
c34,.052, normal_pdf,0.05,0.05;
end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%*** end of estim params priors ****%%%%%%%%%%%%%/
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%%%%%%%%%/
%%%%%%Observed variables%%%%%%%%%/
%%%%%%%%%%%%%%%%%%%%%%%%/
varobs  
data_ystardiffU data_pistarU 
data_RstarU  

 ;



% to avoid java heap space errors 
% and if increasing java memory in Matlab preferences
% does not work:
options_.console_mode=1; %(default: 0)

estimation(order=1,datafile=myOriginalData, mh_replic=1000000, mh_nblocks=1, mh_drop=.9,  mh_jscale=0.2, 
 mode_compute=6,  nodisplay, 
 smoother,
 nodiagnostic,
  mode_check,
 %diffuse_filter, 
 %filter_step_ahead = [1 4 8 12], forecast=20,
 filtered_vars, 
% cova_compute=0, 
% do this only for the final model b/c takes a lot of time
% moments_varendo,
% conditional_variance_decomposition=[1 4 8 20],
% 
graph_format=(eps,pdf))

lmuzplusU  
lmupsiU 

lRstarU 
lystarU lpistarU 
lmuzU 

data_ystardiffU 
data_pistarU 
data_RstarU  

;

stoch_simul(order=1,irf=20, nodisplay, 
graph_format = (pdf))  
lmuzplusU 
lmupsiU 
lRstarU 
lystarU lpistarU 
lmuzU 
data_ystardiffU 
data_pistarU 
data_RstarU  
; 



 