%**************************************************************************************************
%                          HOUSING PRICES, WELFARE AND MONETARY POLICY                            *
%               Berrak Buyukkarabacak and Olena Mykhaylova, University of Richmond                *
%                                         March 2009                                              *
%**************************************************************************************************


periods     20000;

%   All variables are expressed in logs, expect for
%       Lagrange multipliers: Lamh, Lame, Lamc;
%       profits: divr, divc;
%       adjustment costs: adjHp, adjHh, adjHe, adjK;
%       credit constraints: mh, me, mc, mb;
%       bond holdings: bp, bh, be, bc, f.
%   NOTE: pateint household borrowings bp are REDEFINED as savings.

var         cp, ch, ce, ip, ih, itot,                   % Consumption and capital investment
            lp, lh, npe, npc, nhe, nhc, wp, wh,         % Labor market variables
            q, hp, hh, he, h,                           % Housing variables
            bp, bh, be, bc, f,                          % Lending and borrowing
            rp, rh, re, rc, r, rgap,                    % (Net) interest rates
            y, x, ke, kp, rk, markup,                   % Production variables
            prn, pa, pb,                                % Retailer price setting
            wpn, wpb, wpa, whn, whb, wha, wm,           % Consumer wage setting
            Lamh, Lame, Lamc,                           % Lagrange multipliers
            ginf, inf, divr, divc,                      % Inflation, profits
            a, mh, me, mc, mb,                          % Productivity and credit constraints
            mhvar, mevar, mcvar, mbvar,                 % Credit constraints, auxiliary variables
            adjHp, adjHh, adjHe, adjK,                  % Adjustment costs, auxiliary variables
            vp, vh,                                     % Consumer utility
            zeta, ksi, psai, rho, phi, tao, stigma,     % Parameters necessary for steady state calculations
            zbig, aux1, aux2, aux3,
            A1, A11, A2, A21, A3, A31, z,
            rhy, chy, pcy, bcy, ihy, aggc, aggs, gdp,   % Useful steady state ratios
            pc,                                         % Private borrowing
            fcp, fch, fce, fip, fih, fitot,
            flp, flh, fnpe, fnpc, fnhe, fnhc,
            fwp, fwh,
            fq, fhp, fhh, fhe, fh,
            fbp, fbh, fbe, fbc, ff,
            frp, frh, fre, frc, fr,
            frgap, fy, fx, fke, fkp, frk,
            fmarkup, fprn, fpa, fpb, 
            fwpn, fwpb, fwpa, fwhn, fwhb, fwha,
            fLamh, fLame, fLamc,
            fginf, finf, fgdp,
            fadjHp, fadjK, fadjHh, fadjHe;
            
varexo      epsa, epsi,                 % Shocks to productivity and the world interest rate
            epsh, epse, epsc, epsb,     % Shocks to credit constraints
            epspi;                      % Shock to the union-wide inflation


% PARAMETER RESTRICTIONS
% (1/betaP - mbss*(1/(betaP*(1+Rss))-1)+kappa)^(-1) > betaH
% (1/betaP - mbss*(1/(betaP*(1+Rss))-1)+kappa)^(-1) > betaE

parameters  betaP, betaH, betaE, betaC, theta, chi, gama,
            phiP, phiH, phiE, psiK, delta, deltaH,
            mu, alpha, nu, eta, omegaP, omegaW, kappa, rhoI, rhoR, rhoA,
            rhoH, rhoE, rhoC, rhoB, mhss, mess, mcss, mbss, sigma, phiW,
            Rss, rhoPi, rhoGap, rhoQ, tl, gdpss, qss;


            betaP = 0.99;           % Patient household's discount factor
            betaH = 0.95;           % Impatient household's discount factor
            betaE = 0.962;          % Entrepreneur's discount factor
            betaC = 0.962;          % Construction sector discount factor
            theta = 1;              % Consumption risk aversion
            chi = 3;                % Inverse of Frisch labor supply elasticity
            gama = 0.25;            % Weight on housing services for households
            phiP = 0;               % Housing adjustment costs, patient households
            phiH = 0;               % Housing adjustment costs, impatient households
            phiE = 0;               % Housing adjustment costs, entrepreneurs
            psiK = 9;               % Capital adjustment costs
            delta = 0.03;           % Capital depreciation rate
            deltaH = 0.0055;        % Housing depreciation rate
            mu = 0.30;              % Share of capital in production
            alpha = 0.79;           % Income share of patient households
            nu = 0.11;              % Share of housing in production
            eta = 0.10;             % Share of land in construction
            omegaP = 0.75;          % Price stickiness parameter
            omegaW = 0.75;          % Wage stickiness parameter
            kappa = 0.03;           % Bank's overhead costs
            rhoI = 0.63;            % Inertia in the home Taylor rule
            rhoR = 0.77;            % Inertia in the world interest rate
            rhoA = 0.67;            % Inertia in the productivity process
            rhoH = 0.32;            % Inertia in the credit process, impatient
            rhoE = 0.65;            % Inertia in the credit process, entrepreneur
            rhoC = 0.65;            % Inertia in the credit process, construction
            rhoB = 0.92;            % Inertia in the credit process, banks
            mhss = 0.66;            % Fraction of collateral to be used against borrowing, impatient
            mess = 0.27;            % Fraction of collateral to be used against borrowing, entrepreneur
            mcss = 0.27;            % Fraction of collateral to be used against borrowing, construction
            mbss = 0.34;            % Fraction of collateral to be used against borrowing, banks (0.34)
            sigma = 3.9;            % Elasticity in the goods varieties aggregator
            phiW = 4.3;             % Elasticity in the goods varieties aggregator
            Rss = 0.0025;           % Steady state value of the world interest rate
            rhoPi = 2.3;            % Taylor weight on inflation
            rhoGap = 0.39;          % Taylor weight on output gap (0.39)
            rhoQ = 0;               % Taylor weight on housing inflation
            tl = 0;                 % Log of land available for construction
            gdpss = 0.622475;       % Steady state value of (log) gdp
            qss = 0.311787;         % Steady state value of housing price


