% MSDSGE with Rotemberg pricing, 02/07/2015
%  
% 

%----------------------------------------------------------------
% 0. Housekeeping
%----------------------------------------------------------------

close all

%----------------------------------------------------------------
% 1. Defining variables
%----------------------------------------------------------------

// Endogenous Variables
var c, llambda, R, ppi, W, mc, H, Y, Util, Welf, xstar, xcross;

// Exogenous Variables
var A, Ystar, a, zzeta, g, G, d, m, Ycross;

// Shocks
varexo epsa epsg epsd epsm;

// Paramters 
parameters nbeta ttau oomega ggamma gbar ppsipi ppibar ppsiy rrhor rrhoa rrhog
  ssigmaa ssigmag pphip iiota ddelta mmu rrhod rrhom ssigmad ssigmam 
  Y_SS ppi_SS c_SS R_SS a_SS g_SS A_SS H_SS Util_SS Welf_SS zzeta_SS
  llambda_SS W_SS mc_SS G_SS Ystar_SS xstar_SS m_SS d_SS Ycross_SS xcross_SS;

%----------------------------------------------------------------
% 2. Calibration
%----------------------------------------------------------------

// Parameters
nbeta = 0.9987;
ttau = 2.5;
oomega = 0.0039;
ggamma = 0.66;
gbar = 1.25;
ppibar = 1;
ppsipi = 1.5;
ppsiy = 0.5;
rrhor = 0.7;
rrhoa = 0.8;
rrhog = 0.9;
ssigmaa = 0.003;
ssigmag = 0.012;
pphip = 160;
iiota = 0;

ddelta = 0.9999/0.9987;
mmu = 0.1;
rrhod = 0.9;
rrhom = 0.9;
ssigmad = 0.01;
ssigmam = 0.01;


// Analytic Steady State
m_SS = mmu;
ppi_SS = ppibar;
a_SS = 1;
g_SS = gbar;
A_SS = exp(oomega);
mc_SS = 1/(1+m_SS) + m_SS/(1+m_SS) * (1-nbeta)*pphip*(ppi_SS-ppi_SS^iiota)*ppi_SS;
zzeta_SS = (gbar-1)/gbar;
c_SS = (mc_SS)^(1/ttau);
Y_SS = c_SS/(1-pphip/2*(ppi_SS-ppi_SS^iiota)^2 - zzeta_SS);
R_SS = ppibar*exp(oomega)/nbeta;
H_SS = Y_SS;

Util_SS = d_SS*((c_SS^(1-ttau)-1)/(1-ttau) - H_SS);
Welf_SS = Util_SS/(1-nbeta);

llambda_SS = ddelta/mc_SS;
W_SS = mc_SS;
G_SS = ((gbar-1)/gbar)*Y_SS;
Ystar_SS = ((1/(1+m_SS))^(1/ttau))*gbar;
xstar_SS = Y_SS/Ystar_SS;
d_SS = ddelta;
Ycross_SS = 1/(1-zzeta_SS);
xcross_SS = Y_SS/Ycross_SS;



// Initialize Loop
for ii = 1:length(X(1,:))
    ppsipi = X(1,ii);
    ppsiy = X(2,ii);
    rrhor = X(3,ii);
    ppibar = X(4,ii);

// Steady State Adjustment
m_SS = mmu;
ppi_SS = ppibar;
a_SS = 1;
g_SS = gbar;
A_SS = exp(oomega);
mc_SS = 1/(1+m_SS) + m_SS/(1+m_SS) * (1-nbeta)*pphip*(ppi_SS-ppi_SS^iiota)*ppi_SS;
zzeta_SS = (gbar-1)/gbar;
c_SS = (mc_SS)^(1/ttau);
Y_SS = c_SS/(1-pphip/2*(ppi_SS-ppi_SS^iiota)^2 - zzeta_SS);
R_SS = ppibar*exp(oomega)/nbeta;
H_SS = Y_SS;

Util_SS = d_SS*((c_SS^(1-ttau)-1)/(1-ttau) - H_SS);
Welf_SS = Util_SS/(1-nbeta);

