%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The model code is based on Lawrence J. Christiano, Mathias Trabandt and Karl Walentin, 2011. Introducing financial frictions and unemployment into a small open economy model, Journal of Economic Dynamics and Control, Elsevier, vol. 35(12), pages 1999-2041.
% Adapted by Ginters Buss for Latvia as a small open economy in a monetary
% union.
% 
% 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This version implements the financial frictions version of
% the model, i.e. without employment frictions.
% gb2013
% don't use alpha beta gamma psi

% global structure for parameters and steady states 
global M_;
 
var            
%endogenous variables
lmcU lrkbarU lwbarU lRfU lpxU lpixU lkU lRnustarU 
lHU lpiU lmcxU lpimcU lmcmcU lpimxU lmcmxU lpimiU 
lmcmiU lpicU lpmcU lpiiU lpmiU lkbarU liU lpsizplusU 
lcU lpkprimeU luU lpinvestU %lsU 
aU lqU lxU 
lpcU  lRU lyU  lRxU lpmxU lRkU 
%dummy variables for investment cost function
StildeU SprimetildeU

lpitildedpiU lKdU lFdU lphaloU 
lpitildexpiU lKxU lFxU lphaloxU
lpitildemcpiU lKmcU lFmcU lphalomcU lcmU
lpitildemipiU lKmiU lFmiU lphalomiU limU
lpitildemxpiU lKmxU lFmxU lphalomxU lxmU

lpitildewpiU lKwU lFwU lwhaloU lsmallhU
%lwhalohaloU  

aofuU lpiwU

%variables for the unemployment part
%lintensityU lJzplusU lwU   
%lUzplusU lVxzplusU lvU  
%lG1U lG2U lG3U
%lOm1U lOm2U lOm3U lOm4U
%lvtilde0U lvtilde1U lvtilde2U lvtilde3U
%lVzplus0U lVzplus1U lVzplus2U lVzplus3U 
%ll0U ll1U ll2U ll3U lLU lmU lBigQU lpitildewU
%lchi0U lchi1U lchi2U lchi3U lfU lVwU JwU lUnemprateU

%variables for the fin fric part
lomegabarU lnU lspreadU lbankruptcyrateU

%law of motions for exogenous processes 
lmuzplusU lepsilonU lmupsiU lpitargetU 
lUpsilonU lzetacU lRstarU phitildeU ltauyU lzetahU 
lystarU lpistarU lgU lmuzU 
ltaudU ltauxU ltaumcU ltaumiU ltaumxU

%unemp model exog processes
%lsigmamU letaU

%fin fric model exog processes
lsigmaU lgammaU

%data
data_pidU data_RU data_wdiffU data_cdiffU data_idiffU data_qdiffU data_HdiffU 
data_ydiffU data_xdiffU data_impdiffU data_picU data_piiU data_ystardiffU 
data_pistarU data_RstarU data_gdiffU data_spreadU data_ndiffU

lylessdU

nxdivGDPU
lSlevU
data_HhatU
data_spreaddiffU
GDP12q_diffU
GDP4q_diffU
pic4qU

%data_sdevU
lbU
etaaU
;

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

%finfric model shocks
sigma_eps gamma_eps

%unemp model shocks
%sigmam_eps eta_eps

%measurement shocks
data_pidU_eps data_wdiffU_eps data_cdiffU_eps data_idiffU_eps data_qdiffU_eps data_HdiffU_eps 
data_ydiffU_eps data_xdiffU_eps data_impdiffU_eps data_picU_eps data_piiU_eps data_ystardiffU_eps 
data_pistarU_eps data_ndiffU_eps data_unempdiffU_eps data_spreadU_eps data_gdiffU_eps
data_HhatU_eps
data_spreaddiffU_eps
%data_sdevU_eps
lsU_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 pibar Upsilon zetac tauc tauk R tauy tauw zetah
nustar nuf nux etag lambdad lambdax lambdamc lambdami lambdamx 

lambdaw 

wbar pic psizplus H q pimi mcmi mc pix mcmx mcx pimx 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_epsilon sig_upsilon sig_zetac sig_zetah sig_phitilde
sig_epsR sig_pitarget sig_g sig_tauy 

s Stilde Sprimetilde 

xiw 

rkbar Rf Rnustar x Rx
varkappad varkappax varkappamx varkappami varkappamc varkappaw 
thetaw pibreve pitildedpi Kd Fd phalo pitildexpi Kx Fx phalox 
pitildemcpi Kmc Fmc phalomc cm pitildemipi Kmi Fmi phalomi im 
pitildemxpi Kmx Fmx phalomx xm

pitildewpi Kw Fw whalo  smallh 

%whalohalo

aofu 
 
taud taux taumc taumi taumx rhotaud rhotaux rhotaumc rhotaumi rhotaumx piw taub
pitarget etaa


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

%finfric parameters
we mu Fomegabar n omegabar sigma gammaa spread
rhosigma rhogamma

%unemp parameters
%recruitshare bshare rhosigmam rhoeta BigN rho iota sigma_match 
%kappa eta bu Jw intensity Vw BigQ v m sigmam f w 
%Jzplus Vxzplus Uzplus vtilde recruiting L pitildew G1 G2 G3 chi0 chi1 chi2 chi3
%Om1 Om2 Om3 Om4 vtilde0  vtilde1 
%vtilde2 vtilde3 Vzplus0 Vzplus1 Vzplus2 Vzplus3 l0 l1 l2 l3

%recruitsharePercent 

SppScaled muzPercent mupsiPercent wesharePercent
intensity_target iy_target xy_target  nk_target

ylessd hours_target
sigmaLScaled
work_cap_para
%rs
rRstar
;


%tex names for the tex option in the estimation
M_.param_names_tex=M_.param_names;
M_.endo_names_tex=M_.endo_names;
M_.exo_names_tex=M_.exo_names;
% /* need to update (include levels) tex-names if using them.
% options_.TeX_varobs=['data_RU         '
%                      'data_wdiffU     '
%                      'data_cdiffU     '
%                      'data_ydiffU     '
%                      'data_ystardiffU '
%                      'data_pistarU    '
%                      'data_RstarU     '
%                      'data_HdiffU     '
%                      'data_impdiffU   '
%                      'data_piiU       '
%                      'data_pidU       '
%                      'data_picU       '
%                      'data_xdiffU     '
%                      'data_qdiffU     '
%                      'data_idiffU     '
%                      'data_ndiffU     ' 
%                      'data_spreaddiffU'
%                      'data_gdiffU     '];
% */ 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% **** begin of parameter setting ****%%%%%%%%%%%%%%%%%%%%/
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

lambdaw=1.5; % **** ??? motivation for choice of above value ??? 
             % the only restriction is that it can't be too low (e.g. 1.2.)
             % b/c then indeterminancy hits.
xiw=0.75;

% *****************************************************
% ********** begin pasting from CTW.mod, 
% **    CTW_110208_2010MASTER_but_large_pi-MEs2 based on full sample, -2010q3
% *****************************************************

alphaa=.4; %.35;  %to get a Pk*k/y around 8 and reasonable c/y ratio in steady state. (But, the latter fails on both counts at this param vector)
betaa=0.995; % .995 *** so that like realR = 2.25, i.e. Real R Q=1.005573  (but, 3.46-1.50=1.96 in data) 
             %(2% infl; roughly 4.25% nom.Interest)  
            % GB: if real annual interest rate is 3%, then set betaa = 1.03^(-1/4)
            % LV: ~13%, so betaa = 1.13^(-.25);  
            % use the EA data, so impose betaa about .99;
            % 1.0145^(-.25)=.996 

% *** residual, not set anymore: mupsiPercent=0.04; % both tech-trends calculated for the sample 1995Q1-2009q2
muzplus=1.005;    % set to match growth rate of GDP as given in datafile
muzPercent=100*(muzplus-1); % Yields mupsi=1 %  %0.39;   % only investment prices: 0.63  
                   % 0.54 if equal weight on prices and quants

%Steady state inflation rates to be set
pibreve=1.005;
pitarget=1.005;
pistar=1.005;

etaa=0.0;   %Steady state net foreign asset position as a share of GDP
etag=0.202;	%LV:0.19  SE: 0.291 government expenditure share in GDP sample 1995Q1-2010q3

                
phitildea=0.01; 

lambdad     =	1.3	;
lambdamx	=	1.2	;
lambdamc	=	1.3	;
lambdami	=	1.3	;
lambdax     =	1.2 ;   %Note: ALLV lambdax=1 does not work here 
                        %due to non-linear pricing equations
 
% calibrated from IO-tables:
omegac=.45;   %0.45;    %LV(Santa): 0.4-0.5 SE: 0.25 import share in consumption goods 
omegai=0.65;    %LV(Santa): 0.5 or more; SE: 0.43 import share in invertment goods 
omegax=0.55; %.57;   %LV(Santa): 0.57; SE: .35 import share in export goods 

% taxes
tauc=0.18;   % LV: 0.18; SE: .35
tauw=0.33;  % LV: .33; SE: .35
tauy=0.30;   % LV: .3; SE: .3
tauk=0.1;   %LV: 0.1; SE: .25
taub=0.0;     % to get low bond interest rate

%working capital fractions
nuf=1; 
nustar=1; 
nux=1;

% wage indexation:
thetaw=1;   %weight on steady state z-plus growth

% mean values of exogenous shock processes
epsilon=1;
Upsilon=1;
zetac=1;
zetah=1;

%Steady state values of shocks to marginal costs
taux=1; 
taumc=1; 
taumi=1; 
taumx=1; 
taud=1;

%Setting some steady state values as in ALLV
u=1;
Stilde=0;
Sprimetilde=0;

%Persistence of shocks to marginal costs ("markup shocks")
rhotaux=0; 
rhotaumc=0; 
rhotaumi=0; 
rhotaumx=0; 
rhotaud=0; 