model;

    % ***************************************    HOUSEHOLDS    *************************************
    % exp(q) * (exp(hp) - exp(hp(-1))) + adjHp + exp(cp) + bp + exp(ip) + adjK = 
    %    exp(wp) * exp(lp) + (1+rp(-1)) * bp(-1) / (1+inf) + rk * exp(kp(-1)) + divr + divc;
    adjHp = (phiP / 2) * exp(q) * exp(hp(-1)) * (exp(hp)/exp(hp(-1)) - 1)^2;
    adjK = (psiK / 2) * exp(kp(-1)) * (exp(ip)/exp(kp(-1)) - delta)^2;
    -rp = log(betaP) - theta * (cp(+1) - cp) - inf(+1);
    % chi * lp = wp - theta * cp;
    (exp(cp))^(-theta) * exp(q) * (1 + phiP * (exp(hp)/exp(hp(-1)) - 1)) = gama * (exp(hp))^(-theta) +
        betaP * (exp(cp(+1)))^(-theta) * exp(q(+1)) * (1 + (phiP / 2) * ((exp(hp(+1))/exp(hp))^2 -1));
    (exp(cp))^(-theta) * (1 + psiK * (exp(kp)/exp(kp(-1)) - 1)) = betaP * (exp(cp(+1)))^(-theta) *
        (1 - delta + (psiK / 2) * ((exp(kp(+1))/exp(kp))^2 -1) + rk(+1));

    exp(q) * (exp(hh) - exp(hh(-1))) + adjHh + exp(ch) + (1+rh(-1)) * bh(-1) / (1+inf) =
        exp(wh) * exp(lh) + bh;
    adjHh = (phiH / 2) * exp(q) * exp(hh(-1)) * (exp(hh)/exp(hh(-1)) - 1)^2;
    bh = mh * exp(q(+1)) * exp(hh) * (inf(+1)+1) / (rh+1);
    (exp(ch))^(-theta) = betaH * (1+rh) * (exp(ch(+1)))^(-theta) / (1+inf(+1)) + Lamh * (1+rh);
    % chi* lh = wh - theta * ch;
    (exp(ch))^(-theta) * exp(q) * (1 + phiH * (exp(hh)/exp(hh(-1)) - 1)) - Lamh * mh * exp(q(+1)) * (1+inf(+1)) =
        gama * (exp(hh))^(-theta) + betaH * (exp(ch(+1)))^(-theta) * exp(q(+1)) * (1 + (phiH / 2) *
        ((exp(hh(+1))/exp(hh))^2 -1));

    wm = phiW / (phiW - 1);
    exp(wpn * (phiW*chi+1)) = wm * wpb / wpa;
    % exp(wp) = exp(lp * chi - theta * cp);
    wpb = exp(wp * phiW*(1+chi) + lp * (1+chi)) + omegaW * betaP * wpb(+1) * (inf(+1)+1)^(phiW*(1+chi));
    wpa = exp(-theta * cp + phiW * wp + lp) + omegaW * betaP * wpa(+1) * (inf(+1)+1)^(phiW-1);
    exp(wp * (1-phiW)) = (1 - omegaW) * exp(wpn * (1-phiW)) + omegaW * exp(wp(-1) * (1-phiW)) * (inf+1)^(phiW-1);

    % exp(wh) = (whb/wha)^(1/(phiW*chi+1));
    exp(wh) = exp(lh * chi - theta * ch);
    whb = exp(wh * phiW*(1+chi) + lh * (1+chi)) + omegaW * betaH * whb(+1) * (inf(+1)+1)^(phiW*(1+chi));
    wha = exp(-theta * ch + phiW * wh + lh) + omegaW * betaH * wha(+1) * (inf(+1)+1)^(phiW-1);
    exp(wh * (1-phiW)) = (1 - omegaW) * exp(whn * (1-phiW)) + omegaW * exp(wh(-1) * (1-phiW)) * (inf+1)^(phiW-1);
    % wh = whn;

    vp = (cp - lp^(1+chi) / (1 + chi) + gama * hp) + betaP * vp(+1);
    vh = (ch - lh^(1+chi) / (1 + chi) + gama * hh) + betaH * vh(+1);


    % ************************************    FLEX HOUSEHOLDS    ***********************************
    fadjHp = (phiP / 2) * exp(fq) * exp(fhp(-1)) * (exp(fhp)/exp(fhp(-1)) - 1)^2;
    fadjK = (psiK / 2) * exp(fkp(-1)) * (exp(fip)/exp(fkp(-1)) - delta)^2;
    -frp = log(betaP) - theta * (fcp(+1) - fcp) - finf(+1);
    (exp(fcp))^(-theta) * exp(fq) * (1 + phiP * (exp(fhp)/exp(fhp(-1)) - 1)) = gama * (exp(fhp))^(-theta) +
        betaP * (exp(fcp(+1)))^(-theta) * exp(fq(+1)) * (1 + (phiP / 2) * ((exp(fhp(+1))/exp(fhp))^2 -1));
    (exp(fcp))^(-theta) * (1 + psiK * (exp(fkp)/exp(fkp(-1)) - 1)) = betaP * (exp(fcp(+1)))^(-theta) *
        (1 - delta + (psiK / 2) * ((exp(fkp(+1))/exp(fkp))^2 -1) + frk(+1));

    exp(fq) * (exp(fhh) - exp(fhh(-1))) + fadjHh + exp(fch) + (1+frh(-1)) * fbh(-1) / (1+finf) =
        exp(fwh) * exp(flh) + fbh;
    fadjHh = (phiH / 2) * exp(fq) * exp(fhh(-1)) * (exp(fhh)/exp(fhh(-1)) - 1)^2;
    fbh = mh * exp(fq(+1)) * exp(fhh) * (finf(+1)+1) / (frh+1);
    (exp(fch))^(-theta) = betaH * (1+frh) * (exp(fch(+1)))^(-theta) / (1+finf(+1)) + fLamh * (1+frh);
    (exp(fch))^(-theta) * exp(fq) * (1 + phiH * (exp(fhh)/exp(fhh(-1)) - 1)) - fLamh * mh * exp(fq(+1)) * (1+finf(+1)) =
        gama * (exp(fhh))^(-theta) + betaH * (exp(fch(+1)))^(-theta) * exp(fq(+1)) * (1 + (phiH / 2) *
        ((exp(fhh(+1))/exp(fhh))^2 -1));

    exp(fwpn * (phiW*chi+1)) = wm * fwpb / fwpa;
    fwpb = exp(fwp * phiW*(1+chi) + flp * (1+chi));
    fwpa = exp(-theta * fcp + phiW * fwp + flp);
    exp(fwp * (1-phiW)) = exp(fwpn * (1-phiW));

    exp(fwh) = exp(flh * chi - theta * fch);
    fwhb = exp(fwh * phiW*(1+chi) + flh * (1+chi));
    fwha = exp(-theta * fch + phiW * fwh + flh);
    exp(fwh * (1-phiW)) = exp(fwhn * (1-phiW));


    % ************************************    ENTREPRENEURS    *************************************
    exp(x) * exp(y) + be = exp(ce) + exp(q) * (exp(he) - exp(he(-1))) + adjHe + exp(wp) * exp(npe) +
        exp(wh) * exp(nhe) + (1+re(-1)) * be(-1) / (1+inf) + rk * exp(ke(-1));
    adjHe = (phiE / 2) * exp(q) * exp(he(-1)) * (exp(he)/exp(he(-1)) - 1)^2;
    y = a + mu*ke(-1) + alpha * (1 - mu - nu) * npe + (1-alpha) * (1 - mu - nu) * nhe + nu * he(-1);
    be = me * exp(q(+1)) * exp(he) * (inf(+1)+1) / (re+1);
    (exp(ce))^(-theta) = betaE * (1+re) * (exp(ce(+1)))^(-theta) / (1+inf(+1)) + Lame * (1+re);
    wp + npe = log(alpha * (1 - mu - nu)) + x + y;
    wh + nhe = log((1 - alpha) * (1 - mu - nu)) + x + y;
    log(rk) + ke(-1) = log(mu) + x + y;
    (exp(ce))^(-theta) * exp(q) * (1 + phiE * (exp(he)/exp(he(-1)) - 1)) - Lame * me * exp(q(+1)) * (1+inf(+1)) =
        betaE * (exp(ce(+1)))^(-theta) * (nu * exp(x(+1)) * exp(y(+1)) / exp(he) + exp(q(+1)) *
        (1 + (phiE / 2) * ((exp(he(+1))/exp(he))^2 -1)));


    % **********************************    FLEX ENTREPRENEURS    **********************************
    exp(fx) * exp(fy) + fbe = exp(fce) + exp(fq) * (exp(fhe) - exp(fhe(-1))) + fadjHe + exp(fwp) * exp(fnpe) +
        exp(fwh) * exp(fnhe) + (1+fre(-1)) * fbe(-1) / (1+finf) + frk * exp(fke(-1));
    fadjHe = (phiE / 2) * exp(fq) * exp(fhe(-1)) * (exp(fhe)/exp(fhe(-1)) - 1)^2;
    fy = a + mu*fke(-1) + alpha * (1 - mu - nu) * fnpe + (1-alpha) * (1 - mu - nu) * fnhe + nu * fhe(-1);
    fbe = me * exp(fq(+1)) * exp(fhe) * (finf(+1)+1) / (fre+1);
    (exp(fce))^(-theta) = betaE * (1+fre) * (exp(fce(+1)))^(-theta) / (1+finf(+1)) + fLame * (1+fre);
    fwp + fnpe = log(alpha * (1 - mu - nu)) + fx + fy;
    fwh + fnhe = log((1 - alpha) * (1 - mu - nu)) + fx + fy;
    log(frk) + fke(-1) = log(mu) + fx + fy;
    (exp(fce))^(-theta) * exp(fq) * (1 + phiE * (exp(fhe)/exp(fhe(-1)) - 1)) - fLame * me * exp(fq(+1)) * (1+finf(+1)) =
        betaE * (exp(fce(+1)))^(-theta) * (nu * exp(fx(+1)) * exp(fy(+1)) / exp(fhe) + exp(fq(+1)) *
        (1 + (phiE / 2) * ((exp(fhe(+1))/exp(fhe))^2 -1)));


    % **************************************    RETAILERS    ***************************************
    markup = log(sigma / (sigma - 1));
    exp(prn) = exp(markup) * pb / pa;
    pb = omegaP * betaP * pb(+1) + exp(x - cp + y);
    pa = omegaP * betaP * pa(+1) + exp(-cp + y);
    1 = (1-omegaP) * exp((1-sigma) * prn) + omegaP * ((1+inf)^(sigma-1));

    % inf = (1 - omegaP) * (1 - omegaP * betaP) * (x + markup) / omegaP + betaP * inf(+1);
    divr = (1 - exp(x)) * exp(y);


    % ************************************    FLEX RETAILERS    ************************************
    fmarkup = markup;
    exp(fprn) = exp(markup) * fpb / fpa;
    fpb = exp(fx - fcp + fy);
    fpa = exp(-fcp + fy);
    1 = exp((1-sigma) * fprn);


    % *************************************    CONSTRUCTION    *************************************
    ih = alpha*(1-eta) * npc + (1-alpha)*(1-eta) * nhc + eta*tl;
    Lamc = 1 / (1+rc) - betaC;

    % divc = exp(q + ih) + bc - exp(wp + npc) - exp(wh + nhc) - (1+rc(-1)) * bc(-1) / (1 + inf);        % No delay
    % bc = mc * exp(q(+1) + ih) * (inf(+1)+1) / (rc+1);
    % exp(wp) = alpha * (1-eta) * exp(ih - npc) * (Lamc*mc*exp(q(+1))*(1+inf(+1)) + q);
    % exp(wh) = (1-alpha) * (1-eta) * exp(ih - nhc) * (Lamc*mc*exp(q(+1))*(1+inf(+1)) + q);

    % divc = exp(q + ih(-1)) + bc - exp(wp + npc) - exp(wh + nhc) - (1+rc(-1)) * bc(-1) / (1 + inf);    % 1-period delay
    % bc = mc * exp(q(+1) + ih) * (inf(+1)+1) / (rc+1);
    % exp(wp) = alpha * (1-eta) * exp(ih - npc) * exp(q(+1)) * (1+inf(+1)) * (Lamc * mc + betaC);
    % exp(wh) = (1-alpha) * (1-eta) * exp(ih - nhc) * exp(q(+1)) * (1+inf(+1)) * (Lamc * mc + betaC);

    divc = exp(q + ih(-4)) + bc - exp(wp + npc) - exp(wh + nhc) - (1+rc(-1)) * bc(-1) / (1 + inf);      % 4-period delay
    bc = mc * exp(q(+1) + ih(-3)) * (inf(+1)+1) / (rc+1);
    exp(wp) = alpha * (1-eta) * (betaC^3) * exp(ih - npc) * exp(q(+4)) * (1+inf(+4)) * (1+inf(+3)) * (1+inf(+2)) * (1+inf(+1)) * (Lamc(-3) * mc(-3) + betaC);
    exp(wh) = (1-alpha) * (1-eta) * (betaC^3) * exp(ih - nhc) * exp(q(+4)) * (1+inf(+4)) * (1+inf(+3)) * (1+inf(+2)) * (1+inf(+1)) * (Lamc(-3) * mc(-3) + betaC);

    % divc = exp(q + ih(-8)) + bc - exp(wp + npc) - exp(wh + nhc) - (1+rc(-1)) * bc(-1) / (1 + inf);    % 8-period delay
    % bc = mc * exp(q(+1) + ih(-7)) * (inf(+1)+1) / (rc+1);
    % exp(wp) = alpha*(1-eta)*(betaC^7)*exp(ih-npc)*exp(q(+8))*(1+inf(+8))*(1+inf(+7))*(1+inf(+6))*(1+inf(+5))*(1+inf(+4))*(1+inf(+3))*(1+inf(+2))*(1+inf(+1))*(Lamc(-7)*mc(-7)+betaC);
    % exp(wh) = (1-alpha)*(1-eta)*(betaC^7)*exp(ih - nhc)*exp(q(+8))*(1+inf(+8))*(1+inf(+7))*(1+inf(+6))*(1+inf(+5))*(1+inf(+4))*(1+inf(+3))*(1+inf(+2))*(1+inf(+1))*(Lamc(-7)*mc(-7)+betaC);
    

    % ***********************************    FLEX CONSTRUCTION    **********************************
    fih = alpha*(1-eta) * fnpc + (1-alpha)*(1-eta) * fnhc + eta*tl;
    fLamc = 1 / (1+frc) - betaC;

    % fbc = mc * exp(fq(+1) + fih) * (finf(+1)+1) / (frc+1);                                                % No delay
    % exp(fwp) = alpha * (1-eta) * exp(fih - fnpc) * (fLamc*mc*exp(fq(+1))*(1+finf(+1)) + fq);
    % exp(fwh) = (1-alpha) * (1-eta) * exp(fih - fnhc) * (fLamc*mc*exp(fq(+1))*(1+finf(+1)) + fq);

    % fbc = mc * exp(fq(+1) + fih) * (finf(+1)+1) / (frc+1);
    % exp(fwp) = alpha * (1-eta) * exp(fih - fnpc) * exp(fq(+1)) * (1+finf(+1)) * (fLamc * mc + betaC);     % 1-period delay
    % exp(fwh) = (1-alpha) * (1-eta) * exp(fih - fnhc) * exp(fq(+1)) * (1+finf(+1)) * (fLamc * mc + betaC);

    fbc = mc * exp(fq(+1) + fih(-3)) * (finf(+1)+1) / (frc+1);                                              % 4-period delay
    exp(fwp) = alpha * (1-eta) * (betaC^3) * exp(fih - fnpc) * exp(fq(+4)) * (1+finf(+4)) * (1+finf(+3)) * (1+finf(+2)) * (1+finf(+1)) * (fLamc(-3) * mc(-3) + betaC);
    exp(fwh) = (1-alpha) * (1-eta) * (betaC^3) * exp(fih - fnhc) * exp(fq(+4)) * (1+finf(+4)) * (1+finf(+3)) * (1+finf(+2)) * (1+finf(+1)) * (fLamc(-3) * mc(-3) + betaC);

    % fbc = mc * exp(fq(+1) + fih(-7)) * (finf(+1)+1) / (frc+1);                                            % 8-period delay
    % exp(fwp) = alpha*(1-eta)*(betaC^7)*exp(fih-fnpc)*exp(fq(+8))*(1+finf(+8))*(1+finf(+7))*(1+finf(+6))*(1+finf(+5))*(1+finf(+4))*(1+finf(+3))*(1+finf(+2))*(1+finf(+1))*(fLamc(-7)*mc(-7)+betaC);
    % exp(fwh) = (1-alpha)*(1-eta)*(betaC^7)*exp(fih - fnhc)*exp(fq(+8))*(1+finf(+8))*(1+finf(+7))*(1+finf(+6))*(1+finf(+5))*(1+finf(+4))*(1+finf(+3))*(1+finf(+2))*(1+finf(+1))*(fLamc(-7)*mc(-7)+betaC);
    

    % ***************************************    BANKING    ****************************************
    bp + f = bh + be + bc;
    f = mb * (bh(+1) + be(+1) + bc(+1)) * (1+inf(+1)) / (1+r);
    rc = rp - mb(-1) * (rp(-1) - r(-1)) / (1+r(-1)) + kappa;
        % rc = rp - mb(+1) * (rp(+1) - r(+1)) / (1+r(+1)) + kappa;
    re = rc;
    rh = rc;
    rgap = rh - rp;
    

    % *************************************    FLEX BANKING    *************************************
    fbp + ff = fbh + fbe + fbc;
    ff = mb * (fbh(+1) + fbe(+1) + fbc(+1)) * (1+finf(+1)) / (1+r);
    frc = frp - mb(-1) * (frp(-1) - r(-1)) / (1+r(-1)) + kappa;
    fre = frc;
    frh = frc;
    frgap = frh - frp;
    

    % **********************************    MARKET CLEARING    *************************************
    exp(cp) + exp(ch) + exp(ce) + exp(ip) = exp(y) + f - (1+r(-1)) * f(-1) / (1 + inf);
    exp(npe) + exp(npc) = exp(lp);
    exp(nhe) + exp(nhc) = exp(lh);
    kp = ke;
    exp(kp) = (1 - delta) * exp(kp(-1)) + exp(ip);

    % No delay in investment
    % exp(h(-1)) + exp(ih) = exp(hp) + exp(hh) + exp(he);
    % exp(h) = (1 - deltaH) * exp(h(-1)) + exp(ih);

    % One period delay in investment
    % exp(h(-1)) = exp(hp) + exp(hh) + exp(he);
    % exp(h) = (1 - deltaH) * exp(h(-1)) + exp(ih);

    % Four period delay in investment
    exp(h(-1)) = exp(hp) + exp(hh) + exp(he);
    exp(h) = (1 - deltaH) * exp(h(-1)) + exp(ih(-3));

    % Eight period delay in investment
    % exp(h(-1)) = exp(hp) + exp(hh) + exp(he);
    % exp(h) = (1 - deltaH) * exp(h(-1)) + exp(ih(-7));


    % ********************************    FLEX MARKET CLEARING    **********************************
    exp(fcp) + exp(fch) + exp(fce) + exp(fip) = exp(fy) + ff - (1+r(-1)) * ff(-1) / (1 + finf);
    exp(fnpe) + exp(fnpc) = exp(flp);
    exp(fnhe) + exp(fnhc) = exp(flh);
    fkp = fke;
    exp(fkp) = (1 - delta) * exp(fkp(-1)) + exp(fip);

    % No delay in investment
    % exp(fh(-1)) + exp(fih) = exp(fhp) + exp(fhh) + exp(fhe);
    % exp(fh) = (1 - deltaH) * exp(fh(-1)) + exp(fih);

    % One period delay in investment
    % exp(fh(-1)) = exp(fhp) + exp(fhh) + exp(fhe);
    % exp(fh) = (1 - deltaH) * exp(fh(-1)) + exp(fih);

    % Four period delay in investment
    exp(fh(-1)) = exp(fhp) + exp(fhh) + exp(fhe);
    exp(fh) = (1 - deltaH) * exp(fh(-1)) + exp(fih(-3));

    % Eight period delay in investment
    % exp(fh(-1)) = exp(fhp) + exp(fhh) + exp(fhe);
    % exp(fh) = (1 - deltaH) * exp(fh(-1)) + exp(fih(-7));


    % **********************************    MONETARY POLICY    *************************************
    % ginf = (1 - rhoQ) * inf + rhoQ * (q - q(-1));
    ginf = inf + rhoQ * (q - q(-1));
    % rp = -(1 - rhoI) * log(betaP) + rhoI * rp(-1) + (1 - rhoI) * rhoPi * ginf + epspi;
    rp = -(1 - rhoI) * log(betaP) + rhoI * rp(-1) + (1 - rhoI) * (rhoPi * ginf + rhoGap * (gdp - fgdp) ) + epspi;
    % rp = -(1 - rhoI) * log(betaP) + rhoI * rp(-1) + (1 - rhoI) * 1.01 * inf(-1) + epspi;


    % ********************************    FLEX MONETARY POLICY    **********************************
    fginf = finf + rhoQ * (fq - fq(-1));
    frp = -(1 - rhoI) * log(betaP) + rhoI * frp(-1) + (1 - rhoI) * rhoPi * fginf + epspi;


    % ********************************    EXOGENOUS PROCESSES    ***********************************
    r = (1 - rhoR) * Rss + rhoR * r(-1) + epsi;
    fr = r;
    a  = rhoA * a(-1) + epsa;
    mh  = mhss * exp(mhvar);
        mhvar = rhoH * mhvar(-1) + epsh;
    me  = mess * exp(mevar);
        mevar = rhoE * mevar(-1) + epse;
    mc  = mcss * exp(mcvar);
        mcvar = rhoC * mcvar(-1) + epsc;
    mb  = mbss * exp(mbvar);
        mbvar = rhoB * mbvar(-1) + epsb;


    % ***********************************    STEADY STATE    **************************************
    zeta = (1 - mu - nu) * rk / mu;
    ksi = (1 - eta) * (Lamc * mcss + betaC) / zeta;
    psai = 1/(1 - betaH - mhss*( (1/(rh+1)) - betaH));
    rho = (1 + (rh/(rh+1)) * mhss * gama * psai)^(1/(1+chi));
    phi = 1/(1/betaE - (1/(betaE*(re+1))-1)*mess - 1);
    tao = 1 / ( 1/delta + ksi*nu*rk*phi/mu );
    stigma = delta + rk * (nu - sigma/(sigma-1)) / mu;
    aux1 = -r * mb / ((1+r) * (1+rh));
    aux2 = (1 - alpha) * gama * zeta / (rho^(1+chi));
    aux3 = alpha * gama * zeta / (1 - betaP);
        A11 = aux1 * (1 - tao * (1/delta + mcss/mess)) - ksi*tao*stigma/mess;
    A1 = ((re/(re+1)) * (1 - tao/delta) + A11) * aux3 * mess + alpha * zeta;
        A21 = (re/(re+1)) - aux1 * (mhss/mess - 1);
    A2 = aux2/gama + stigma + aux2*psai*mess * A21;
        A31 = aux1 * (mess/delta + mcss) / ksi + stigma + (re/(re+1)) * mess / (delta * ksi);
    A3 = ksi * phi * tao * (nu*rk/mu + aux2) * A31;
    z = rho / exp(lp);
    zbig = z^((1-alpha)*(1-eta));


    % ******************************    SOME STEADY STATE RATIOS   ********************************
    gdp = log(exp(qss + ih) + exp(y));
    itot = log(exp(ip) + exp(q + ih));                  % Total investment
    rhy = (exp(q + hp) + exp(q + hh)) / (4 * gdp);      % Ratio of residential housing to GDP
    chy = exp(q + he) / (4 * gdp);                      % Ratio of commercial real estate to GDP
    pcy = (bh + be + bc) / gdp;                         % Ratio of private credit to GDP
    pc = bh + be + bc;
    bcy = (be + bc) / gdp;                              % Ratio of business credit to GDP
    ihy = deltaH * (exp(q+hh) + exp(q+hp)) / gdp;       % Ratio of housing investment to GDP
    aggc = log(exp(cp) + exp(ch) + exp(ce));
    aggs = bp + f;

    fgdp = log(exp(qss + fih) + exp(fy));
    fitot = log(exp(fip) + exp(fq + fih));


