Handling Portfolio Choices: Two-asset HANK (Auclert et al, 2021)

In this example, we solve the two-asset HANK model as presented in Auclert, Bardóczy, Rognlie, and Straub (2021, Appendix B.3 and Appendix E.1). In this model, each heterogeneous household selects an optimal portfolio consisting of liquid and illiquid assets, taking into account a convex portfolio adjustment cost. Additionally, New Keynesian sticky-price firms determine intermediate goods prices, considering Rotemberg-type price adjustment costs, and make dynamic investment decisions with quadratic capital adjustment costs. Labor unions control wages to maximize the average welfare of their members.

This example is employed to showcase the toolbox’s capability in solving complex decision problems with multi-dimensional endogenous state and control variables, as well as a large equilibrium system of equations.

We adhere to the notation, parametrization, and calibration procedure outlined in their original paper. For the detailed derivation of the equilibrium system of equations, we direct readers to the original paper.

The Model

Households

Households are heterogeneous in their Markov labor efficiency ($e$), and the initial holding liquid assets ($b$) and illiquid assets ($a$). They solve the Bellman equation:

\[\begin{split} V_t\left(e_{it}, b_{it-1}, a_{it-1}\right) =\max_{c_{it}, b_{it}, a_{it}}\left\{\frac{c_{it}^{1-\sigma}}{1-\sigma} - \varphi \frac{N_{t}^{1+\nu}}{1+\nu} +\beta \mathbb{E}_t V_{t+1}\left(e_{it+1}, b_{it}, a_{it}\right)\right\} \\ \quad s.t. \quad c_{it} + a_{it} + b_{it} = \left(1-\tau_{t}\right)w_t e_{it} N_{t} + \left(1+r_{t}^{a}\right) a_{it-1} \\ +\left(1+r_t^b\right) b_{it-1}-\Phi_{t}\left(a_{it}, a_{it-1}\right) \\ a_{it} \geq 0, \quad b_{it} \geq \underline{b}, \end{split}\]

where $\tau_t$ denotes the flat tax rate, $w_t$ the real wage rate, $N_t$ the labor hours, $r_t^b$ and $r_t^{a}$ the real interest rates for liquid assets and illiquid assets, and $\Phi_t$ the adjust cost function specified as $$ \Phi_{t}\left(a_{it}, a_{it-1}\right)=\frac{\chi_1}{\chi_2}\left|\frac{a_{it}-\left(1+r_t^a\right) a_{it-1}}{\left(1+r_t^a\right) a_{it-1}+\chi_0}\right|^{\chi_2}\left[\left(1+r_t^a\right) a_{it-1}+\chi_0\right], $$ with $\chi_0 > 0$, $\chi_1 > 0$, and $\chi_2 > 1$.

Households choose consumption ($c_{it}$) and a portfolio $\left(b_{it}, a_{it}\right)$, subject to a borrowing constraint $a_{it} \geq 0$ and $b_{it} \geq \underline{b}$ , and supply a uniform labor ($N_{t}$) determined by the Labor Unions to accommodate labor demand from the firm. Without loss of generality, the mean of the labor efficiency is normalized to equal to $1$. Hence $N$ represents both average labor hours and the total efficiency unit of the labor supply.

Denote the measure that represents the distribution over households states by $\Gamma_t$, and policy functions by $g_{c,t},g_{b,t},g_{a,t}$.

The Financial Intermediary

Households deposit both liquid and illiquid asset holdings into a representative financial intermediary. The financial intermediary invests the liquid deposit ($\int{b_{it}di}$) to illiquid government bonds at a per unit liquidity transformation cost $\omega$, so that $r_t^b = r_t - \omega$ with $r_t$ denoting realized real interest rate of government bonds at $t$. The financial intermediary also invests illiquid deposit ($\int{a_{it}di}$) into a mutual fund consisting of government bonds and firm equity which implies

\[ 1 + r_t^a = \Theta_{pt-1} \frac{D_t + p_t}{p_{t-1}} + (1 - \Theta_{pt-1}) \left(1 + r_t\right), \]

where the mutual fund equity share $\Theta_{pt-1} = \frac{p_{t-1}}{p_{t-1} + B^{g}_{t-1} - B^{h}_{t-1}}$ is determined by equity share price ($p_{t-1}$), aggregate government bond $\left(B^{g}_{t-1}\right)$, and aggregate liquid asset holdings $\left(B^{h}_{t-1} = \int{b_{it-1}di}\right)$ at $t-1$. In a perfect foresight equilibrium, the no arbitrage condition implies the equalization of future returns for $t \geq 1$

\[ 1 + r_t^a = \frac{D_t + p_t}{p_{t-1}} = 1 + r_t, \forall t \geq 1, \]

However, the returns of assets at $t = 0$ may differ due to surprise shocks at $t = 0$.

Firms

Intermediate goods firms operate a Cobb-Douglas production technology with capital and labor to produce differentiated goods, which are aggregated by competitive final goods firms using a CES aggregator with an elasticity of substitution $\mu_{p}/\left(\mu_p - 1\right)$. Intermediate goods firms engage in monopolistic competition on the intermediate goods market, make dynamic investment decision subject to a quadratic adjustment cost, and face price adjustment cost of a quadratic form in logs.

