
% 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
aofuU lpiwU

%variables for the unemployment part
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

% 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 lsigma_aU

%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_ndiffU data_unempdiffU data_spreadU data_gdiffU
lavg_NashU lylessdU

% Endog. separation variables
lepscall0U lepscall1U lepscall2U lepscall3U
lfcall0U lfcall1U lfcall2U lfcall3U
lintensity0U lintensity1U lintensity2U lintensity3U
lmcall0U lmcall1U lmcall2U lmcall3U     
la0U la1U la2U la3U
lgcallprime0U lgcallprime1U lgcallprime2U lgcallprime3U
lfcallprime0U lfcallprime1U lfcallprime2U lfcallprime3U
lJtildezplus0U lJtildezplus1U lJtildezplus2U lJtildezplus3U

% Additional observed variables
data_HhatU data_unemphatU
nxdivGDPU
lSlevU
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
sigma_a_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_unemphatU_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 

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

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 
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
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   Vw BigQ v m sigmam f w 
Jzplus Vxzplus Uzplus   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 work_cap_para wesharePercent
hours_target iy_target xy_target  nk_target 
avg_Nash ylessd

sigma_a 

epscall0  epscall1  epscall2  epscall3 fcall0  fcall1  fcall2  fcall3 
intensity0  intensity1  intensity2  intensity3 
mcall0  mcall1  mcall2  mcall3 a0  a1  a2  a3 

SW SE   % Switches for surplus: INDICATOR variables SW & SE for worker and employer surplus respectively
gcallprime0 gcallprime1 gcallprime2 gcallprime3
fcallprime0 fcallprime1 fcallprime2 fcallprime3
Jtildezplus0 Jtildezplus1 Jtildezplus2 Jtildezplus3

rhosigma_a sig_sigma_a

BigF BigFPercent ab 
sigmaLScaled
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;
options_.TeX_varobs=['data_RU         '
                     'data_wdiffU     '
                     'data_cdiffU     '
                     'data_ydiffU     '
                     'data_ystardiffU '
                     'data_pistarU    '
                     'data_RstarU     '
                     'data_HhatU      '
                     'data_impdiffU   '
                     'data_piiU       '
                     'data_pidU       '
                     'data_picU       '
                     'data_xdiffU     '
                     'data_qdiffU     '
                     'data_idiffU     '
                     'data_ndiffU     ' 
                     'data_spreaddiffU'
                     'data_unempdiffU ' 
                     'data_gdiffU     '];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% **** begin of parameter setting ****%%%%%%%%%%%%%%%%%%%%/
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Switches for surplus measure:
SW=0;
SE=1;  


alphaa=.4;  %to get a Pk*k/y around 8 
betaa=0.995; % *** so that real rate = 2.25, i.e. Real R Q=1.005573 

% 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); % .97 % Yields mupsi=1   
% *** residual, not set anymore: mupsiPercent=0.04;                  

%Steady state inflation rates 
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;	% Avg in 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=0.45; %.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.3; %.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;

% Some steady state values
u=1;
Stilde=0;
Sprimetilde=0;

% Persistence of shocks to marginal costs ("markup shocks")
rhotaux=0; 
rhotaumc=0; 
rhotaumi=0; 
rhotaumx=0; 
rhotaud=0; 

% The unused arguments in the Taylor rule (i.e. coefficient = 0)
rdeltapi=	0	;
rq      =	0	;
rdeltay	=	0	;

% Financial friction params
wesharePercent  =   0.1;   % we/ybar expressed in percent
Fomegabar       =   0.02;  % Based on UC (Credit registry) data, sample avg 95-09q2 

% Labor params
BigN            =   4;
sigma_match     =   0.5; % crucial parameter t oavoid fata errors; % 0.5
L               =   1-0.137; % .117% 0.078 sample avg 1995Q1-2010q3
iota            =	1; 

sigma_a         =   0.1; % dummy-setting, gets real value in ss-file  %standard deviation of idiosyncratic productivity of workers in steady state
                                    
rho             =   0.97; %.972  %  rho-BigF=0.9695 implies 1/f=3 quarters (if U=0.078)
sigmam          =   0.4;%.5%.57   % such that quarterly Q is roughly 0.9, implying monthly Q roughly 0.6 as in Swedish data

% shocks
sig_muz     =	0.15;
sig_epsilon	=	0.5	;
sig_upsilon	=	0.5	;
sig_zetac	=	0.15; 
sig_phitilde=	0.15;
sig_epsR	=	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_ystar	=	0.5	;
sig_pistar	=	0.5	;
sig_Rstar	=	.15; %gb: originally 0.5; % different scaling

sig_zetah	=	0.15;

%Following shocks are deactivated
sig_eta     =	0	;
sig_sigma_a =   0   ;
sig_sigmam  =	0; %0	;
sig_tauy	=	0	;
sig_pitarget=	0	;
sig_mupsi	=	0;
% sig_mupsi set to zero again below;
sig_sigma	=	0;
% sig_sigma set to zero again below

phitilde=0;  %OBS! Keep to zero! Always! 
rhopi=0.975;    % unused 


% ******************************************************************
% Below are the estimated parameters, and these are set to 
% their prior mean here.
% Param value 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  =	.75; % 0.2; % 0.75	;
b           =	0.65;
SppScaled	=	0.5	;
sigmaa      =	0.2	;

% *** Taylor rule ***			
rhoR	=	0.8	;
rpi     =	1.7	;
ry      =	0.125	;
rRstar  =   0.6; 

% *** open econ. stuff	
etai	=	 1.5 ; % fixed b/c estimated below 1
etaf	=	1.5	;
etac	=	1.5	;
etax	=	1.5	;
phitildes=	0; %1.25	; 

% *** financial frictions parameters ***	
mu              =	0.3;   

% *** labor params
recruitsharePercent	= 0.1;   % denoted hshare in the paper
bshare              = 0.75;  % unemployment benefits
BigFPercent         = 0.25;  % steady state endogenous separation rate, can't have any wage dispersion if BigF>0 

% *** 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	; 


% steady state targets to be matched
hours_target=.24; %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;    % 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 **** %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%**** 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;

% **************************************%%%%%%%%%%%%%%%%%%%%%
% ***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))*
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 eeq4 Jzplus - Stationary value of a job for the firm
-exp(lJzplusU) +(((exp(lwbarU)*exp(lepscall0U)/(1-exp(lfcall0U))-exp(lwU)*exp(lwbarU))*exp(lintensity0U)-(kappa/2)*(exp(lvtilde0U)^2)))*(1-exp(lfcall0U))  
+betaa*  (exp(lpsizplusU(+1))/exp(lpsizplusU))*((exp(lwbarU(+1))*exp(lepscall1U(+1))/(1-exp(lfcall1U(+1)))-exp(lG1U(+1))*exp(lwU)*exp(lwbarU))*exp(lintensity1U(+1))-(kappa/2)*(exp(lvtilde1U(+1))^2))*exp(lOm1U(+1))
+betaa^2*(exp(lpsizplusU(+2))/exp(lpsizplusU))*((exp(lwbarU(+2))*exp(lepscall2U(+2))/(1-exp(lfcall2U(+2)))-exp(lG2U(+2))*exp(lwU)*exp(lwbarU))*exp(lintensity2U(+2))-(kappa/2)*(exp(lvtilde2U(+2))^2))*exp(lOm2U(+2))
+betaa^3*(exp(lpsizplusU(+3))/exp(lpsizplusU))*((exp(lwbarU(+3))*exp(lepscall3U(+3))/(1-exp(lfcall3U(+3)))-exp(lG3U(+3))*exp(lwU)*exp(lwbarU))*exp(lintensity3U(+3))-(kappa/2)*(exp(lvtilde3U(+3))^2))*exp(lOm3U(+3))
+betaa^4*(exp(lpsizplusU(+4))/exp(lpsizplusU))*exp(lJzplusU(+4))*exp(lOm4U(+4))/(1-exp(lfcall0U(+4))); 

