# Mendoza (2010): Sudden Stops with Asset Price Deflation¶

Mendoza (2010) studies an open economy model with a borrowing constraint tied to the capital price. A negative shock that lowers the capital price tightens the collateral constraint, leading to deleveraging and reductions in financing for both investment and working capital, so the effect of the shock gets amplified. The model generates financial crises with current account reversals resembling those experienced by many emerging economies. The collateral constraints are occasionally binding in the ergodic set. The equilibrium policy and state transition functions are highly nonlinear. The model has two endogenous states, capital and bond, and can be solved using the toolbox very easily (in 200 lines of GDSGE code and within a minute on a regular laptop).

Though the model can be directly solved over the natural state space (capital, bond) using the toolbox, practically it is sometimes difficult to determine the feasible state space over which an equilibrium exists. For example, in the current context, given capital stock, an economy too much in debt may not be able to find feasible allocations. This suggests that the lower bound of bond depends on the capital stock and the feasible region of the state space is not a simple box (rectangle). We demonstrate how such a model can solved with state transformation and implicit state transition functions, without the hassle of determining the feasible region ex-ante.

## The Model¶

There is a single consumption good. The inter-temporal expected utility is given by

$\mathbb E_{0} \sum_{t=0}^{\infty} \beta^t[ u\left(c_{t}-N\left(L_{t}\right)\right)]$

where $$\beta$$ is the discount factor, $$c_t$$ is consumption, $$L_t$$ is labor supply. $$u(\cdot)$$ is the utility function.

The production function is

$A_tF(k_t,L_t,v_t),$

where $$A_t$$ is productivity which follows an exogenous Markov process. $$F$$ is a constant-return-to-scale production function. $$k_t$$ is capital, $$L_t$$ is labor, and $$v_t$$ is intermediate goods used.

$$\phi$$ fraction of the cost of hiring $$L_t$$ and $$v_t$$ needs to be financed externally, at a (world) gross interest rate $$R_t$$ that follows an exogenous process. The economy can trade a state non-contingent bond $$b_{t+1}$$ (who negative values imply borrowing) with the world, but borrowing is subject to the following constraint

$q_{t}^{b} b_{t+1}-\phi R_{t}\left(w_{t} L_{t}+p_{t} v_{t}\right) \geq-\kappa q_{t} k_{t+1},$

which says that the borrowing plus the finance for working capital cannot exceed $$\kappa$$ fraction of the value of the next period capital $$k_{t+1}$$, where $$p_t$$ is intermediate goods price following an exogenous process, $$w_t$$ is wage and $$q_t$$ is the endogenous asset price of the next period capital determined below. It should be clear that the binding collateral constraint will affect the financing of working capital.

Investment is subject to a convex adjustment cost

$i_{t}=k_{t+1} - (1-\delta)k_t+\frac{1}{2} a \frac{(k_{t+1}-k_t)^{2}}{k_t},$

where $$a$$ is a parameter determining the level of adjustment cost. The next period capital price is competitively determined as

$q_t=\frac{\partial i_t}{\partial k_{t+1}}=1+a\frac{k_{t+1}-k_t}{k_t}.$

We also write down the shadow price of current capital

$\begin{split}q_t+d_t & \equiv A_tF_1(k_t,L_t,v_t)-\frac{\partial i_t}{\partial k_{t}}\\ &=A_tF_1(k_t,L_t,v_t)- (1-\delta) -a\frac{k_{t+1}-k_t}{k_t}- \frac{a}{2}\frac{(k_{t+1}-k_t)^{2}}{k_t^2}.\end{split}$

Under the competitive equilibrium, the agent maximizes the inter-temporal expected utility subject to

$\begin{split}(1+\tau)c_{t}+i_{t}=A_t F\left(k_{t}, L_{t}, v_{t}\right)-p_{t} v_{t}-\phi\left(R_{t}-1\right)\left(w_{t} L_{t}+p_{t} v_{t}\right)-q_{t}^{b} b_{t+1}+b_{t}, \\ q_{t}^{b} b_{t+1}-\phi R_{t}\left(w_{t} L_{t}+p_{t} v_{t}\right) \geq-\kappa q_{t} k_{t+1}, \\ i_{t}=k_{t+1} - (1-\delta)k_t+\frac{1}{2} a^2 \frac{(k_{t+1}-k_t)^{2}}{k_t},\end{split}$