The demand, technology and adjustment costs for an intermediate goods firm $j$ are defined by the following equations:

Product Demand Function:

\[ y_{jt} = \left(\frac{p_{jt}}{P_{t}}\right)^{-\frac{\mu_{p}}{\mu_p - 1}} Y_{t} \]

Production Function:

\[ y_{jt} = Z_t k_{jt-1}^\alpha n_{jt}^{1-\alpha} \]

Capital Adjustment Cost Function:

\[ \zeta \left(k_{jt},k_{jt-1}\right) k_{jt-1} = \frac{1}{2 \delta \epsilon_I}\left(\frac{k_{jt}}{k_{jt-1}} - 1\right)^2 k_{jt-1} \]

Price Adjustment Cost Function:

\[ \psi^p_t \left(p_{jt},p_{jt-1}\right)=\frac{\mu_p}{\mu_p-1} \frac{1}{2 \kappa_p}\left[\log \left(\frac{p_{jt}}{p_{jt-1}}\right)\right]^2 Y_t \]

Capital Accumulation Equation:

\[ i_{jt} = k_{jt} - \left(1 - \delta\right) k_{jt-1} + \zeta \left(k_{jt},k_{jt-1}\right) k_{jt-1} \]

Dividend Equation:

\[ d_{jt} = y_{jt} - w_{t} n_{jt} - i_{jt} - \psi^p_t \]

The firm $j$ maximizes the discounted value of the dividend flow. In a symmetric equilibrium, firms’ optimal solution can be characterized by the following equations:

The Price Phillips Curve:

\[ \log(1+\pi_t) = \kappa_p \left(mc_t - \frac{1}{\mu_p} \right) + \frac{1}{1+r_{t+1}} \frac{Y_{t+1}}{Y_t} \log(1+\pi_{t+1}) \]

Investment Equation:

\[ Q_t = 1 + \frac{1}{\delta \epsilon_I}\left(\frac{K_t-K_{t-1}}{K_{t-1}}\right) \]

Valuation Equation

\[ (1+r_{t+1})Q_{t} = \alpha \frac{Y_{t+1}}{K_t} mc_{t+1} - \left[\frac{K_{t+1}}{K_t} - (1-\delta) + \frac{1}{2\delta \epsilon_I}\left(\frac{K_{t+1} - K_t}{K_t}\right)^2\right] + \frac{K_{t+1}}{K_t}Q_{t+1} \]

where $mc_t$ is a short-hand notation for

Marginal Cost Function:

\[ mc_t = \frac{w_t }{(1-\alpha)\frac{Y_t}{N_t}} \]

Labor Unions

A competitive labor recruiter aggregates differentiated labor provided by a continuum of monopolistically competitive labor unions. Each labor union $k$ faces a labor demand function with an elasticity of substitution $\mu_{w}/\left(\mu_{w} - 1\right)$, and sets nominal wage rates ($W_{kt}$) to maximize the average utility of its members, subjective to a quadratic adjustment cost in utils

\[ \psi_t^w \left(W_{kt},W_{kt-1}\right)=\frac{\mu_w}{\mu_w - 1} \frac{1}{2 \kappa_w}\left[\log \left(\frac{W_{kt}}{W_{kt-1}}\right)\right]^2 \]

In a symmetric equilibrium, union’s optimal solution leads to the Wage Phillips Curve:

\[ \log \left(1+\pi_t^w\right)=\kappa_w\left[\varphi N_t^{1+v}-\frac{\left(1-\tau_t\right) w_t N_t}{\mu_w} \int{e_{it}c_{it}^{-\sigma}di}\right]+\beta \log \left(1+\pi_{t+1}^w\right), \]

where $1 + \pi_t^w = W_t /W_{t-1} = (1 + \pi_t) w_t/w_{t-1}$ is the wage inflation rate at $t$.

Monetary and Fiscal Policies

The monetary authority implements a Taylor rule with a zero steady state target inflation rate:

\[ i_t = r_{\ast} + \phi \pi_t + \phi_y \left(Y_t - Y_{\ast}\right) + m_t, \]

where $m_t$ is a monetary policy shock, and $r_{\ast}$ and $Y_{\ast}$ are steady state real interest rate and GDP, respectively.

The nominal, real interest rates and inflation rate are linked through the Fisher equation:

\[ 1+r_t=\frac{1+i_{t-1}}{1+\pi_{t}}. \]

Government runs a balanced budget to maintain a constant debt in real value ($B^g$), and sets taxes $\tau_t$ to finance exogenous purchases ($G_t$) and interest payments on bonds:

\[ \tau_t =\frac{r_{t}B^{g} + G_t}{w_t N_t}. \]

Market clearing conditions

Asset market clearing:

\[ p_{t} + B^{g}=\int{\left(a_{it} + b_{it}\right)di} \]

goods market clearing:

\[ Y_t = \int{c_{it}di} + \int{\Phi_{t}\left(a_{it}, a_{it-1}\right)di} + \omega \int{b_{it}di} + I_t + \psi^p \left(P_{t},P_{t-1}\right) + G_t. \]

