%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file implements the model of ’Introducing Financial    %
% Frictions and Unemployment into a Small Open Economy Model’ %
% by Lawrence J. Christiano, Mathias Trabandt and Karl       %
% Walentin, Journal of Economic Dynamics and Control         %
% Implemented for Dynare 4.1.0.                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% gb2013 SVAR only
% don't use alpha beta gamma psi

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

lmuzplusU %lepsilonU 
lmupsiU %lpitargetU 
%lUpsilonU lzetacU 
lRstarU %phitildeU ltauyU lzetahU 
lystarU lpistarU %lgU 
lmuzU 

data_ystardiffU 
data_pistarU 
data_RstarU 

;


varexo 
mupsi_eps muz_eps %epsilon_eps pitarget_eps
%Upsilon_eps zetac_eps 
Rstar_eps %phitilde_eps tauy_eps 
%zetah_eps epsR_eps 
ystar_eps pistar_eps %g_eps 
%taud_eps taux_eps taumc_eps taumi_eps taumx_eps
data_ystardiffU_eps 
data_pistarU_eps 

;


parameters  
alphaa %kappad 
betaa %rhopi xid xix kappax ximc etaf
%ximx kappamc kappamx kappami ximi omegac etac
%omegai etai Spp 
muzplus mupsi %delta b sigmab sigmaa
%phitildea a phitildes rhoR rpi ry rq rdeltapi rdeltay etax omegax phi 
rhomupsi %epsilon 

muz %ybar
%u kappaw AL sigmaL Rk varphi 
ystar
%mcmc px pimc pmc pmi kbar i c pc pkprime pinvest pmx pii 
pistar
Rstar

%further parameters
rhomuz  %rhoepsilon rhoUpsilon rhozetac  
%rhozetah rhog rhotauy rhophitilde phitilde 

sig_muz sig_mupsi %sig_epsilon sig_upsilon sig_zetac sig_zetah sig_phitilde

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

%SppScaled 
muzPercent %mupsiPercent

;


% *** 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; %0.35; %  SE:0.2; to get same prior k/y=8, as the prior value in CTW
            % doesnot find ss with 0.35
            % Aleksejs has about k/y=1.5
            % alpha=0.35 gives about the same

betaa=0.99;  % SE: 0.99943
              % 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;      %LV: 1.0137 SE: 1.005 set to match growth rate of GDP as given in datafile

muzPercent=100*(muzplus-1);    % if  = 100*(muzplus-1), then mupsi=1 %
                               % if < 100*(muzplus-1), mupsi>1 (Swedish case)
                               % 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;
% 400*log(Rstar)
ystar=4;

% 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/100+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/100+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); 
% lepsilonU =log(epsilon);  
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, .15,        inv_gamma_pdf,0.15,inf;  
 %%% stderr mupsi_eps,         inv_gamma_pdf,0.15,inf; 
stderr ystar_eps, .5,      inv_gamma_pdf,0.5,inf;  
stderr pistar_eps, .5,     inv_gamma_pdf,0.5,inf;  
stderr Rstar_eps,  .15,     inv_gamma_pdf,.15,inf; % scaled differently, by 1000


% *** Exog. AR1 processes ***
rhomuz, .5,beta_pdf,0.5, 0.075; 
% *** rhomupsi,    beta_pdf,0.5, 0.065; ***

% *** Foreign VAR ***

% Std setup:
a11,.8, normal_pdf,0.5,0.5;  
a22,0, normal_pdf,0.0,0.5; 
a33,.8, normal_pdf,0.5,0.5; 

a12,0, normal_pdf,0.0,0.5; 
a13,0, normal_pdf,0.0,0.5;
a21,0, normal_pdf,0.0,0.5;
a23,0, normal_pdf,0.0,0.5;
a24,0, normal_pdf,0.0,0.5;
a31,0, normal_pdf,0.0,0.5;
a32,0,normal_pdf,0.0,0.5; 
a34,0, normal_pdf,0.0,0.5;

c21,0, normal_pdf,0.0,0.5;
c31,0, normal_pdf,0.0,0.5;
c32,0, normal_pdf,0.0,0.5;
c24,0, normal_pdf,0.0,0.5;
c34,0, normal_pdf,0.0,0.5;
end;


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