% Finfric params
wesharePercent  =   0.1;   % we/ybar expressed in percent
Fomegabar       =   0.02;  % Based on UC (Credit registry) data, sample avg 95-09q2 *** Using SCB data 0.0063 is sample avg 95-09q2.
                           % gb: bancruptcy rate
                           % gb: calibrate to s.t. forecasts AND 
                           % spread decomp unbiased
% /*                    
% % Labor params
% BigN=4;
% sigma_match=.5;
% L = 1-0.088;            % 0.078 sample avg 1995Q1-2010q3
%                         %LV about 0.1 or higher
% iota                =	1; 
% 
% sigma_a             =  .1564; % dummy-setting, gets real value in ss-file  %standard deviation of idiosyncratic productivity of workers in steady state
% BigFPercent         =  .25; % steady state endogenous separation rate, can't have any wage dispersion if BigF>0
%                                     
% rho=0.972; % old 0.9735     %  rho-BigF=0.9695 implies 1/f=3 quarters (if U=0.078)
% sigmam=0.57; % old 0.548   % such that quarterly Q is roughly 0.9 
%                  %(note: monthly Q roughly 0.6 (0.9=0.6+0.6*(1-0.6)+0.6*(1-0.6)^2 )in Swedish data)
% */

% shocks
sig_muz     =	0.15;
sig_mupsi	=	0;
% *** 0.15	; *** also set below;
sig_epsilon	=	0.5	;
sig_upsilon	=	0.5	;
sig_zetac	=	0.15; 
sig_phitilde=	0.15;
sig_epsR	=	0; %0.15;
sig_g       =	0.5	;    
sig_taud	=	0.15;
sig_taux	=	0.15;
sig_taumc	=	0.15;
sig_taumi	=	0.15;
sig_taumx	=	0.15;
sig_gamma	=	0.5	;
sig_sigma	=	0;
% *** 0.5	; *** also set below

sig_ystar	=	0.5	;
sig_pistar	=	0.5	;
sig_Rstar	=	.15; %gb: revised from 0.5	;   % different scaling

sig_zetah	=	0.15;

%Following shocks are deactivated
sig_eta     =	0	;
sig_sigma_a =   0   ;
sig_sigmam  =	0	;

sig_tauy	=	0	;
sig_pitarget=	0	;


%you better do not touch these:
phitilde=0;     %OBS! Keep to zero! Always! 
rhopi=0.975;    % unused 


% ******************************************************************
% Below SHOULD be estimated parameters, and these SHOULD be set to 
% their prior mean here.
% (param values only needed for the case that param is NOT estimated. 
% O/w redundant and overwritten in estimation
%  later.)
% ******************************************************************

% *** stickiness ****			
xid     =	0.75	;
xix     =	0.75	;
ximc	=	0.75	;
ximi	=	0.75	;
ximx	=	0.75	;

kappad	=	0.5	;
kappax	=	0.5	;
kappamc	=	0.5	;
kappami	=	0.5	;
kappamx	=	0.5	;

kappaw=0;   %weight on pibreve=1-kappaw

work_cap_para=0.5;

% *** Domestic prefs and tech ***			
sigmaLScaled      =	0.3; % 0.75	;
b           =	0.65;
SppScaled	=	0.5	;
sigmaa      =	0.2	;

% *** Taylor rule ***			
rhoR	=	0.8	;
rpi     =	1.7; %1; %1.7	;
rdeltapi=	0	;
rq      =	0	;
ry      =	.1; %0.1; %0.125	;
rdeltay	=	0	;
%rs      =   1.2;
rRstar  =   .6;   
% *** open econ. stuff	
etai	=	1.1; %1.5	;
etaf	=	1.5	;
etac	=	1.5	;
etax	=	1.5	;
phitildes=	0; %1.25	; 

% *** new financial frictions parameters ***	
mu              =	0.3;    % such that external finance premium roughly 1.5 percent annually (0.375 quarterly), as 1995Q1-2010q3


% *** Exog. AR1 processes ***	
rhomuz      =	0.5	;
rhomupsi	=	0.5	;
rhoepsilon	=	0.85	;
rhoUpsilon	=	0.85	;
rhozetac	=	0.85	;
rhophitilde	=	0.85	;
rhog        =	0.85	;
rhosigma	=	0.85	;
rhogamma	=	0.85	;
rhoeta      =	0.85	;
rhosigma_a  =   0.85    ;
   
% *** 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	;

% needed although unused
rhotauy     =	0.85	;   
rhozetah	=	0.85	;
rhosigmam	=	0.85	; 


%some steady state targets to be matched
hours_target=0.27;     %SE:0.25  = L*varsigma(hours per employee) AL to be set accordingly    
iy_target=0.2549;      %SE: 0.18 avg over sample;  delta to be set accordingly 
xy_target=0.4618;      %SE: 0.44 x/GDP avg 95-10q3. varphi to be set accordingly
nk_target=0.6;         %SE: 0.5   gamma to be set accordingly
                       % LV might be ~0.14

delta=0.03;    %SE: 0.015 Only initial, b/c calculated in ss-file
varphi=2;  % Only initial, b/c calculated in ss-file


load foreign_svar50_oo_.mat;
mypars=foreign_svar50_oo_.posterior_mean.parameters;
myparsmat = cell2mat(struct2cell(mypars));
myshocks=foreign_svar50_oo_.posterior_mean.shocks_std;
myshocksmat=cell2mat(struct2cell(myshocks));
% % *** Pasting param values from posterior

       rhomuz            = myparsmat(1)           ; 

         a11            = myparsmat(2)            ; 
         a22            = myparsmat(3)            ; 
         a33            = myparsmat(4)            ; 
         a12            = myparsmat(5)            ; 
         a13            = myparsmat(6)            ; 
         a21            = myparsmat(7)           ; 
         a23            = myparsmat(8)            ; 
         a24            = myparsmat(9)            ; 
         a31            = myparsmat(10)            ; 
         a32            = myparsmat(11)            ; 
         a34            = myparsmat(12)           ; 
         c21            = myparsmat(13)           ; 
         c31            = myparsmat(14)            ; 
         c32            = myparsmat(15)           ; 
         c24            = myparsmat(16)           ; 
         c34            = myparsmat(17)           ;  
 

sig_muz=myshocksmat(1); %(sig_muz)^2;
 sig_ystar=myshocksmat(2);  %(sig_ystar)^2;
 sig_pistar=myshocksmat(3) ; %(sig_pistar)^2;
 sig_Rstar=myshocksmat(4) ; % (sig_Rstar)^2;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% **** end of parameter setting **** %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/

%calculate no-indexation parameters 
varkappad=1-kappad;
varkappax=1-kappax;
varkappamx=1-kappamx;
varkappami=1-kappami;
varkappamc=1-kappamc;

varkappaw=1-kappaw;

intensity_target=hours_target;

%steady state get initval values
%finfric_steadystate;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%**** begin of model *****%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model;
%1 eq1 - mc: domestic homog. good marginal costs
exp(lmcU)=exp(ltaudU)*(1/(1-alphaa))^(1-alphaa)*
         (1/alphaa)^(alphaa)*(exp(lrkbarU))^(alphaa)*
         (exp(lwbarU)*exp(lRfU))^(1-alphaa)/exp(lepsilonU);

%2 eq2 - 2nd definition of mc
exp(lmcU)=exp(ltaudU)*exp(lmupsiU)^(alphaa)*(exp(lwbarU)*exp(lRfU))/
 (exp(lepsilonU)*(1-alphaa)*((exp(lkU)/exp(lmuzplusU)/exp(lHU)))^(alphaa));


%Non-linear pricing equations for domestic intermediate goods producer
%3 eq40a
exp(lFdU)-exp(lpsizplusU+lyU)-betaa*xid*exp(lpitildedpiU(+1)/(1-lambdad)+lFdU(+1));

%4 eq40b
exp(lKdU)-lambdad*exp(lpsizplusU+lyU+lmcU)-betaa*xid*exp(lpitildedpiU(+1)*lambdad/(1-lambdad)+lKdU(+1));

%5 eq40c
exp(lphaloU*lambdad/(1-lambdad))-(1-xid)*((1-xid*exp(lpitildedpiU/(1-lambdad)))/(1-xid))^lambdad
   -xid*(exp(lpitildedpiU+lphaloU(-1)))^(lambdad/(1-lambdad));

%6 eq40d
exp(lKdU-lFdU)=((1-xid*exp(lpitildedpiU/(1-lambdad)))/(1-xid))^(1-lambdad);

%7 eq40e
exp(lpitildedpiU)=exp(-lpiU+kappad*lpiU(-1)+(1-kappad-varkappad)*lpitargetU+varkappad*log(pibreve));


%Non-linear pricing equations for export goods producer
% 8 eq4a
exp(lFxU)-exp(lpsizplusU+lqU+lpcU+lpxU+lxU)-betaa*xix*exp(lpitildexpiU(+1)/(1-lambdax)+lFxU(+1));

%9 eq4b
exp(lKxU)-lambdax*exp(lpsizplusU+lqU+lpcU+lpxU+lxU+lmcxU)-betaa*xix*exp(lpitildexpiU(+1)*lambdax/(1-lambdax)+lKxU(+1));

%10 eq4c
exp(lphaloxU*lambdax/(1-lambdax))-(1-xix)*((1-xix*exp(lpitildexpiU/(1-lambdax)))/(1-xix))^lambdax
   -xix*(exp(lpitildexpiU+lphaloxU(-1)))^(lambdax/(1-lambdax));

%11 eq4d
exp(lKxU-lFxU)=((1-xix*exp(lpitildexpiU/(1-lambdax)))/(1-xix))^(1-lambdax);