Notice that we have implicitly imposed labor market clearing conditions by using $N_t$ to denote both labor demand in the firms’ problem and labor supply in the Wage Phillips Curve.

The Equilibrium System of Equations: Transition Path

Given an initial distribution over households states, $\Gamma_0$, a perfect foresight competitive equilibrium is a sequence of (1) distributions over households states $\{\Gamma_t\}_{t=0}^{\infty}$, (2) households’ value and policy functions $\{V_t,g_{c,t},g_{b,t},g_{a,t}\}_{t=0}^{\infty}$, and (3) aggregate quantities and prices $\{Y_t,K_t,N_t,I_t,\psi^{p}_{t},D_t,B_t^h,w_t,mc_t,p_t,Q_t,\pi_t,\pi_{t}^{w},i_t,r_t,r_t^a,r_t^b,\tau_t\}_{t=0}^{\infty}$ such that

\[\begin{split} Y_t = Z_t K_{t-1}^\alpha N_{t}^{1-\alpha} \\ I_t = K_{t} - \left(1 - \delta\right) K_{t-1} + \zeta \left(K_{t},K_{t-1}\right) K_{t-1} \\ \psi^p_t =\frac{\mu_p}{\mu_p-1} \frac{1}{2 \kappa_p}\left[\log \left(1 + \pi_t\right)\right]^2 Y_t \\ D_t = Y_t - w_t N_t - I_t - \psi^{p}_{t} \\ w_t = (1-\alpha)\frac{Y_t}{N_t} mc_t \\ \log(1+\pi_t) = \kappa_p \left(mc_t - \frac{1}{\mu_p} \right) + \frac{1}{1+r_{t+1}} \frac{Y_{t+1}}{Y_t} \log(1+\pi_{t+1}) \\ Q_t = 1 + \frac{1}{\delta \epsilon_I}\left(\frac{K_t-K_{t-1}}{K_{t-1}}\right) \\ (1+r_{t+1})Q_{t} = \alpha \frac{Y_{t+1}}{K_t} mc_{t+1} - \left[\frac{K_{t+1}}{K_t} - (1-\delta) + \frac{1}{2\delta \epsilon_I}\left(\frac{K_{t+1} - K_t}{K_t}\right)^2\right] + \frac{K_{t+1}}{K_t}Q_{t+1} \\ 1 + \pi_t^w = (1 + \pi_t) w_t/w_{t-1} \\ \log \left(1+\pi_t^w\right)=\kappa_w\left[\varphi N_t^{1+v}-\frac{\left(1-\tau_t\right) w_t N_t}{\mu_w} \int{e g_{c,t}^{-\sigma}d \Gamma_t}\right]+\beta \log \left(1+\pi_{t+1}^w\right) \\ i_t = r_{\ast} + \phi \pi_t + \phi_y \left(Y_t - Y_{\ast}\right) + m_t \\ 1+r_t=\frac{1+i_{t-1}}{1+\pi_{t}} \\ r_t^b = r_t - \omega \\ \frac{D_{t+1} + p_{t+1}}{p_{t}} = 1 + r_{t+1} \\ \Theta_{pt-1} = \frac{p_{t-1}}{p_{t-1} + B^{g} - B^{h}_{t-1}} \\ 1 + r_t^a = \Theta_{pt-1} \frac{D_t + p_t}{p_{t-1}} + (1 - \Theta_{pt-1}) \left(1 + r_t\right) \\ \tau_t =\frac{r_{t}B^{g} + G_t}{w_t N_t} \\ B^{h}_{t} = \int{g_{b,t} d \Gamma_t} \\ p_{t} + B^{g}=\int{\left(g_{a,t} + g_{b,t}\right)d \Gamma_t} \end{split}\]

and goods market clearing implied by Walras’s law. The system can be further simplified. See the hmod file below that defines the simplified system of equations.

The Equilibrium System of Equations: Stationary Equilibrium

In a stationary equilibrium with a constant $\left(Z_t, \pi_t, G_t\right) = \left(z_{\ast}, 0, G_{\ast}\right) $, the equilibrium system of equations, in the order as presented above, reduce to

\[\begin{split} Y_{\ast} = z_{\ast} K_{\ast}^\alpha N_{\ast}^{1-\alpha} \\ I_{\ast} = \delta K_{\ast} \\ \psi^p_{\ast} = 0 \\ D_{\ast} = Y_{\ast} - w_{\ast} N_{\ast} - \delta K_{\ast} \\ w_{\ast} = (1-\alpha)\frac{Y_{\ast}}{N_{\ast}} mc_{\ast} \\ mc_{\ast} = \frac{1}{\mu_p} \\ Q_{\ast} = 1 \\ r_{\ast} + \delta = \alpha \frac{Y_{\ast}}{K_{\ast}} mc_{\ast} \\ \pi_{\ast}^w = 0 \\ \varphi N_{\ast}^{v} = \frac{\left(1-\tau_{\ast}\right) w_{\ast}}{\mu_w} \int{e g_{c,{\ast}}^{-\sigma}d \Gamma_{\ast}} \\ i_{\ast} = r_{\ast} \\ i_{\ast} = r_{\ast} \\ r_{\ast}^b = r_{\ast} - \omega \\ D_{\ast} = r_{\ast} p_{\ast} \\ \Theta_{p{\ast}} = \frac{p_{\ast}}{p_{\ast} + B^{g} - B^{h}_{\ast}} \\ r_{\ast}^a = r_{\ast} \\ \tau_{\ast} =\frac{r_{\ast}B^{g} + G_{\ast}}{w_{\ast} N_{\ast}} \\ B^{h}_{\ast} = \int{g_{b,{\ast}} d \Gamma_{\ast}} \\ p_{\ast} + B^{g}=\int{\left(g_{a,{\ast}} + g_{b,{\ast}}\right)d \Gamma_{\ast}} \end{split}\]

