// Ai, Croce, Li (2011)
// Toward a Quantitative General Equilibrium Asset Pricing Model
// with Intangible Capital
//
// Prepared for Dynare by Jaroslav Borovicka and Lars Peter Hansen as an
// illustration of the methods developed in "Examining Macroeconomic Models
// through the Lens of Asset Pricing". This Dynare model file should be
// used in conjuction with the elasticity toolbox developed for Dynare.
// See "Shock Elasticities in Dynare" by the above two authors,
// downloadable from
// http://home.uchicago.edu/~borovicka/software.html
//
// We thank Hengjie Ai, Max Croce and Kai Li for providing us with details
// of their model. All errors are our own.
//
// This model is the extended version of the model, with adjustment costs
// in the accumulation of intangible capital.

// ========================================================================
// definition of variables
// households
var v ev c;
// producers
var k s i j m h hj;
// prices
var qex qexstar qb q qstar rk rkstar rf rkrf rkstarrf;

// auxiliary variables
var imins mmins jmins mmini;

// exogenous processes
var x x1 x2 da omega;

// exogenous shocks, together with the perturbations of the model used for
// the calculation of partial shock elasticities
// eps_da1 and eps_x1 correspond to the neutral technology channel
// eps_da2 and eps_x2 correspond to the investment-specific channel
varexo eps_da eps_x eps_da1 eps_x1 eps_da2 eps_x2;

// INCREMENTS IN ADDITIVE FUNCTIONAL FOR THE ELASTICITY TOOLBOX
var sdf_household;
var cf_k cf_kstar cf_q cf_qstar;
var cf_rkrf_cum cf_rkstarrf_cum cf_consumption;

// ========================================================================
// definition of parameters

parameters BETA RHO GAMMA;
parameters NU LAMBDA LAMBDASTAR PHI ETA;
parameters MU_da SIGMA_da RHO_x SIGMA_x;
parameters A1 A2 XI;
// auxiliary parameters
parameters QB_BAR IS_BAR MS_BAR JS_BAR MI_BAR QK_BAR K_BAR M_BAR S_BAR JK_BAR;

parameters SCALE;

// preferences (original parameter values following)
BETA = 0.971;         // discount rate
GAMMA = 10;           // GAMMA = RA in the recursive utility
RHO = 0.5;            // RHO^(-1) = IES in the recursive utility

// firms
NU = 0.30;            // capital share
LAMBDA = 0.11;        // depreciation of physical capital
LAMBDASTAR = 0.11;    // depreciation of intangible capital
PHI = 0.88;           // weight in G(I,S)
ETA = 2.5;            // elasticity of substitution in G(I,S)

// parameters for the exogeneous processes
MU_da = 0.02;         // trend growth rate
SIGMA_da = 0.0508;    // volatility of the permament shock 
RHO_x = 0.925;        // persistence of the long-run risk shock
SIGMA_x = 0.008636;   // volatility of the long-run risk shock

// adjustment cost on intangible capital accumulation
QB_BAR = log(BETA) - RHO*(MU_da);
IS_BAR = ETA*log((exp(-QB_BAR)-1+LAMBDASTAR)*PHI/(1-PHI));
MS_BAR = 1/(1-1/ETA)*log(PHI*exp((1-1/ETA)*IS_BAR) + (1-PHI));
JS_BAR = log(exp(MU_da) - (1-exp(MS_BAR))*(1-LAMBDASTAR));
MI_BAR = 1/(1-1/ETA)*log(PHI + (1-PHI)*exp((1-1/ETA)*(-IS_BAR)));
QK_BAR = log(1 - LAMBDASTAR + 1/PHI*exp(-(1/ETA)*MI_BAR));
K_BAR = 1/(NU-1)*log((exp(QK_BAR)*(1/exp(QB_BAR)-1+LAMBDA))/NU);
M_BAR = K_BAR + log(exp(MU_da) - 1 + LAMBDA);
S_BAR = M_BAR - MS_BAR;
JK_BAR = JS_BAR + S_BAR - K_BAR;

XI = 5;
A1 = exp(JK_BAR)^(1/XI);
A2 = -(1/(XI-1))*exp(JK_BAR);

// scaling of the small orthogonal perturbations
SCALE = 0.001;