%12 eq4e
exp(lpitildexpiU)=exp(-lpixU+kappax*lpixU(-1)+(1-kappax-varkappax)*log(pix)+varkappax*log(pibreve));


%Non-linear pricing equations for consumption importer
%13 eq5ac
exp(lFmcU)-exp(lpsizplusU+lpmcU+lcmU)-betaa*ximc*exp(lpitildemcpiU(+1)/(1-lambdamc)+lFmcU(+1));

%14 eq5bc
exp(lKmcU)-lambdamc*exp(lpsizplusU+lpmcU+lmcmcU+lcmU)-betaa*ximc*exp(lpitildemcpiU(+1)*lambdamc/(1-lambdamc)+lKmcU(+1));

%15 eq5cc
exp(lphalomcU*lambdamc/(1-lambdamc))-(1-ximc)*((1-ximc*exp(lpitildemcpiU/(1-lambdamc)))/(1-ximc))^lambdamc
   -ximc*(exp(lpitildemcpiU+lphalomcU(-1)))^(lambdamc/(1-lambdamc));

%16 eq5dc
exp(lKmcU-lFmcU)=((1-ximc*exp(lpitildemcpiU/(1-lambdamc)))/(1-ximc))^(1-lambdamc);

%17 eq5ec
exp(lpitildemcpiU)=exp(-lpimcU+kappamc*lpimcU(-1)+(1-kappamc-varkappamc)*lpitargetU+varkappamc*log(pibreve));


%Non-linear pricing equations for export inputs importer
%18 eq5ax
exp(lFmxU)-exp(lpsizplusU+lpmxU+lxmU)-betaa*ximx*exp(lpitildemxpiU(+1)/(1-lambdamx)+lFmxU(+1));

%19 eq5bx
exp(lKmxU)-lambdamx*exp(lpsizplusU+lpmxU+lmcmxU+lxmU)-betaa*ximx*exp(lpitildemxpiU(+1)*lambdamx/(1-lambdamx)+lKmxU(+1));

%20 eq5cx
exp(lphalomxU*lambdamx/(1-lambdamx))-(1-ximx)*((1-ximx*exp(lpitildemxpiU/(1-lambdamx)))/(1-ximx))^lambdamx
   -ximx*(exp(lpitildemxpiU+lphalomxU(-1)))^(lambdamx/(1-lambdamx));

%21 eq5dx
exp(lKmxU-lFmxU)=((1-ximx*exp(lpitildemxpiU/(1-lambdamx)))/(1-ximx))^(1-lambdamx);

%22 eq5ex
exp(lpitildemxpiU)=exp(-lpimxU+kappamx*lpimxU(-1)+(1-kappamx-varkappamx)*lpitargetU+varkappamx*log(pibreve));

%Non-linear pricing equations for investment importer
%23 eq5ai
exp(lFmiU)-exp(lpsizplusU+lpmiU+limU)-betaa*ximi*exp(lpitildemipiU(+1)/(1-lambdami)+lFmiU(+1));

%24 eq5bi
exp(lKmiU)-lambdami*exp(lpsizplusU+lpmiU+lmcmiU+limU)-betaa*ximi*exp(lpitildemipiU(+1)*lambdami/(1-lambdami)+lKmiU(+1));

%25 eq5ci
exp(lphalomiU*lambdami/(1-lambdami))-(1-ximi)*((1-ximi*exp(lpitildemipiU/(1-lambdami)))/(1-ximi))^lambdami
   -ximi*(exp(lpitildemipiU+lphalomiU(-1)))^(lambdami/(1-lambdami));

%26 eq5di
exp(lKmiU-lFmiU)=((1-ximi*exp(lpitildemipiU/(1-lambdami)))/(1-ximi))^(1-lambdami);

%27 eq5ei
exp(lpitildemipiU)=exp(-lpimiU+kappami*lpimiU(-1)+(1-kappami-varkappami)*lpitargetU+varkappami*log(pibreve));

%28 eq6 - Domestic consumption inflation
exp(lpicU)=exp(lpiU)*
((1-omegac+omegac*(exp(lpmcU))^(1-etac))/
(1-omegac+omegac*(exp(lpmcU(-1)))^(1-etac)))^(1/(1-etac));


%29 eq7 - Domestic investment inflation
exp(lpiiU)=exp(lpiU)/exp(lmupsiU)*
((1-omegai+omegai*(exp(lpmiU))^(1-etai))/
(1-omegai+omegai*(exp(lpmiU(-1)))^(1-etai)))^(1/(1-etai));

%30 eq8 - Law of motion for physical capital
exp(lkbarU)=(1-delta)/(exp(lmuzplusU)*exp(lmupsiU))*
    exp(lkbarU(-1))+exp(lUpsilonU)*(1-StildeU)*exp(liU);

%31 eq9 - Household consumption FOC
exp(lzetacU)/(exp(lcU)-b*exp(lcU(-1))/exp(lmuzplusU))-
betaa*b*(exp(lzetacU(+1))/
(exp(lcU(+1))*exp(lmuzplusU(+1))-b*exp(lcU)))
-exp(lpsizplusU)*exp(lpcU)*(1+tauc)=0;


%Commented out when finanical frictions part of the model
%32 eq10 - Household capital FOC
%-exp(lpsizplusU)+betaa*exp(lpsizplusU(+1))/
%(exp(lpiU(+1))*exp(lmuzplusU(+1)))*exp(lRkU(+1));


% **************************************%%%%%%%%%%%%%%%%%%%%%
% ***The financial frictions block******%%%%%%%%%%%%%%%%%%%%%
% **************************************%%%%%%%%%%%%%%%%%%%%%
%32 eq34 law of motion for net worth
-exp(lnU)+exp(lgammaU-lpiU-lmuzplusU)* % old: log((exp(lmupsiU))^(alphaa/(1-alphaa))*exp(lmuzU))*
(exp(lRkU+lpkprimeU(-1)+lkbarU(-1)) 
-exp(lRU(-1))*(exp(lpkprimeU(-1)+lkbarU(-1))-exp(lnU(-1)))
-mu*normcdf((lomegabarU +exp(lsigmaU(-1))^2/2)/exp(lsigmaU(-1)) -  exp(lsigmaU(-1)),0,1)   
*exp(lRkU+lpkprimeU(-1)+lkbarU(-1))) + we;

%33 eq36 first order condition associated with contract
(1-(normcdf((lomegabarU(+1)+exp(lsigmaU)^2/2)/exp(lsigmaU) - exp(lsigmaU),0,1)   
+exp(lomegabarU(+1))*(1-normcdf((lomegabarU(+1)+exp(lsigmaU)^2/2)/exp(lsigmaU),0,1)) ))*
( (exp(lRkU(+1)))/(exp(lRU)) )
+( (1-( normcdf((lomegabarU(+1)+exp(lsigmaU)^2/2)/exp(lsigmaU),0,1) ))
 /(1-( normcdf((lomegabarU(+1)+exp(lsigmaU)^2/2)/exp(lsigmaU),0,1) )
-mu*exp(lomegabarU(+1))*( exp(-0.5*(((lomegabarU(+1)+exp(lsigmaU)^2/2)/exp(lsigmaU))/1)^2)/(1*sqrt(2*3.141592653589793))/exp(lomegabarU(+1))/exp(lsigmaU) )) )
*(( ((exp(lRkU(+1))))/(exp(lRU)) )*( (1-mu)*normcdf((lomegabarU(+1)+exp(lsigmaU)^2/2)/exp(lsigmaU) - exp(lsigmaU),0,1)
+exp(lomegabarU(+1))*(1-normcdf((lomegabarU(+1)+exp(lsigmaU)^2/2)/exp(lsigmaU),0,1)) )-1); 

%34 eq35 zero profit condition
(exp(lpkprimeU(-1)+lkbarU(-1)+lRkU)*( (1-mu)*normcdf((lomegabarU +exp(lsigmaU(-1))^2/2)/exp(lsigmaU(-1)) -  exp(lsigmaU(-1)),0,1)
+exp(lomegabarU)*(1-normcdf((lomegabarU+exp(lsigmaU(-1))^2/2)/exp(lsigmaU(-1)),0,1)) )
 /(exp(lnU(-1)+lRU(-1))))-(exp(lpkprimeU(-1)+lkbarU(-1))/exp(lnU(-1)))+1; 
% ***********************************************%%%%%%%%%%%%%%%%%%
% *******end of financial frictions block *******%%%%%%%%%%%%%%%%%%
% ***********************************************%%%%%%%%%%%%%%%%%%

%35 eq11 - Household investment FOC
-exp(lpsizplusU)*exp(lpinvestU)+exp(lpsizplusU)*exp(lpkprimeU)*
exp(lUpsilonU)*(1-StildeU-SprimetildeU*exp(lmuzplusU)*
exp(lmupsiU)*exp(liU)/exp(liU(-1)))+betaa*exp(lpsizplusU(+1))*exp(lpkprimeU(+1))*
exp(lUpsilonU(+1))*SprimetildeU(+1)*(exp(liU(+1))/exp(liU))^2*exp(lmuzplusU(+1))*
exp(lmupsiU(+1))=0;

%36 eq12 - Household capital util. FOC
exp(lrkbarU)=exp(lpinvestU)*(sigmab*sigmaa*
   (exp(luU))+sigmab*(1-sigmaa));

%37 eq13 - Household foreign bonds FOC
%Note: net asset position variable <a> linearized and not log-lin. 
%since a=0 in steady state might be possible
-exp(lpsizplusU)+betaa*exp(lpsizplusU(+1))/(exp(lmuzplusU(+1))*exp(lpiU(+1)))
*(exp(lsU_eps(+1))*exp(lRstarU)*
exp(-phitildea*(aU-a)-phitildes*(exp(lRstarU)-exp(lRU)-(Rstar-R))+phitildeU)
-taub*(exp(lsU_eps(+1))*(exp(lRstarU))*
% *** old mistake: exp(-phitildea*(aU-a)-phitildes*(exp(lRstarU-lRU)-(Rstar-R))+phitildeU)
exp(-phitildea*(aU-a)-phitildes*(exp(lRstarU)-exp(lRU)-(Rstar-R))+phitildeU)
-exp(lpiU(+1))))=0;

