Handling Non-convex Adjustment Costs: Khan and Thomas (2008)

This is a real business cycle model with heterogeneous firms. Firms face fixed costs when making capital adjustments, which generates an \((S,s)\) style investment rule. The toolbox solves the steady state and the linear/non-linear transition path after an unexpected aggregate shock to productivity. Following Boppart, Krusell and Mitman (2018) and Auclert et al. (2021), the deterministic transition path characterizes the first order effect of the aggregate shock.

We use the parameterization in Winberry (2018) which adjusts the parameter values in Khan and Thomas (2008) for no trend growth. We follow the notation in Khan and Thomas (2008).

Model

Heterogeneous firms differ in their idiosyncratic productivity, \(\varepsilon\), which follows an exogenous Markov chain, and capital holding \(k\). When making investment decisions, the firm can incur a fixed cost of \(\xi\) in units of labor, so that its choice of investment becomes unconstrained; otherwise, its choice of investment needs to be constrained in an interval \([-ak,ak]\), leaving next period capital choice constrained in the interval \([(1-\delta+a)k,(1-\delta+b)k]\), with \(\delta>0\) the depreciation rate. \(\xi\) is drawn i.i.d. from the uniform distribution over \([0,\bar{\xi}]\), with the distribution denoted by \(G\).

At time \(t\), a firm with realized \((\varepsilon,k,\xi)\) solves the following problem:

\[\begin{split} v_t^1(\varepsilon,k,\xi)=\max_{n,k^*,k^C} \Big[z_t \varepsilon k^{\alpha}n^{\nu}-\omega_t n+(1-\delta)k+\max\{-\xi \omega_t+r_t(\varepsilon,k^*),r_t(\varepsilon,k^C)\}\Big] \\ s.t. \quad n\geq 0, k^* \geq 0 \\ k^C\in [(1-\delta-a)k,(1-\delta+a)k]. \end{split}\]

