% financial frictions in small open economy
% this Dynare code is for small open economy with financial frictions
% written on 2010/10/8, revised: 2010/10/30


// endogenous variables, 'lnx' means log(x)
var lnpc lnpmc lnq lnR lnpibar lnpic lnpx lnN lnpstar af 
    lnPhi lns lnx lnc lny lnKp lnFp  lnimp lnk lni
    lnPk lnmc lnrk lnw lnn lnRk omega ZR
// exogenous variables
lnRf lnyf phi lng lnA tau lnpif sigma gamm
// additional variables for easier computation
d Sadj Sadjprime Fprime Ft Gt Gammat;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% pc: scaled price of consumptiong good, Pc/P
% pmc: scaled price of imported good, Pmc/P
% q: real exchange rate, S*Pf/Pc
% R: risk free nominal rate of interest
% pibar: price inflation
% pic: consumer price inflation
% px: terms of trade, Px/Pf
% N: labor supply
% pstar: 'efficiency distortion'
% af: beginning-of-period t net stock of scaled real, domestic value of foreign assets, S(t)Af(t+1)/P(t)A(t)
% Phi: premium on foreign asset returns, over foreign risk free rate, Rf
% s: s(t) = S(t)/S(t-1)
% x: export, X/A
% c: scaled consumption, C/A
% y: scaled output, Y/A
% Kp, Fp: variables from Calvo pricing
% imp: import
% k: capital
% i: investment
% Pk: price of capital goods, Q = Pk x P
% mc: marginal cost
% rk: return from using capital
% w: scaled real wage, W/PxA
% n: net worth
% Rk: average return on capital
% omega: cut-off shock level
% ZR: risk premium in financial market, interest rate paid by entrepreneurs (Z) - riskfree rate (R)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

varexo lnRf_eps lnyf_eps phi_eps lng_eps lnA_eps tau_eps lnpif_eps money_eps sigma_eps gamm_eps;

parameters theta epsil beta nu varphi psi rho_R r_pi r_y r_q r_ZR omega_c eta_c psi_f eta_f nu_x eta_g 
    phi_a phi_s eta_a rho_Rf rho_yf rho_phi rho_g rho_A rho_tau rho_pif rho_sigma rho_gamm alpha delta kappa mu we
    lnpcSS lnpmcSS lnqSS lnRSS lnpibarSS lnpicSS lnpxSS lnNSS lnpstarSS afSS 
    lnPhiSS lnsSS lnxSS lncSS lnKpSS lnFpSS lnySS lnimpSS lnkSS lniSS
    lnPkSS lnmcSS lnrkSS lnwSS lnnSS lnRkSS omegaSS ZRSS
    lnRfSS lnyfSS lngSS lnASS tauSS lnpifSS sigmaSS gammSS;

theta = 0.75;            % degree of stickiness in prices
epsil = 6;               % markup papameter
beta = 1.03^(-0.25);    % time discount
nu = 1/epsil;           % subsidy, set it to eliminate the monopoly distortion
varphi = 1;             % elasticity of labor supply
psi = 0.1;                % fraction of wage bill that must be financed in advance in domestic good
nu_x = 0.1;               % fraction of wage bill financed in advance by exporters
psi_f = 0.1  ;              % fraction of wage bill that must be financed in advance in foreign good
rho_R = 0.9;            % smoothing parameter in Taylor rule
r_pi = 1.5;             % weight on expected inflation in Taylor rule
r_y = 0.15;             % weight on output in Taylor rule
r_q = 0.00;             % weight on exchange rate in Taylor rule
r_ZR = 0.00;            % weight on financial risk premium in Taylor rule
omega_c = 0.4;          % approximate share of domestic goods in consumption
eta_c = 1.5;              % elasticity of substitution b/w domestic and foreign goods
eta_f = 1.5;            % elasticity of demand for exports
eta_g = 0.3;            % ratio fo gov't spending 
phi_a = 0.03;           % weight on net foreign assets in risk term
phi_s = 1.5;            % weight on interest rate differential in risk term
eta_a = 0.0;              % steady state foreign assets
rho_Rf = 0.95;          % AR(1) coefficient on shocks
rho_yf = 0.95;
rho_phi = 0.95;
rho_g = 0.95;
rho_A = 0.95;
rho_tau = 0.95;
rho_pif = 0.95;
rho_sigma = 0.95;
rho_gamm = 0.95;   
alpha = 0.33;           % capital income share
delta = 0.025;            % depreciation     
kappa = 0.03;           % parameter on investment adjustment
mu = 0.12;               % monitoring cost as a share of total assets
we = 0.03;                 % this is obtained while getting initial values