%38 eq14 - Fisher equation 
-exp(lpsizplusU)+betaa*exp(lpsizplusU(+1))/exp(lmuzplusU(+1))
*(exp(lRU)-taub*(exp(lRU)-exp(lpiU(+1))))
/exp(lpiU(+1))=0;


 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% *********The labor market block**************%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%39 Jzplus - Stationary value of a job for the firm
%-exp(lJzplusU) +((exp(lwbarU)-exp(lwU)*exp(lwbarU))*exp(lintensityU)-(kappa/2)*(exp(lvtilde0U)^2))  % was wrong
%+betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU))*((exp(lwbarU(+1))-exp(lG1U(+1))*exp(lwU)*exp(lwbarU))*exp(lintensityU(+1))-(kappa/2)*(exp(lvtilde1U(+1))^2))*exp(lOm1U)
%+betaa^2*(exp(lpsizplusU(+2))/exp(lpsizplusU))*((exp(lwbarU(+2))-exp(lG2U(+2))*exp(lwU)*exp(lwbarU))*exp(lintensityU(+2))-(kappa/2)*(exp(lvtilde2U(+2))^2))*exp(lOm2U(+1))
%+betaa^3*(exp(lpsizplusU(+3))/exp(lpsizplusU))*((exp(lwbarU(+3))-exp(lG3U(+3))*exp(lwU)*exp(lwbarU))*exp(lintensityU(+3))-(kappa/2)*(exp(lvtilde3U(+3))^2))*exp(lOm3U(+2))
%+betaa^4*(exp(lpsizplusU(+4))/exp(lpsizplusU))*exp(lJzplusU(+4))*exp(lOm4U(+3));

%40-42 G updater: GjU(t)=G_t,j, dated in Dynare at time t+j, which is the date that the value of G_t,j is realized
%exp(lG1U)=exp(lpitildewU)/exp(lpiU)/(exp(lmuzplusU));
%exp(lG2U)=exp(lG1U(-1))*(exp(lpitildewU))/exp(lpiU)/(exp(lmuzplusU));
%exp(lG3U)=exp(lG2U(-1))*(exp(lpitildewU))/exp(lpiU)/(exp(lmuzplusU));

%43 Pitildew Inflation updater
%exp(lpitildewU)=(exp(lpicU(-1)))^(kappaw)*(exp(lpitargetU))^(1-kappaw);
%exp(lpitildewU)=exp(kappaw*lpicU(-1)+(1-kappaw-varkappaw)*lpitargetU+varkappaw*log(pibreve)+thetaw*log(muzplus)); 

%44 L-total employment
%exp(lLU)=exp(ll0U) + exp(ll1U) + exp(ll2U) + exp(ll3U);

%45-48 Omega updater
%exp(lOm1U)=exp(lchi0U)+rho;
%exp(lOm2U)=exp(lOm1U(-1))*(exp(lchi1U)+rho);
%exp(lOm3U)=exp(lOm2U(-1))*(exp(lchi2U)+rho);
%exp(lOm4U)=exp(lOm3U(-1))*(exp(lchi3U)+rho);

%49-52 Hiring rates
%exp(lchi0U)=exp(lBigQU)^(1-iota)*exp(lvtilde0U);
%exp(lchi1U)=exp(lBigQU)^(1-iota)*exp(lvtilde1U);
%exp(lchi2U)=exp(lBigQU)^(1-iota)*exp(lvtilde2U);
%exp(lchi3U)=exp(lBigQU)^(1-iota)*exp(lvtilde3U);

%53-56 F.O.C. of J wrt vtilde for different cohorts at a fixed date
%cohort 0
%kappa*exp(lvtilde0U)/(exp(lBigQU))^(1-iota)-betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU))*((exp(lwbarU(+1))-exp(lG1U(+1))*exp(lwU)*exp(lwbarU))*exp(lintensityU(+1))
%+(kappa/2)*(exp(lvtilde1U(+1))^2)+rho*kappa*exp(lvtilde1U(+1))/(exp(lBigQU(+1)))^(1-iota));

%cohort 1
%kappa*exp(lvtilde1U)/(exp(lBigQU))^(1-iota)-betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU))*((exp(lwbarU(+1))-exp(lG2U(+1))*exp(lwU(-1))*exp(lwbarU(-1)))*exp(lintensityU(+1))
%+(kappa/2)*(exp(lvtilde2U(+1))^2)+rho*kappa*exp(lvtilde2U(+1))/(exp(lBigQU(+1)))^(1-iota));

%cohort 2 (up to BigN-2)
%kappa*exp(lvtilde2U)/(exp(lBigQU))^(1-iota)-betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU))*((exp(lwbarU(+1))-exp(lG3U(+1))*exp(lwU(-2))*exp(lwbarU(-2)))*exp(lintensityU(+1))
%+(kappa/2)*(exp(lvtilde3U(+1))^2)+rho*kappa*exp(lvtilde3U(+1))/(exp(lBigQU(+1)))^(1-iota));

%same F.O.C. (with slight difference , for the BigN-1 cohort)
%kappa*exp(lvtilde3U)/(exp(lBigQU))^(1-iota)-betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU))*((exp(lwbarU(+1))-exp(lwU(+1))*exp(lwbarU(+1)))*exp(lintensityU(+1))
%+(kappa/2)*(exp(lvtilde0U(+1))^2)+rho*kappa*exp(lvtilde0U(+1))/(exp(lBigQU(+1)))^(1-iota));


%57-60 Vzplus for cohorts - value of a job for a worker
%exp(lVzplus0U)-exp(lwU)*exp(lwbarU)*((1-exp(ltauyU))/(1+tauw))*exp(lintensityU)+exp(lzetahU)*AL*exp(lintensityU)^(1+sigmaL)/((1+sigmaL)*exp(lpsizplusU))-betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU))*(rho*exp(lVzplus1U(+1))+(1-rho)*exp(lUzplusU(+1)));
%exp(lVzplus1U)-exp(lG1U)*exp(lwU(-1))*exp(lwbarU(-1))*((1-exp(ltauyU))/(1+tauw))*exp(lintensityU) +exp(lzetahU)*AL*exp(lintensityU)^(1+sigmaL)/((1+sigmaL)*exp(lpsizplusU))-betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU))*(rho*exp(lVzplus2U(+1))+(1-rho)*exp(lUzplusU(+1)));
%exp(lVzplus2U)-exp(lG2U)*exp(lwU(-2))*exp(lwbarU(-2))*((1-exp(ltauyU))/(1+tauw))*exp(lintensityU) +exp(lzetahU)*AL*exp(lintensityU)^(1+sigmaL)/((1+sigmaL)*exp(lpsizplusU))-betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU))*(rho*exp(lVzplus3U(+1))+(1-rho)*exp(lUzplusU(+1)));
%exp(lVzplus3U)-exp(lG3U)*exp(lwU(-3))*exp(lwbarU(-3))*((1-exp(ltauyU))/(1+tauw))*exp(lintensityU) +exp(lzetahU)*AL*exp(lintensityU)^(1+sigmaL)/((1+sigmaL)*exp(lpsizplusU))-betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU))*(rho*exp(lVzplus0U(+1))+(1-rho)*exp(lUzplusU(+1)));

%61 Uzplus : value of being unemployed for the worker      
%exp(lUzplusU)-bu*(1-exp(ltauyU))
%-betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU))*(exp(lfU)*exp(lVxzplusU(+1))+(1-exp(lfU))*exp(lUzplusU(+1)));

%62 job finding rate f
%exp(lfU)-exp(lmU)/(1-exp(lLU));

%63 Vxzplus : average value of finding a job 
%exp(lVxzplusU)-1/exp(lmU(-1))*(
%exp(lchi0U(-1))*exp(ll0U(-1))*exp(lVzplus1U)
%+exp(lchi1U(-1))*exp(ll1U(-1))*exp(lVzplus2U)
%+exp(lchi2U(-1))*exp(ll2U(-1))*exp(lVzplus3U)
%+exp(lchi3U(-1))*exp(ll3U(-1))*exp(lVzplus0U)
%);

%64 BigQ Probability of filling a vacancy
%exp(lBigQU)=exp(lmU)/exp(lvU);

%65 M total new matches      
%exp(lmU)= exp(lchi0U)*exp(ll0U)+exp(lchi1U)*exp(ll1U)
%+exp(lchi2U)*exp(ll2U) +exp(lchi3U)*exp(ll3U);

%66 Matching function
%exp(lmU)=(exp(lsigmamU)*((1-exp(lLU))^sigma_match)*(exp(lvU)^(1-sigma_match)));

%67 Aggregate leisure/intensity decision
%((1-exp(ltauyU))/(1+tauw))*exp(lwbarU)=exp(lzetahU)*AL*(exp(lintensityU)^sigmaL)/exp(lpsizplusU);

%68-71 L.O.M. for employees in each cohort 
%exp(ll0U)-(exp(lchi3U(-1))+rho)*exp(ll3U(-1));
%exp(ll1U)-(exp(lchi0U(-1))+rho)*exp(ll0U(-1));
%exp(ll2U)-(exp(lchi1U(-1))+rho)*exp(ll1U(-1));
%exp(ll3U)-(exp(lchi2U(-1))+rho)*exp(ll2U(-1));