\(z_t\) is the aggregate productivity; \(\omega_t\) is the wage rate; \(r_t(\varepsilon,k')\) is the continuation value associated with a future capital level \(k'\), defined below. Firms choose labor demand \(n\), unconstrained capital adjustment level \(k^*\), and constrained capital level \(k^C\). The \(\max\) operator in the objective compares the continuation value under unconstrained or constrained adjustment, accounting for the fixed cost \(\xi \omega_t\) that needs to be incurred to unlock the unconstrained adjustment. The continuation value is given by

\[ r_t(\varepsilon,k') = -k' + d_t\cdot\mathbb{E}[v^0_{t+1}(\varepsilon',k')|\varepsilon], \]

where \(d_t\) is the discount factor inherited from the households, and \(v_{t+1}^0\) is the expected value integrating over future adjustment cost \(\xi'\):

\[ v_{t+1}^0(\varepsilon',k')=\int_0^{\bar{\xi}} v_{t+1}^1(\varepsilon',k',\xi') G (d\xi'). \]

Firms’ decision for whether to incur the adjustment cost can be characterized by a threshold rule: a firm will incur the adjustment cost if and only if the benefit of unconstrained adjusting outweighs the fixed cost, i.e., with a draw of \(\xi\) that satisfies

\[ \xi \leq \xi^*_t(\varepsilon,k)\equiv \max\{\min\{\frac{r_t(\varepsilon,g_{k^*,t}(\varepsilon,k))-r_t(\varepsilon,g_{k^C,t}(\varepsilon,k))}{\omega_t},\bar{\xi}\},0\}, \]

which applies that \(\xi\) is drawn from a uniform distribution over \([0,\bar{\xi}]\). \(g_{n,t},g_{k^*,t},g_{k^C,t}\) are policy functions for \(n,k^*,k^C\) respectively.

Denote \(\Gamma_t\) the measure that represents the distribution over firms endogenous states \((\varepsilon,k)\).

Representative households value consumption and leisure, and maximize the lifetime utility

\[ \sum_{t=0}^{\infty} \beta^t \Big[\frac{C_t^{1-\sigma}-1}{1-\sigma}-\chi N_t\Big]. \]

Households’ optimality gives the conditions:

\[\begin{split} d_t = \beta \Big[\frac{C_{t+1}}{C_t}\Big]^{-\sigma} \\ \omega_t C_t^{-\sigma}=\chi. \end{split}\]

Market clearing conditions for labor:

\[ N_t = \underbrace{\int g_{n,t}(\varepsilon,k) d \Gamma_t}_{\text{labor used in production}}+ \underbrace{\iint_0^{\xi_t^*(\varepsilon,k)} \xi G(d\xi)d\Gamma_t}_{\text{labor used in capital adjustment cost}}, \]

and for goods:

\[ C_t=\int \{z_t \varepsilon k^{\alpha}[g_{n,t}(\varepsilon,k)]^{\nu} + (1-\delta)k\} d \Gamma_t - \int_0^{\xi_t^*(\varepsilon.k)} g_{k^*,t}(\varepsilon,k) d\Gamma_t-\int_{\xi_t^*(\varepsilon.k)}^{\bar{\xi}} g_{k^C,t}(\varepsilon,k) d\Gamma_t. \]

A sequential competitive equilibrium, given an initial distribution over firms \(\Gamma_0\), is a sequence of (1) firms’ value and policy functions \(\{v_t^1,v_t^0,r_t,\xi_t^*,g_{k^*,t},g_{k^C,t}\}\); (2) aggregate quantities and prices \(\{C_t,N_t,\omega_t,d_t\}\); (3) distribution over firms states \((\varepsilon,k)\), \(\{\Gamma_t\}\), such that

  • value and policy functions solve firms’ Bellman equation

  • \(\{C_t,N_t,\omega_t,d_t\}\) satisfy the four set of equations that characterize households’ optimality and market clearing conditions, listed above

  • \(\{\Gamma_t\}\) are consistent with the exogenous transition of \(\varepsilon\) and policy functions \(\{\xi_t^*,g_{k^*,t},g_{k^C,t}\}\). We write out the transition of \(\Gamma_t\) to highlight the stochastic transition of states due to idiosyncratic draws of \(\xi\):

\[ \Gamma_{t+1}(\mathcal{E}',\mathcal{K}')=\int \left( \int_0^{\xi_t^*(\varepsilon,k)}\mathcal{I}(g_{k^*,t}(\varepsilon,k)\in \mathcal{K}')G(d\xi)+ \int_{\xi_t^*(\varepsilon,k)}^{\bar{\xi}}\mathcal{I}(g_{k^C,t}(\varepsilon,k)\in \mathcal{K}')G(d\xi) \right)P(\mathcal{E}'|\varepsilon)\Gamma_t(d \varepsilon, dk), \forall \text{ Boreal set } (\mathcal{E}',\mathcal{K}') \]

where \(P(\mathcal{E}'|\varepsilon)\) is the transition measure of \(\varepsilon\).

Characterizations

Several simplifications of firms’ problem in order. First, the choice of labor \(n\) is independent of the investment decision, and the constrained and unconstrained investment problems are independent of each other. The decision problem thus can be decoupled and written as

\[ v_t^1(\varepsilon,k,\xi)=(1-\delta)k+\max_{n} (z_t \varepsilon k^{\alpha}n^{\nu} -\omega_t n) +\max\Big\{-\xi \omega_t+ \max_{k^*\geq0}r_t(\varepsilon,k^*), \max_{k^C\in [(1-\delta-a)k,(1-\delta+a)k]} r_t(\varepsilon,k^C)\Big\}. \]

Next, with the assumption that \(\xi\) is drawn from the uniform distribution, we can integrate over the expected value under the cutoff rule:

\[\begin{split} v_t^0(\varepsilon,k)=\int_0^{\bar{\xi}} v_t^1(\varepsilon,k,\xi)G(d\xi) \\ =\pi_t(\varepsilon,k)+\int_0^{\xi_t^*(\varepsilon,k)} [R_t^a(\varepsilon,k)-\xi \omega_t]\frac{1}{\bar{\xi}} d\xi+\int_{\xi_t^*(\varepsilon,k)}^{\bar{\xi}}R_t^c(\varepsilon,k)\frac{1}{\bar{\xi}} d\xi \\ =\pi_t(\varepsilon,k)+\frac{\xi_t^*(\varepsilon,k)}{\overline{\xi}}R_t^a(\varepsilon,k)-\frac{\omega_t}{2\bar{\xi}}[\xi_t^*(\varepsilon,k)]^2+ \left(1-\frac{\xi_t^*(\varepsilon,k)}{\overline{\xi}}\right)R_t^c(\varepsilon,k) \end{split}\]

where \(\pi_t(\varepsilon,k)\equiv(1-\delta)k+\max_{n} (z_t \varepsilon k^{\alpha}n^{\nu} -\omega_t n)\) is the maximized profit, \(R_t^a(\varepsilon,k)\equiv\max_{k^*\geq0}r_t(\varepsilon,k^*)\), and \(R_t^c(\varepsilon,k)\equiv\max_{k^C\in [(1-\delta-a)k,(1-\delta+a)k]}r_t(\varepsilon,k^C)\). The threshold \(\xi_t^*(\varepsilon,k)\) satisfies

\[ \xi^*_t(\varepsilon,k) = \max\{\min\{\frac{R_t^a(\varepsilon,k)-R_t^c(\varepsilon,k)}{\omega_t},\bar{\xi}\},0\}. \]

With these characterizations ready, we are able to input the model into the toolbox script file.

The hmod File

The model can be represented using KT2008.hmod, listed below.

 1parameters beta sigma alpha nu delta xi_max a w z d;
 2% The parameters are taken from Khan and Thomas (2008)
 3beta = 0.961; % subjective discount factor 
 4sigma = 1.0;  % inverse of EIS
 5alpha = 0.256; % capital elasticity of production function
 6nu = 0.64;  % labor share
 7delta = 0.085; % capital depreciation rate
 8xi_max = 0.0083; % upper bound on adjustment cost
 9a = 0.011; % no-adjustment cost region
10w = 1.0;    % wage rate
11z = 1.0;    % aggregate productivity
12d = beta;   % discount factor
13
14var_shock e;
15% idiosyncratic productivity log(e') = 0.859*log(e) + N(0,0.022^2)
16% discretized as a 15 State Markov Chain using Rouwenhorst Method 
17setup_shock_process;
18
19var_state k;
20k_min = 1e-3;
21k_max = 4;
22nkgrid = 200;
23k = exp(linspace(log(k_min), log(k_max + k_min), nkgrid)) - k_min;
24
25var_policy ka kc;
26initial ka k;
27initial kc (1-delta)*k;
28
29var_aux p nd kp y c;
30
31vfi;
32    % labor demand from accounting profit max
33    nd_input = (nu*exp(z)*exp(e)*(k^alpha)/w)^(1/(1-nu));
34    y = exp(z)*exp(e)*(k^alpha)*(nd_input^nu);
35    
36    % fixed part of Bellman
37    Tv_fixed = y - w*nd_input + (1-delta)*k; 
38
39    Ra = -ka + d*EXPECT(v(ka));   % objective for unconstrained firm
40    Rc = -kc + d*EXPECT(v(kc));   % objective for constrained firm
41    objective = Ra + Rc;          % multi objective optimimzation 
42
43    ka >= 0.0;
44    kc >= (1 - delta - a)*k;    % adjustment constrained
45    kc <= (1 - delta + a)*k;
46
47    xi_tilde = (Ra - Rc)/w; % adjustment threshold before considering domain
48    xi_star = min(max(xi_tilde,0.0),xi_max); % adjusting for domain
49    p = xi_star/xi_max; % adjustment probability
50    nd_adjust = xi_star*xi_star/(2*xi_max); % adjustment cost in labor unit
51    
52    % Bellman equation update
53    Tv = Tv_fixed + p*Ra - w*nd_adjust + (1-p)*Rc;
54    
55    kp = p*ka + (1-p)*kc;       % average capital holding for the next period
56    nd = nd_input + nd_adjust;  % firm's gross labor demand
57    c = y + (1-delta)*k - kp;   % firm's sale to consumers
58    
59    % Stochastic transition
60    k(+1) = {ka, kc} {p, 1-p};
61end;
62
63var_agg C;
64C = 0.4135;        % initial
65chi = 2.32531;
66n_star = 1/3;   % calibration target
67
68var_agg_shock z;
69z = 0.0;        % steady state value
70
71model_cali(C, chi);
72    % update parameters that enter vfi
73    w = chi / C^(-sigma);       % implied by FOC for labor
74    
75    % equations
76    C + kp - (1-delta)*k == y;  % goods demand = goods supply
77    nd == n_star;               % labor demand = cali labor supply
78end;
79
80model;
81    % update parameters that enter vfi
82    d = beta*(C(+1)/C)^(-sigma);
83    w = chi / C^(-sigma);       % implied by FOC for labor
84    
85    % equations
86    C + kp - (1-delta)*k == y;  % goods demand = goods supply
87    
88    % post evaluation
89    N = nd; % aggregate labor demand
90end;

We highlight several new features of the toolbox demonstrated by this example. First, the toolbox handles multi-object optimization by defining the “objective” as the sums of each problem:

31vfi;
32    % labor demand from accounting profit max
33    nd_input = (nu*exp(z)*exp(e)*(k^alpha)/w)^(1/(1-nu));
34    y = exp(z)*exp(e)*(k^alpha)*(nd_input^nu);
35    
36    % fixed part of Bellman
37    Tv_fixed = y - w*nd_input + (1-delta)*k; 
38
39    Ra = -ka + d*EXPECT(v(ka));   % objective for unconstrained firm
40    Rc = -kc + d*EXPECT(v(kc));   % objective for constrained firm
41    objective = Ra + Rc;          % multi objective optimimzation 
42
43    ka >= 0.0;
44    kc >= (1 - delta - a)*k;    % adjustment constrained
45    kc <= (1 - delta + a)*k;
46
47    xi_tilde = (Ra - Rc)/w; % adjustment threshold before considering domain
48    xi_star = min(max(xi_tilde,0.0),xi_max); % adjusting for domain
49    p = xi_star/xi_max; % adjustment probability
50    nd_adjust = xi_star*xi_star/(2*xi_max); % adjustment cost in labor unit
51    
52    % Bellman equation update
53    Tv = Tv_fixed + p*Ra - w*nd_adjust + (1-p)*Rc;
54    
55    kp = p*ka + (1-p)*kc;       % average capital holding for the next period
56    nd = nd_input + nd_adjust;  % firm's gross labor demand
57    c = y + (1-delta)*k - kp;   % firm's sale to consumers
58    
59    % Stochastic transition
60    k(+1) = {ka, kc} {p, 1-p};
61end;

Here, the problems under constrained and unconstrained adjustment are indepenent, with objectives specified in Ra and Rc respectively.

The updated value in the Bellman equation, “Tv”, needs not agree with “objective”.

Then, we can specify the stochastic transition of endogenous states, here, due to idiosyncratic investment cost draws.

60    k(+1) = {ka, kc} {p, 1-p};

The syntax is as follows: var_state_name(+1) = {value_1, value_2, …} {probability_1, probability_2, …}.

Next, we can have a different system of equations than the main model block, enclosed in a model_modname(var_agg_1, var_agg_2,…) block:

71model_cali(C, chi);
72    % update parameters that enter vfi
73    w = chi / C^(-sigma);       % implied by FOC for labor
74    
75    % equations
76    C + kp - (1-delta)*k == y;  % goods demand = goods supply
77    nd == n_star;               % labor demand = cali labor supply
78end;

Here, we define a system of equations for calibration, which adds from the equilibrium condition one extra unknown \(\chi\) and the associated extra equation, total labor demand = targetd labor supply. The toolbox will generate a solve_modname.m script for each alternative model block, which can be called for solving each alternative system defined.

Use the Toolbox

After parsing the script file, the toolbox generates MATLAB files: solve_vfi.m, solve_ss.m, solve_trans_linear.m, solve_trans_nonlinear.m, solve_cali.m and other functions that can be called to solve the steady state, stationary distribution, transition paths, and alternative system (here, calibration) of the model. The usages of the generated files are illustrated below.

Untitled
Locate KT2008.hmod, call hans in MATLAB to compile the script file.
hans KT2008
Parsing... Symbolic parser failed. Switching to AD parser...ok! Compiling dpopt_vfi Building with 'MinGW64 Compiler (C++)'. MEX completed successfully.
ans = struct with fields:
parameters: {'beta' 'sigma' 'alpha' 'nu' 'delta' 'xi_max' 'a' 'w' 'z' 'd'} var_shock: 'e' var_state: 'k' var_pre_vfi: {1×0 cell} var_policy: {'ka' 'kc'} var_aux: {'p' 'nd' 'kp' 'y' 'c'} var_agg: 'C' var_agg_shock: 'z' var_agg_params: {'chi' 'delta' 'sigma' 'beta'} var_agg_assign: {'d' 'w'} var_agg_in_ind_params: {'d' 'w' 'z'} var_ind_in_eq: {'y' 'nd' 'kp' 'k'}
We first solve the calibration block.
tic;
cali_rslt = solve_cali;
Evaluating at var_agg: 0.4135, 2.32531, VFI converged (metric_v) 9.85215e-09, (metric_pol) 4.29862e-07 in 390 iterations Range of k: 0, 3.68 Equilbirium residual: 0.000197203, -0.000149756, Iteration Line Search Func-count ||f(x)||^2 fnorm(x) Norm of Step Step Size Algorithm 0 0 1 6.131594e-08 2.476205e-04 Evaluating at var_agg: 0.4136, 2.32531, VFI converged (metric_v) 9.64107e-09, (metric_pol) 2.48728e-07 in 227 iterations Range of k: 0, 3.68 Equilbirium residual: 0.000959078, -0.000766295, Evaluating at var_agg: 0.4135, 2.32554, VFI converged (metric_v) 9.86429e-09, (metric_pol) 3.76106e-07 in 213 iterations Range of k: 0, 3.68 Equilbirium residual: 0.000470949, -0.000404774, Evaluating at var_agg: 0.413464, 2.32538, VFI converged (metric_v) 9.81415e-09, (metric_pol) 7.01103e-08 in 216 iterations Range of k: 0, 3.68 Equilbirium residual: 3.96032e-08, -1.34074e-09, 1 1 4 1.570212e-15 3.962590e-08 7.735502e-05 1.000000e+00 Newton Evaluating at var_agg: 0.413464, 2.32538, VFI result reused Range of k: 0, 3.68 Equilbirium residual: 3.96032e-08, -1.34074e-09,
toc;
Elapsed time is 1.464644 seconds.
cali_rslt.var_agg
ans = struct with fields:
C: 0.4135 chi: 2.3254 z: 0 w: 0.9615
We then solve the steady state, by passing in the calibrated parameters.
ss_rslt = solve_ss(cali_rslt);
Evaluating at var_agg: 0.413464, VFI result reused Range of k: 0, 3.68 Equilbirium residual: 3.97218e-08, Iteration Line Search Func-count ||f(x)||^2 fnorm(x) Norm of Step Step Size Algorithm 0 0 1 1.577820e-15 3.972178e-08 Evaluating at var_agg: 0.413464, VFI result reused Range of k: 0, 3.68 Equilbirium residual: 3.97218e-08,
Let's inspect the solved aggregate variables at the steady state
ss_rslt.var_agg
ans = struct with fields:
C: 0.4135 z: 0 d: 0.9610 w: 0.9615 N: 0.3333
Let's inspect the policy functions. The probability for incuring the adjustment cost:
mesh(ss_rslt.vfi_rslt.var_state.k, ss_rslt.vfi_rslt.var_shock.e, ss_rslt.vfi_rslt.var_aux.p);
ylabel('log(varepsilon)');
xlabel('k');
title('Probability of Incurring the Adjustment Cost')
As shown, the investment decision exhibits an (S,s) rule: the adjustment cost is incurred when capital level is either too low or too high compared to the desired level.
The stationary distribution
plot(ss_rslt.vfi_rslt.var_state.k,ss_rslt.dist); title('Histogram of Capital');
Note that the zig-zag shape of the histogram is a feature. One can replicate such patterns by increasing the number of capital grids in the code for solving steady states in Winberry (2018).
We next generate a time series of aggregate productivity shock.
rho_z = 0.9;
size_z_shock = 0.01;
T = 200;
z_t = size_z_shock*rho_z.^[0:T-1];
plot(z_t); title('z'); xlabel('T')
We feed the aggregate shock into the model and solve the transition path.
Note solve_trans_nonlinear takes three arguments, the first for the initial steady state and the second for the final steady state.
options = struct;
options.z_t = z_t;
tic;
trans_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 3.386038e-07 5.818968e-04 1 1 2 2.281655e-09 4.776667e-05 5.418998e-05 1.000000e+00 Broyden
toc;
Elapsed time is 10.657985 seconds.
Let's inspect the transition path.
plot(trans_rslt.var_agg_t.p); title('Probability of Adjustment');
plot(trans_rslt.var_agg_t.y); title('Aggregate Output');
As shown, the probability of incurring the fixed adjustment cost rises after a positive productivity shock, showing that an important share of investment adjustment is at the extensive margin.
Does non-convex adjustment cost generate aggregate non-linearity?
Let's increase the size of the shock, and resolve the transition path.
options.z_t = 10*z_t;
trans_large_shock = 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 3.372703e-03 5.807498e-02 1 1 2 1.824057e-04 1.350577e-02 5.670312e-03 1.000000e+00 Broyden 2 1 3 3.383330e-06 1.839383e-03 4.916697e-04 1.000000e+00 Broyden 3 1 4 9.411725e-08 3.067854e-04 2.803021e-05 1.000000e+00 Broyden 4 1 5 1.220428e-10 1.104730e-05 3.834561e-06 1.000000e+00 Broyden
plot(trans_large_shock.var_agg_t.p); title('Probability of Adjustment');
plot(trans_large_shock.var_agg_t.y); title('Aggregate Output');
As shown, increasing the size of the shock by 10 times leads to responses of endogenous variables to increase by roughly the same factor.
This shows that the model does not exhibit aggregate nonlinearity despite micro-level non-covex adjustment cost, echoing findings by Khan and Thomas (2008).

What’s Next?

Behind the scene: understand the details of the Algorithm.

Check the Toolbox API reference.

More examples: