%   Keinsley (2017) Biased Indexation Project
%
%

%-----
% ENDOGENOUS VARIABLES
%-----

    % Household
    var c infl b w infl1 tau l d r lam3;
    
    % Retailer
    var y1 y2 y infl2;
    
    % Poduction
    var z a1 a2 b1 b2 infl1_star infl2_star mc Lam;
    
    % Government
    var g;
    
    % Market Clearing
    var v1 v2 d1 d2;
    
%-----
% EXOGENOUS VARIABLES
%-----

    % Production
    varexo e_z;

    % Government
    varexo e_g e_r;

%-----
% PARAMETERS
%-----

    % Household
    parameters sigc sigl eta taxinc;
    parameters cs infls bs ws infl1s taus ls ds rs lam3s;

    % Retailer
    parameters delta theta;
    parameters y1s y2s ys infl2s;

    % Production
    parameters ep phi1 phi2 rhoz sigz;
    parameters zs a1s a2s b1s b2s infl1_stars infl2_stars mcs Lams;

    % Government
    parameters thetat gammat rhog rhor rhop sigg sigr;
    parameters gs;

    % Market Clearing
    parameters v1s v2s d1s d2s;

%-----
% CALIBRATION
%-----

    % Structural/Dynamic Parameters
    delta = 0.20;
    theta = 5;
    phi1  = 0.30;
    phi2  = 0.70;
    beta  = 0.95;
    ep    = 6;
    thetat= 0.901;
    gammat= 0.134;
    sigc  = 1.39;
    sigl  = 1.2;
    eta   = 4.358479527824188;
    rhoz  = 0.60;
    rhor  = 0.20;
    rhog  = 0.40;
    rhop  = 1.50;
    sigz  = 0.01;
    sigr  = 0.0025;
    sigg  = 0.01;

    % Steady State Values
    cs          = 0.8;
    infls       = 1.02;
    bs          = 6.142513339476156;
    ws          = 3.839865087795592;
    infl1s      = 1.03;
    taus        = 0.102561694052862;
    ls          = 1/3;
    ds          = -0.226794285398488;
    rs          = 1.04;
    lam3s       = -12.910604908776458;
    y1s         = 0.539804876821161;
    y2s         = 0.509274647340516;
    ys          = 1;
    infl2s      = 1.017574580483186;
    zs          = 3.161785531426725;
    a1s         = -9.953495192634451;
    a2s         = -8.945548265470215;
    b1s         = -11.596305078797418;
    b2s         = -10.296272059679081;
    infl1_stars = 1.03;
    infl2_stars = 1.015859662091025;
    mcs         = 1.214460958730142;
    Lams        = 0.95;
    gs          = 0.20;
    v1s         = 0.976212475776451;
    v2s         = 1.034734907754879;
    d1s         = -0.093959263882966;
    d2s         = -0.132228272883066;
    taxinc      = ws*ls+(rs-1)*bs/infls+ds;


%-----
% MODEL EQUATIONS
%-----

model (linear);

    % Household
cs*c+bs*b = (1-taus)*ws*ls*(w+l)+(1-taus)*ds*d + (1-taus*(rs-1))/infls*bs*b(-1)
            -bs/infls^2*(1-taus*(rs-1))*infl - taxinc*tau-bs/infls*taus*r(-1);

tau = thetat*gammat/(taxinc*infl1s)*(ws*ls*(w+l)+bs/infls*r(-1)+(rs-1)*bs/infls*b(-1)
      -(rs-1)*bs/infls^2*infl+ds*d) + thetat*gammat/infl1s^(gammat+1)*infl1;

sigl*l = w+1/(cs^(-sigc)*(1+taus)+lam3s*gammat*(1-taus)*taxinc^(-1))*
         (-cs^(-sigc)*sigc*(1+taus)*c+(-cs^(-sigc)-lam3s*gammat/taxinc)*tau
         +lam3s*(1-taus)*gammat/taxinc*(lam3-ws*ls*(w+l)/taxinc-(rs-1)*bs/(taxinc*infls)*b(-1)
         -bs/(taxinc*infls)*r(-1)+bs*(rs-1)/(taxinc*infls^2)*infl-ds/taxinc*d));

-sigc*c = sigc/infls*(rs*(1-taus)-taus)*c(+1) + 1/infls^2*(taus*(rs-1)-rs+lam3s*gammat*(1-taus)*(rs-1)*(ws*ls+ds)/(c^(-sigc)*taxinc^2))*infl(+1)
          +(1-taus)/infls*(1-lam3s*gammat/taxinc*cs^sigc*(1-(rs-1)*bs/infls/taxinc/infls))*r
          -(rs-1)/infls*(1-lam3s*gammat/taxinc*cs^sigc)*tau(+1)
          +lam3s*gammat*(1-taus)*(rs-1)/taxinc*cs^sigc/infls*(-lam3(+1)+ws*ls/taxinc*(w(+1)+l(+1))
          +(rs-1)*bs/infls/taxinc*b + ds/taxinc*d(+1));