llambda_SS = ddelta/mc_SS;
W_SS = mc_SS;
G_SS = ((gbar-1)/gbar)*Y_SS;
Ystar_SS = ((1/(1+m_SS))^(1/ttau))*gbar;
xstar_SS = Y_SS/Ystar_SS;
d_SS = ddelta;
Ycross_SS = 1/(1-zzeta_SS);
xcross_SS = Y_SS/Ycross_SS;

%----------------------------------------------------------------
% 3. Model   
%----------------------------------------------------------------

// model in logs 
model;

Util = exp(d)*((exp(c)^(1-ttau)-1)/(1-ttau) - exp(H));
Welf = Util + nbeta*Welf(+1);

exp(d)/exp(c)^ttau - exp(llambda);
exp(llambda) - nbeta*(exp(llambda(+1))/A(+1))*(exp(R)/exp(ppi(+1)));
exp(d) - exp(llambda)*exp(W);
exp(W) - exp(mc);

pphip*(exp(ppi)-ppi_SS^iiota)*exp(ppi)*exp(Y) + (1/exp(m))*exp(Y) - ((1+exp(m))/exp(m))*exp(mc)*exp(Y) - nbeta*(exp(llambda(+1))/exp(llambda))*pphip*(exp(ppi(+1))-ppi_SS^iiota)*exp(ppi(+1))*exp(Y(+1));

exp(G) - exp(zzeta)*exp(Y);
exp(Y) - exp(c) - exp(G) -(pphip/2)*((exp(ppi) - ppi_SS^iiota)^2)*exp(Y);
exp(H) - exp(Y);
exp(Ystar) - (1/(1+exp(m)))^(1/ttau) * 1/(1-exp(zzeta));

log(A) - oomega - a;
a - rrhoa*a(-1) - ssigmaa*epsa;
exp(g) - 1/(1-exp(zzeta));
g - (1-rrhog)*log(g_SS) - rrhog*g(-1)  - ssigmag*epsg;
d - (1-rrhod)*log(d_SS) - rrhod*d(-1)  - ssigmad*epsd;
m - (1-rrhom)*log(m_SS) - rrhom*m(-1)  - ssigmam*epsm;

exp(xstar) - exp(Y)/exp(Ystar);
exp(Ycross) - 1/(1-exp(zzeta));
exp(xcross) - exp(Y)/exp(Ycross);

// Monetary Policy Rule
(exp(R(-1))/R_SS)^rrhor*((exp(ppi)/ppi_SS)^ppsipi*(exp(Y)/Y_SS)^ppsiy)^(1-rrhor) - exp(R)/R_SS;



end;

%----------------------------------------------------------------
% 4. Computation
%----------------------------------------------------------------

// Analytic Steady State
initval;
  Util = Util_SS;
  Welf = Welf_SS;
  H = log(H_SS);  
  Y = log(Y_SS);
  ppi = log(ppi_SS);
  c = log(c_SS);
  R = log(R_SS);
  a = log(a_SS);
  g = log(g_SS);
  A = A_SS;
  zzeta = log(zzeta_SS);
  llambda = log(llambda_SS);
  W = log(W_SS);
  mc = log(mc_SS);
  Ystar = log(Ystar_SS);
  xstar = log(xstar_SS);
  d = log(d_SS);
  m = log(m_SS);
  G = log(G_SS);
  Ycross = log(Ycross_SS);
  xcross = log(xcross_SS);
end;


resid;
check;


// Setting Variances of Shocks
shocks;
    var epsa = 1;
    var epsg = 1;
    var epsd = 1;
    var epsm = 1;

end;

// Stochastic Simulation

// End for-loop

      try
          stoch_simul(order=2,nocorr,nofunctions,irf=0,nograph);
          X(5,ii) = oo_.mean(10);
          X(6,ii) = oo_.mean(3)- 2*sqrt(oo_.var(3,3));
      catch
          X(5,ii) = -Inf;
      end

      ii

end