%72 marginal value of a job for the worker wrt. wage
%exp(lVwU)=((exp(lintensityU)*((1-exp(ltauyU))/(1+tauw))
%     +(rho*betaa)*(exp(lpsizplusU(+1))/exp(lpsizplusU)*exp(lG1U(+1))*exp(lintensityU(+1))*((1-exp(ltauyU(+1)))/(1+tauw)))
%     +(rho*betaa)^2*(exp(lpsizplusU(+2))/exp(lpsizplusU)*exp(lG2U(+2))*exp(lintensityU(+2))*((1-exp(ltauyU(+2)))/(1+tauw)))
 %    +(rho*betaa)^3*(exp(lpsizplusU(+3))/exp(lpsizplusU)*exp(lG3U(+3))*exp(lintensityU(+3))*((1-exp(ltauyU(+3)))/(1+tauw)))));

%73 Nash bargaining foc      
%exp(letaU)*exp(lVwU)*exp(lJzplusU)+(1-exp(letaU))*(exp(lVzplus0U)-exp(lUzplusU))*JwU;


%74 Jw marginal value of job match for a firm for a change in the wage 
%JwU=((-1)*(exp(lintensityU)
 %    +betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU)*exp(lG1U(+1))*exp(lOm1U)*exp(lintensityU(+1)))
 %    +betaa^2*(exp(lpsizplusU(+2))/exp(lpsizplusU)*exp(lG2U(+2))*exp(lOm2U(+1))*exp(lintensityU(+2))) 
%     +betaa^3*(exp(lpsizplusU(+3))/exp(lpsizplusU)*exp(lG3U(+3))*exp(lOm3U(+2))*exp(lintensityU(+3)))));

%75 Total hours
%exp(lHU)=exp(lintensityU)*exp(lLU);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% *** end of labor block *************%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
%START OLD CALVO STICKY WAGE EQUATIONS NOT APPLICABLE UNDER UNEMPLOYMENT%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
%Non-linear pricing equations for wage setting
%39 eq15a
exp(lFwU)-1/lambdaw*exp(lpsizplusU-lwhaloU*lambdaw/(1-lambdaw)+lsmallhU)*((1-exp(ltauyU))/(1+tauw))
  - betaa*xiw*exp(lwbarU(+1)-lwbarU+(1/(1-lambdaw))*lpitildewpiU(+1)+lFwU(+1));
%40 eq15b
exp(lKwU)-exp(lzetahU+(-lwhaloU*lambdaw/(1-lambdaw)+lsmallhU)*(1+sigmaL))-betaa*xiw*exp(lpitildewpiU(+1)*lambdaw*(1+sigmaL)/(1-lambdaw)+lKwU(+1));
%41 eq15c
exp(lKwU-lFwU)=exp(lwbarU)*(((1-xiw*exp(lpitildewpiU/(1-lambdaw)))/(1-xiw))^(1-lambdaw*(1+sigmaL)))/AL;
%42 eq15d
exp(lwhaloU*lambdaw/(1-lambdaw))-(1-xiw)*((1-xiw*exp(lpitildewpiU/(1-lambdaw)))/(1-xiw))^lambdaw
   -xiw*(exp(lpitildewpiU+lwhaloU(-1)))^(lambdaw/(1-lambdaw));
%43 eq15e
exp(lpitildewpiU)=exp(-lpiwU+kappaw*lpicU(-1)+(1-kappaw-varkappaw)*lpitargetU+varkappaw*log(pibreve)+thetaw*log(muzplus));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
%END OLD CALVO STICKY WAGE EQUATIONS NOT APPLICABLE UNDER UNEMPLOYMENT%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/


%76 eq16 - Taylor Rule
% comment this if you want to calculate Ramsey optimal policy
% *** w. measured GDP ***
%lRU-log(R)=rhoR*(lRU(-1)-log(R))+(1-rhoR)*(lpitargetU-log(pitarget)+
% rpi*(lpicU(-1)-lpitargetU)+
%  ry*(lylessdU(-1)-log(ylessd))+rq*(lqU(-1)-log(q)))    
%  +rdeltapi*(lpicU-lpicU(-1))+rdeltay*(lylessdU-lylessdU(-1))+epsR_eps/100;

% ** Below Taylor as in closed ec paper: leading inflation, current output**
% *** BUT contemp inflation instead of leading
% lRU-log(R)=rhoR*(lRU(-1)-log(R))+(1-rhoR)*(lpitargetU-log(pitarget)+
%   rpi*(lpicU -lpitargetU)+
%   ry*(lylessdU-log(ylessd)) +
%   rRstar*(lRstarU - lRstarU(-1)))     %rs*data_sdevU
%    + epsR_eps/10; 
   
   %gb: originally epsR_eps/100
%+rq*(lqU-log(q)))    
%  +rdeltapi*(lpicU-lpicU(-1))+rdeltay*(lylessdU-lylessdU(-1));

%77 eq17 - Adjusted Resource Constraint
exp(lyU)=exp(lgU)+(1-omegac)*(exp(lpcU))^(etac)*exp(lcU)+(1-omegai)*(exp(lpinvestU))^(etai)*
(exp(liU)+aofuU*exp(lkbarU(-1))/(exp(lmupsiU)*exp(lmuzplusU))) % *** old mistake (without utilization), current line was just: exp(liU) 
+(omegax*exp(lpmxU)^(1-etax)+1-omegax)^(etax/(1-etax))*     %old factor eliminated: exp(lRxU)^(etax)
(1-omegax)*exp(lphaloxU*lambdax/(1-lambdax))*exp(lpxU)^(-etaf)*exp(lystarU)
+mu*normcdf((lomegabarU +exp(lsigmaU(-1))^2/2)/exp(lsigmaU(-1)) -  exp(lsigmaU(-1)),0,1)
*exp(lRkU+lpkprimeU(-1)+lkbarU(-1)-lpiU-lmuzplusU);

%+((kappa/2)*(exp(lvtilde0U)^2*exp(ll0U)+exp(lvtilde1U)^2*exp(ll1U)+exp(lvtilde2U)^2*exp(ll2U) 
%+ exp(lvtilde3U)^2*exp(ll3U)));
%

%78 eq38 - Current Account
aU+exp(lqU+lpcU+lRnustarU)*( exp(lcmU+(lambdamc/(1-lambdamc))*lphalomcU) 
  +exp(limU+(lambdami/(1-lambdami))*lphalomiU) 
  +exp(lxmU+(lambdamx/(1-lambdamx))*lphalomxU) )
  -exp(lqU+lpcU+lpxU+lxU)-exp(lRstarU(-1)-phitildea*(aU(-1)-a)
  -phitildes*(exp(lRstarU(-1))-exp(lRU(-1))-(Rstar-R))
  +phitildeU(-1))*exp(lsU_eps-lpiU-lmuzplusU)*aU(-1); 

%79 eq19 - marginal costs of final export goods
exp(lmcxU)=exp(ltauxU)*exp(lRxU)/(exp(lqU)*exp(lpcU)*exp(lpxU))*
(omegax*exp(lpmxU)^(1-etax)+1-omegax)^(1/(1-etax));

%80 eq20a - marginal costs of consumption importers
exp(lmcmcU)=exp(ltaumcU)*exp(lqU)*exp(lpcU)*exp(lRnustarU)/exp(lpmcU);

%81 eq20b - marginal costs of investment importers
exp(lmcmiU)=exp(ltaumiU)*exp(lqU)*exp(lpcU)*exp(lRnustarU)/exp(lpmiU);

%82 eq20c - marginal costs of export importers
exp(lmcmxU)=exp(ltaumxU)*exp(lqU)*exp(lpcU)*exp(lRnustarU)/exp(lpmxU);

%83 eq21 - Real GDP, defined from the production side
exp(lyU)=exp(lphaloU)^(lambdad/(lambdad-1))*
        %((exp(lsmallhU-(lambdaw/(1-lambdaw))*lwhaloU))^(1-alphaa)*exp(lepsilonU)*
        ((exp(lHU))^(1-alphaa)*exp(lepsilonU)*
        (exp(lkU)/(exp(lmuzplusU)*exp(lmupsiU)))^(alphaa)-phi);

%84 eq41 Capacity utilization
aofuU = 0.5*sigmab*sigmaa*(exp(luU))^2+sigmab*(1-sigmaa)*exp(luU)
        +sigmab*((sigmaa/2)-1);

%85 eq22 - Definition of Rf
exp(lRfU)=nuf*(exp(lRU))+1-nuf;

%86 eq23 - Definition of Rnustar
exp(lRnustarU)=nustar*(exp(lRstarU))+1-nustar;

%87 eq24 - Total foreign export demand
exp(lxU)=(exp(lpxU))^(-etaf)*exp(lystarU);

%88 eq25 - Relative price of final consumption good
exp(lpcU)=(1-omegac+omegac*(exp(lpmcU))^(1-etac))^(1/(1-etac));

%89 eq26 - Relative price of final investment good
exp(lpinvestU)=(1-omegai+omegai*(exp(lpmiU))^(1-etai))^(1/(1-etai));

%90 eq27 - Definition of Rx
exp(lRxU)=nux*(exp(lRU))+1-nux;

%91 eq28 - Definition of capital utilitzation
exp(lkU)=exp(lkbarU(-1)+luU);

%92 eq29 - Dynamics of pmx
exp(lpmxU)/exp(lpmxU(-1))=exp(lpimxU)/exp(lpiU);

%93 eq30 - Dynamics of pmc
exp(lpmcU)/exp(lpmcU(-1))=exp(lpimcU)/exp(lpiU);

%94 eq31 - Dynamics of pmi
exp(lpmiU)/exp(lpmiU(-1))=exp(lpimiU)/exp(lpiU);

%95 eq32 - Dynamics of px
exp(lpxU)/exp(lpxU(-1))=exp(lpixU)/exp(lpistarU);

%96 eq33 - Dynamics of the real exchange rate
exp(lqU)/exp(lqU(-1))=exp(lsU_eps)*exp(lpistarU)/exp(lpicU);

%97 eq42 Dummy variable for Investment cost function
StildeU=0.5*(exp(sqrt(Spp)*(exp(liU)/exp(liU(-1))*exp(lmuzplusU)*exp(lmupsiU)-muzplus*mupsi))+
exp(-sqrt(Spp)*(exp(liU)/exp(liU(-1))*exp(lmuzplusU)*exp(lmupsiU)-muzplus*mupsi))-2);

%98 eq42b Dummy variable for Investment cost function derivative
SprimetildeU=0.5*sqrt(Spp)*(exp(sqrt(Spp)*(exp(liU)/exp(liU(-1))*exp(lmuzplusU)*exp(lmupsiU)-muzplus*mupsi))
-exp(-sqrt(Spp)*(exp(liU)/exp(liU(-1))*exp(lmuzplusU)*exp(lmupsiU)-muzplus*mupsi)));

%99 eq44 definition of rate of return on capital
exp(lRkU)=((1-tauk)*(exp(luU)*exp(lpinvestU)*(sigmab*sigmaa*(exp(luU))
 +sigmab*(1-sigmaa))
 -exp(lpinvestU)*aofuU) +(1-delta)*exp(lpkprimeU) 
 +delta*tauk*exp(lmupsiU-lpiU+lpkprimeU(-1)))
  /(exp(lpkprimeU(-1)+lmupsiU-lpiU));

%100 eq 45 Relation between smallh and  H 
exp(lsmallhU)=exp(lwhaloU)^(lambdaw/(1-lambdaw))*exp(lHU);

%101 eq 46 Imported consumption
exp(lcmU)=exp(lcU)*omegac*(exp(lpcU)/exp(lpmcU))^etac;

%102 eq 47 Imported investment
exp(limU)=omegai*(exp(lpinvestU)/exp(lpmiU))^etai
*(exp(liU)+aofuU*exp(lkbarU(-1))/(exp(lmupsiU)*exp(lmuzplusU)));

%103 eq49 Imported export inputs
exp(lxmU)-omegax*(((omegax*exp(lpmxU*(1-etax))
  +1-omegax)^(1/(1-etax))*exp(-lpmxU))^etax)*
  exp(-(lambdax/(lambdax-1))*lphaloxU-etaf*lpxU+lystarU);

%104 eq15f
lpiwU=lwbarU+lmuzplusU+lpiU-lwbarU(-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
% equations for the law of motions of exogenous processes %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/

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

%106-127 
ltaudU -log(taud) = rhotaud*(ltaudU(-1)-log(taud)) + taud_eps; %gb: originally taud_eps/10
ltauxU -log(taux) = rhotaux*(ltauxU(-1)-log(taux)) + taux_eps; %gb: originally taux_eps/10
ltaumcU -log(taumc) = rhotaumc*(ltaumcU(-1)-log(taumc)) + taumc_eps; %gb: originally taumc_eps/10
ltaumiU -log(taumi) = rhotaumi*(ltaumiU(-1)-log(taumi)) + taumi_eps; %gb: originally taumi_eps/10
ltaumxU -log(taumx) = rhotaumx*(ltaumxU(-1)-log(taumx)) + taumx_eps; %gb: originally taumx_eps/10 

% *** All below shocks enter by some factor in order to obtain same order of 
% *** magnitude for the optimization

lepsilonU =(1-rhoepsilon)*log(epsilon)+rhoepsilon*lepsilonU(-1)+epsilon_eps/10; 
lUpsilonU =(1-rhoUpsilon)*log(Upsilon)+rhoUpsilon*lUpsilonU(-1)+Upsilon_eps; %gb: originally Upsilon_eps/10 
 
lzetacU =(1-rhozetac)*log(zetac)+rhozetac*lzetacU(-1)+zetac_eps; % gb: originally zetac_eps/10  
lzetahU =(1-rhozetah)*log(zetah)+rhozetah*lzetahU(-1)+zetah_eps; % gb: originally zetah_eps/10  

   
%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; 
                     %gb: originally Rstar_eps/1000

%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;

%Fiscal shocks
ltauyU =(1-rhotauy)*log(tauy)+rhotauy*ltauyU(-1)+tauy_eps/100; 
lgU =(1-rhog)*log(etag*ybar)+rhog*lgU(-1)+g_eps/10; %gb: originally g_eps/100

%Obs: zero mean
phitildeU =(1-rhophitilde)*phitilde+rhophitilde*phitildeU(-1)+phitilde_eps/100; 
  
% inflation target shock
lpitargetU =(1-rhopi)*log(pitarget)+rhopi*lpitargetU(-1)+pitarget_eps/1000; 

% Fin fric shocks
lsigmaU=(1-rhosigma)*log(sigma) + rhosigma*lsigmaU(-1)+sigma_eps/100;
lgammaU=(1-rhogamma)*log(gammaa) + rhogamma*lgammaU(-1)+gamma_eps/10; %gb: originally gamma_eps/100

% Labor market shocks
%lsigmamU=(1-rhosigmam)*log(sigmam) + rhosigmam*lsigmamU(-1)+sigmam_eps/100;
%letaU=(1-rhoeta)*log(eta) + rhoeta*letaU(-1)+eta_eps/100;


%Measurement equations 128-142
%linking data and model variables
data_pidU-data_pidU_eps+ 400*log(pibar)=400*log(exp(lpiU)); 
data_piiU-data_piiU_eps+ 400*log(pii)=400*log(exp(lpiiU));   
data_picU-data_picU_eps+400*log(pic)=400*log(exp(lpicU)); 
data_RU=400*(exp(lRU)-1);

data_ystardiffU-data_ystardiffU_eps+100*log(muzplus)=100*(lmuzplusU+(lystarU - lystarU(-1)));
data_pistarU-data_pistarU_eps=400*log(exp(lpistarU));
data_RstarU=400*(exp(lRstarU)-1);

data_HdiffU-data_HdiffU_eps=100*(lsmallhU-lsmallhU(-1));

data_qdiffU-data_qdiffU_eps=100*(lqU-lqU(-1));
data_cdiffU-data_cdiffU_eps+100*log(muzplus)=100*(lmuzplusU+(lcU - lcU(-1)));

%  exclude utilization costs, recruitment costs and monitoring costs from GDP,
% done by using lylessdU instead of lyU in GDP measurement eq.
data_ydiffU-data_ydiffU_eps+100*log(muzplus)=100*(lmuzplusU+log(exp(lylessdU) ) -
                                                            log(exp(lylessdU(-1)) ) );

% *** version w/o cap. util in meas. eq.
data_idiffU-data_idiffU_eps+100*(log(muzplus)+log(mupsi))=100*(lmuzplusU+lmupsiU+ liU - liU(-1));

data_wdiffU-data_wdiffU_eps+100*log(muzplus)=100*(lmuzplusU+(lwbarU-lwbarU(-1)));

data_xdiffU-data_xdiffU_eps+100*log(muzplus)=100*(lmuzplusU+(lxU - lxU(-1)));
data_impdiffU-data_impdiffU_eps+100*log(muzplus)=100*(lmuzplusU
             +log(exp(lcmU+(lambdamc/(1-lambdamc))*lphalomcU) 
             +exp(limU+(lambdami/(1-lambdami))*lphalomiU)
             +exp(lxmU+(lambdamx/(1-lambdamx))*lphalomxU)) 
             - log(exp(lcmU(-1)+(lambdamc/(1-lambdamc))*lphalomcU(-1))
             + exp(limU(-1)+(lambdami/(1-lambdami))*lphalomiU(-1)) 
             +exp(lxmU(-1)+(lambdamx/(1-lambdamx))*lphalomxU(-1))));  
 
data_ndiffU-data_ndiffU_eps+100*log(muzplus)=100*(lmuzplusU+(lnU - lnU(-1)));
data_spreaddiffU-data_spreaddiffU_eps=100*(lspreadU-lspreadU(-1));

data_spreadU-data_spreadU_eps=400*exp(lspreadU);
%gb2013 demeaned
% data_spreadU-data_spreadU_eps=400*(exp(lspreadU)-spread);

data_gdiffU-data_gdiffU_eps+100*log(muzplus)=100*(lmuzplusU+(lgU - lgU(-1)));

data_HhatU-data_HhatU_eps=100*(exp(lsmallhU)-smallh)/smallh; % *** added KW

%exp(lUnemprateU)=1-exp(lLU);
exp(lspreadU)=exp(lomegabarU(+1))*exp(lRkU(+1))/(1-exp(lnU)/(exp(lpkprimeU)*exp(lkbarU)))-exp(lRU);
exp(lbankruptcyrateU)=normcdf((lomegabarU(+1)+exp(lsigmaU)^2/2)/exp(lsigmaU),0,1);

% GDP w/o utilization costs, monitoring costs (or recruitment costs):
exp(lylessdU)=exp(lyU)-(mu*normcdf((lomegabarU +exp(lsigmaU(-1))^2/2)/exp(lsigmaU(-1)) -  exp(lsigmaU(-1)),0,1)
*exp(lRkU+lpkprimeU(-1)+lkbarU(-1)-lpiU-lmuzplusU) )
-exp(lpinvestU)*aofuU*exp(lkbarU(-1))/(exp(lmupsiU)*exp(lmuzplusU));

% net exports/Y: (value of exports - value of imports)/Y
nxdivGDPU= ( exp(lqU+lpcU+lpxU+lxU) - ( exp(lqU+lpcU+lRnustarU)*( exp(lcmU+(lambdamc/(1-lambdamc))*lphalomcU) 
  +exp(limU+(lambdami/(1-lambdami))*lphalomiU) 
  +exp(lxmU+(lambdamx/(1-lambdamx))*lphalomxU) ) ) ) / exp(lyU);

% The nominal x-rate in levels
exp(lSlevU)/exp(lSlevU(-1))=exp(lsU_eps);

% gb2013: observed deviation of nominal x-rate
%exp(lSlevU)/exp(lSlevU(-1))=exp(data_sdevU)/exp(data_sdevU(-1));
%data_sdevU - data_sdevU(-1) - data_sdevU_eps= lsU;

% Cumulated GDP-growth, up to 12 Q
GDP12q_diffU=12*100*log(muzplus)+data_ydiffU+data_ydiffU(-1)+data_ydiffU(-2)+data_ydiffU(-3)+data_ydiffU(-4)+data_ydiffU(-5)+data_ydiffU(-6)
            +data_ydiffU(-7)+data_ydiffU(-8)+data_ydiffU(-9)+data_ydiffU(-10)+data_ydiffU(-11);
GDP4q_diffU =4*100*log(muzplus)+data_ydiffU+data_ydiffU(-1)+data_ydiffU(-2)+data_ydiffU(-3);
pic4qU      =(data_picU+data_picU(-1)+data_picU(-2)+data_picU(-3))/4;

% credit amount in levels, gb2014
exp(lbU) = exp(lpkprimeU+lkbarU)-exp(lnU);
% net assets to gdp ratio, gb2014
etaaU = aU/exp(lylessdU);
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%**** end of model *****%%%%%%%%%%%%%%%%%%%%%%/
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


initval;
lmcU=log(mc); 
lrkbarU=log(rkbar); 
lwbarU=log(wbar);
lRfU=log(Rf);
lpxU=log(px);
lpixU=log(pix);
lkU=log(kbar);
lRnustarU=log(Rnustar); 
lHU=log(H); 
lpiU=log(pibar); 
lmcxU=log(mcx);
lpimcU=log(pimc);
lmcmcU=log(mcmc);
lpimxU=log(pimx);
lmcmxU=log(mcmx);
lpimiU=log(pimi); 
lmcmiU=log(mcmi); 
lpicU=log(pic); 
lpmcU=log(pmc); 
lpiiU=log(pii);
lpmiU=log(pmi);
lkbarU=log(kbar);
liU=log(i);
lpsizplusU=log(psizplus); 
lcU=log(c); 
lpkprimeU=log(pkprime);
luU=log(1);
lpinvestU=log(pinvest);
%lsU=log(s); 
aU=a; 
lqU=log(q);
lxU=log(x); 
lpcU=log(pc);
lRU=log(R); 
lyU=log(ybar);  
lRxU=log(Rx); 
lpmxU=log(pmx);
lRkU=log(Rk);
StildeU=0; 
SprimetildeU=0;
lmuzplusU=log(muzplus); 
lepsilonU =log(epsilon);  
lmupsiU =log(mupsi);  
lpitargetU =log(pibar);  
lUpsilonU =log(Upsilon);  
lzetacU =log(zetac);    
lRstarU =log(R);  
phitildeU =0;  
ltauyU =log(tauy);    
lzetahU =log(zetah);  
lystarU =log(ystar);   
lpistarU =log(pibar);  
lgU =log(etag*ybar); 
lmuzU=log(muz);


% Finfric variables:
lgammaU=log(gammaa);
lsigmaU=log(sigma);
lnU=log(n);
lomegabarU=log(omegabar);

data_pidU=0; %400*log(pibar);
data_RU=400*log(R);
data_wdiffU=0;
data_cdiffU=0;
data_idiffU=0;
data_qdiffU=0;    
data_HdiffU=0;
data_ydiffU=0;
data_xdiffU=0;
data_impdiffU=0;
data_picU=0; %400*log(pic);
data_piiU=0; %400*log(pii);
data_ystardiffU=0;
data_pistarU=400*log(pistar);
data_RstarU=400*log(Rstar);
data_ndiffU=0;
%data_unempdiffU=0;
data_spreadU=400*spread;
data_gdiffU=0;

lpitildedpiU=log(pitildedpi);
lKdU=log(Kd);
lFdU=log(Fd);
lphaloU=log(phalo); 
lpitildexpiU=log(pitildexpi);
lKxU=log(Kx);
lFxU=log(Fx);
lphaloxU=log(phalox);
lpitildemcpiU=log(pitildemcpi);
lKmcU=log(Kmc);
lFmcU=log(Fmc);
lphalomcU=log(phalomc);
lcmU=log(cm);
lpitildemipiU=log(pitildemipi);
lKmiU=log(Kmi);
lFmiU=log(Fmi);
lphalomiU=log(phalomi);
limU=log(im);
lpitildemxpiU=log(pitildemxpi);
lKmxU=log(Kmx);
lFmxU=log(Fmx);
lphalomxU=log(phalomx);
lxmU=log(xm);

lpitildewpiU=log(pitildewpi);
lKwU=log(Kw);
lFwU=log(Fw);
lwhaloU=log(whalo);
%lwhalohaloU=log(whalohalo);
lsmallhU=log(smallh);

aofuU=aofu;


ltaudU=log(taud);
ltauxU=log(taux);
ltaumcU=log(taumc);
ltaumiU=log(taumi);
ltaumxU=log(taumx);
lpiwU=log(piw);

lspreadU=log(spread);
%lUnemprateU=log(1-L);
lbankruptcyrateU=log(Fomegabar);
lylessdU=log(ylessd);
%data_sdevU=0;
lbU=log(pkprime*kbar-n);
etaaU=a/ylessd;
end;

resid(1); 

steady;

check; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
%BEGINNING OF ESTIMATION%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/


shocks;
var muz_eps=sig_muz^2;
var mupsi_eps = 0; % sig_mupsi^2;     
var epsilon_eps =sig_epsilon^2;
var Upsilon_eps =sig_upsilon^2;
var zetac_eps =sig_zetac^2;
var zetah_eps =sig_zetah^2;
var phitilde_eps =sig_phitilde^2;
var pitarget_eps=sig_pitarget^2;  
var epsR_eps =sig_epsR^2;
var tauy_eps =sig_tauy^2;       
var g_eps=sig_g^2;
var ystar_eps= sig_ystar^2;
var pistar_eps=sig_pistar^2;
var Rstar_eps= sig_Rstar^2;

%fin fric innovations  
var gamma_eps=sig_gamma^2;  
var sigma_eps=0; % sig_sigma^2; 

%unemployment innovations
%var sigmam_eps=sig_sigmam^2;
%var eta_eps=sig_eta^2;
%var sigma_a_eps=sig_sigma_a^2;
 
% Markup variances:  
var taud_eps = sig_taud^2; 
var taux_eps = sig_taux^2;
var taumc_eps =sig_taumc^2;
var taumi_eps =sig_taumi^2;
var taumx_eps =sig_taumx^2;

% measurement errors. Updated 10q3

% 10 percent of variance of time series
var data_ydiffU_eps= 	0.5354; %0.106;
var data_pidU_eps= 	 7.0402; %0.374;     % released by GB2013
var data_wdiffU_eps= 	 0.5534; %0.071;
var data_cdiffU_eps= 	0.8045; % b/c outlier deleted 0.9611; % 0.061;
var data_idiffU_eps= 	26.6239; % 0.639;
var data_qdiffU_eps= 	0.6311; % 0.706;
var data_HdiffU_eps= 	0.4819; % 0.035;
var data_xdiffU_eps= 	1.1653; % 0.620;
var data_impdiffU_eps= 	3.9705; % 0.731;
var data_picU_eps= 	3.9595; % 0.203;      % released by GB2013
var data_piiU_eps= 	264.7161; % 0.304;     % released by GB2013
var data_ndiffU_eps= 	11.1925; % 11.170;
var data_spreaddiffU_eps= 247.7939; % 18.520;
var data_spreadU_eps= 0.50664; %2.0266; % b/c multiplied by 2 not 4 %8.1063;  

var data_gdiffU_eps= 	2.9773; % 0.053;
var data_HhatU_eps= 	7.1971;
%var data_sdevU_eps= 	6.5e-06;

%var data_unempdiffU_eps= 7.3042; %	1.777;
%var data_unemphatU_eps= 	46.223;

%25 percent of variance
%var data_pidU_eps= 31.4531; %  0.9355;
%var data_picU_eps= 9.8987; % 0.5069;
%var data_piiU_eps= 182.7017; % 0.7603;

end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
% *** estim params priors begin ****%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/

estimated_params; 
% ********** begin pasting from CTW.mod, 
% **    CTW_110208_2010MASTER_but_large_pi-MEs2 based on full sample, -2010q3
% *****************************************************


%shocks
% init-val, prior shape, mean, std;
% *** all shock priors re-scaled by some factor; see the above equations
% stderr muz_eps,      .44,inv_gamma_pdf,0.15,inf;  
% *** stderr mupsi_eps,    inv_gamma_pdf,0.15,inf; ***
stderr epsilon_eps,  .123,inv_gamma_pdf,.15,inf; % rescaled
stderr Upsilon_eps,  .15,inv_gamma_pdf,0.15,inf;  
stderr zetac_eps,    .21,inv_gamma_pdf,0.15,inf;    
stderr zetah_eps,    .483,inv_gamma_pdf,0.5,inf;        
stderr phitilde_eps, .539,inv_gamma_pdf,0.5,inf;     
%stderr epsR_eps,     .22,inv_gamma_pdf,0.15,inf;
stderr g_eps,        .464,   inv_gamma_pdf,0.5,inf;
stderr taud_eps,     .313,inv_gamma_pdf,0.5,inf; % old: 0.26136; % 4 df obtained using [s,nu]=inverse_gamma_specification(0.5,0.26136,1)
stderr taux_eps,     .95,inv_gamma_pdf,0.5,inf; % old:0.26136; 
stderr taumc_eps,    .599,inv_gamma_pdf,0.5,inf; % old: 0.26136; 
stderr taumi_eps,    .389,inv_gamma_pdf,0.5,inf; % old:0.26136; 
stderr taumx_eps,    1.1,inv_gamma_pdf,0.5,inf; % old:0.26136;  
stderr gamma_eps,    .282,inv_gamma_pdf,0.5,inf; 
% stderr sigma_eps,    inv_gamma_pdf,0.5,inf; 
% stderr eta_eps,         inv_gamma_pdf,0.5,inf;
% stderr ystar_eps, .14,inv_gamma_pdf,0.5,inf;  
% stderr pistar_eps, .3,inv_gamma_pdf,0.5,inf;  
% stderr Rstar_eps, .22,inv_gamma_pdf,0.15,inf; % *** REVISED from 0.5 % scaled differently, by 1000
                                         % gb: revised back to 0.15   


% *** stickiness ****
xid,     .797, beta_pdf,0.75,0.075;
xix,     .87,beta_pdf,0.75,0.075;
ximc,    .79,beta_pdf,0.75,0.075;
ximi,    .39,beta_pdf,0.65,0.075;
ximx,    .57,beta_pdf,0.65,0.1; 

kappad,     .132,beta_pdf,0.4,0.15;
kappax,     .29,beta_pdf,0.4,0.15;
kappamc,    .37,beta_pdf,0.4,0.15;
kappami,    .21,beta_pdf,0.4,0.15;
kappamx,    .302,beta_pdf,0.4,0.15;

kappaw,     .2,beta_pdf,0.4,0.15;

work_cap_para, .4,beta_pdf,0.5,0.25;

% preference and technology 
sigmaLScaled,.16, gamma_pdf,  0.3, 0.15;  % CTW: 0.75,   0.2; 
b,           .9,beta_pdf,   0.65,   0.15;
SppScaled,   .15,gamma_pdf,  0.5,    0.15;
sigmaa,   .7,gamma_pdf,  0.2,    0.075;  % Prior moments, but not function, from Smets and Wouters, JEEA 2003 (they use normal(0.2;0.075). old: (0.5; 0.25);  , then (0.2,0.15)


% *** Taylor rule ***
% changed to identical priors to SW 2003.
% rhoR,    .73,beta_pdf,  0.8,     0.1;    % CTW closed: (0.8;0.1)
% rpi,     1.8,normal_pdf,1.7,     0.15;   % CTW closed: (1.6;0.15)  
% % rdeltapi, normal_pdf,0.3,     0.1;
% ry,      .005,beta_pdf,0.1,   0.1;     % SW 2003: 0.125,  0.05; CTW closed: (0.2;0.15)
% % rdeltay, gamma_pdf, 0.0625,   0.05;   % SW 2003: 0.0625,  0.05;
% %rs,         normal_pdf,1.2,   0.1;
% rRstar, .55,beta_pdf, .6, .1; 
%
% *** open econ. stuff
etax,         1.48,  gamma_pdf,1.5,0.25; % all etas trunc'ed below
etac,         1.26,  gamma_pdf,1.5,0.25;
%etai,         1.06,  gamma_pdf,1.5,0.25;
etaf,          1.57, gamma_pdf,1.5,0.25;
%phitildes,     .95, gamma_pdf,1.25,0.1;  % old: 1.5,0.15 


% *** financial frictions parameters ***
mu,  .27,beta_pdf, 0.3,  0.075;   
% wesharePercent,gamma_pdf,0.1,0.05;


% /*
% % *** labor market parameters ***
% recruitsharePercent, gamma_pdf,     0.1, 0.075;   
% bshare,           beta_pdf,      0.75,0.075;  
% % iota,       beta_pdf,         0.5, 0.25;    
% 
% BigFPercent,      beta_pdf,      0.25,0.05; 
% */

% *** Exog. AR1 processes ***
%rhomuz,    .58,  beta_pdf,0.5, 0.075; % REVISED from 0.5, 0.2. 
% *** rhomupsi,    beta_pdf,0.5, 0.1; ***
rhoepsilon,  .858,beta_pdf,0.85,0.075;
rhoUpsilon,  .5,beta_pdf,0.85,0.075;
rhozetac,    .82,beta_pdf,0.85,0.075; 
rhozetah,    .85,beta_pdf,0.85,0.075; 
rhophitilde, .947,beta_pdf,0.85,0.075; 
rhog,        .83,beta_pdf,0.85,0.075;
% rhosigma,    beta_pdf,0.85,0.075; 
rhogamma,    .823,beta_pdf,0.85,0.075;
% rhoeta,         beta_pdf,0.85,0.075;

% *** 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;
etax, 1.01,20;
etac, 1.01,20;
%etai, 1.01,20;
etaf, 1.01,20; 
end;

 





%%%%%%%%%%%%%%%%%%%%%%%%/
%%%%%%Observed variables%%%%%%%%%/
%%%%%%%%%%%%%%%%%%%%%%%%/

varobs  data_RU data_wdiffU data_cdiffU    
data_ydiffU  data_ystardiffU data_pistarU 
data_RstarU data_HdiffU % *** data_HdiffU ***
data_impdiffU data_piiU 
data_pidU data_picU data_xdiffU data_qdiffU data_idiffU
data_ndiffU data_spreadU %data_spreaddiffU
%data_unempdiffU 
data_gdiffU
 %data_sdevU
 ; % GDP12q_diffU GDP4q_diffU pic4qU; % gb2013 
 

%finfric_oo_irfs=oo_.irfs;
%save finfric_oo_irfs finfric_oo_irfs  

options_.console_mode=1; %(default: 0)

estimation(order=1,datafile=myOriginalData, mh_replic=250000, mh_nblocks=2, mh_drop=.8,  mh_jscale=0.2, 
 mode_compute=6, 
 filter_step_ahead = [1 4 8 12],forecast=20,
 smoother, filtered_vars, diffuse_filter,
 % get the IRF conf bands as otherwise stoch_simul does not provide the
 % bands:
 mh_conf_sig=.95,bayesian_irf,irf=20,
nodisplay, 
mode_check,
% 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))
% 
%lmcU lrkbarU lwbarU lRfU lpxU lpixU lkU lRnustarU 
%lHU lpiU lmcxU lpimcU lmcmcU lpimxU lmcmxU lpimiU 
%lmcmiU lpicU lpmcU lpiiU lpmiU lkbarU liU lpsizplusU 
%lcU lpkprimeU luU lpinvestU %lsU 
%aU %lqU lxU 
%lpcU  lRU lyU  lRxU lpmxU lRkU 
%lylessdU 
%lomegabarU lnU lspreadU 
lbankruptcyrateU

%law of motions for exogenous processes 
%lmuzplusU lepsilonU lmupsiU %lpitargetU 
%lUpsilonU lzetacU lRstarU phitildeU ltauyU lzetahU 
%lystarU lpistarU lgU lmuzU 
%ltaudU ltauxU ltaumcU ltaumiU ltaumxU

nxdivGDPU % 2
data_RU 
data_wdiffU 
data_cdiffU    
data_ydiffU  % 4-th
data_ystardiffU 
data_pistarU 
data_RstarU  
data_HdiffU 
data_impdiffU 
data_piiU 
data_pidU 
data_picU % 73-rd
data_xdiffU 
data_qdiffU 
data_idiffU 
data_gdiffU

data_ndiffU 
data_spreadU
%data_sdevU
lbU
etaaU
;


stoch_simul(order=1,irf=20, nodisplay, 
%order=2, pruning,periods=100,simul_replic=5,
graph_format = (eps, pdf))
%lmcU lrkbarU lwbarU lRfU lpxU lpixU lkU lRnustarU 
%lHU lpiU lmcxU lpimcU lmcmcU lpimxU lmcmxU lpimiU 
%lmcmiU lpicU lpmcU lpiiU lpmiU lkbarU liU lpsizplusU 
%lcU lpkprimeU luU lpinvestU %lsU 
%aU %lqU lxU 
%lpcU  lRU lyU  lRxU lpmxU lRkU 
%lylessdU 
%lomegabarU lnU lspreadU 
%lbankruptcyrateU

%law of motions for exogenous processes 
%lmuzplusU lepsilonU lmupsiU %lpitargetU 
%lUpsilonU lzetacU lRstarU phitildeU ltauyU lzetahU 
%lystarU lpistarU lgU lmuzU 
%ltaudU ltauxU ltaumcU ltaumiU ltaumxU

nxdivGDPU % 61
data_RU 
data_wdiffU 
data_cdiffU    
data_ydiffU  % 65-th
data_ystardiffU 
data_pistarU 
data_RstarU  
data_HdiffU 
data_impdiffU 
data_piiU 
data_pidU 
data_picU % 73-rd
data_xdiffU 
data_qdiffU 
data_idiffU 
data_gdiffU

data_ndiffU 
data_spreadU
%data_sdevU
%lbU
etaaU
; 

% 
shock_decomposition(parameter_set=posterior_mean)
data_ydiffU 
%data_picU
%data_spreadU 
%data_sdevU
;



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


/*
load('finfric125_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','ndiff',...
'spread'};

pid_ss=oo_.steady_state(95);
R_ss=oo_.steady_state(96);
Rstar_ss=oo_.steady_state(109); 
pic_ss=oo_.steady_state(105);
pii_ss=oo_.steady_state(106);
pistar_ss=oo_.steady_state(108);
spread_ss=oo_.steady_state(111);
mylevel = [0,pid_ss,R_ss,0,0,0,0,0,0,0,pic_ss,pii_ss,0,Rstar_ss,pistar_ss,0,0,spread_ss];


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','sigma','omegabar','Rk','kbar','wbar',...
    'H','ystar','rkbar','AL','varphi','gammaa'};
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
myParamVals
disp('Capital to output ratio:')
exp(oo_.steady_state(7)-oo_.steady_state(34))


*/