% estimated_params_bounds;
% a33, -.999, .999; 
% end;



%%%%%%%%%%%%%%%%%%%%%%%%/
%%%%%%Observed variables%%%%%%%%%/
%%%%%%%%%%%%%%%%%%%%%%%%/
varobs  %data_RU data_wdiffU data_cdiffU    
%data_ydiffU  
data_ystardiffU data_pistarU 
data_RstarU  %data_HdiffU %data_HhatU %replaced by gb2013

 ;



% 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=600000, mh_nblocks=1, mh_drop=.9,  mh_jscale=0.2, 
 mode_compute=6,  nodisplay, 
 %nodiagnostic,
 %smoother,
  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 %lepsilonU 
lmupsiU %lpitargetU 
%lUpsilonU lzetacU 
lRstarU %phitildeU ltauyU lzetahU 
lystarU lpistarU %lgU 
lmuzU 

data_ystardiffU 
data_pistarU 
data_RstarU  

;

% due to time constraints, do this only for the final calibration
% ,moments_varendo, conditional_variance_decomposition=[1 4 8]
% for recursive forecasts use `nobs=[]' integer but do this maybe with calibrated vars

% stoch_simul(order=1,irf=20, nodisplay, 
% %order=2, pruning,periods=100,simul_replic=5,
% graph_format = (pdf)) % eps, 
% lmuzplusU %lepsilonU 
% lmupsiU %lpitargetU 
% %lUpsilonU lzetacU 
% lRstarU %phitildeU ltauyU lzetahU 
% lystarU lpistarU %lgU 
% lmuzU 
% 
% data_ystardiffU 
% data_pistarU 
% data_RstarU  
% 
% ; 

% shock_decomposition(parameter_set=posterior_mean)
% data_ydiffU 
% data_picU
% %data_spreaddiffU 
% %data_unempdiffU
% ;

% do this only for the final model b/c time-consuming

%identification(graph_format = (eps, pdf));

%dynare_sensitivity(identification=1,morris=1,
%graph_format = (eps, pdf), conf_sig=0.90,
%datafile=ctwBaselineData);


% dynare_sensitivity(identification=1, morris=2); % Iskrev

%model_diagnostics(M_,options_,oo_) %run it when some error for some further info


/*
load('baseline200_results.mat')
whichmoment='Mean'; % 'Median' or 'Mean'

smoothedvars= eval(['oo_.SmoothedVariables.',whichmoment]);
updatedvars = eval(['oo_.UpdatedVariables.',whichmoment]);
filteredvars = eval(['oo_.FilteredVariables.',whichmoment]);

myvarlist = {'ydiff','pid','R','wdiff','cdiff','idiff','qdiff','Hdiff',...
'xdiff','impdiff','pic','pii','gdiff','Rstar','pistar','ystardiff'};

R_ss=oo_.steady_state(90);
Rstar_ss=oo_.steady_state(103);
pid_ss=oo_.steady_state(89); 
pic_ss=oo_.steady_state(99);
pii_ss=oo_.steady_state(100);
pistar_ss=oo_.steady_state(102);

mylevel = [0,pid_ss,R_ss,0,0,0,0,0,0,0,pic_ss,pii_ss,0,Rstar_ss,pistar_ss,0];

for myj = 1:length(myvarlist)
 smoothedvar = eval(['smoothedvars.data_',myvarlist{myj},'U']);
 updatedvar = eval(['updatedvars.data_',myvarlist{myj},'U']);
 filteredvar = eval(['filteredvars.data_',myvarlist{myj},'U']);
 
    figure(myj)
    plot(smoothedvar,'blue','linewidth',2)
    hold on
    plot(filteredvar +mylevel(myj),'--red', 'linewidth',2)
    legend([myvarlist{myj},' original'], [myvarlist{myj},' filtered'])  
    hold off
end


findWhatParams={'delta','varphi','AL'};
myParamVals=zeros(length(findWhatParams),1);
myparamcell=cellstr(M_.param_names);
for myI=1:length(findWhatParams)
    myIndex=find(strcmp(myparamcell,findWhatParams(myI)));
    myParamVals(myI)=M_.params(myIndex);
end
disp('matched from moments:');
myParamVals

*/