We will use these equatios for both calibration and the computation of stationary equilibria.

The hmod File

The model can be represented using an hmod file listed below. To maintain consistency, the hmod file predominantly employs the same variable notations and equation representations as those in the model presentation. The file should be self-explanatory, aided by comments within the code.

% The fixed (non-calibrated) parameter values are from Auclert, Bardóczy, Rognlie, and Straub (2021, Appendix B.3)

parameters beta sigma varphi nu chi0 chi1 omega ra r w tau N;

beta = 0.976; % Calibrated Parameter: HANS final calibrated value is 0.977599682193789 
sigma = 2.0;  % Fixed Parameter
% Note: a typo of varphi in Table B.III: varphi = 2.073 in the table, 
% but the actual varphi value in the SSJ code is around 1.713  
varphi = 1.713; % Calibrated Parameter: HANS final calibrated value is 1.681157908161363
nu = 1.0;   % Fixed Parameter

chi0 = 0.25; % Fixed Parameter
chi1 = 6.416; % Calibrated Parameter: HANS final calibrated value is 6.518415164643450
omega = 0.005; % Fixed Parameter

var_shock e;
shock_trans = [
   0.966289   0.033422   0.000289
   0.016711   0.966578   0.016711
   0.000289   0.033422   0.966289
   ];  % need this to be the full transition matrix

e = [0.18315644, 0.67277917, 2.47128522];

% State Variable grid
var_state b a;
% For replication, the following grids were used in SSJ
% b = 10.^(linspace(log10(0.25), log10(50+0.25), 50)) - 0.25;
% a = 10.^(linspace(log10(0.25), log10(4000+0.25), 70)) - 0.25;
b = 10.^(linspace(log10(0.25), log10(40+0.25), 50)) - 0.25;
a = 10.^(linspace(log10(0.25), log10(120+0.25), 70)) - 0.25;

% Precompute cash on hand, as it does not change over vfi
var_pre_vfi coh;
coh = (1+ra)*a + (1+r-omega)*b + (1-tau)*w*e*N;

% Declare three choice variables in vfi and provide initial guess
var_policy bp ap c;
initial bp 0.0;
initial ap (1+ra)*a;
initial c (1+r-omega)*b + (1-tau)*w*e*N;

% Compute additional variables from vfi for later use: 
% (a) chi is used later in Goods Market Clearing Condition
% (b) uce is used later in the Wage Phillips Curve
var_aux chi uce;

vfi;
    % Define constraint equations
    chi = (chi1/2)*(ap-(1+ra)*a)^2/((1+ra)*a+chi0);    
    %c + ap + bp == (1+ra)*a + (1+r-omega)*b + (1-tau)*w*e*N - chi;
    c + ap + bp == coh - chi;

    % Define the Bellman Equation: note two endogenous state variables (bp,ap)
    u = c^(1.0 - sigma)/(1.0 - sigma) - varphi*N^(1.0+nu)/(1.0+nu);
    Tv = u + beta*EXPECT(v(bp,ap));
    
    % Define Bound Constraints
    bp >= 0.0;
    ap >= 0.0;
    c  >= 1e-8;

    % post evaluation after value function iteration
    uce = e*(c^(-sigma)); % efficiency unit weighted marginal utility used in Wage Phillips Curve
end;

% Endogenous Aggregate Variables
var_agg r pii Y K N w p d Bh;

% Fixed Parameters
delta = 0.02; % depreciation
epsI = 4.0; % capital adjustment cost parameter
kappap = 0.1; % price adjustment cost parameter
muw = 1.1; % wage market power parameter
kappaw = 0.1; % wage adjustment cost parameter

% Calibration Targets: Six Variables
r = 0.0125; % calibration target
Y = 1.0;   % calibration target
K = 10.0;  % calibration target
N = 1.00;  % calibration target
p = 11.2;  % calibration target
Bh = 1.04; % calibration target

% Parameters implied by Calibration Targets
d = 0.14;  % d = p*r = 11.2*0.0125 = 0.14
ra = r; % steady state no-arbitrage condition
Bg = 2.8; % Debt/GDP ratio
G = 0.2; % G/GDP ratio

% Calibrated Parameters: 3 preference parameters above and 3 production parameters below
Z = 0.468; % HANKS Final Calibrated Parameters: Z = 0.467789814531232 
alpha = 0.33; % HANKS Final Calibrated Parameters: alpha = 0.329949238578680 
mup = 1.015; % HANKS Final Calibrated Parameters: mup = 1.015228426395939;

