处理带有非凸性调整成本的模型：Khan and Thomas (2008)

本文是一个带有异质性厂商的实际商业周期模型。在这个模型中，厂商调整资本规模时会面临固定成本，因此投资决策会服从 $(S,s)$ 规则。我们的工具包可以用以求解此类问题的稳态，以及宏观冲击后的线性/非线性转移路径。参照 Boppart, Krusell & Mitman (2018) 和 Auclert et al. (2021)，我们求解的确定性转移路径刻画了宏观冲击的一阶影响，这与采用Winberry (2018)中局部扰动法得到的结果是一致的。

我们参照 Winberry (2018) 对 Khan & Thomas (2008) 中的参数进行了调整，以剔除趋势性。

模型

异质性厂商 模型中的厂商在个体生产率 $\varepsilon$ 以及资本 $k$ 上有异质性。其中生产率的变动服从一个外生的马尔科夫链。当进行投资决策时，如果厂商的投资超出 $[-ak,ak]$，就需要支付一定的固定成本，即 $\xi$ 单位的劳动。如果不支付这个固定成本，厂商下一期的资本规模必须落在 $[(1-\delta+a)k,(1-\delta+b)k]$ 的区间内，其中 $\delta>0$ 表示资本的折旧率。投资的固定成本 $\xi$ 是一个 i.i.d. 的随机变量，服从 $[0,\bar{\xi}]$ 的均匀分布，它的分布函数我们用 $G$ 来表示。

在每一个时间点 $t$，厂商的状态变量包括 $(\varepsilon,k,\xi)$。厂商的最优化问题可以表示为：