// ========================================================================
model;

    // household
    // recursive utility continuation value
    // GAMMA = RA, RHO^(-1) = IES
    exp((1-GAMMA)*ev) = exp((1-GAMMA)*(v(+1)+da(+1)));
    exp(v) = ( (1-BETA)*exp((1-RHO)*c) + BETA*exp((1-RHO)*ev) )^(1/(1-RHO));
    sdf_household = log(BETA) - RHO*(c-c(-1)+da) + (RHO-GAMMA)*(v+da - ev(-1));

    // aggregate feasibility
    exp(c) = exp(NU*k) - exp(i) - exp(j);

    // laws of motion
    exp(k+da) = exp(k(-1))*(1-LAMBDA) + exp(omega+m(-1));
    exp(m) = ( PHI*exp((1-1/ETA)*i) + (1-PHI)*exp((1-1/ETA)*s))^(1/(1-1/ETA));   // G(I,KSTAR)
    exp(h) = exp(k)*(A2 + A1/(1-1/XI)*exp((1-1/XI)*(j-k)));                      // H(ISTAR,K)
    exp(hj) = A1*exp(-1/XI*(j-k));                                               // H_ISTAR(ISTAR,K)
    exp(s+da) = (exp(s(-1)) - exp(m(-1)))*(1-LAMBDASTAR) + exp(h(-1));

    // Euler equations
    // q and qstar correspond to pk and ps in Ai, Croce and Li (2010)
    // qex and qexstar correspond to qk and qs in Ai, Croce and Li (2010)
    1/PHI*exp((1/ETA)*(i-m)) + (1-LAMBDASTAR)*exp(qex) = exp(sdf_household(+1)+omega(+1)+q(+1));
    exp(qex) = 1/exp(hj);
    exp(qexstar) = exp(sdf_household(+1)+q(+1));
    exp(qex) = exp(sdf_household(+1)+qstar(+1));
    exp(q) = NU*exp((NU-1)*k) + exp(qex)*(exp(h-k)-exp(hj+j-k)) + (1-LAMBDA)*exp(qexstar);
    exp(qstar) = ((1-PHI)/PHI*exp(1/ETA*(i-s)) + (1-LAMBDASTAR)*exp(qex));

    // pricing
    exp(qb) = exp(sdf_household(+1));
    rk = q - qexstar(-1);
    rkstar = qstar - qex(-1);
    rf = - qb;
    rkrf = rk - rf(-1);
    rkstarrf = rkstar - rf(-1);
    
    // exogenous processes, including small perturbations
    x = RHO_x * x(-1) + SIGMA_x * eps_x;
    x1 = RHO_x * x1(-1) + SIGMA_x * eps_x1;
    x2 = RHO_x * x2(-1) + SIGMA_x * eps_x2;

    da = MU_da + x(-1) + SIGMA_da * eps_da + SCALE*(x1(-1) + SIGMA_da * eps_da1);

    omega = -(1-NU)/NU*(x(-1) + SIGMA_da*eps_da + SCALE*(x2(-1) + SIGMA_da * eps_da2));

    // auxiliary variables
    imins = i - s;
    mmins = m - s;
    jmins = j - s;
    mmini = m - i;

    // INCREMENTS IN ADDITIVE FUNCTIONAL FOR THE ELASTICITY TOOLBOX
    cf_k = k - k(-1) + da;            // physical capital stock
    cf_kstar = s - s(-1) + da;        // intangible capital stock
    cf_q = q - q(-1);                 // price of physical capital
    cf_qstar = qstar - qstar(-1);     // price of intangible capital
    cf_rkrf_cum = rkrf;               // cumulative excess return on physical capital
    cf_rkstarrf_cum = rkstarrf;           // cumulative excess return on physical capital
    cf_consumption = c - c(-1) + da;  // consumption

end;

// ========================================================================
// initial values for steady state calculation

initval;

    x = 0;
    x1 = 0;
    x2 = 0;
    da = MU_da;
    omega = 0;

    sdf_household = log(BETA) - RHO*(da);
    qb = sdf_household;

    imins = ETA*log((exp(-sdf_household)-1+LAMBDASTAR)*PHI/(1-PHI));
    mmins = 1/(1-1/ETA)*log(PHI*exp((1-1/ETA)*imins) + (1-PHI));
    jmins = log(exp(da) - (1-exp(mmins))*(1-LAMBDASTAR));
    mmini = 1/(1-1/ETA)*log(PHI + (1-PHI)*exp((1-1/ETA)*(-imins)));

    qexstar = log(1 - LAMBDASTAR + 1/PHI*exp(-(1/ETA)*mmini));
    k = 1/(NU-1)*log((exp(qexstar)*(1/exp(sdf_household)-1+LAMBDA))/NU);
    q = log(NU*exp((NU-1)*k) + (1-LAMBDA)*exp(qexstar));
    m = k + log(exp(da) - 1 + LAMBDA);
    s = m - mmins;
    i = imins + s;
    j = jmins + s;
    c = log(exp(NU*k) - exp(i) - exp(j));

    qex = 0;    
    h = k + log(A2 + A1/(1-1/XI)*exp((1-1/XI)*(j-k)));
    hj = log(A1) - 1/XI*(j-k);

    v = 1/(1-RHO)*log( ((1-BETA)*(exp((1-RHO)*c)))/(1-BETA*exp((1-RHO)*da)) );
    ev = v + da;
    
    rk = q - qexstar;
    qstar = log((1-PHI)/PHI*exp(1/ETA*(i-s)) + 1 - LAMBDASTAR);
    rkstar = qstar;
    rf = - qb;
    rkrf = rk - rf;
    rkstarrf = rkstar - rf;

    // INCREMENTS IN ADDITIVE FUNCTIONAL FOR THE ELASTICITY TOOLBOX
    cf_k = da;
    cf_kstar = da;
    cf_q = 0;
    cf_qstar = 0;
    cf_rkrf_cum = rkrf;
    cf_rkstarrf_cum = rkstarrf;
    cf_consumption = da;

end;

// ========================================================================

shocks;
    var eps_da;
    stderr 1;
    var eps_x;
    stderr 1;
    var eps_da1;
    stderr 1;
    var eps_x1;
    stderr 1;
    var eps_da2;
    stderr 1;
    var eps_x2;
    stderr 1;
end;

options_.pruning = 0;

stoch_simul(nograph,noprint,order=2);