% Parameters implied by Calibration Targets and Calibrated Parameters
mc = 0.9850; % mu = 1/mup
w = 0.66; % w = (1-alpha)*Y/N * mc

% Taylor Rule Parameters
pii = 0.0000;  % inflation
phi = 1.5;  % Taylor rule
phi_y = 0.0; % Taylor rule
pi_star = 0.0; % Taylor rule
rstar = 0.0125; % Taylor rule
Y_star = 1.0; % Taylor rule

var_agg_shock m_shock;
m_shock = 0.0; % placeholder for monetary shocks

% Calibration Block: Six Calibration Targets (r,Y,K,N,p,Bh) with Six Parameters (beta,varphi,chi1,Z,alpha,mup),
% among which four parameters (varphi,Z,alpha,mup)  can be solved analytically.
% Hence, calibration involves an equation solver with two variables and two equilibrium conditions

model_cali(beta, chi1);
    % The calibration of mc follows from the following analytical conditions:
    % d = rp and d = Y- w N - delta K implies
    % rp = Y- w N - delta K so that r (p - K) = Y- w N - ( r + delta) K
    % Now use ( r + delta) = mc*alpha*Y/K and w = mc * (1 - alpha) * Y / N to get
    % r (p - K) = Y - mc Y so that mc = 1 - r (p - K) / Y;
    mc = 1 - r * (p - K) / Y;
    
    mup = 1 / mc; % steady state Price Phillips Curve
    alpha = (r + delta) * K / Y / mc; % stead state Valuation Equation: r + delta = mc*alpha*Y/K
    Z = Y/(K^alpha*N^(1-alpha)); % production function: Y = Z K^alpha N^(1-alpha)

    w = mc * (1 - alpha) * Y / N; % marginal cost function: mc = w/((1-alpha)*Y/N)
    ra = r; % no-arbitrage condition
    tau = (r*Bg + G)/(w*N); % government budget

    % equilibrium conditions
    ap + bp == p + Bg; % asset demand = asset supply
    Bh == bp; % liquid asset supply = liquid asset demand

    % Bound Constraints for Calibrated Parameters
    beta <= 1 - 1e-10;
    beta >= 1e-10;
    chi1 >= 1e-10;
    
    % post evaluation
    varphi =((1-tau)*w*uce/muw)/(N.^nu); % Wage Phillips Curve in Steady State
    Y_star = Y; % update Y_star in Taylor Rule
end;

% Stationary Equilibrium: An equation solver with two variables and two equilibrium conditions

model_ss(r, Y);
    ra = r; % no-arbitrage condition
    mc = 1/mup; % steady state Price Phillips Curve
    % Infer K from the stead state Valuation Equation: r + delta = mc*alpha*Y/K;
    K = mc*alpha*Y/(r + delta);
    % Infer N from the production function: Y = z K^alpha N^(1-alpha)
    N = (Y/Z/(K^alpha))^(1/(1-alpha));
    % Infer w from the marginal cost function: mc = w/((1-alpha)*Y/N)
    w = mc*(1-alpha)*Y/N;
    % Infer d from I = delta*K and zero steady state price adjustment cost
    d = Y - w*N - delta*K;
    % Infer p from asset pricing equation 
    p = d/r;
    % Infer tau from government budget
    tau = (r*Bg + G)/(w*N);

    % equilibrium conditions
    ap + bp == p + Bg; % asset demand = asset supply
    N == ((1-tau)*w*uce/(varphi*muw)).^(1/nu); % Wage Phillips Curve

    % Bound Constraints for Equilibrium Variables
    r <= (1/beta - 1  - 1e-10);
    Y >= 0.5;

    % post evaluation
    rstar = r; % update rstar in Taylor Rule
    Y_star = Y; % update Y_star in Taylor Rule
    pii = pi_star; % exogenously fixed steady state inflation
    Bh = bp;    % liquid asset 
end;

% Transition Path Equilibrium Block: Nine Equations with Nine Unknowns
% The Equation System can be further simplified by utilizing the analytical structure. 
% We keep it in the current form for ease of presentation.

model;
    ii = rstar + phi*pii(-1) + phi_y*(Y(-1) - Y_star) + m_shock(-1); % Taylor Rule: recall ii was determined from last period

    pshare = p(-1)/(p(-1)+Bg-Bh(-1)); % equity share of mutual fund portfolio
    ra = pshare*(d+p)/p(-1) + (1- pshare)*(1+r) - 1; % the return of illiquid asset

    I = K - (1-delta)*K(-1) + ((K/K(-1) - 1).^2)*K(-1)/(2*delta*epsI); % capital accumulation equation
    psip = ((log(1+pii))^2)*Y*mup/((mup - 1)*2*kappap); % price adjustment cost
    Q = 1 + (K/K(-1) - 1)/(delta*epsI); % Tobin's Q from the investment equation
    
    tau = (r*Bg + G)/(w*N); % government budget constraint

    % equilibrium conditions: 9 equations with 9 unknowns
    r == (1.0 + ii)/(1.0 + pii) - 1.0; % Fisher Equation      
    ap + bp == p + Bg; % asset demand = asset supply
    log(1 + pii) + log(w/w(-1)) ==  kappaw*(varphi*N^(1+nu) - (1-tau)*w*N*uce/muw) + beta*(log(1+pii(+1)) + log(w(+1)/w)); % Wage Phillips Curve
    log(1 + pii) ==  kappap*(w*N/((1-alpha)*Y) - 1/mup) + log(1+pii(+1))*Y(+1)/Y/(1+r(+1)); % Price Phillips Curve
    % valuation equation
    (1 + r(+1))*(1 + (K/K(-1) - 1)/(delta*epsI)) == alpha*w(+1)*N(+1)/((1-alpha)*K) - (K(+1)/K - (1-delta) + (K(+1)/K - 1).^2/(2*delta*epsI)) + (1 + (K(+1)/K - 1)/(delta*epsI))*K(+1)/K;
    Y == Z*(K(-1)^alpha)*(N^(1-alpha)); % production function
    d == Y - w*N - I - psip; % dividend equation
    0 == (d(+1) + p(+1)) - p*(1.0 + r(+1)); % no-arbitrage equation
    bp == Bh; % liquid asset supply = demand
    
    % Bound Constraints for Unknowns
    r <= (1/beta - 1  - 1e-10);
    pii >= beta - 1.0 + 1e-10;   % 1+r <= 1/beta implies 1+pii >=beta*(1+ii)>= beta
    N >= 0.5; % tau <=1 implies N > = (r*Bg + G)/w = 0.3556 when r = 0.0125

    % post evaluation
    goods_mkt = c + I + G + psip + chi + omega*bp - Y; % residual = aggregate demand - aggregate supply;
end;

Use the Toolbox

Parse .hmod Files

Suppose you have hank2_ssj.hmod in the working directory. Call in MATLAB

clear_hans
hans hank2_ssj
Parsing...ok!
Compiling dpopt_vfi
Building with 'MinGW64 Compiler (C++)'.
MEX completed successfully.
ans = struct with fields:
               parameters: {'beta'  'sigma'  'varphi'  'nu'  'chi0'  'chi1'  'omega'  'ra'  'r'  'w'  'tau'  'N'}
                var_shock: 'e'
                var_state: {'b'  'a'}
              var_pre_vfi: 'coh'
               var_policy: {'bp'  'ap'  'c'}
                  var_aux: {'chi'  'uce'}
                  var_agg: {'r'  'pii'  'Y'  'K'  'N'  'w'  'p'  'd'  'Bh'}
            var_agg_shock: 'm_shock'
           var_agg_params: {'G'  'Y_star'  'muw'  'kappap'  'kappaw'  'phi_y'  'beta'  'delta'  'omega'  'alpha'  'Bg'  'rstar'  'mup'  'phi'  'varphi'  'nu'  'epsI'  'Z'}
           var_agg_assign: {'ii'  'pshare'  'ra'  'I'  'psip'  'Q'  'tau'}
    var_agg_in_ind_params: {'r'  'N'  'w'  'ra'  'tau'}
            var_ind_in_eq: {'c'  'chi'  'uce'  'ap'  'bp'}

As shown, HANS parases the hmod file, compiles the mex files, and returns a structure that describes the model.

clear_hans in the first line clears any cached files generated before.

Calibration and The Stationary Equilibrium

Call the generated code "solve_cali.m" to calibrate model parameters to satisfy data targets in stationary equilibria. By construction, the result of the calibration also returns the corresponding stationary equilibrium at the calibrated parameters.

tic;
cali_rslt = solve_cali;
Evaluating at var_agg: 0.976, 6.416, 
VFI converged (metric_v) 5.93717e-08, (metric_pol) 0 in 745 iterations
Range of b: 0, 21.3544
Range of a: 0, 76.6147
Equilbirium residual: -2.11182, -0.137683, 

     Iteration    Line Search     Func-count    ||f(x)||^2     fnorm(x)      Norm of Step      Step Size      Algorithm

             0              0              1   4.478758e+00   2.116308e+00
Evaluating at var_agg: 0.9761, 6.416, 
VFI converged (metric_v) 1.76167e-08, (metric_pol) 0 in 566 iterations
Range of b: 0, 21.3544
Range of a: 0, 76.6147
Equilbirium residual: -1.99684, -0.128231, 
Evaluating at var_agg: 0.976, 6.41664, 
VFI converged (metric_v) 2.4074e-08, (metric_pol) 0 in 551 iterations
Range of b: 0, 21.3544
Range of a: 0, 76.6147
Equilbirium residual: -2.11176, -0.137825, 
Evaluating at var_agg: 0.977823, 6.57271, 
VFI converged (metric_v) 5.67098e-08, (metric_pol) 0 in 686 iterations
Range of b: 0, 21.3544
Range of a: 0, 91.6827
Equilbirium residual: 0.352755, 0.0130091, 

             1              1              4   1.246056e-01   3.529952e-01   1.567168e-01   1.000000e+00         Newton
