% Used in the forecast comparison with (Bayesian)SVAR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If these codes are useful, cite: 
% Ginters Buss (2015) "Search-and-matching frictions and labor market dynamics in Latvia", Working Papers 2015/04, Latvijas Banka.
% 
% 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 %lepsilonU 
lmupsiU %lpitargetU 
%lUpsilonU lzetacU 
lRstarU %phitildeU ltauyU lzetahU 
lystarU lpistarU %lgU 
lmuzU 

data_ydiffU 
data_picU 
data_RU 

;


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_ydiffU_eps 
data_picU_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.995;  %  %nominal i-rate ~ 4.4% 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=.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/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_ydiffU-data_ydiffU_eps+100*log(muzplus)=100*(lmuzplusU+(lystarU - lystarU(-1))); % foreign output qoq growth %
data_picU-data_picU_eps=400*log(exp(lpistarU)); % foreign inflation qoq annualized %
data_RU=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_ydiffU=0;
data_picU=400*log(pistar);
data_RU=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, .9,        inv_gamma_pdf,1,inf;  
 %%% stderr mupsi_eps,         inv_gamma_pdf,0.15,inf; 
stderr ystar_eps, 1.5,      inv_gamma_pdf,1,inf;  
stderr pistar_eps, 1.09,     inv_gamma_pdf,1,inf;  % rescaled, dividing by 1000
stderr Rstar_eps,  .438,     inv_gamma_pdf,.5,inf; 


% *** Exog. AR1 processes ***
rhomuz, .77,beta_pdf,0.7, 0.1; 
% *** rhomupsi,    beta_pdf,0.5, 0.065; ***

% *** Foreign VAR ***

% Std setup:
a11,.94, beta_pdf,0.9,0.05; % std .22 = var .05  = alpha,beta=2 in beta_pdf, see wiki
a22,.45, beta_pdf,0.5,0.05; 
a33,.88, beta_pdf,0.9,0.05; 

a12,.02, normal_pdf,-.05,0.05; 
a13,-.12, normal_pdf,-.2,0.1;
a21,.17, normal_pdf,.15,0.1;
a23,.006, normal_pdf,-.05,0.05;
a24,.09, normal_pdf,.1,0.1;
a31,.1, normal_pdf,.05,0.05;
a32,-.09,normal_pdf,-.1,0.1; 
a34,-.06, normal_pdf,.05,0.05;

c21,.083, normal_pdf,0.08,0.05;
c31,.066, normal_pdf,0.07,0.05;
c32,.018, normal_pdf,0.02,0.02;
c24,.039, normal_pdf,0.04,0.04;
c34,.01, normal_pdf,0.01,0.01;
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_ydiffU data_picU 
data_RU  %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=0, 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 %lepsilonU 
lmupsiU %lpitargetU 
%lUpsilonU lzetacU 
lRstarU %phitildeU ltauyU lzetahU 
lystarU lpistarU %lgU 
lmuzU 

data_ydiffU 
data_picU 
data_RU  

;

% 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_ydiffU 
data_picU 
data_RU  

; 

% 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('foreign_svar36_results.mat')
whichmoment='Mean'; % 'Median' or 'Mean'

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

myvarlist = {'Rstar','pistar','ystardiff'};

Rstar_ss=oo_.steady_state(9);
pistar_ss=oo_.steady_state(8);

mylevel = [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

disp('posterior mean of params:')
oo_.posterior_mean.parameters
disp('posterior std of params:')
oo_.posterior_std.parameters

disp('posterior mean of shock std:')
oo_.posterior_mean.shocks_std
disp('posterior std of shock std:')
oo_.posterior_std.shocks_std

% conditional variance decomposition
oo_.conditional_variance_decomposition

%impulse responses for variables:

*/

 