%40-42 eeq17 (or is it eeq17a?) 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 equation "up", Pitildew Inflation updater
exp(lpitildewU)=exp(kappaw*lpicU(-1)+(1-kappaw-varkappaw)*lpitargetU+varkappaw*log(pibreve)+thetaw*log(muzplus)); 

%44 eeq20  L-total employment
exp(lLU)=(1-exp(lfcall0U))*exp(ll0U) + (1-exp(lfcall1U))*exp(ll1U) + (1-exp(lfcall2U))*exp(ll2U) + (1-exp(lfcall3U))*exp(ll3U);

%45-48 eeq18 Omega updater
exp(lOm1U)=(1-exp(lfcall1U))*(exp(lchi0U(-1))+rho)*(1-exp(lfcall0U(-1))) ;
exp(lOm2U)=(1-exp(lfcall2U))*(exp(lchi1U(-1))+rho)*exp(lOm1U(-1));
exp(lOm3U)=(1-exp(lfcall3U))*(exp(lchi2U(-1))+rho)*exp(lOm2U(-1));
exp(lOm4U)=(1-exp(lfcall0U))*(exp(lchi3U(-1))+rho)*exp(lOm3U(-1));

%49-52 eeq1 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 eeq2 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)+lepscall1U(+1))-exp(lG1U(+1))*exp(lwU)*exp(lwbarU)*(1-exp(lfcall1U(+1))))*exp(lintensity1U(+1))
+kappa*(1-exp(lfcall1U(+1)))*( 1/2*(exp(lvtilde1U(+1))^2)+rho*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)+lepscall2U(+1))-exp(lG2U(+1))*exp(lwU(-1))*exp(lwbarU(-1))*(1-exp(lfcall2U(+1))))*exp(lintensity2U(+1))
+kappa*(1-exp(lfcall2U(+1)))*( 1/2*(exp(lvtilde2U(+1))^2)+rho*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)+lepscall3U(+1))-exp(lG3U(+1))*exp(lwU(-2))*exp(lwbarU(-2))*(1-exp(lfcall3U(+1))))*exp(lintensity3U(+1))
+kappa*(1-exp(lfcall3U(+1)))*( 1/2*(exp(lvtilde3U(+1))^2)+rho*exp(lvtilde3U(+1))/(exp(lBigQU(+1)))^(1-iota)) );

% eeq2a 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)+lepscall0U(+1))-exp(lwU(+1))*exp(lwbarU(+1))*(1-exp(lfcall0U(+1))))*exp(lintensity0U(+1))
+kappa*(1-exp(lfcall0U(+1)))*( 1/2*(exp(lvtilde0U(+1))^2)+rho*exp(lvtilde0U(+1))/(exp(lBigQU(+1)))^(1-iota)) );


%57-60 eeq6 Vzplus for cohorts - value of a job for a worker
exp(lVzplus0U)-exp(lwU)*exp(lwbarU)*((1-exp(ltauyU))/(1+tauw))*exp(lintensity0U)+exp(lzetahU)*AL*exp(lintensity0U)^(1+sigmaL)/((1+sigmaL)*exp(lpsizplusU))-betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU))*(rho*(1-exp(lfcall1U(+1)))*exp(lVzplus1U(+1))+(1-rho*(1-exp(lfcall1U(+1))))*exp(lUzplusU(+1)));
exp(lVzplus1U)-exp(lG1U)*exp(lwU(-1))*exp(lwbarU(-1))*((1-exp(ltauyU))/(1+tauw))*exp(lintensity1U) +exp(lzetahU)*AL*exp(lintensity1U)^(1+sigmaL)/((1+sigmaL)*exp(lpsizplusU))-betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU))*(rho*(1-exp(lfcall2U(+1)))*exp(lVzplus2U(+1))+(1-rho*(1-exp(lfcall2U(+1))))*exp(lUzplusU(+1)));
exp(lVzplus2U)-exp(lG2U)*exp(lwU(-2))*exp(lwbarU(-2))*((1-exp(ltauyU))/(1+tauw))*exp(lintensity2U) +exp(lzetahU)*AL*exp(lintensity2U)^(1+sigmaL)/((1+sigmaL)*exp(lpsizplusU))-betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU))*(rho*(1-exp(lfcall3U(+1)))*exp(lVzplus3U(+1))+(1-rho*(1-exp(lfcall3U(+1))))*exp(lUzplusU(+1)));
exp(lVzplus3U)-exp(lG3U)*exp(lwU(-3))*exp(lwbarU(-3))*((1-exp(ltauyU))/(1+tauw))*exp(lintensity3U) +exp(lzetahU)*AL*exp(lintensity3U)^(1+sigmaL)/((1+sigmaL)*exp(lpsizplusU))-betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU))*(rho*(1-exp(lfcall0U(+1)))*exp(lVzplus0U(+1))+(1-rho*(1-exp(lfcall0U(+1))))*exp(lUzplusU(+1)));


%61 eeq7 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 eeq12 job finding rate f
exp(lfU)-exp(lmU)/(1-exp(lLU));

%63 eeq8 Vxzplus : average value of finding a job 
exp(lVxzplusU)-1/exp(lmU(-1))*(
 exp(lchi0U(-1))*(1-exp(lfcall0U(-1)))*exp(ll0U(-1))*(exp(lfcall1U)*exp(lUzplusU)  + (1-exp(lfcall1U))*exp(lVzplus1U))
+exp(lchi1U(-1))*(1-exp(lfcall1U(-1)))*exp(ll1U(-1))*(exp(lfcall2U)*exp(lUzplusU)  + (1-exp(lfcall2U))*exp(lVzplus2U))   
+exp(lchi2U(-1))*(1-exp(lfcall2U(-1)))*exp(ll2U(-1))*(exp(lfcall3U)*exp(lUzplusU)  + (1-exp(lfcall3U))*exp(lVzplus3U))   
+exp(lchi3U(-1))*(1-exp(lfcall3U(-1)))*exp(ll3U(-1))*(exp(lfcall0U)*exp(lUzplusU)  + (1-exp(lfcall0U))*exp(lVzplus0U)) );

%64 eeq21 BigQ Probability of filling a vacancy
exp(lBigQU)=exp(lmU)/exp(lvU);

%65 eeq19  m total new matches      
exp(lmU)= exp(lchi0U)*(1-exp(lfcall0U))*exp(ll0U)
+exp(lchi1U)*(1-exp(lfcall1U))*exp(ll1U)
+exp(lchi2U)*(1-exp(lfcall2U))*exp(ll2U)
+exp(lchi3U)*(1-exp(lfcall3U))*exp(ll3U);


%66 eeq10 Matching function
exp(lmU)=(exp(lsigmamU)*((1-exp(lLU))^sigma_match)*(exp(lvU)^(1-sigma_match)));

%67-70 eeq15 leisure/intensity decision, per cohort
lintensity0U = ( lpsizplusU + lwbarU + log((1-exp(ltauyU))/(1+tauw)) + lepscall0U - log(1-exp(lfcall0U)) - lzetahU - log(AL) )/sigmaL;
lintensity1U = ( lpsizplusU + lwbarU + log((1-exp(ltauyU))/(1+tauw)) + lepscall1U - log(1-exp(lfcall1U)) - lzetahU - log(AL) )/sigmaL;
lintensity2U = ( lpsizplusU + lwbarU + log((1-exp(ltauyU))/(1+tauw)) + lepscall2U - log(1-exp(lfcall2U)) - lzetahU - log(AL) )/sigmaL;
lintensity3U = ( lpsizplusU + lwbarU + log((1-exp(ltauyU))/(1+tauw)) + lepscall3U - log(1-exp(lfcall3U)) - lzetahU - log(AL) )/sigmaL;