y = delta^(1/theta)*(y1s/ys)^((theta-1)/theta)*y1 + (1-delta)^(1/theta)*(y2s/ys)^((theta-1)/theta)*y2;

infl = delta*(infl1s/infls)^(-theta)*infl1 + (1-delta)*(infl2s/infls)^(-theta)*infl2;

z = rhoz*z(-1)+e_z;

mc = w - z;

Lam = sigc*(c-c(+1));

infl1_star = a1-b1;

infl2_star = a2-b2;

a1 = (1/infls+(2-ep)/a1s*(infl1s/infls)^(ep-2)*mcs*y1s)*infl 
     +mcs*(ep-2)/a1s*(infl1s/infls)^(ep-3)*y1s*infl1
     +infls/a1s*(infl1s/infls)^(ep-2)*mcs*y1s*(mc+y1)
     +infls^ep*Lams*(Lam(+1)+a1(+1))+(ep-1)*Lams*infls^(ep-1)*infl(+1);

b1 = (ep-1)/b1s*(infl1s/infls)^(ep-1)*y1s/infls*(infl1s/infls*infl1-infl+infls/(ep-1)*y1)
     +Lams*infls^(ep-1)*(Lam(+1)+(ep-1)/infls*infl(+1) + b1(+1));

a2 = (1/infls+(2-ep)/a2s*(infl2s/infls)^(ep-2)*mcs*y2s)*infl 
     +mcs*(ep-2)/a2s*(infl2s/infls)^(ep-3)*y2s*infl2
     +infls/a2s*(infl2s/infls)^(ep-2)*mcs*y2s*(mc+y2)
     +infls^ep*Lams*(Lam(+1)+a2(+1))+(ep-1)*Lams*infls^(ep-1)*infl(+1);

b2 = (ep-1)/b2s*(infl2s/infls)^(ep-1)*y2s/infls*(infl2s/infls*infl2-infl+infls/(ep-1)*y2)
     +Lams*infls^(ep-1)*(Lam(+1)+(ep-1)/infls*infl(+1) + b2(+1));

infl1 = (1-phi1)*(infl1_stars/infl1s)^(-ep)*infl1_star;

infl2 = (1-phi2)*(infl2_stars/infl2s)^(-ep)*infl2_star;

r = rhor*r(-1)+(1-rhor)*rhop*infl+e_r;

gs*g = taxinc*tau+taus*(ws*ls*(w+l)+ds*d)+bs*(b+(taus*(rs-1)-rs)/infls*b(-1))
       +bs/(gs*infls)*((taus*(rs-1)-rs)/infls*infl-(1-taus)*r(-1));

g = rhog*g(-1)+e_g;

y1 = z+l-v1;

y2 = z+l-v2;

v1 = ep/v1s*((1-phi1)*(infl1_stars/infl1s)^(-ep)/infl1_stars*infl1_star
        +1/infl1s*((1-phi1)*(infl1_stars/infl1s)^(-ep)+phi1*(infl1s/infls)^(-ep)*v1s)*infl1)
      +phi1*(infl1s/infls)^(-ep)*(ep/infls*infl(-1)-ep/infl1s*infl1(-1)+v1(-1));

v2 = ep/v2s*((1-phi2)*(infl2_stars/infl2s)^(-ep)/infl2_stars*infl2_star
        +1/infl2s*((1-phi2)*(infl2_stars/infl2s)^(-ep)+phi2*(infl2s/infls)^(-ep)*v2s)*infl2)
      +phi2*(infl2s/infls)^(-ep)*(ep/infls*infl(-1)-ep/infl2s*infl2(-1)+v2(-1));

ds*infls*d = d1s*infl1s*d1 + d2s*infl2s*d2 + d1s*infl1 + d2s*infl2 - (d1s*infl1s+d2s*infl2s)/infls*infl;

d1 = y1 - mcs*v1s*infls/infl1s/(1-mcs*v1s*infls/infl1s)*(mc+v1+infl/infls-infl1/infl1s);

d2 = y2 - mcs*v2s*infls/infl2s/(1-mcs*v2s*infls/infl2s)*(mc+v2+infl/infls-infl2/infl2s);

y = cs/ys*c+gs/ys*g;

end;

check;

shocks;
var e_z = sigz^2; 
var e_g = sigg^2;
var e_r = sigr^2;
end;


%-----
% SOLVING THE RATIONAL EXPECTATIONS MODEL
%-----
stoch_simul(order=1, irf=20) y infl infl1 infl2 b tau r;