Evaluating at var_agg: 0.977615, 6.52444, 
VFI converged (metric_v) 2.97529e-08, (metric_pol) 0 in 616 iterations
Range of b: 0, 21.3544
Range of a: 0, 83.8117
Equilbirium residual: 0.0235029, 0.000380076, 

             2              1              5   5.525299e-04   2.350595e-02   4.826829e-02   1.000000e+00        Broyden
Evaluating at var_agg: 0.977603, 6.51925, 
VFI converged (metric_v) 1.15759e-08, (metric_pol) 0 in 532 iterations
Range of b: 0, 21.3544
Range of a: 0, 83.8117
Equilbirium residual: 0.00450506, 0.000139501, 

             3              1              6   2.031503e-05   4.507220e-03   5.188706e-03   1.000000e+00        Broyden
Evaluating at var_agg: 0.977599, 6.51834, 
VFI converged (metric_v) 9.93413e-09, (metric_pol) 3.90935e-08 in 479 iterations
Range of b: 0, 21.3544
Range of a: 0, 83.8117
Equilbirium residual: -0.000365811, -1.01377e-05, 

             4              1              7   1.339202e-07   3.659510e-04   9.066127e-04   1.000000e+00        Broyden
Evaluating at var_agg: 0.9776, 6.51842, 
VFI converged (metric_v) 9.86051e-09, (metric_pol) 4.01353e-08 in 365 iterations
Range of b: 0, 21.3544
Range of a: 0, 83.8117
Equilbirium residual: -7.92847e-06, -3.685e-07, 

             5              1              8   6.299649e-11   7.937033e-06   7.206637e-05   1.000000e+00        Broyden
Evaluating at var_agg: 0.9776, 6.51842, 
VFI result reused
Range of b: 0, 21.3544
Range of a: 0, 83.8117
Equilbirium residual: -7.92847e-06, -3.685e-07, 
fprintf('Time used for Calibration:\n');
Time used for Calibration:
toc;
Elapsed time is 151.809596 seconds.

We can inspect the values of cali_rslt

cali_rslt
cali_rslt = struct with fields:
                vfi_rslt: [1×1 struct]
                   stats: [1×1 struct]
                 var_agg: [1×1 struct]
                    dist: [3×50×70 double]
              parameters: [1×1 struct]
          var_agg_params: [1×1 struct]
             shock_trans: [3×3 double]
    shock_invariant_dist: [3×1 double]
                   resid: [2×1 double]
               exit_flag: 1
               rslt_type: 'ss'
           solver_output: [1×1 struct]
           ind_rslt_dict: [1×1 ArrayKeyMap]

With all parameters and the corresponding stationary equilibrium calculated, we can directly proceed to the transition path computation. Nevertheless, as an extra check we can resolve the stationary equilibrium by passing the returned structure "cali_rslt" to "solve_ss.m" function:

tic;
ss_rslt = solve_ss(cali_rslt);
Evaluating at var_agg: 0.0125, 1, 
VFI converged (metric_v) 0.00771125, (metric_pol) 0 in 33 iterations
Range of b: 0, 21.3544
Range of a: 0, 83.8117
Equilbirium residual: -7.92235e-06, 1.57066e-10, 

     Iteration    Line Search     Func-count    ||f(x)||^2     fnorm(x)      Norm of Step      Step Size      Algorithm

             0              0              1   6.276356e-11   7.922346e-06
Evaluating at var_agg: 0.0125, 1, 
VFI result reused
Range of b: 0, 21.3544
Range of a: 0, 83.8117
Equilbirium residual: -7.92235e-06, 1.57066e-10, 
fprintf('Time used for solving stationary equilibrium:\n');
Time used for solving stationary equilibrium:
toc;
Elapsed time is 0.435768 seconds.

As expected, the initial condition has already satisfied all conditions for the stationary equilibrium. As a result, the computaton stops right after the initial iteration.

All stationary equilibrium information is returned in the structure ss_rslt.

ss_rslt
ss_rslt = struct with fields:
                vfi_rslt: [1×1 struct]
                   stats: [1×1 struct]
                 var_agg: [1×1 struct]
                    dist: [3×50×70 double]
              parameters: [1×1 struct]
          var_agg_params: [1×1 struct]
             shock_trans: [3×3 double]
    shock_invariant_dist: [3×1 double]
                   resid: [2×1 double]
               exit_flag: 1
               rslt_type: 'ss'
           solver_output: [1×1 struct]
           ind_rslt_dict: [1×1 ArrayKeyMap]

We can inspect the values of equilibrium variables and verify that all calibration targets are satisfied:

ss_rslt.var_agg
ans = struct with fields:
          r: 0.0125
          Y: 1
    m_shock: 0
         ra: 0.0125
         mc: 0.9850
          K: 10
          N: 1
          w: 0.6600
          d: 0.1400
          p: 11.2000
        tau: 0.3561
      rstar: 0.0125
     Y_star: 1
        pii: 0
         Bh: 1.0400

To visualize the results, let's plot the policy functions for a given idiosyncratic shock in a 3-Dimensional graph:

% Pick an idiosyncratic shock and the corresponding (b,a) mesh in 2-D

shock_idx = 2;

b_mesh_2d = squeeze(ss_rslt.vfi_rslt.var_state_mesh.b_mesh(shock_idx,:,:));