end;

initval;

    % Exogenous variables
    mhvar = 0;
    mevar = 0;
    mcvar = 0;
    mbvar = 0;
    mh = mhss;
    me = mess;
    mc = mcss;
    mb = mbss;
    a = 0;
    ginf = 0;
    inf = 0;
    adjHp = 0;
    adjHh = 0;
    adjHe = 0;
    adjK = 0;
    markup = log(sigma / (sigma - 1));
    wm = phiW/(phiW-1);
    x = -markup;

    % Interest rates
    r = Rss;
    rp = -log(betaP);
    rk = 1/betaP - 1 + delta;
    rc = rp - mbss * (rp - r) / (1 + r) + kappa;
    re = rc;
    rh = rc;
    rgap = rh - rp;
    Lamc = 1 / (1+rc) - betaC;

    % Auxiliary variables
    zeta = (1 - mu - nu) * rk / mu;
    ksi = (1 - eta) * (Lamc * mcss + betaC) / zeta;
    psai = 1/(1 - betaH - mhss*( (1/(rh+1)) - betaH));
    rho = (1 + (rh/(rh+1)) * mhss * gama * psai)^(1/(1+chi));
    phi = 1/(1/betaE - (1/(betaE*(re+1))-1)*mess - 1);
    tao = 1 / ( 1/delta + ksi*nu*rk*phi/mu );
    stigma = delta + rk * (nu - sigma/(sigma-1)) / mu;
    aux1 = -r * mb / ((1+r) * (1+rh));
    aux2 = (1 - alpha) * gama * zeta / (rho^(1+chi));
    aux3 = alpha * gama * zeta / (1 - betaP);
        A11 = aux1 * (1 - tao * (1/delta + mcss/mess)) - ksi*tao*stigma/mess;
    A1 = ((re/(re+1)) * (1 - tao/delta) + A11) * aux3 * mess + alpha * zeta;
        A21 = (re/(re+1)) - aux1 * (mhss/mess - 1);
    A2 = aux2/gama + stigma + aux2*psai*mess * A21;
        A31 = aux1 * (mess/delta + mcss) / ksi + stigma + (re/(re+1)) * mess / (delta * ksi);
    A3 = ksi * phi * tao * (nu*rk/mu + aux2) * A31;

    % Rest of the model
    lp = (1/(1+chi)) * log(A1 / (A3 - A2));
    npc = lp + log( ksi*phi*tao*(nu*rk/mu + aux2) + aux3*ksi*tao/exp((chi+1)*lp) );
    he = -eta*npc + (1-alpha)*(1-eta)*log(rho - lp) + log( exp(npc)/delta - aux2*ksi*psai*exp(lp) - aux3*ksi/exp(chi*lp) );
    npe = log( exp(lp) - exp(npc) );
    z = rho / exp(lp);
    zbig = z^((1-alpha)*(1-eta));
    hh = -eta*npc + log((1 - alpha)*gama*zeta*ksi*psai) - (1+chi)*log(z) + log(zbig) - chi*lp;
    q = ( log(mu) + x - (1-mu)*log(ksi) - log(rk) + (1-alpha)*(eta*(1-mu)-nu)*log(z) - nu*(npe-he) + eta*(1-mu)*npc ) / (1-mu);
    Lamh = (1/(1+rh) - betaH) * z^(chi+1) * exp(chi*lp + eta*npc) / (zbig*(1-alpha)*zeta*ksi*exp(q));
    Lame = (1/(1+re) - betaE) * (( ksi*nu*rk*zbig*(exp(q+npe - eta*npc))/mu - (re/(re+1))*mess*exp(q+he) )^(-theta));
    ke = log(zbig) + q + npe - eta*npc + log(ksi);
    h = log(zbig) - log(delta) + (1-eta)*npc;
    hp = log( exp(h) - exp(hh) - exp(he));
    cp = ( log(alpha*zeta) + ke - npe - chi*lp ) / theta;
    ch = ( theta*cp - log(alpha) + log(1-alpha) - (1+chi)*log(z) ) / theta;
    ce = log( nu*rk*exp(ke)/mu - (re/(re+1))*mess*exp(q+he) );
    nhe = log(z) + npe;
    nhc = log(z) + npc;
    y = log(rk/mu) + ke - x;
    ih = alpha*(1-eta)*npc + (1-alpha)*(1-eta)*nhc;
    gdp = log(exp(q + ih) + exp(y));
    wp = log(alpha * zeta) + ke - npe;
    wh = log((1-alpha) * zeta) + ke - nhe;
    f = (mbss * exp(q) / ((1+r)*(1+rh))) * (mhss*exp(hh) + mess*exp(he) + mcss * exp(ih));
    bh = mhss * exp(q+hh) / (1+rh);
    be = mess * exp(q+he) / (1+re);
    divr = (1 - exp(x))*exp(y);
    bc = mcss * exp(q+ih) / (1+rc);
    divc = exp(q + ih) + bc - exp(wp+npc) - exp(wh + nhc) - (rc+1)*bc;
    bp = bh + be + bc - f;
    lh = log( exp(nhe) + exp(nhc));
    kp = ke;
    ip = log(delta) + kp;
    itot = log(exp(ip) + exp(q + ih));
    rhy = (exp(q + hp) + exp(q + hh)) / (4 * exp(y));
    chy = exp(q + he) / (4 * exp(y));
    pcy = (bh + be + bc) / exp(y);
    pc = bh + be + bc;
    bcy = (be + bc) / exp(y);
    ihy = deltaH * (exp(q + hh) + exp(q + hp)) / exp(y);
    aggc = log(exp(cp) + exp(ch) + exp(ce));
    aggs = bp + f;

    vp = (cp - lp^(1+chi) / (1 + chi) + gama * hp) / (1 - betaP);
    vh = (ch - lh^(1+chi) / (1 + chi) + gama * hh) / (1 - betaH);
    
    pb = exp(-cp + x + y) / (1 - omegaP * betaP);
    pa = exp(-cp + y) / (1 - omegaP * betaP);
    prn = 0;
    wpb = exp(wp * phiW*(1+chi) + lp * (1+chi)) / (1 - omegaW * betaP);
    wpa = exp(-theta * cp + wp * phiW + lp) / (1 - omegaW * betaP);
    wpn = wp;
    whb = exp(wh * phiW*(1+chi) + lh * (1+chi)) / (1 - omegaW * betaH);
    wha = exp(-theta * ch + wh * phiW + lh) / (1 - omegaW * betaH);
    whn = wh;

    fr = r;
    frgap = rgap;
    fmarkup = markup;
    fcp = cp;
    fch = ch;
    fce = ce;
    fip = ip;
    fih = ih;
    flp = lp;
    flh = lh;
    fnpe = npe;
    fnpc = npc;
    fnhe = nhe;
    fnhc = nhc;
    fwp = wp;
    fwh = wh;
    fq = q;
    fhp = hp;
    fhh = hh;
    fhe = he;
    ff = f;
    fh = h;
    fbp = bp;
    fbh = bh;
    fbe = be;
    fbc = bc;
    frp = rp;
    frh = rh;
    fre = re;
    frc = rc;
    fy = y;
    fx = x;
    fke = ke;
    fkp = kp;
    frk = rk;
    fprn = prn;
    fwpn = wpn;
    fwhn = whn;
    fpb = pb;
    fpa = pa;
    fwpb = wpb;
    fwpa = wpa;
    fwhb = whb;
    fwha = wha;
    fLamh = Lamh;
    fLame = Lame;
    fLamc = Lamc;
    fginf = ginf;
    finf = inf;
    fgdp = gdp;
    fadjHp = adjHp;
    fadjK = adjK;
    fadjHh = adjHh;
    fadjHe = adjHe;
    fitot = itot;
    