% 71-74 eeq16 L.O.M. for employees in each cohort 
exp(ll1U)-(exp(lchi0U(-1))+rho)*(1-exp(lfcall0U(-1)))*exp(ll0U(-1));
exp(ll2U)-(exp(lchi1U(-1))+rho)*(1-exp(lfcall1U(-1)))*exp(ll1U(-1));
exp(ll3U)-(exp(lchi2U(-1))+rho)*(1-exp(lfcall2U(-1)))*exp(ll2U(-1));
exp(ll0U)-(exp(lchi3U(-1))+rho)*(1-exp(lfcall3U(-1)))*exp(ll3U(-1));

% 75 eeq5  marginal value of a job for the worker wrt. wage
exp(lVwU)=exp(lmcall0U)*exp(lintensity0U)*((1-exp(ltauyU))/(1+tauw))
     +(rho*betaa)*exp(lmcall1U(+1))*(exp(lpsizplusU(+1))/exp(lpsizplusU)*exp(lG1U(+1))*exp(lintensity1U(+1))*((1-exp(ltauyU(+1)))/(1+tauw)))
     +(rho*betaa)^2*exp(lmcall2U(+2))*(exp(lpsizplusU(+2))/exp(lpsizplusU)*exp(lG2U(+2))*exp(lintensity2U(+2))*((1-exp(ltauyU(+2)))/(1+tauw)))
     +(rho*betaa)^3*exp(lmcall3U(+3))*(exp(lpsizplusU(+3))/exp(lpsizplusU)*exp(lG3U(+3))*exp(lintensity3U(+3))*((1-exp(ltauyU(+3)))/(1+tauw)));

% 76 eeq9 Nash bargaining foc      
exp(letaU)*exp(lVwU)*exp(lJzplusU)+(1-exp(letaU))*(exp(lVzplus0U)-exp(lUzplusU))*JwU;

% 77 eeq3 Jw marginal value of job match for a firm for a change in the wage 
%uncomment the following in case you want union bargaining. In that case, you must comment out the following batch of code
% /*
% JwU=((-1)*(exp(lintensity0U)*(1-exp(lfcall0U))
%      +betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU)*exp(lG1U(+1))*exp(lOm1U(+1))*exp(lintensity1U(+1)))
%      +betaa^2*(exp(lpsizplusU(+2))/exp(lpsizplusU)*exp(lG2U(+2))*exp(lOm2U(+2))*exp(lintensity2U(+2))) 
%      +betaa^3*(exp(lpsizplusU(+3))/exp(lpsizplusU)*exp(lG3U(+3))*exp(lOm3U(+3))*exp(lintensity3U(+3)))));
% */

% Below cdt holds for atomistic bargaining (connected to steady state file!)
JwU=((-1)*(exp(lintensity0U)*(1-exp(lfcall0U))
       +betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU)*exp(lG1U(+1))*(1-exp(lfcall1U(+1)))*rho*(1-exp(lfcall0U))*exp(lintensity1U(+1)))
     +betaa^2*(exp(lpsizplusU(+2))/exp(lpsizplusU)*exp(lG2U(+2))*(1-exp(lfcall2U(+2)))*rho*(1-exp(lfcall1U(+1)))*rho*(1-exp(lfcall0U))*exp(lintensity2U(+2))) 
     +betaa^3*(exp(lpsizplusU(+3))/exp(lpsizplusU)*exp(lG3U(+3))*(1-exp(lfcall3U(+3)))*rho*(1-exp(lfcall2U(+2)))*rho*(1-exp(lfcall1U(+1)))*rho*(1-exp(lfcall0U))*exp(lintensity3U(+3)))));

% 78 eeq11 Total EFFECTIVE hours (NOT to be used in measurement equation)
exp(lHU)=exp(lintensity0U)*exp(lepscall0U)*exp(ll0U)
        +exp(lintensity1U)*exp(lepscall1U)*exp(ll1U)
        +exp(lintensity2U)*exp(lepscall2U)*exp(ll2U)
        +exp(lintensity3U)*exp(lepscall3U)*exp(ll3U);


% *** Below large block of equations controls endogenous separations:
% 79-82 M-caligraphic product-updater that summarize probability of future endogen. separation
exp(lmcall0U)=(1-exp(lfcall0U));
exp(lmcall1U)=(1-exp(lfcall0U(-1)))*(1-exp(lfcall1U));
exp(lmcall2U)=(1-exp(lfcall0U(-2)))*(1-exp(lfcall1U(-1)))*(1-exp(lfcall2U));
exp(lmcall3U)=(1-exp(lfcall0U(-3)))*(1-exp(lfcall1U(-2)))*(1-exp(lfcall2U(-1)))*(1-exp(lfcall3U));

% 83-86 Definition of lfcall0U, probability of endog. separation., j=0,...,N-1
exp(lfcall0U)=normcdf((la0U + 1/2*exp(lsigma_aU)^2)/exp(lsigma_aU),0,1);
exp(lfcall1U)=normcdf((la1U + 1/2*exp(lsigma_aU)^2)/exp(lsigma_aU),0,1);
exp(lfcall2U)=normcdf((la2U + 1/2*exp(lsigma_aU)^2)/exp(lsigma_aU),0,1);
exp(lfcall3U)=normcdf((la3U + 1/2*exp(lsigma_aU)^2)/exp(lsigma_aU),0,1);

% 87 - 90 Definition of lepscall0U, density of idio. productivity, j=0,...,N-1
exp(lepscall0U)=1-normcdf((la0U + 1/2*exp(lsigma_aU)^2)/exp(lsigma_aU) - exp(lsigma_aU),0,1);
exp(lepscall1U)=1-normcdf((la1U + 1/2*exp(lsigma_aU)^2)/exp(lsigma_aU) - exp(lsigma_aU),0,1);
exp(lepscall2U)=1-normcdf((la2U + 1/2*exp(lsigma_aU)^2)/exp(lsigma_aU) - exp(lsigma_aU),0,1);
exp(lepscall3U)=1-normcdf((la3U + 1/2*exp(lsigma_aU)^2)/exp(lsigma_aU) - exp(lsigma_aU),0,1);

% **** Begin endogenous separation cut-off equations
% Recall, INDICATOR variables SW & SE for worker and employer surplus respectively

% 91 eeq90, cut-off for cohort 0:
SW*(exp(lwU)*exp(lwbarU)*(1-exp(ltauyU))/(1+tauw)-AL*exp(lzetahU)*exp(lintensity0U)^sigmaL/exp(lpsizplusU))*
1/sigmaL*exp(lintensity0U)^(1-sigmaL)*exp(lwU)*exp(lwbarU)*exp(lpsizplusU)*(1-exp(ltauyU))/(1+tauw)*exp(lgcallprime0U)/(exp(lzetahU)*AL)+  % this line only "intensity0-prime"
SE*exp(lwbarU)*((exp(lepscall0U)/(1-exp(lfcall0U))-exp(lwU))*
1/sigmaL*exp(lintensity0U)^(1-sigmaL)*exp(lwU)*exp(lwbarU)*exp(lpsizplusU)*(1-exp(ltauyU))/(1+tauw)*exp(lgcallprime0U)/(exp(lzetahU)*AL)+  % this line only "intensity0-prime"
exp(lgcallprime0U)*exp(lintensity0U) ) =
( SW*(exp(lVzplus0U) - exp(lUzplusU)) + SE*exp(lJtildezplus0U) )*exp(lfcallprime0U) / (1-exp(lfcall0U));