%%%%%%%% getting the steady state values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% start from setting pibarSS
pibarSS = 1.005;
% given or assumed
pifSS = pibarSS;
tauSS = 1.0;
ASS = 1;
gammSS = 0.97;
PkSS = 1;

picSS = pibarSS;
pstarSS = ((1-theta*pibarSS^epsil)/(1-theta))/((1-theta*pibarSS^(epsil-1))/(1-theta))^(epsil/(epsil-1));
FpSS = 1/(1-beta*theta*pibarSS^(epsil-1));
KpSS = FpSS*(1-theta*pibarSS^(epsil-1))/(1-theta);
RSS = picSS/beta;
RfSS = RSS*pifSS/picSS;

mcSS = (1-beta*theta*pibarSS)*((epsil-1)/epsil)*KpSS;
sSS = picSS/pifSS;
PhiSS = 1;

% set pcSS*qSS = varphi_w = 1 for convenience
varphi_w = 1;           
pmcSS = varphi_w*(1-psi_f+psi_f*RfSS);
pxSS = (1-nu_x+nu_x*RSS)/varphi_w;
pcSS = (1-omega_c+omega_c*(pmcSS^(1-eta_c)))^(1/(1-eta_c));
qSS = varphi_w/pcSS;

% fix rk and adjust until (A6) is satisfied
rkSS = 0.05;
RkSS = pibarSS*(rkSS+1-delta);

FSS = 0.03;         % fix the bankruptcy rate, F(omegabar), a priori, then proceed to get omegabar and sigma

sigmaSS = 0.5;omegaSS=0.8;
%[omegaSS,sigmaSS] = findomega_sigma(RkSS,RSS,sigmaSS,mu);
%sigmaSS = 0.1;
%omegaSS = findomega(RkSS,RSS,sigmaSS,mu);    % need findomega.m and getomega.com

z_std = (log(omegaSS)+0.5*sigmaSS^2)/sigmaSS;

GSS = normcdf(z_std-sigmaSS,0,1);
GammaSS = omegaSS*(1-FSS)+GSS;

knSS = -1/((GammaSS-mu*GSS)*(RkSS/RSS)-1); 
nSS = we/(1-gammSS*(1-GammaSS)*RkSS*knSS/pibarSS);
kSS = knSS*nSS;   

iSS = delta*kSS;
wSS = (mcSS/((1-nu)*tauSS*((1/(1-alpha))^(1-alpha))*((1/alpha)^alpha)*(rkSS^alpha)*((1-psi+psi*RSS)^(1-alpha))))^(1/(1-alpha));
NSS = kSS*((1-nu)*tauSS*wSS*(1-psi+psi*RSS)/((1-alpha)*mcSS))^(-1/alpha);
ySS = pstarSS*(NSS^(1-alpha))*(kSS^alpha);
gSS = eta_g*ySS;
afSS = eta_a*ySS;

cSS = ((1-eta_g)*ySS+(afSS-sSS*RfSS*afSS/pibarSS)/pcSS*qSS*pxSS -iSS-mu*GSS*RkSS*kSS/pibarSS)/((1-omega_c)*(pcSS^eta_c)+pmcSS*omega_c*((pcSS/pmcSS)^eta_c)/pcSS*qSS*pxSS);
xSS = (afSS+pmcSS*omega_c*((pcSS/pmcSS)^eta_c)*cSS-sSS*RfSS*afSS/pibarSS)/(pcSS*qSS*pxSS);

yfSS = xSS*pxSS^eta_f;

impSS = omega_c*((pcSS/pmcSS)^eta_c)*cSS;

ZRSS = omegaSS*RkSS/(1-(nSS/PkSS*kSS))-RSS;

% taking logs
lnpcSS = log(pcSS);
lnpmcSS = log(pmcSS);
lnqSS = log(qSS);
lnRSS = log(RSS);
lnpibarSS = log(pibarSS);
lnpicSS = log(picSS);
lnpxSS = log(pxSS);
lnNSS = log(NSS);
lnpstarSS = log(pstarSS);
lnPhiSS = log(PhiSS);
lnsSS = log(sSS);
lnxSS = log(xSS);
lncSS = log(cSS);
lnKpSS = log(KpSS);
lnFpSS = log(FpSS);
lnySS = log(ySS);
lnRfSS = log(RfSS);
lnyfSS = log(yfSS);
lngSS = log(gSS);
lnASS = log(ASS);
lnpifSS = log(pifSS);
lnimpSS = log(impSS);
lnkSS = log(kSS);
lniSS = log(iSS);
lnPkSS = log(PkSS);
lnmcSS = log(mcSS);
lnrkSS = log(rkSS);
lnwSS = log(wSS);
lnnSS = log(nSS);
lnRkSS = log(RkSS);
omegaSS = omegaSS;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

model;

// optimal consumption/portfolio/labor supply
// equation (1)
log(beta)+lns(+1)+lnRf+lnPhi+lnc-lnpic(+1)-lnc(+1)-lnA(+1)+lnA;
// equation (2)   
log(beta)+lnR+lnc-lnpic(+1)-lnc(+1)-lnA(+1)+lnA; 
// equation (3) 
log(tau)+varphi*lnN+lnc = lnw;

// price distortion
// equation (4)                                   
exp(lnFp) = 1 + beta*theta*exp(lnFp(+1)+(epsil-1)*lnpibar(+1));      
// equation (5)              
exp(lnKp-lnFp) = ((1-theta*exp(lnpibar*(epsil-1)))/(1-theta))^(1/(1-epsil));       
// equation (6)
(epsil/(epsil-1))*exp(lnmc)+beta*theta*exp(epsil*lnpibar(+1))*exp(lnKp(+1))-exp(lnKp);
// equation (7)
(((1-theta)*(1-theta*exp(lnpibar)^(epsil-1))/(1-theta))^(epsil/(epsil-1))
    +theta*exp(lnpibar)^epsil/exp(lnpstar(-1)))*exp(lnpstar);

// equation (8): import
lnimp - lnc = log(omega_c)+eta_c*(lnpc-lnpmc);                          
                   
// equation (9): price of consumption good
exp(lnpc) = (1-omega_c+omega_c*exp(lnpmc)^(1-eta_c))^(1/(1-eta_c));   

// equation (10): inflation of consumption good
exp(lnpic) = exp(lnpibar)*((1-omega_c+omega_c*exp(lnpmc)^(1-eta_c))
    /(1-omega_c+omega_c*exp(lnpmc(-1))^(1-eta_c)))^(1/(1-eta_c));

// equation (11): demand for export
lnx = -eta_f*lnpx+lnyf;

// equation (12)
lnpmc = lnpc+lnq+log(1-psi_f+psi_f*exp(lnRf));  

// equation (13)
lnq+lnpx+lnpc = log(1-nu_x+nu_x*exp(lnR));

// equation (14)
lnq-lnq(-1)=lns+lnpif-lnpic;

// Taylor rule
// equation (15)
// when NOT responding to real exchange rates and risk premium: r_q = 0, r_ZR = 0
lnR - lnRSS = rho_R*(lnR(-1)-lnRSS)+(1-rho_R)*(r_pi*(lnpic(+1)-lnpicSS)           
    +r_y*(lny(+1)-lnySS)+r_q*(lnq(+1)-lnqSS))+r_ZR*(log(ZR(+1)+exp(lnR))-lnR)+money_eps;
%lnR - lnRSS = rho_R*(lnR(-1)-lnRSS)+(1-rho_R)*(r_pi*(lnpic-lnpicSS)           
%    +r_y*(lny-lnySS)+r_q*(lnq-lnqSS))+r_ZR*(log(ZR+exp(lnR(-1)))-lnR(-1))+money_eps;

// aggregate condition (2 equations): without financial frictions
// equation (16) and (17)
%lny = lnpstar+lnN;
%exp(lny) = (1-omega_c)*exp(lnpc*eta_c+lnc)+exp(lnx)+exp(lng);

// equation (18): current account
af+exp(lnpmc+log(omega_c)+eta_c*(lnpc-lnpmc)+lnc) = exp(lnq+lnpc+lnpx+lnx)
    +exp(log(af(-1)+lns+lnRf(-1)+lnPhi(-1)-lnpibar-lnA+lnA(-1))); // previosly dropped lnc in LHS by mistake

// equation (19): risk term
lnPhi = -phi_a*(af-afSS)-phi_s*(exp(lnRf)-exp(lnR)-exp(lnRfSS)+exp(lnRSS))+ phi;

// equations for marginal cost
// equation (20)
lnmc = log(1-nu)+log(tau)-(1-alpha)*log(1-alpha)-alpha*log(alpha)+alpha*lnrk+(1-alpha)*(lnw+log(1-psi+psi*exp(lnR)));
%mc = (1-nu)*tau*(exp(lnN*varphi+lnc))*(1-psi+psi*exp(lnR));
// equation (21)
lnmc = log(1-nu)+log(tau)+lnw+log(1-psi+psi*exp(lnR))-log(1-alpha)-alpha*(lnk(-1)-lnA+lnA(-1)-lnN);

// equation (22): production function
lny = lnpstar+(1-alpha)*lnN+alpha*(lnk(-1)-lnA+lnA(-1));

// equation (23): law of motion for capital
Sadj = 0.5*kappa*(exp(lnA-lnA(-1)+lni-lni(-1))-1)^2;    
Sadjprime = kappa*(exp(lnA-lnA(-1)+lni-lni(-1))-1);
exp(lnk) = (1-delta)*exp(lnk(-1))*exp(lnA(-1)-lnA) + (1-Sadj)*exp(lni);

// equation (24): optimal investment
(1-Sadj-Sadjprime*exp(lnA-lnA(-1)+lni-lni(-1)))*exp(lnPk) = 1-beta*exp(lnc-lnc(+1)-lnA(+1)+lnA)*exp(lnPk(+1))*Sadjprime(+1)*(exp(lnA(+1)-lnA+lni(+1)-lni)^2);

// equation (25), check t or t+1
exp(lnRk) = exp(lnpibar)*(exp(lnrk)+(1-delta)*exp(lnPk))/exp(lnPk(-1));


// for functional forms, see the footnote in 'financial frictions' section
// Ft means F(omega(t))
Fprime = (1/sqrt(2*3.141592))*exp(-0.5*((log(omega)+sigma(-1)^2/2)/sigma(-1))^2)/(omega*sigma(-1));
Gt = normcdf((log(omega)+sigma(-1)^2/2)/sigma(-1)-sigma(-1),0,1);
Ft = normcdf((log(omega)+sigma(-1)^2/2)/sigma(-1),0,1);
Gammat = omega*(1-Ft)+Gt;


// n (net worth)  should be lagged like k
// equation (26): zero profit condition
exp(lnPk(-1)+lnk(-1)-lnn(-1))*(Gammat-mu*Gt)*exp(lnRk-lnR(-1)) = exp(lnPk(-1)+lnk(-1)-lnn(-1))-1;

// equation (27): optimality condition for entrepreneurial contract
(1-Gammat)*exp(lnRk-lnR(-1)) + ((1-Ft)/(1-Ft-mu*omega*Fprime))
    *(exp(lnRk-lnR(-1))*(Gammat-mu*Gt)-1);

// equation (28): law of motion for aggregate wealth
exp(lnn) = gamm*(1-Gammat)*exp(lnRk+lnPk(-1)+lnk(-1)-lnpibar-lnA+lnA(-1))+we;

// equation (17) or (29): aggregate condition
exp(lny) = (1-omega_c)*exp(lnpc*eta_c+lnc)+exp(lnx)+exp(lng)+exp(lni)+d;
// need to define d
d = mu*Gt*omega*exp(lnPk(-1)+lnRk+lnk(-1)-lnpibar-lnA+lnA(-1));

// risk premium in financial market
ZR = omega*exp(lnRk)/(1-exp(lnn(-1)-lnPk(-1)-lnk(-1)))-exp(lnR(-1));


// shocks
lnRf = (1-rho_Rf)*lnRfSS + rho_Rf*lnRf(-1) + lnRf_eps; 
lnyf = (1-rho_yf)*lnyfSS + rho_yf*lnyf(-1) + lnyf_eps;
phi = rho_phi*phi(-1) + phi_eps; 
lng = (1-rho_g)*lngSS+rho_g*lng(-1)+lng_eps;
lnA = (1-rho_A)*lnASS + rho_A*lnA(-1) + lnA_eps;
tau = (1-rho_tau)*tauSS + rho_tau*tau(-1) + tau_eps;
lnpif = (1-rho_pif)*lnpifSS + rho_pif*lnpif(-1) + lnpif_eps;
log(sigma/sigmaSS) = rho_sigma*log(sigma(-1)/sigmaSS) + sigma_eps;
log(gamm/gammSS) = rho_gamm*log(gamm(-1)/gammSS) + gamm_eps;

end;

initval;
lnpc = lnpcSS;
lnpmc = lnpmcSS;
lnq = lnqSS;
lnR = lnRSS;
lnpibar = lnpibarSS;
lnpic = lnpicSS;
lnpx = lnpxSS;
lnN = lnNSS;
lnpstar = lnpstarSS;
af = afSS;
lnPhi = lnPhiSS;
lns = lnsSS;
lnx = lnxSS;
lnc = lncSS;
lnKp = lnKpSS;
lnFp = lnFpSS;
lny = lnySS;
lnimp = lnimpSS;
phi = 0;
lnRf = lnRfSS;
lnyf = lnyfSS;
lng = lngSS;
lnA = lnASS;
tau = tauSS;        
lnpif = lnpifSS;
lnk = lnkSS;
lni = lniSS;
lnPk = lnPkSS;
lnmc = lnmcSS;
lnrk = lnrkSS;
lnw = lnwSS;
lnn = lnnSS;
lnRk = lnRkSS;
omega = omegaSS;
sigma = sigmaSS;
gamm = gammSS;
end;

shocks;
var lnRf_eps;
stderr 0.01;
var lnyf_eps;
stderr 0.01;
var phi_eps;
stderr 0.01;
var money_eps;
stderr 0.01;
var lng_eps;
stderr 0.01;
var lnA_eps;
stderr 0.01;
var tau_eps;
stderr 0.01;
var lnpif_eps;
stderr 0.01;
var sigma_eps;
stderr 0.01;
var gamm_eps;
stderr 0.01;
end;

steady;
check;

stoch_simul(periods=5000,simul_seed=1,irf=20,order=1,nograph) 
lny lnc lni lnpic lnR ZR lnn lnRk omega;

IRFplot(lny_money_eps,lnc_money_eps,lni_money_eps,lnpic_money_eps, lnRSS, lnR_money_eps, ZR_money_eps, lnn_money_eps, lnRk_money_eps, omega_money_eps);
%hout = suptitle('response to monetary policy shock');       % use suptitle.m written by CTW