taking exogenous processes $$A_t,R_t,p_t$$ and $$q_t^b(=1/R_t)$$, and equilibrium prices $$w_t$$ and $$q_t$$ as given. Denote $$\lambda_t$$ as the multiplier for the budget constraint and $$\mu_t$$ as the multiplier for the collateral constraint, the first order conditions and complementarity slackness conditions read:

$\begin{split}& \{b_{t+1}\}: -q_t^b\lambda_t+\mu_tq_t^b+\beta \mathbb{E}_t \lambda_{t+1}=0 \\ & \{k_{t+1}\}: \beta \mathbb{E}_t \lambda_{t+1}A_{t+1}F_1(k_{t+1},L_{t+1},v_{t+1}) - \beta \mathbb{E}_t\lambda_{t+1} \frac{\partial i_{t+1}}{\partial k_{t+1}} - \lambda_t q_t +\mu_t \kappa q_t=0 \\ & \{\mu_t\}: [q_{t}^{b} b_{t+1}-\phi R_{t}\left(w_{t} L_{t}+p_{t} v_{t}\right) +\kappa q_{t} k_{t+1}]{\mu}_t=0 \\ & \{L_t\}: A_{t}F_2(k_t,L_t,v_t)=w_t (1+\phi(R_t-1)+\frac{\mu_t}{\lambda_t} \phi R_t); \ w_t=N'(L_{t})(1+\tau) \\ & \{v_t\}: A_{t}F_3(k_t,L_t,v_t)=p_{t}(1+\phi(R_t-1)+\frac{\mu_t}{\lambda_t} \phi R_t) \\ & \{c_t\}: u'(c_t-N(L_t)) = \lambda_t(1+\tau)\end{split}$

The exogenous states are $$z=(A,R,p)$$. The natural endogenous states are $$(k,b)$$. A recursive equilibrium is functions $$b'(z,k,b), k',\mu,L,v,c,\lambda,q,w$$ that satisfy the above equations. The system can be further simplified. Before doing so, we introduce the state transformation to ensure the feasible state space is square.

## State Space Transformation¶

If the economy is too indebted, it may never be able to finance positive consumption subject to the current collateral constraint, and the equilibrium does not exist. Therefore, the lower bound of the bond holding of the feasible region in the state space depends on the exogenous shocks and capital, rendering the feasible state space non-rectangle. This non-rectangle feasible region also cannot be determined ex-ante, since the collateral constraint depends on the endogenous capital price. This is a well-known challenge encountered by many researchers studying models with a borrowing constraint tied to an asset price. The challenge can be solved by noticing the exact boundary condition of the initial problem: the boundary of the feasible region is such that consumption (more precisely, the GHH consumption-labor bundle) is zero. Therefore, if we transform the state space to consumption and capital, then the feasible region will be a simple rectangle.

Define the GHH consumption-labor bundle as

$\tilde{c}_t \equiv c_t-N(L_t).$

With such state transformation, we write down the minimum recursive system that represents the equilibrium. A recursive equilibrium is functions $$k'(z,\tilde{c},k)$$, $$b'(z,\tilde{c},k)$$, $$L(z,\tilde{c},k)$$, $$v(z,\tilde{c},k)$$, $$\mu(z,\tilde{c},k)$$, $$\tilde{c}'(z';z,\tilde{c},k)$$ such that (omitting the dependence on $$(z,\tilde{c},k)$$ when there is no confusion)

$\begin{split}& -q^b\lambda+\mu q^b+\beta \mathbb{E} [\lambda'(z',\tilde{c}'(z'),k')|z]=0 \\ & \beta \mathbb{E}[ \lambda'(z',\tilde{c}'(z'),k')\{d'(z',\tilde{c}'(z'),k')+q'(z',\tilde{c}'(z'),k')\}|z ]- \lambda q +\mu \kappa q=0 \\ & AF_2(k,L,v)=w (1+\phi(R-1)+\frac{\mu}{\lambda} \phi R) \\ & AF_3(k,L,v)=p(1+\phi(R-1)+\frac{\mu}{\lambda} \phi R) \\ & (1+\tau)\tilde{c}'(z')=B'(z',\tilde{c}'(z'),k') +b', \forall z'\end{split}$

where $$\lambda,w,q,\lambda'$$ are intermediate variables and can be evaluated according to

$\begin{split}& \lambda = u'(\tilde{c}) \\ & w= N'(L) \\ & q=1+a\frac{k'-k}{k} \\ & \lambda'= u'(\tilde{c}'(z'))\end{split}$

The last equation $$(1+\tau)\tilde{c}'(z')=B'(z',\tilde{c}'(z'),k') + b'$$ specifies the consistency equation for $$\tilde{c}_{t+1}$$, which is simply the period-$$(t+1)$$ budget constraint:

$(1+\tau)\underbrace{[c_{t+1}-N(L_{t+1})]}_{\tilde{c}_{t+1}}=B_{t+1} + b_t,$

where

$\begin{split}B_{t+1}&=&A_{t+1} F\left(k_{t+1}, L_{t+1}, v_{t+1}\right)-p_{t+1} v_{t+1}-\phi\left(R_{t+1}-1\right)\left(w_{t+1} L_{t+1}+p_{t+1} v_{t+1}\right)\\ &&-q_{t+1}^{b} b_{t+2}-i_{t+1}-(1+\tau)N(L_{t+1}).\end{split}$

Notice, how the budget constraint is decreased by $$(1+\tau)N(L_{t+1})$$ so that this becomes a consistency equation for $$\tilde{c}_{t+1}$$.

The recursive system is solved by taking $$d'(z,\tilde{c},k),q'(z,\tilde{c},k),B'(z,\tilde{c},k)$$ as known at each time iteration step, and updating these objects following their definitions.

## The gmod File¶

We use the exact specification and parameterization as in Mendoza (2010). In particular, for preference

$u(c) = \frac{c^{1-\sigma}}{1-\sigma}$
$N(L) = L^{\omega}/\omega.$

The production function is Cobb-Douglas:

$F(k,L,v) = k^{\gamma}L^{\alpha}v^{\eta}.$

For parameters, see the gmod file (mendoza2010) listed below

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 % Parameters parameters gamma sigma omega beta alpha eta delta a phi kappa tau; parameters flow_scale; beta = 0.92; % Discount factor sigma = 2; % CRRA omega = 1.846; % Labor elasticity gamma = 0.306; % Capital share alpha = 0.592; % Labor share eta = 0.102; % Intermediate share delta = 0.088; % Depreciation rate a = 2.75; % Capital adjustment cost phi = 0.26; % Working capital exernal finance share kappa = 0.2; % Tightness of collateral constraint tau = 0.17; % Consumption tax flow_scale = 1e4; % Budget normalization factor % Toolbox options PrintFreq = 10; SaveFreq = 100; SIMU_RESOLVE=0; SIMU_INTERP=1; % Use interp for fast simulation % Shocks var_shock A R p; shock_num = 8; % Markov process of shocks (from Mendoza AER (2010): % Transition probability matrix shock_trans = [0.51363285, 0.16145114, 0.07773228, 0.02443373, 0.16145114, 0.03201937, 0.02443373, 0.00484575;... 0.03201937, 0.64306463, 0.00484575, 0.09732026, 0.16145114, 0.03201937, 0.02443373, 0.00484575;... 0.07773228, 0.02443373, 0.51363285, 0.16145114, 0.02443373, 0.00484575, 0.16145114, 0.03201937;... 0.00484575, 0.09732026, 0.03201937, 0.64306463, 0.02443373, 0.00484575, 0.16145114, 0.03201937;... 0.03201937, 0.16145114, 0.00484575, 0.02443373, 0.64306463, 0.03201937, 0.09732026, 0.00484575;... 0.03201937, 0.16145114, 0.00484575, 0.02443373, 0.16145114, 0.51363285, 0.02443373, 0.07773228;... 0.00484575, 0.02443373, 0.03201937, 0.16145114, 0.09732026, 0.00484575, 0.64306463, 0.03201937;... 0.00484575, 0.02443373, 0.03201937, 0.16145114, 0.02443373, 0.07773228, 0.16145114, 0.51363285]; % Values of the stochastic process AMean = 7 % Scale to match average output with SCU preference in Mendoza (2010) A = AMean*[0.986595988257; 1.013404011625; 0.986595988257; 1.013404011625; 0.986595988257; 1.013404011625; 0.986595988257; 1.013404011625]; R = [1.064449652804; 1.064449652804; 1.064449652804; 1.064449652804; 1.106970405389; 1.106970405389; 1.106970405389; 1.106970405389]; p = [0.993455002946; 0.993455002946; 1.062225686321; 1.062225686321; 0.993455002946; 0.993455002946; 1.062225686321; 1.062225686321]; % State var_state cTilde k; cTilde_min = 100; cTilde_max = 300; cTilde_shift = 10; cTilde_pts = 30; cTilde = linspace(cTilde_min,cTilde_max,cTilde_pts); k_min = 600; k_max = 900; k_shift = 10; k_pts = 30; k = linspace(k_min,k_max,k_pts); % Initial period var_policy_init L v; inbound_init L 0 1000; inbound_init v 0 1000; var_aux_init flow q d; model_init; % Some trivial equations w = L^(omega-1)*(1+tau); Y = A* k^gamma * L^alpha * v^eta; F_1 = gamma*Y/k; F_2 = alpha*Y/L; F_3 = eta*Y/v; NL = L^omega/omega; flow = Y - p*v - phi*(R-1)*(w*L+p*v) - NL*(1+tau); lambda = cTilde^(-sigma); q = 1; d = F_1-delta; equations; F_2 - w*(1+phi*(R-1)); F_3 - p*(1+phi*(R-1)); end; end; % Interpolation var_interp flow_interp q_plus_d_interp; % Initial initial flow_interp flow/flow_scale; initial q_plus_d_interp q + d; % Transition flow_interp = flow/flow_scale; q_plus_d_interp = q + d; % Policy ADAPTIVE_FACTOR = 1.5; var_policy L v mu nbNext kNext cTildeNext[8]; inbound L 0 1000 adaptive(ADAPTIVE_FACTOR); inbound v 0 1000 adaptive(ADAPTIVE_FACTOR); inbound mu 0 1; inbound nbNext 0 2000 adaptive(ADAPTIVE_FACTOR); inbound kNext 0 1000 adaptive(ADAPTIVE_FACTOR); inbound cTildeNext 0.0 500 adaptive(ADAPTIVE_FACTOR); % Extra variables returned var_aux lambda flow q d c Y bNext b inv nx gdp b_gdp nx_gdp lev wkcptl; var_output c Y cTildeNext q mu inv kNext bNext b v nx gdp b_gdp nx_gdp lev wkcptl; model; % Output and marginal product Y = A* k^gamma * L^alpha * v^eta; F_1 = gamma*Y/k; F_2 = alpha*Y/L; F_3 = eta*Y/v; w = L^(omega-1)*(1+tau); % Some calculations lambda = cTilde^(-sigma)/(1+tau); inv = kNext - k*(1-delta) + a/2*(kNext-k)^2/k; q = 1+a*(kNext-k)/k; d = F_1 - delta + a/2*(kNext-k)^2/(k^2); qb = 1/R; % Transform back bNext bNext = (nbNext + phi*R*(w*L+p*v) - kappa*q*kNext)/qb; % Interpolation [flow_future',q_plus_d_future'] = GDSGE_INTERP_VEC'(cTildeNext',kNext); flow_future' = flow_future'*flow_scale; % some named equations lambda_future' = cTildeNext'^(-sigma)/(1+tau); euler_bond_residual = -1 + mu + beta*GDSGE_EXPECT{lambda_future'}/(qb*lambda); euler_capital_residual = -1 + kappa*mu + beta*GDSGE_EXPECT{lambda_future'*q_plus_d_future'}/(q*lambda); % consistency equation cTilde_consis' = flow_future' + bNext - cTildeNext'*(1+tau); % Extra variables returned NL = L^omega/omega; flow = Y - p*v - phi*(R-1)*(w*L+p*v) - qb*bNext - inv - NL*(1+tau); c = cTilde + NL; b = cTilde*(1+tau) - flow; gdp = Y-p*v; b_gdp = b / gdp; nx = qb*bNext - b + phi*(R-1)*(w*L+p*v); nx_gdp = nx / gdp; lev = (qb*bNext-phi*R*(w*L+p*v)) / (q*kNext); wkcptl = phi*(w*L+p*v); equations; euler_bond_residual; nbNext*mu; euler_capital_residual; F_2 - w*(1+phi*(R-1)+mu*phi*R); F_3 - p*(1+phi*(R-1)+mu*phi*R); cTilde_consis'; end; end; simulate; num_periods = 50000; num_samples = 24; initial cTilde 200; initial k 800; initial shock ceil(shock_num/2); var_simu c Y q inv b mu bNext v nx gdp b_gdp nx_gdp lev wkcptl; cTilde' = cTildeNext'; k' = kNext; end; 

Most of the gmod codes are straightforward and should be self-explanatory. We do want to highlight that the state transformation and implicit state transition are implemented by including $$\tilde{c}'(z')$$ as unknowns in Line 97

 97 var_policy L v mu nbNext kNext cTildeNext[8]; 

defining the consistency equations in Line 136

 136  cTilde_consis' = flow_future' + bNext - cTildeNext'*(1+tau); 

and including them in the system of equations in Line 156

 156  cTilde_consis'; 

We also discuss the necessary transformation when solving the system of equations with potentially non-box constraints (i.e., the collateral constraint which specifies that borrowing cannot exceed a fraction of the value of capital that depends on the endogenous capital price). This should be familiar if readers already came from the Bianchi (2011) or Heaton and Lucas (1996) examples. This is done by defining a variable and including it as unknown

$nb_{t+1}=q_{t}^{b} b_{t+1}-\phi R_{t}\left(w_{t} L_{t}+p_{t} v_{t}\right) +\kappa q_{t} k_{t+1}$

so the collateral constraint becomes $$nb_{t+1} \geq 0$$. And transform $$nb_{t+1}$$ back to $$b_{t+1}$$ in Line

 125  bNext = (nbNext + phi*R*(w*L+p*v) - kappa*q*kNext)/qb; 

Another transformation to ensure numerical stability is to normalize each of the equations in the system so they take value in similar scales (the equation solver handles this automatically to some extent but researchers can do better with prior information). For example, the Euler equation for bond is written as

$-1 +\mu_t / \lambda_t +\beta \mathbb{E}_t \lambda_{t+1} / (\lambda_t q_t^b) = 0$

so the equation takes value between [0,1]. Accordingly, we have defined the multiplier for the collateral constraint to be normalized by the marginal utility $$\hat{\mu}_t=\mu_t/\lambda_t$$ ($$\mu$$ in the code actually corresponds to $$\hat{\mu}$$), so from the bond Euler equation we know $$\hat{\mu}$$ must lie in [0,1], which is specified as the bound for $$\mu$$ entering the equation solver. These procedures are not necessary but are simple and effective ways to solve the system faster and more accurately.

Now let’s run the policy iterations by calling iter_mendoza2010.m returned from the online GDSGE compiler, which produces output

>> IterRslt = iter_mendoza2010;

Iter:10, Metric:0.00819545, maxF:3.83276e-09
Elapsed time is 5.885684 seconds.
Iter:20, Metric:0.00329338, maxF:1.72695e-09
Elapsed time is 2.855950 seconds.
Iter:30, Metric:0.00139174, maxF:1.72276e-09
Elapsed time is 2.474352 seconds.

...

Iter:117, Metric:9.38205e-07, maxF:1.53505e-09
Elapsed time is 1.252342 seconds.


Notice so far, relatively crude grids for the state space are specified in the gmod file. We can refine the solution on finer grids starting from a previously converged solutions, by simply specifying the converged solution as a WarmUp. Here we enlarge the state space to have 80 grid points along both $$\tilde{c}$$ and $$k$$ space, increasing the number of (exogenous and endogenous) states to $$8\times 80\times 80 = 51200$$.

>> options = struct;
cTilde_pts_fine = 80;
options.cTilde = linspace(IterRslt.var_state.cTilde(1),IterRslt.var_state.cTilde(end),...
cTilde_pts_fine);
k_pts_fine = 80;
options.k = linspace(IterRslt.var_state.k(1),IterRslt.var_state.k(end),k_pts_fine);
options.WarmUp = IterRslt;
IterRslt = iter_mendoza2010(options);

Iter:120, Metric:0.0018683, maxF:8.36063e-09
Elapsed time is 6.374272 seconds.
Iter:130, Metric:6.5268e-05, maxF:8.48359e-09
Elapsed time is 16.458581 seconds.
Iter:140, Metric:1.57919e-05, maxF:8.50337e-09
Elapsed time is 13.727288 seconds.
Iter:150, Metric:4.02936e-06, maxF:8.41128e-09
Elapsed time is 12.966158 seconds.
Iter:160, Metric:1.0877e-06, maxF:7.73927e-09
Elapsed time is 11.709191 seconds.
Iter:161, Metric:9.56597e-07, maxF:7.23566e-09
Elapsed time is 1.117845 seconds.


The refined iterations continue to converge in another minute. It’s always a good convention to start with some cruder grid to assess the properties (such as the range of the ergodic set) of the solution of a model, and refine it using the crude solution as a WarmUp.

Now we inspect the policy functions. We first come back to our motivation for state space transformation: the feasible region in the natural state space $$(k,b)$$ is non-rectangle. This can be verified by projecting the state variable $$\tilde{c}$$ onto policy $$b$$ and state variable $$k$$.

>> figure;
[cTilde_mesh,k_mesh] = ndgrid(IterRslt.var_state.cTilde,IterRslt.var_state.k);
surf(squeeze(IterRslt.var_aux.b(1,:,:)),k_mesh,cTilde_mesh);
xlabel('b');
ylabel('k');
title('Feasible Region of Shock Index 1');
view(-0.06,90);


As shown in the figure, the level of bond that gives a same lower bound of consumption (here $$\tilde{c}$$ equal to 100) is declining in capital, suggesting that a higher capital allows the current debt to be higher. However, determining this lower bound of bond level is difficult; but determining the boundary condition for consumption is much easier ($$\tilde{c}=0$$).

We can inspect the policy functions by projecting them to the natural state space as well. For example, the capital price:

>> figure;
[cTilde_mesh,k_mesh] = ndgrid(IterRslt.var_state.cTilde,IterRslt.var_state.k);
surf(squeeze(IterRslt.var_aux.b(1,:,:)),k_mesh,squeeze(IterRslt.var_aux.q(1,:,:)));
xlabel('b');
ylabel('k');
title('Capital Price, Shock Index 1');


And the multiplier (normalized by marginal utility) for collateral constraint:

>> figure;
[cTilde_mesh,k_mesh] = ndgrid(IterRslt.var_state.cTilde,IterRslt.var_state.k);
surf(squeeze(IterRslt.var_aux.b(1,:,:)),k_mesh,squeeze(IterRslt.var_policy.mu(1,:,:)));
xlabel('b');
ylabel('k');
title('Multiplier for the Collateral Constraint, Shock Index 1');
view(-16,26);


These two policy functions show that the collateral constraint tends to bind when bond is low (borrowing is high). And capital price falls rapidly in borrowing when the collateral constraint binds.

With the converged policy and state transition functions, we can simulate the economy’s ergodic set. To remind readers, here we directly interpolate the state transition and policy functions to simulate, which is faster than resolving the system (the simulation with 24 paths of 50,000 periods completes in less than half a minute on a regular laptop). This is done by specifying in the gmod file

 22 SIMU_RESOLVE=0; SIMU_INTERP=1; % Use interp for fast simulation 

In this case, the variables to be recorded in the simulation need to be explicitly specified in var_output

 107 var_output c Y cTildeNext q mu inv kNext bNext b v nx gdp b_gdp nx_gdp lev wkcptl; 

and can then be included in var_simu

 167  var_simu c Y q inv b mu bNext v nx gdp b_gdp nx_gdp lev wkcptl; 

With implicit state transition functions for $$\tilde{c}$$ solved as part of the policy functions, the transitions need to be specified in the simulate block

 168  cTilde' = cTildeNext'; 

which means that the next period cTilde should pick policy cTildeNext at a sampled exogenous shock index.

Calling simulate_mnendoza2010 produces the following

>> SimuRslt = simulate_mendoza2010(IterRslt);

Periods: 1000
shock  cTilde       k       c       Y       q     inv       b   bNext       v      nx     gdp   b_gdp  nx_gdp     lev  wkcptl
1   166.9   747.1   282.5     429   1.048   78.98   109.7   86.46   43.32  -23.58   385.9  0.2843 -0.06110.0002417   76.13

...

Periods: 50000
shock  cTilde       k       c       Y       q     inv       b   bNext       v      nx     gdp   b_gdp  nx_gdp     lev  wkcptl
5   149.8   741.1   262.8   423.8  0.9844   61.04  -33.48  -31.26   42.34   13.19   381.7-0.08769 0.03455 -0.1525    74.4


We can inspect the histogram of the natural state variables $$(k,b)$$ by calling in MATLAB:

>> figure;
hist3([SimuRslt.b(:),reshape(SimuRslt.k(:,1:end-1),[],1)],[50,50],'CDataMode','auto','FaceColor','interp');
view([15,30]);
xlabel('Bond');
ylabel('Capital');
title('Histogram of Capital and Bond');
print('figures/histogram_b_k.png');


We are interested in the fraction of time that the economy spends in a crisis, i.e., states with binding collateral constraint:

>> muGreaterThanZero = SimuRslt.mu > 1e-7;
fractionCrisis = sum(muGreaterThanZero(:)) / numel(muGreaterThanZero(:));
fprintf('The fraction in crisis: %g\n',fractionCrisis);
The fraction in crisis: 0.131997


We can calculate the unconditional moments:

>> statsList = {'mean','std'};
varList = {'gdp','c','inv','nx_gdp','k','b_gdp','q','lev','v','wkcptl'};
ratioList = {'nx_gdp','b_gdp','lev'};
meanStats = struct;
stdStats = struct;
corrStats = struct;
autoCorrStats = struct;
startPeriod = 1001;
endPeriod = 50000;
for i_var = 1:length(varList)
varName = varList{i_var};
meanStats.(varName) = mean(reshape(SimuRslt.(varName)(:,startPeriod:endPeriod),[],1));
if max(strcmp(varName,ratioList))==1
stdStats.(varName) = std(reshape(SimuRslt.(varName)(:,startPeriod:endPeriod),[],1));
else
stdStats.(varName) = std(reshape(SimuRslt.(varName)(:,startPeriod:endPeriod),[],1)) ...
/ meanStats.(varName);
end
corrStats.(varName) = corr(reshape(SimuRslt.(varName)(:,startPeriod:endPeriod),[],1), ...
reshape(SimuRslt.gdp(:,startPeriod:endPeriod),[],1));
autoCorrStats.(varName) = corr(reshape(SimuRslt.(varName)(1,startPeriod+1:endPeriod),[],1), ...
reshape(SimuRslt.(varName)(1,startPeriod:endPeriod-1),[],1));
end

 GDSGE Mendoza (2010), Table 3, Panel C Mean gdp 400.646 388.339 c 277.349 267.857 inv 68.838 65.802 nx/gdp 0.018 0.024 k 780.598 747.709 b/gdp -0.031 -0.104 q 1.000 1.000 leverage -0.124 -0.159 v 43.268 41.949 wkcptl 78.555 75.455 Std in percent gdp 3.88% 3.85 c 3.68% 3.69 inv 13.30% 13.45 nx/gdp 2.67% 2.58 k 4.32% 4.31 b/gdp 13.97% 8.9 q 3.19% 3.23 leverage 6.53% 4.07 v 5.85% 5.84 wkcptl 4.30% 4.26 Corr w/ gdp gdp 1.000 1.000 c 0.898 0.931 inv 0.649 0.641 nx/gdp -0.131 -0.184 k 0.753 0.744 b/gdp -0.187 -0.298 q 0.408 0.406 leverage -0.174 -0.258 v 0.830 0.823 wkcptl 0.995 0.987 Auto corr gdp 0.820 0.815 c 0.796 0.766 inv 0.506 0.483 nx/gdp 0.514 0.447 k 0.963 0.963 b/gdp 0.985 0.087 q 0.452 0.428 leverage 0.990 0.040 v 0.773 0.764 wkcptl 0.798 0.777

Most statistics are similar to Mendoza (2010) except that Mendoza (2010) uses a Stationary Cardinal Utility (SCU) utility function and solves the model by solving a quasi social planner’s problem which produces (1) less persistent leverage and bond, and (2) higher leverage and more negative bond/gdp ratio. These differences are also noted in Mendoza and Villalvazo (2020).

## What’s Next?¶

This example illustrates how to solve a model with two endogenous state variables of which the feasible region may be non-rectangle and hard to determine ex-ante. The crucial step is to transform the boundary condition to that of an endogenous state variable which makes the feasible region a simple rectangle. In the current example, the feasible debt level is high so the model can be solved directly over the natural state space—bond and capital. Such transformation becomes essential for more complex problems, and sometimes the boundary condition needs to be solved using a different equilibrium system. The toolbox accommodates this by allowing a different system of equations for problems at the boundary. See Cao and Nie (2017) for such an example, in which consumption is zero and the Inada condition is violated at the boundary.

The policy functions in the current example are highly nonlinear. We learned how to start with a crude grid and refine the solution with finer grids. A more effective approach is to let the toolbox automatically adapt the grids to increase points in regions that are nonlinear. The example in Bianchi (2011) shows how to use adaptive grids for a model with a single-dimension endogenous state. The current model can also be solved with adaptive grids by simply specifying one option in the gmod file, which we leave as an exercise.

Or you can directly proceed to Toolbox API.