% 92-94 eeq91, cut-off for cohort j=1,2,3 (only the cohort-digit changes, nothing else, no discounting differences)
% j=1
SW*(exp(lG1U)*exp(lwU(-1))*exp(lwbarU(-1))*(1-exp(ltauyU))/(1+tauw)-AL*exp(lzetahU)*exp(lintensity1U)^sigmaL/exp(lpsizplusU))*
1/sigmaL*exp(lintensity1U)^(1-sigmaL)*exp(lwU)*exp(lwbarU)*exp(lpsizplusU)*(1-exp(ltauyU))/(1+tauw)*exp(lgcallprime1U)/(exp(lzetahU)*AL)+  % this line only "intensity1-prime"
SE*( (exp(lwbarU)*exp(lepscall1U)/(1-exp(lfcall1U))-exp(lG1U)*exp(lwU(-1))*exp(lwbarU(-1)))*
1/sigmaL*exp(lintensity1U)^(1-sigmaL)*exp(lwU)*exp(lwbarU)*exp(lpsizplusU)*(1-exp(ltauyU))/(1+tauw)*exp(lgcallprime1U)/(exp(lzetahU)*AL)+  % this line only "intensity1-prime"
exp(lwbarU)*exp(lgcallprime1U)*exp(lintensity1U) ) =
( SW*(exp(lVzplus1U) - exp(lUzplusU)) + SE*exp(lJtildezplus1U) )*exp(lfcallprime1U) / (1-exp(lfcall1U));

% j=2
SW*(exp(lG2U)*exp(lwU(-2))*exp(lwbarU(-2))*(1-exp(ltauyU))/(1+tauw)-AL*exp(lzetahU)*exp(lintensity2U)^sigmaL/exp(lpsizplusU))*
1/sigmaL*exp(lintensity2U)^(1-sigmaL)*exp(lwU)*exp(lwbarU)*exp(lpsizplusU)*(1-exp(ltauyU))/(1+tauw)*exp(lgcallprime2U)/(exp(lzetahU)*AL)+  % this line only "intensity1-prime"
SE*( (exp(lwbarU)*exp(lepscall2U)/(1-exp(lfcall2U))-exp(lG2U)*exp(lwU(-2))*exp(lwbarU(-2)))*
1/sigmaL*exp(lintensity2U)^(1-sigmaL)*exp(lwU)*exp(lwbarU)*exp(lpsizplusU)*(1-exp(ltauyU))/(1+tauw)*exp(lgcallprime2U)/(exp(lzetahU)*AL)+  % this line only "intensity1-prime"
exp(lwbarU)*exp(lgcallprime2U)*exp(lintensity2U) ) =
( SW*(exp(lVzplus2U) - exp(lUzplusU)) + SE*exp(lJtildezplus2U) )*exp(lfcallprime2U) / (1-exp(lfcall2U));

% j=3
SW*(exp(lG3U)*exp(lwU(-3))*exp(lwbarU(-3))*(1-exp(ltauyU))/(1+tauw)-AL*exp(lzetahU)*exp(lintensity3U)^sigmaL/exp(lpsizplusU))*
1/sigmaL*exp(lintensity3U)^(1-sigmaL)*exp(lwU)*exp(lwbarU)*exp(lpsizplusU)*(1-exp(ltauyU))/(1+tauw)*exp(lgcallprime3U)/(exp(lzetahU)*AL)+  % this line only "intensity1-prime"
SE*( (exp(lwbarU)*exp(lepscall3U)/(1-exp(lfcall3U))-exp(lG3U)*exp(lwU(-3))*exp(lwbarU(-3)))*
1/sigmaL*exp(lintensity3U)^(1-sigmaL)*exp(lwU)*exp(lwbarU)*exp(lpsizplusU)*(1-exp(ltauyU))/(1+tauw)*exp(lgcallprime3U)/(exp(lzetahU)*AL)+  % this line only "intensity1-prime"
exp(lwbarU)*exp(lgcallprime3U)*exp(lintensity3U) ) =
( SW*(exp(lVzplus3U) - exp(lUzplusU)) + SE*exp(lJtildezplus3U) )*exp(lfcallprime3U) / (1-exp(lfcall3U));


% 95-98 eeq92, definition of Jtildezplus, employers value of an employee AFTER endog. separation has taken place
% j=0
exp(lJtildezplus0U)=(exp(lwbarU)*exp(lepscall0U)/(1-exp(lfcall0U))-exp(lwU)*exp(lwbarU))*exp(lintensity0U)-(kappa/2)*(exp(lvtilde0U)^2)  
+betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU))*(exp(lchi0U)+rho)*(
((exp(lwbarU(+1))*exp(lepscall1U(+1))/(1-exp(lfcall1U(+1)))-exp(lG1U(+1))*exp(lwU)*exp(lwbarU))*exp(lintensity1U(+1))-(kappa/2)*(exp(lvtilde1U(+1))^2))*(1-exp(lfcall1U(+1)))
+betaa*(exp(lpsizplusU(+2))/exp(lpsizplusU(+1)))*
((exp(lwbarU(+2))*exp(lepscall2U(+2))/(1-exp(lfcall2U(+2)))-exp(lG2U(+2))*exp(lwU)*exp(lwbarU))*exp(lintensity2U(+2))-(kappa/2)*(exp(lvtilde2U(+2))^2))*(1-exp(lfcall2U(+2)))*(exp(lchi1U(+1))+rho)*(1-exp(lfcall1U(+1)))
+betaa^2*(exp(lpsizplusU(+3))/exp(lpsizplusU(+1)))*
((exp(lwbarU(+3))*exp(lepscall3U(+3))/(1-exp(lfcall3U(+3)))-exp(lG3U(+3))*exp(lwU)*exp(lwbarU))*exp(lintensity3U(+3))-(kappa/2)*(exp(lvtilde3U(+3))^2))*(1-exp(lfcall3U(+3)))*(exp(lchi2U(+2))+rho)*(1-exp(lfcall2U(+2)))*(exp(lchi1U(+1))+rho)*(1-exp(lfcall1U(+1)))
+betaa^3*(exp(lpsizplusU(+4))/exp(lpsizplusU(+1)))*
exp(lJzplusU(+4))*(exp(lchi3U(+3))+rho)*(1-exp(lfcall3U(+3)))*(exp(lchi2U(+2))+rho)*(1-exp(lfcall2U(+2)))*(exp(lchi1U(+1))+rho)*(1-exp(lfcall1U(+1))) ); 

% 96 for j=1 
exp(lJtildezplus1U)=(exp(lwbarU)*exp(lepscall1U)/(1-exp(lfcall1U))-exp(lG1U)*exp(lwU(-1))*exp(lwbarU(-1)))*exp(lintensity1U)-(kappa/2)*(exp(lvtilde1U)^2)  
+betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU))*(exp(lchi1U)+rho)*(
((exp(lwbarU(+1))*exp(lepscall2U(+1))/(1-exp(lfcall2U(+1)))-exp(lG2U(+1))*exp(lwU(-1))*exp(lwbarU(-1)))*exp(lintensity2U(+1))-(kappa/2)*(exp(lvtilde2U(+1))^2))*(1-exp(lfcall2U(+1)))
+betaa*(exp(lpsizplusU(+2))/exp(lpsizplusU(+1)))*
((exp(lwbarU(+2))*exp(lepscall3U(+2))/(1-exp(lfcall3U(+2)))-exp(lG3U(+2))*exp(lwU(-1))*exp(lwbarU(-1)))*exp(lintensity3U(+2))-(kappa/2)*(exp(lvtilde3U(+2))^2))*(1-exp(lfcall3U(+2)))*(exp(lchi2U(+1))+rho)*(1-exp(lfcall2U(+1)))
+betaa^2*(exp(lpsizplusU(+3))/exp(lpsizplusU(+1)))*
exp(lJzplusU(+3))*(exp(lchi3U(+2))+rho)*(1-exp(lfcall3U(+2)))*(exp(lchi2U(+1))+rho)*(1-exp(lfcall2U(+1))) ); 

% 97 for j=2:
exp(lJtildezplus2U)=(exp(lwbarU)*exp(lepscall2U)/(1-exp(lfcall2U))-exp(lG2U)*exp(lwU(-2))*exp(lwbarU(-2)))*exp(lintensity2U)-(kappa/2)*(exp(lvtilde2U)^2)  
+betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU))*(exp(lchi2U)+rho)*(
((exp(lwbarU(+1))*exp(lepscall3U(+1))/(1-exp(lfcall3U(+1)))-exp(lG3U(+1))*exp(lwU(-2))*exp(lwbarU(-2)))*exp(lintensity3U(+1))-(kappa/2)*(exp(lvtilde3U(+1))^2))*(1-exp(lfcall3U(+1)))
+betaa*(exp(lpsizplusU(+2))/exp(lpsizplusU(+1)))*
exp(lJzplusU(+2))*(exp(lchi3U(+1))+rho)*(1-exp(lfcall3U(+1))) ); 

% 98 for j=3:
exp(lJtildezplus3U)=(exp(lwbarU)*exp(lepscall3U)/(1-exp(lfcall3U))-exp(lG3U)*exp(lwU(-3))*exp(lwbarU(-3)))*exp(lintensity3U)-(kappa/2)*(exp(lvtilde3U)^2)  
+betaa*(exp(lpsizplusU(+1))/exp(lpsizplusU))*(exp(lchi3U)+rho)*(
exp(lJzplusU(+1)) );     % this line implicitly contains betaa^0*exp(lpsizplusU(+1))/exp(lpsizplusU(+1))
                         % Expositional problem: counter for chi and fcall in eeq92 goes wrong way: supposed to fall from N-1 to j+1. 
                         % Here that means rising from 3 to 4. Interpreted as empty set.

% 99-102 eq. "diffG_full", j=0,...,N-1
% j=0:
exp(lgcallprime0U)=(-exp(la0U)*exp(lfcallprime0U)*(1-exp(lfcall0U))+exp(lepscall0U)*exp(lfcallprime0U)) / (1-exp(lfcall0U))^2;

% j=1:
exp(lgcallprime1U)=(-exp(la1U)*exp(lfcallprime1U)*(1-exp(lfcall1U))+exp(lepscall1U)*exp(lfcallprime1U)) / (1-exp(lfcall1U))^2;

% j=2:
exp(lgcallprime2U)=(-exp(la2U)*exp(lfcallprime2U)*(1-exp(lfcall2U))+exp(lepscall2U)*exp(lfcallprime2U)) / (1-exp(lfcall2U))^2;

% j=3:
exp(lgcallprime3U)=(-exp(la3U)*exp(lfcallprime3U)*(1-exp(lfcall3U))+exp(lepscall3U)*exp(lfcallprime3U)) / (1-exp(lfcall3U))^2;


% 103-106 eq. "diffF_full" (just an implementation of a normal PDF with variance exp(lsigma_aU)^2 and mean -exp(lsigma_aU)^2/2 evaluated at exp(lsigma_aU))
% j=0
exp(lfcallprime0U)=1/(exp(la0U)*exp(lsigma_aU)*sqrt(2*3.141592653589793))*exp(-(la0U+0.5*exp(lsigma_aU)^2)^2/(2*exp(lsigma_aU)^2));

% j=1
exp(lfcallprime1U)=1/(exp(la1U)*exp(lsigma_aU)*sqrt(2*3.141592653589793))*exp(-(la1U+0.5*exp(lsigma_aU)^2)^2/(2*exp(lsigma_aU)^2));

% j=2
exp(lfcallprime2U)=1/(exp(la2U)*exp(lsigma_aU)*sqrt(2*3.141592653589793))*exp(-(la2U+0.5*exp(lsigma_aU)^2)^2/(2*exp(lsigma_aU)^2));

% j=3
exp(lfcallprime3U)=1/(exp(la3U)*exp(lsigma_aU)*sqrt(2*3.141592653589793))*exp(-(la3U+0.5*exp(lsigma_aU)^2)^2/(2*exp(lsigma_aU)^2));

% *** end endogenous separation cut-off equations block ***


% /* 
% %the following three sets of code should be uncommented in case you want fixed separations. In that case, comment out the previous three sets of code
% lfcall0U=log(fcall0);
% lfcall1U=log(fcall1);
% lfcall2U=log(fcall2);
% lfcall3U=log(fcall3);
% 
% lepscall0U=log(epscall0);
% lepscall1U=log(epscall1);
% lepscall2U=log(epscall2);
% lepscall3U=log(epscall3);
% 
% la0U=log(a0);
% la1U=log(a1);
% la2U=log(a2);
% la3U=log(a3);
% */

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% *** end of labor block *************%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%113 eq16 - Taylor Rule
% comment this if you want to calculate Ramsey optimal policy
% lRU-log(R)=rhoR*(lRU(-1)-log(R))+(1-rhoR)*(lpitargetU-log(pitarget)+
%   rpi*(lpicU -lpitargetU)+
%   ry*(lylessdU-log(ylessd))+rRstar*(lRstarU - lRstarU(-1)) ) + epsR_eps/10; %gb:originally epsR_eps/100


%114 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))) 
+(omegax*exp(lpmxU)^(1-etax)+1-omegax)^(etax/(1-etax))*     
(1-omegax)*exp(lphaloxU*lambdax/(1-lambdax))*exp(lpxU)^(-etaf)*exp(lystarU)
+((kappa/2)*(exp(lvtilde0U)^2*exp(ll0U)*(1-exp(lfcall0U))+exp(lvtilde1U)^2*exp(ll1U)*(1-exp(lfcall1U))+exp(lvtilde2U)^2*exp(ll2U)*(1-exp(lfcall2U)) 
+ exp(lvtilde3U)^2*exp(ll3U)*(1-exp(lfcall3U))))
+mu*normcdf((lomegabarU +exp(lsigmaU(-1))^2/2)/exp(lsigmaU(-1)) -  exp(lsigmaU(-1)),0,1)
*exp(lRkU+lpkprimeU(-1)+lkbarU(-1)-lpiU-lmuzplusU);


%115 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); 

%116 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));

%117 eq20a - marginal costs of consumption importers
exp(lmcmcU)=exp(ltaumcU)*exp(lqU)*exp(lpcU)*exp(lRnustarU)/exp(lpmcU);

%118 eq20b - marginal costs of investment importers
exp(lmcmiU)=exp(ltaumiU)*exp(lqU)*exp(lpcU)*exp(lRnustarU)/exp(lpmiU);

%119 eq20c - marginal costs of export importers
exp(lmcmxU)=exp(ltaumxU)*exp(lqU)*exp(lpcU)*exp(lRnustarU)/exp(lpmxU);

%120 eq21 - Real GDP, defined from the production side
exp(lyU)=exp(lphaloU)^(lambdad/(lambdad-1))*
        % *** commented out: ((exp(lsmallhU-(lambdaw/(1-lambdaw))*lwhaloU))^(1-alphaa)*exp(lepsilonU)* 
        ((exp(lHU))^(1-alphaa)*exp(lepsilonU)*
        (exp(lkU)/(exp(lmuzplusU)*exp(lmupsiU)))^(alphaa)-phi);

%121 eq41 Capacity utilization
aofuU = 0.5*sigmab*sigmaa*(exp(luU))^2+sigmab*(1-sigmaa)*exp(luU)
        +sigmab*((sigmaa/2)-1);

%122 eq22 - Definition of Rf
exp(lRfU)=nuf*(exp(lRU))+1-nuf;

%123 eq23 - Definition of Rnustar
exp(lRnustarU)=nustar*(exp(lRstarU))+1-nustar;

%124 eq24 - Total foreign export demand
exp(lxU)=(exp(lpxU))^(-etaf)*exp(lystarU);

%125 eq25 - Relative price of final consumption good
exp(lpcU)=(1-omegac+omegac*(exp(lpmcU))^(1-etac))^(1/(1-etac));

%126 eq26 - Relative price of final investment good
exp(lpinvestU)=(1-omegai+omegai*(exp(lpmiU))^(1-etai))^(1/(1-etai));

%127 eq27 - Definition of Rx
exp(lRxU)=nux*(exp(lRU))+1-nux;

%128 eq28 - Definition of capital utilitzation
exp(lkU)=exp(lkbarU(-1)+luU);

%129 eq29 - Dynamics of pmx
exp(lpmxU)/exp(lpmxU(-1))=exp(lpimxU)/exp(lpiU);

%130 eq30 - Dynamics of pmc
exp(lpmcU)/exp(lpmcU(-1))=exp(lpimcU)/exp(lpiU);

%131 eq31 - Dynamics of pmi
exp(lpmiU)/exp(lpmiU(-1))=exp(lpimiU)/exp(lpiU);

%132 eq32 - Dynamics of px
exp(lpxU)/exp(lpxU(-1))=exp(lpixU)/exp(lpistarU);

%133 eq33 - Dynamics of the real exchange rate
exp(lqU)/exp(lqU(-1))=exp(lsU_eps)*exp(lpistarU)/exp(lpicU);

%134 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);

%135 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)));

%136 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));

%137 eq 45 Relation between smallh and  H, commented out b/c only relevant for s.s wage dispersion.
%exp(lsmallhU)=exp(lwhaloU)^(lambdaw/(1-lambdaw))*exp(lHU);

%138 eq 46 Imported consumtion
exp(lcmU)=exp(lcU)*omegac*(exp(lpcU)/exp(lpmcU))^etac;

%139 eq 47 Imported investment
exp(limU)=omegai*(exp(lpinvestU)/exp(lpmiU))^etai
*(exp(liU)+aofuU*exp(lkbarU(-1))/(exp(lmupsiU)*exp(lmuzplusU)));

%140 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);

%141 eq15f
lpiwU=lwbarU+lmuzplusU+lpiU-lwbarU(-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
% equations for the law of motions of exogenous processes %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/

%142 Composite technology growth
exp(lmuzplusU)=(exp(lmupsiU))^(alphaa/(1-alphaa))*exp(lmuzU);

%143-165 
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

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;
lsigma_aU=(1-rhosigma_a)*log(sigma_a) + rhosigma_a*lsigma_aU(-1)+sigma_a_eps/100;

%Measurement equations 166-189
%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*( log (exp(lintensity0U)*(1-exp(lfcall0U))*exp(ll0U)
                                 +exp(lintensity1U)*(1-exp(lfcall1U))*exp(ll1U)
                                 +exp(lintensity2U)*(1-exp(lfcall2U))*exp(ll2U)
                                 +exp(lintensity3U)*(1-exp(lfcall3U))*exp(ll3U)) 
                               -log( exp(lintensity0U(-1))*(1-exp(lfcall0U(-1)))*exp(ll0U(-1))
                                 +exp(lintensity1U(-1))*(1-exp(lfcall1U(-1)))*exp(ll1U(-1))
                                 +exp(lintensity2U(-1))*(1-exp(lfcall2U(-1)))*exp(ll2U(-1))
                                 +exp(lintensity3U(-1))*(1-exp(lfcall3U(-1)))*exp(ll3U(-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 measured 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)) ) );

% Investment w/o capacity util in meas. eq.
data_idiffU-data_idiffU_eps+100*(log(muzplus)+log(mupsi))=100*(lmuzplusU+lmupsiU+ liU - liU(-1));

% We match the weighted AVERAGE wage 
data_wdiffU-data_wdiffU_eps+100*log(muzplus)=100*(lmuzplusU+(lavg_NashU - lavg_NashU(-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_unempdiffU-data_unempdiffU_eps=100*(log((1-exp(lLU)))-log((1-exp(lLU(-1)))));
%data_spreaddiffU-data_spreaddiffU_eps=100*(lspreadU-lspreadU(-1));
data_spreadU-data_spreadU_eps=400*exp(lspreadU);
data_gdiffU-data_gdiffU_eps+100*log(muzplus)=100*(lmuzplusU+(lgU - lgU(-1)));

data_unemphatU-data_unemphatU_eps=100*(1-exp(lLU) - (1-L) ) / (1-L);
% Note ss of Hmeas=intensity*(1-Fcall)*4*l0=intensity*L (or a sum of this over N cohorts if ss wage dispersion)
data_HhatU-data_HhatU_eps=100*( (exp(lintensity0U)*(1-exp(lfcall0U))*exp(ll0U)
             +exp(lintensity1U)*(1-exp(lfcall1U))*exp(ll1U)
             +exp(lintensity2U)*(1-exp(lfcall2U))*exp(ll2U)
             +exp(lintensity3U)*(1-exp(lfcall3U))*exp(ll3U)) - intensity0*L ) / (intensity0*L);

% *** some further variables of interest ***
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);

exp(lavg_NashU)=1/exp(lLU)*(exp(ll0U)*exp(lwU)*exp(lwbarU)*(1-exp(lfcall0U)) +
                            exp(ll1U)*exp(lG1U)*exp(lwU(-1))*exp(lwbarU(-1))*(1-exp(lfcall1U))  +
                            exp(ll2U)*exp(lG2U)*exp(lwU(-2))*exp(lwbarU(-2))*(1-exp(lfcall2U))  +
                            exp(ll3U)*exp(lG3U)*exp(lwU(-3))*exp(lwbarU(-3))*(1-exp(lfcall3U)) );
% 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) )
-((kappa/2)*(exp(lvtilde0U)^2*exp(ll0U)*(1-exp(lfcall0U))+exp(lvtilde1U)^2*exp(ll1U)*(1-exp(lfcall1U))+exp(lvtilde2U)^2*exp(ll2U)*(1-exp(lfcall2U)) 
+ exp(lvtilde3U)^2*exp(ll3U)*(1-exp(lfcall3U))))
-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);
% credit amount in levels
exp(lbU)=exp(lpkprimeU+lkbarU)-exp(lnU);
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
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_spreaddiffU=0;
data_spreadU=400*spread;
data_gdiffU=0;

data_unemphatU=0;
data_HhatU=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);
aofuU=aofu;
ltaudU=log(taud);
ltauxU=log(taux);
ltaumcU=log(taumc);
ltaumiU=log(taumi);
ltaumxU=log(taumx);
lpiwU=log(piw);
  
lJzplusU=log(Jzplus);    
lwU=log(w);   
lUzplusU=log(Uzplus); 
lVxzplusU=log(Vxzplus);   
lvU=log(v);  
lG1U=log(G1); 
lG2U=log(G2); 
lG3U=log(G3);
lOm1U=log(Om1); 
lOm2U=log(Om2); 
lOm3U=log(Om3); 
lOm4U=log(Om4);
lvtilde0U=log(vtilde0); 
lvtilde1U=log(vtilde1); 
lvtilde2U=log(vtilde2); 
lvtilde3U=log(vtilde3);
lVzplus0U=log(Vzplus0); 
lVzplus1U=log(Vzplus1); 
lVzplus2U=log(Vzplus2); 
lVzplus3U=log(Vzplus3); 
ll0U=log(l0); 
ll1U=log(l1); 
ll2U=log(l2); 
ll3U=log(l3);
lLU=log(L); 
lmU=log(m);
lBigQU=log(BigQ);
lpitildewU=log(pitildew);
lchi0U=log(chi0); 
lchi1U=log(chi1); 
lchi2U=log(chi2); 
lchi3U=log(chi3);
lfU=log(f);
JwU=Jw;
lVwU=log(Vw);
lsigmamU=log(sigmam); 
letaU=log(eta);

lspreadU=log(spread);
lUnemprateU=log(1-L);
lbankruptcyrateU=log(Fomegabar);
lavg_NashU=log(avg_Nash);
lylessdU=log(ylessd);

lsigma_aU=log(sigma_a);

lepscall0U=log(epscall0);
lepscall1U=log(epscall1); 
lepscall2U=log(epscall2);
lepscall3U=log(epscall3);
lfcall0U=log(fcall0); 
lfcall1U=log(fcall1); 
lfcall2U=log(fcall2); 
lfcall3U=log(fcall3);
lintensity0U=log(intensity0);
lintensity1U=log(intensity1); 
lintensity2U=log(intensity2); 
lintensity3U=log(intensity3);

lmcall0U=log(mcall0 );
lmcall1U=log(mcall1); 
lmcall2U=log(mcall2 );
lmcall3U=log(mcall3);
la0U=log(a0);
la1U=log(a1);
la2U=log(a2);
la3U=log(a3);

lJtildezplus0U=log(Jtildezplus0); 
lJtildezplus1U=log(Jtildezplus1); 
lJtildezplus2U=log(Jtildezplus2); 
lJtildezplus3U=log(Jtildezplus3); 

lbU=log(pkprime*kbar-n);
etaaU=a/ylessd;
end;

resid(1);

steady;

check; 


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;  % !GB2013 original finfric model has zero!
                            % but not my final finfric model        
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

% 10 percent of variance of time series 1995q1-2010q3
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; % 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= .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; %0.286;
var data_unempdiffU_eps= 10.34; %7.3042; %	1.777; % originally suspended; consider opening
%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; 


%shocks
% *** all shock priors re-scaled by some factor; see the above equations
% init-val, prior shape, mean, std
stderr muz_eps,    sig_muz,  inv_gamma_pdf,0.15,inf;  % .4
% *** stderr mupsi_eps,    inv_gamma_pdf,0.15,inf; ***
stderr epsilon_eps, .131,inv_gamma_pdf,0.15,inf; % rescaled
stderr Upsilon_eps,  .163,inv_gamma_pdf,0.15,inf;  
stderr zetac_eps,    .189,inv_gamma_pdf,0.15,inf;    
stderr zetah_eps,    .262,inv_gamma_pdf,0.15,inf;        
stderr phitilde_eps, .544,inv_gamma_pdf,0.5,inf;     
% stderr epsR_eps,     .24,inv_gamma_pdf,0.15,inf;
stderr g_eps,         .463,  inv_gamma_pdf,0.5,inf;
stderr taud_eps,     .397, inv_gamma_pdf,0.5,inf; 
stderr taux_eps,     1.17,inv_gamma_pdf,0.5,inf; 
stderr taumc_eps,    .66,inv_gamma_pdf,0.5,inf; 
stderr taumi_eps,    .39,inv_gamma_pdf,0.5,inf; 
stderr taumx_eps,   1.55, inv_gamma_pdf,0.5,inf; 
stderr gamma_eps,    .253,inv_gamma_pdf,0.5,inf; 
% stderr sigmam_eps,  .0461,  inv_gamma_pdf,0.1,inf; % gb2014, labor block 
% stderr sigma_eps,    inv_gamma_pdf,0.5,inf; 
% stderr eta_eps,         inv_gamma_pdf,0.5,inf;
stderr ystar_eps, sig_ystar,inv_gamma_pdf,0.5,inf;  % .25
stderr pistar_eps, sig_pistar,inv_gamma_pdf,0.5,inf;  % .26
stderr Rstar_eps, sig_Rstar,inv_gamma_pdf,0.15,inf; % .15% scaled differently, by 1000
                                          % gb: originally 1.5  
% % measurement error std 
% stderr data_ydiffU_eps, 1.25,inv_gamma_pdf,0.5354^.5,inf; 
% stderr data_pidU_eps, 	6.49,inv_gamma_pdf,7.0402^.5,inf; 
% stderr data_wdiffU_eps, 2.02,inv_gamma_pdf,0.5534^.5,inf; 
% stderr data_cdiffU_eps, 1.21,inv_gamma_pdf,0.8045^.5,inf;
% stderr data_idiffU_eps, 13.62,inv_gamma_pdf,26.6239^.5,inf; 
% stderr data_qdiffU_eps, 1.84,inv_gamma_pdf,0.6311^.5,inf;
% stderr data_HdiffU_eps, 1.59,inv_gamma_pdf,0.4819^.5,inf; 
% stderr data_xdiffU_eps, .619,inv_gamma_pdf,1.1653^.5,inf;
% stderr data_impdiffU_eps,4.084,inv_gamma_pdf,3.9705^.5,inf; 
% stderr data_picU_eps, 	.975,inv_gamma_pdf,3.9595^.5,inf; 
% stderr data_piiU_eps, 	49.75,inv_gamma_pdf,264.7161^.5,inf;
% stderr data_ndiffU_eps, 8.78,inv_gamma_pdf,11.1925^.5,inf; 
% stderr data_spreadU_eps, 1.45,inv_gamma_pdf,.50664^.5,inf;  
% stderr data_gdiffU_eps, .774,inv_gamma_pdf,2.9773^.5,inf;
% stderr data_unempdiffU_eps, 7.38,inv_gamma_pdf,10.34^.5,inf; 
                                          
                                                                                    
% *** stickiness ****
xid,     .808,beta_pdf,0.75,0.075;
xix,     .877,beta_pdf,0.75,0.075;
ximc,    .82,beta_pdf,0.75,0.075;
ximi,    .38,beta_pdf,0.75,0.075;
ximx,    .598,beta_pdf,0.66,0.10; 

kappad,     .29,beta_pdf,0.5,0.15;
kappax,     .372,beta_pdf,0.5,0.15;
kappamc,    .476,beta_pdf,0.5,0.15;
kappami,    .278,beta_pdf,0.5,0.15;
kappamx,    .307,beta_pdf,0.5,0.15;

kappaw,     .203,beta_pdf,0.5,0.15;

work_cap_para, .5,beta_pdf,0.5,0.25; % .25

% preference and technology 
sigmaLScaled,.906,gamma_pdf,  .75,   0.2; % .75, .2
b,           .881, beta_pdf,   0.65,   0.15;
SppScaled,   .173,gamma_pdf,  0.5,    0.15;
sigmaa,   .36,gamma_pdf,  0.2,    0.075;  

% *** Taylor rule ***
% identical priors to SW 2003.
% rhoR,    .7,beta_pdf,  0.8,     0.1;    
% rpi,     1.7,normal_pdf,1.7,     0.15;   
% ry,      .02,normal_pdf,0.125,   0.05; 
% rRstar,      beta_pdf,.6,     0.1; 

% *** not estimated *** rdeltapi,,,,normal_pdf,0.3,     0.1;
% *** not estimated *** rdeltay,,,, gamma_pdf, 0.0625,   0.05;   

% *** open econ. params
etax,           1.67,gamma_pdf,1.5,0.25; % all etas trunc'ed below
etac,           1.18,gamma_pdf,1.5,0.25;
etai,           1.06,gamma_pdf,1.5,0.25;
etaf,           1.52,gamma_pdf,1.5,0.25;
%phitildes,      1.09,gamma_pdf,1.25,0.1;  

% *** financial frictions parameters ***
mu, .256,beta_pdf, 0.3,  0.075;   

% *** labor market parameters ***
recruitsharePercent, .41,gamma_pdf,     0.3, 0.15;   
bshare,           .804,beta_pdf,      0.75,0.075;  
BigFPercent,      .38,beta_pdf,      0.25,0.05; 

% *** Exog. AR1 processes ***
% rhomuz,      .6,beta_pdf,0.5, 0.075; 
% *** rhomupsi,    beta_pdf,0.5, 0.075; ***
rhoepsilon,  .866,beta_pdf,0.85,0.075;
rhoUpsilon,  .535,beta_pdf,0.85,0.075;
rhozetac,    .866,beta_pdf,0.85,0.075; 
rhozetah,    .964,beta_pdf,0.85,0.075; 
rhophitilde, .903,beta_pdf,0.85,0.075; 
rhog,        .771,beta_pdf,0.85,0.075;
% rhosigma,    beta_pdf,0.85,0.075; 
rhogamma,    .817,beta_pdf,0.85,0.075;
% rhoeta,         beta_pdf,0.85,0.075;
% rhosigmam,   .75,beta_pdf,0.85,0.075;

% *** Foreign VAR ***

% a11, .92,normal_pdf,0.5,0.5;  
% a22, .07,normal_pdf,0.0,0.5; 
% a33, .9, normal_pdf,0.5,0.5; 
% 
% a12,0,normal_pdf,0.0,0.5; 
% a13,0,normal_pdf,0.0,0.5;
% a21,.1,normal_pdf,0.0,0.5;
% a23,.07,normal_pdf,0.0,0.5;
% a24,0,normal_pdf,0.0,0.5;
% a31,.14,normal_pdf,0.0,0.5;
% a32,-.1,normal_pdf,0.0,0.5; 
% a34,.4,normal_pdf,0.0,0.5;
% 
% c21,-.4,normal_pdf,0.0,0.5;
% c31,0,normal_pdf,0.0,0.5;
% c32,0,normal_pdf,0.0,0.5;
% c24,.05,normal_pdf,0.0,0.5;
% c34,.1,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;
%a33, -.999, .999;
%ry, 0, 2;
end;


%%%%%%%%%%%%%%%%%%%%%%%%/
%%%%%%Observed variables%%%%%%%%%/
%%%%%%%%%%%%%%%%%%%%%%%%/
varobs  data_RU data_wdiffU data_cdiffU    
data_ydiffU  data_ystardiffU data_pistarU 
data_RstarU   data_HdiffU    %GB:  data_HhatU 
data_impdiffU data_piiU 
data_pidU data_picU data_xdiffU data_qdiffU data_idiffU
data_ndiffU data_spreadU data_unempdiffU data_gdiffU;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
%BEGINNING OF ESTIMATION%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
options_.console_mode=1; %(default: 0)

estimation(order=1,datafile=ctwDataRealTime,mode_file=ctw_employer_surplus100forecast1_mode, mh_replic=0, mh_nblocks=1, mh_drop=.8,  mh_jscale=0.2, 
 mode_compute=4, 
  filter_step_ahead = [1 2 3 4 5 6 7 8], forecast=40,
 smoother, filtered_vars, diffuse_filter,
nodisplay, 
 mode_check, 
% selected_variables_only, % doesn't work with laplace only
% cova_compute=0, %boosts speed if var-covar unnecessary
% do only for the final model
% moments_varendo,
% conditional_variance_decomposition=[1 4 8 20],
graph_format=(eps,pdf))
%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
%aofuU lpiwU

%variables for the unemployment part
%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

% 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 lsigma_aU

%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_ndiffU data_unempdiffU data_spreadU data_gdiffU
%lavg_NashU %lylessdU

% Endog. separation variables
%lepscall0U lepscall1U lepscall2U lepscall3U
%lfcall0U lfcall1U lfcall2U lfcall3U
%lintensity0U % lintensity1U lintensity2U lintensity3U
%lmcall0U lmcall1U lmcall2U lmcall3U     
%la0U la1U la2U la3U
%lgcallprime0U lgcallprime1U lgcallprime2U lgcallprime3U
%lfcallprime0U lfcallprime1U lfcallprime2U lfcallprime3U
%lJtildezplus0U lJtildezplus1U lJtildezplus2U lJtildezplus3U

% Additional observed variables
%data_HhatU data_unemphatU
%nxdivGDPU
%lSlevU
%lbU
%etaaU
%lsigmamU
; 


% 
% stoch_simul(order=1,irf=20, nodisplay, graph_format = (eps, pdf))
% %endogenous variables
% %lmcU lrkbarU 
%  %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
% %aofuU lpiwU
% 
% %variables for the unemployment part
% %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
% 
% % 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 lsigma_aU
% 
% %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_ndiffU data_unempdiffU data_spreadU data_gdiffU
% %lavg_NashU %lylessdU
% 
% % Endog. separation variables
% %lepscall0U lepscall1U lepscall2U lepscall3U
% %lfcall0U lfcall1U lfcall2U lfcall3U
% lintensity0U %lintensity1U lintensity2U lintensity3U
% %lmcall0U lmcall1U lmcall2U lmcall3U     
% %la0U la1U la2U la3U
% %lgcallprime0U lgcallprime1U lgcallprime2U lgcallprime3U
% %lfcallprime0U lfcallprime1U lfcallprime2U lfcallprime3U
% %lJtildezplus0U lJtildezplus1U lJtildezplus2U lJtildezplus3U
% 
% % Additional observed variables
% %data_HhatU data_unemphatU
% nxdivGDPU
% %lSlevU
% etaaU
% lwbarU; 

% shock_decomposition //(parameter_set=posterior_mean)
% data_ydiffU 
% %data_picU
% %data_spreadU 
% %data_unempdiffU
% ;



%ctw_oo_irfs=oo_.irfs;
 %save ctw_oo_irfs ctw_oo_irfs L

%model_diagnostics(M_,options_,oo_) %run it when some error for some further info

%load original vars
load('myOriginalData.mat');

%load undemeaned vars
load('myUndemeanedData.mat');

% mean(myUndemeanedData.dy_obs)
% take long run qoq 0.8% from the macromodel assumptions (Igors K.)
figure(1)
% plot(data_ydiffU) % the original demeaned var
hold on
plot(.8+oo_.SmoothedVariables.data_ydiffU,'blue') % smoothed var 
plot(.8+[oo_.SmoothedVariables.data_ydiffU;oo_.forecast.Mean.data_ydiffU],':red') % smoothed var + forecast
hold off

% print gdp + forecast
.8+[oo_.SmoothedVariables.data_ydiffU;oo_.forecast.Mean.data_ydiffU]

% plot filter step ahead
% data_ydiffU is about #136; no: 35
myparamcell=cellstr(M_.endo_names);
myIndex=find(strcmp(myparamcell,'data_ydiffU'));
filterstepahead=oo_.FilteredVariablesKStepAhead;
%for myjj=1:222
subfilterstetahead=squeeze(filterstepahead(:,35,:));  % myjj
subfilterstetahead(subfilterstetahead==0)=nan;
mybigformatrix=nan(size(subfilterstetahead,1)+size(oo_.SmoothedVariables.data_ydiffU,1));
for myi=1:size(oo_.SmoothedVariables.data_ydiffU,1)
    mybigformatrix(myi:myi+size(subfilterstetahead,1)-1,myi)=subfilterstetahead(:,myi);
end
figure(myjj+1)
hold on
plot(oo_.SmoothedVariables.data_ydiffU,'blue') % smoothed var
for myi=1:size(mybigformatrix,1)
plot(mybigformatrix(:,myi),':black') 
end
hold off
%end




/*
load('ctw_employer_surplus60_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','unempdiff', ...
'ndiff','spread'};

pid_ss=oo_.steady_state(129);
R_ss=oo_.steady_state(130);
Rstar_ss=oo_.steady_state(143); 
pic_ss=oo_.steady_state(139);
pii_ss=oo_.steady_state(140);
pistar_ss=oo_.steady_state(142);
spread_ss=oo_.steady_state(146);

mylevel = [0,pid_ss,R_ss,0,0,0,0,0,0,0,pic_ss,pii_ss,0,Rstar_ss,pistar_ss,0,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','varphi','v','sigma_a','psizplus','Uzplus','gammaa','omegabar','rkbar','a0','sigma'};%,'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
myParamVals

findWhatParams={'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
myParamVals

disp('K/Y ratio')
exp(oo_.steady_state(7)-oo_.steady_state(34))

*/