end;

steady;
check;

Sigma_e = [0.000030 0 0 0 0 0 0; 0.000132 0 0 0 0 0; 0.00 0 0 0 0; 0.00 0 0 0; 0.00 0 0; 0.003592 0; 0.000036];
    %     epsa                  epsi                epsh          epse        epsc      epsb        epspi
% Sigma_e = [0.000000 0 0 0 0 0 0; 0.000000 0 0 0 0 0; 0.00 0 0 0 0; 0.00 0 0 0; 0.00 0 0; 0.000000 0; 0.00];
%     %     epsa                  epsi                epsh          epse        epsc      epsb        epspi


% Value function components
% stoch_simul(IRF=0,hp_filter=1600,ORDER=2) vp, vh;

% Behavior of key variables
stoch_simul(IRF=0,hp_filter=1600,ORDER=1) gdp, aggc, ip, ih, q, h, rp, inf, bh;
% stoch_simul(IRF=0,hp_filter=1600,ORDER=1) gdp, ch, ih, q, inf, rp, bh, lh, hh, cp, lp, hp;
% stoch_simul(IRF=20,hp_filter=1600,ORDER=1) gdp, y, npe, npc, nhe, nhc, wp, wh;
% stoch_simul(IRF=20,hp_filter=1600,ORDER=1) gdp, rgap, bp, bh, be, bc, f, rh;
% stoch_simul(IRF=20,hp_filter=1600,ORDER=1) Lamh, Lame, Lamc, boundh;
% stoch_simul(IRF=20,hp_filter=1600,ORDER=1) y, fy, cp, fcp, q, fq;