\[\begin{split} v_t^1(\varepsilon,k,\xi)=\max_{n,k^*,k^C} \Big[z_t \varepsilon k^{\alpha}n^{\nu}-w_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$ 代表总体的生产率水平，$w_t$ 代表工资率，$r_t(\varepsilon,k')$ 表示在下一期拥有资本 $k'$ 的厂商的存续价值。每一个厂商会选择最优的劳动需求、资本规模。在决定最优资本规模时，厂商会比较支付投资调整成本和不支付时的公司价值，并选择是否支付投资调整成本以最大化公司的存续价值。公司的存续价值可以表示为：

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

其中，$d_t$ 是与家庭最优化决策相一致的内生贴现率；$v_{t+1}^0$ 是考虑过未来调整成本之后公司的预期价值。

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

公司是否支付投资调整成本取决于调整成本的相对大小：当自由选择资本规模的净回报高于需要支付的调整成本时，公司会选择支付调整成本：

\[ \xi \leq \xi^*\equiv \frac{r_t(\varepsilon,k^*)-r_t(\varepsilon,k^C)}{\omega_t} \]

因此，公司愿意支付调整成本的概率可以表示为 $g_{a,t}(\varepsilon,k)$。我们用 $\Gamma_t$ 来表示公司在状态变量 $(\varepsilon,k)$ 上的分布。

代表性家庭 的效用函数可以表示为 $$ \sum_{t=0}^{\infty} \beta^t \Big[\frac{C_t^{1-\sigma}}{1-\sigma}-\chi N_t\Big] $$

家庭的最优化条件可以推导出贴现因子和工资率：

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

市场出清条件

劳动力市场出清条件如下：

\[ N_t = \underbrace{\int g_{n,t}(\varepsilon,k) d \Gamma_t}_{\text{用于生产的劳动力}}+ \underbrace{\iint_0^{\xi_t^*(\varepsilon,k)} \xi G(d\xi)d\Gamma_t}_{\text{资本调整成本中的劳动力}} \]

商品市场出清条件如下：

\[ 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 \]

序贯竞争性均衡 给定厂商的初始分布 $\Gamma_0$, 序贯竞争性均衡由以下条件刻画： (1) 厂商的价值函数和决策规则 $\{v_t^1,v_t^0,r_t,\xi_t^*,g_{k^*,t},g_{k^C,t}\}$；(2) 总体数量与价格 $\{C_t,N_t,\omega_t,d_t\}$；(3) 厂商个体状态变量的联合分布 $(\varepsilon,k)$, $\{\Gamma_t\}$，其中：

厂商的价值函数和决策规则满足贝尔曼方程；
$\{C_t,N_t,\omega_t,d_t\}$ 满足家庭单位最优化条件和市场出清条件；
$\{\Gamma_t\}$ 符合 $\varepsilon$的外生转移矩阵和决策规则$\{\xi_t^*,g_{k^*,t},g_{k^C,t}\}$. $\Gamma_t$的迭代规则可以由如下的方程刻画：

\[ \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{博雷尔集} (\mathcal{E}',\mathcal{K}') \]

其中，$P(\mathcal{E}'|\varepsilon)$ 是 $\varepsilon$ 的条件转移概率。

问题的求解

厂商的最优化问题可以做如下的简化。首先，企业的劳动力决策不依赖于投资的决策，因此我们可以单独写出劳动力的决策过程：

\[ 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\} \]

其次，由于我们假设 $\xi$ 的抽样来自于独立的均匀分布，我们可以利用决策的临界规则来算出期望价值：

\[\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}\]

其中，$\pi_t(\varepsilon,k)\equiv(1-\delta)k+\max_{n} (z_t \varepsilon k^{\alpha}n^{\nu} -\omega_t n)$ 是最大化的利润，$R_t^a(\varepsilon,k)\equiv\max_{k^*\geq0}r_t(\varepsilon,k^*)$。$R_t^c(\varepsilon,k)\equiv\max_{k^C\in [(1-\delta-a)k,(1-\delta+a)k]}r_t(\varepsilon,k^C)$。阈值$\xi_t^*(\varepsilon,k)$ 满足：

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

给定上述的求解过程，我们可以将模型输入工具包进行求解。

hmod 文件

原文的模型可以由 KT2008.hmod表示，具体如下：

parameters beta sigma alpha nu delta xi_max a w z d;
% The parameters are taken from Khan and Thomas (2008)
beta = 0.961; % subjective discount factor 
sigma = 1.0;  % inverse of EIS
alpha = 0.256; % capital elasticity of production function
nu = 0.64;  % labor share
delta = 0.085; % capital depreciation rate
xi_max = 0.0083; % upper bound on adjustment cost
a = 0.011; % no-adjustment cost region
w = 1.0;    % wage rate
z = 1.0;    % aggregate productivity
d = beta;   % discount factor

var_shock e;
% idiosyncratic productivity log(e') = 0.859*log(e) + N(0,0.022^2)
% discretized as a 15 State Markov Chain using Rouwenhorst Method 
setup_shock_process;

var_state k;
k_min = 1e-3;
k_max = 4;
nkgrid = 200;
k = exp(linspace(log(k_min), log(k_max + k_min), nkgrid)) - k_min;

var_policy ka kc;
initial ka k;
initial kc (1-delta)*k;

var_aux p nd kp y c;

vfi;
    % labor demand from accounting profit max
    nd_input = (nu*exp(z)*exp(e)*(k^alpha)/w)^(1/(1-nu));
    y = exp(z)*exp(e)*(k^alpha)*(nd_input^nu);
    
    % fixed part of Bellman
    Tv_fixed = y - w*nd_input + (1-delta)*k; 

    Ra = -ka + d*EXPECT(v(ka));   % objective for unconstrained firm
    Rc = -kc + d*EXPECT(v(kc));   % objective for constrained firm
    objective = Ra + Rc;          % multi objective optimimzation 

    ka >= 0.0;
    kc >= (1 - delta - a)*k;    % adjustment constrained
    kc <= (1 - delta + a)*k;

    xi_tilde = (Ra - Rc)/w; % adjustment threshold before considering domain
    xi_star = min(max(xi_tilde,0.0),xi_max); % adjusting for domain
    p = xi_star/xi_max; % adjustment probability
    nd_adjust = xi_star*xi_star/(2*xi_max); % adjustment cost in labor unit
    
    % Bellman equation update
    Tv = Tv_fixed + p*Ra - w*nd_adjust + (1-p)*Rc;
    
    kp = p*ka + (1-p)*kc;       % average capital holding for the next period
    nd = nd_input + nd_adjust;  % firm's gross labor demand
    c = y + (1-delta)*k - kp;   % firm's sale to consumers
    
    % Stochastic transition
    k(+1) = {ka, kc} {p, 1-p};
end;

var_agg C;
C = 0.4135;        % initial
chi = 2.32531;
n_star = 1/3;   % calibration target

var_agg_shock z;
z = 0.0;        % steady state value

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

model;
    % update parameters that enter vfi
    d = beta*(C(+1)/C)^(-sigma);
    w = chi / C^(-sigma);       % implied by FOC for labor
    
    % equations
    C + kp - (1-delta)*k == y;  % goods demand = goods supply
    
    % post evaluation
    N = nd; % aggregate labor demand
end;

我们结合这个案例，介绍工具包的几个功能。首先，我们的工具包可以用于同时求解多个最优化问题，只需把目标函数 “objective” 定义为每一个最优化问题之和即可：

vfi;
    % labor demand from accounting profit max
    nd_input = (nu*exp(z)*exp(e)*(k^alpha)/w)^(1/(1-nu));
    y = exp(z)*exp(e)*(k^alpha)*(nd_input^nu);
    
    % fixed part of Bellman
    Tv_fixed = y - w*nd_input + (1-delta)*k; 

    Ra = -ka + d*EXPECT(v(ka));   % objective for unconstrained firm
    Rc = -kc + d*EXPECT(v(kc));   % objective for constrained firm
    objective = Ra + Rc;          % multi objective optimimzation 

    ka >= 0.0;
    kc >= (1 - delta - a)*k;    % adjustment constrained
    kc <= (1 - delta + a)*k;

    xi_tilde = (Ra - Rc)/w; % adjustment threshold before considering domain
    xi_star = min(max(xi_tilde,0.0),xi_max); % adjusting for domain
    p = xi_star/xi_max; % adjustment probability
    nd_adjust = xi_star*xi_star/(2*xi_max); % adjustment cost in labor unit
    
    % Bellman equation update
    Tv = Tv_fixed + p*Ra - w*nd_adjust + (1-p)*Rc;
    
    kp = p*ka + (1-p)*kc;       % average capital holding for the next period
    nd = nd_input + nd_adjust;  % firm's gross labor demand
    c = y + (1-delta)*k - kp;   % firm's sale to consumers
    
    % Stochastic transition
    k(+1) = {ka, kc} {p, 1-p};
end;

在这里，受约束和不受约束（支付调整成本）的问题是相互独立的，因此我们把这两个问题的目标函数 “objective” 分别定义为 Ra 和 Rc。

贝尔曼方程的迭代值 “Tv” 不一定要和定义的目标函数 “objective” 保持一致。

接下来，我们输入内生状态变量的转移过程。

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

这里的语法是: var_state_name(+1) = {value_1, value_2, …} {probability_1, probability_2, …}.

接下来，我们可以在 model_modname (var_agg_1, var_agg_2,…)的模块中，引入与 model 模块不完全一样的方程组。

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

在这里，我们额外定义了一些用于校准的方程。例如，我们在方程组中引入了一个新的变量$\chi$以及相应的方程，使得劳动总需求等于总供给。我们的工具包可以产生一个叫做 solve_modname.m 的文件，里面包含了这些用于校准的方程组，并且会在求解时调用。

使用工具箱

在设置好输入的文本后，我们的工具包会自动产生Matlab文件，包括solve_vfi.m, solve_ss.m, solve_trans_linear.m, solve_trans_nonlinear.m，solve_cali.m，以及其他用于求解稳态、稳定分布、转移矩阵以及其他功能（例如校准）的方程。这些生成的文件包括：

Untitled

定位于 KT2008.hmod，在 MATLAB 中调用 hans 编译脚本文件。

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'}

我们首先求解校准模块。

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

我们随后通过传递校准后的参数求解稳态。

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, 

让我们检查求解出的稳态总体变量。

ss_rslt.var_agg
ans = struct with fields:
    C: 0.4135
    z: 0
    d: 0.9610
    w: 0.9615
    N: 0.3333

让我们检查策略函数。承担调整成本的概率：

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')

如上所示，投资决策呈现出 (S, s) 规则：当资本水平与预期水平相比过低或过高时，就会产生调整成本。

平稳分布：

plot(ss_rslt.vfi_rslt.var_state.k,ss_rslt.dist); title('Histogram of Capital');

注意，直方图的锯齿状形状是一个特征。可以通过在 Winberry (2018) 求解稳态的代码中增加资本格点数量来复现这样的模式。

我们接下来生成总体生产率冲击的时间序列。

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')

我们将总体冲击传入模型并且求解转移路径。

注意，solve_trans_nonlinear 包含三个输入参数，第一个是起始稳态，第二个是最终稳态。

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.

让我们检查转移路径。

plot(trans_rslt.var_agg_t.p); title('Probability of Adjustment');

plot(trans_rslt.var_agg_t.y); title('Aggregate Output');

如上所示，在一个正向的生产率冲击发生后，承担固定调整成本的概率上升，这表明调整投资的重要份额处于广度边际上。

非凸调整成本能够生成总体非线性吗？

让我们提高冲击的规模，再次求解转移路径。

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');

如图所示，当冲击规模增加10倍时，内生变量的反应也会以大致相同的因子增加。

这表明尽管存在微观层面上的非凸调整成本，但模型并未呈现出总体非线性特征，这与 Khan & Thomas（2008）的研究结果一致。

下一步

理解算法的细节。

熟悉 Toolbox API reference。

更多的案例:

Krusell and Smith (1998)：一个最基础的模型。
McKay, Nakamura & Steinsson (2016) 一个带有异质性个体的新凯恩斯模型，需要求解复杂的高度非线性的均衡。
一个带有两种资产的HANK模型，需要求解资产配置决策。