a_mesh_2d = squeeze(ss_rslt.vfi_rslt.var_state_mesh.a_mesh(shock_idx,:,:));

% Plot the VFI Policy in 3-D Mesh

surf(b_mesh_2d,a_mesh_2d,squeeze(ss_rslt.vfi_rslt.var_policy.bp(shock_idx,:,:)))

xlabel('b');

ylabel('a');

title('Policy functions for b');

surf(b_mesh_2d,a_mesh_2d,squeeze(ss_rslt.vfi_rslt.var_policy.ap(shock_idx,:,:)))

xlabel('b');

ylabel('a');

title('Policy functions for a');

surf(b_mesh_2d,a_mesh_2d,squeeze(ss_rslt.vfi_rslt.var_policy.c(shock_idx,:,:)));

xlabel('b');

ylabel('a');

title('Policy functions for c');

We can also plot the joint stationary distribution over (b, a), as well as the marginal distribution over b or a.

% Extract (b,a) mesh in 2-D

b_mesh_2d = squeeze(ss_rslt.vfi_rslt.var_state_mesh.b_mesh(1,:,:));

a_mesh_2d = squeeze(ss_rslt.vfi_rslt.var_state_mesh.a_mesh(1,:,:));

bgrid = ss_rslt.vfi_rslt.var_state.b(:); % b grid

agrid = ss_rslt.vfi_rslt.var_state.a(:); % a grid

% Plot the histogram of the Stationary Distribution in 3-D Mesh

dist_2d = squeeze(sum(ss_rslt.dist)); % histogram on (b,a) grid

dist_b = sum(dist_2d,2); % histogram on b grid

dist_a = reshape(sum(dist_2d,1),[],1); % histogram on a grid

surf(b_mesh_2d,a_mesh_2d,dist_2d);

xlabel('b');

ylabel('a');

title('Histogram of Stationary Distribution in (b,a)');

plot(bgrid,dist_b);

xlabel('b');

ylabel('prob');

title('Histogram of Stationary Distribution in b');

plot(agrid,dist_a);

xlabel('a');

ylabel('prob');

title('Histogram of Stationary Distribution in a');

Solve the Transition Path After A Temporay Shock

Solve the linearized transition path

We first construct the sequence of monetary policy shocks to the Taylor rule:

T = 300; % Transition periods

rho_r = 0.6; % persistence of monetary shock

sig_r = -0.01/4; % a negative interest rate shock of 25 basis points

m_shock_t = sig_r * rho_r.^(0:(T-1)); % the path of interest rate shock

figure;

plot(m_shock_t)

We then pass the monetary shock and the solved ss_rslt structure into "solve_trans_linear.m" to solve the linearized transition path.

options = struct;
options.T = T;
options.m_shock_t = m_shock_t;
tic;
trans_linear_rslt = solve_trans_linear(ss_rslt, options);
fprintf('Time used for solving the linearized transition path:\n');toc;
Time used for solving the linearized transition path:
Elapsed time is 64.028450 seconds.
trans_linear_rslt
trans_linear_rslt = struct with fields:
             var_agg_t: [1×1 struct]
       var_agg_shock_t: [1×1 struct]
      jacs_ind_wrt_agg: {5×5 cell}
    jac_eqs_of_var_agg: [2700×2700 double]
               irf_ssj: [2700×300 double]

We can plot the transition path (note that K is the next period capital).

figure;

plot(trans_linear_rslt.var_agg_t.K)

title('K');

xlabel('time');

Solve the nonlinear transition path

Now we pass the solved ss_rslt to solve_nonlinear_trans(init_ss, final_ss, options) to solve the nonlinear transition path.

Note that we pass ss_rslt structure as both init_ss and final_ss.

tic;
trans_nonlinear_rslt = solve_trans_nonlinear(ss_rslt, ss_rslt, options);
     Iteration    Line Search     Func-count    ||f(x)||^2     fnorm(x)      Norm of Step      Step Size      Algorithm

             0              0              1   9.999633e-08   3.162220e-04

             1              1              2   2.404074e-10   1.550508e-05   1.405179e-03   1.000000e+00        Broyden
fprintf('Time used for solving the nonlinear transition path:\n');toc;
Time used for solving the nonlinear transition path:
Elapsed time is 113.535368 seconds.
trans_nonlinear_rslt
trans_nonlinear_rslt = struct with fields:
          var_agg_t: [1×1 struct]
    var_agg_shock_t: [1×1 struct]
                  T: 300
           vfi_rslt: [1×1 struct]
        shock_trans: [3×3 double]
           exitflag: 1
             HANS_x: [2700×1 double]
          eqs_resid: [2700×1 double]
          rslt_type: 'nonlinear_trans'
      solver_output: [1×1 struct]

We compare the solution of the nonlinear transition path and the linear transition path.

figure;

plot(trans_linear_rslt.var_agg_t.K);

hold on;

plot(trans_nonlinear_rslt.var_agg_t.K, '--');

xlabel('time');

title('K');

legend({'Linear','Nonlinear'});

We can also inspect other defined equilibrium varialbes along the transition path.

figure;

plot(trans_nonlinear_rslt.var_agg_t.w);

title('w');

xlabel('time');