值函数迭代模拟弹性劳动供给模型遇RuntimeWarning求助

尝试用值函数迭代模拟弹性劳动供给模型时，出现RuntimeWarning: invalid value encountered in subtract错误，根源在于公式错误、数值边界处理不当以及变量存储逻辑问题，以下是具体排查和修正方案：

核心问题分析 #

• 生产函数公式错误：消费计算中的产出部分误用线性组合替代柯布道格拉斯函数，导致产出值异常，消费c频繁为负，进而值函数中出现大量-inf，与初始值0相减时触发数值警告。
• 值函数边界未处理：V_new最后一个元素始终为0，未对应最大资本点的终端条件，迭代中与更新后的V元素相减时产生无效数值。
• 最优劳动供给存储错误：仅保留最后一个资本点的最优劳动供给，丢失其他点的最优值，导致政策函数完全错误。
• 资本积累方程逻辑错误：资本更新公式不符合宏观模型的资本积累规则，无法正确模拟资本路径。

修正后的代码 #

import numpy as np
import matplotlib.pyplot as plt

# 参数设定
beta = 0.95       # 贴现因子
rho = 0.1         # 闲暇的效用权重
alpha = 0.3       # 资本份额
delta = 0.05      # 折旧率（符合标准宏观模型设定）
k_min = 0.1       # 最小资本（避免k=0导致产出为0）
k_max = 10        # 最大资本
n_min = 0.01      # 最小劳动供给（避免log(1-n)无意义）
n_max = 0.99      # 最大劳动供给
tol = 1e-6        # 收敛容忍度
num_points = 100  # 网格点数

# 资本网格
k_grid = np.linspace(k_min, k_max, num_points)

# 初始化值函数和最优劳动供给数组
V = np.zeros(num_points)
opt_n_array = np.zeros(num_points)  # 存储每个资本点的最优劳动供给

# 值函数迭代
while True:
    V_new = np.zeros(num_points)
    
    for i, k in enumerate(k_grid):
        max_value = -np.inf
        best_n = n_min
        
        # 遍历所有可能的未来资本k'
        for j, k_next in enumerate(k_grid):
            # 遍历劳动供给网格
            for n in np.linspace(n_min, n_max, num_points):
                # 标准柯布道格拉斯产出
                y = k**alpha * n**(1 - alpha)
                # 消费=产出+未折旧资本-未来资本
                c = y + (1 - delta)*k - k_next
                
                if c <= 0 or n <= 0 or n >= 1:
                    continue
                
                current_value = np.log(c) + rho * np.log(1 - n) + beta * V[j]
                if current_value > max_value:
                    max_value = current_value
                    best_n = n
        
        V_new[i] = max_value if max_value != -np.inf else 0
        opt_n_array[i] = best_n
    
    # 收敛判断
    if np.max(np.abs(V - V_new)) < tol:
        break
    V = V_new.copy()  # 使用copy避免数组引用问题

# 计算政策函数
c_policy = np.zeros(num_points)
k_policy = np.zeros(num_points)
for i, k in enumerate(k_grid):
    n = opt_n_array[i]
    y = k**alpha * n**(1 - alpha)
    
    # 找到对应最优未来资本
    max_val = -np.inf
    best_k = k_min
    for j, k_next in enumerate(k_grid):
        c = y + (1-delta)*k - k_next
        if c <= 0:
            continue
        val = np.log(c) + rho*np.log(1-n) + beta*V[j]
        if val > max_val:
            max_val = val
            best_k = k_next
    
    k_policy[i] = best_k
    c_policy[i] = y + (1-delta)*k - best_k

# 绘制政策函数
plt.figure(figsize=(10,6))
plt.plot(k_grid, k_grid, 'k--', label='45度线')
plt.plot(k_grid, c_policy, label='消费政策')
plt.plot(k_grid, opt_n_array, label='劳动供给政策')
plt.xlabel('资本存量')
plt.ylabel('政策函数值')
plt.legend()
plt.title('弹性劳动供给模型政策函数')
plt.show()

# 模拟资本路径
k_path = [5.0]
for t in range(20):
    current_k = k_path[-1]
    # 匹配当前资本对应的网格索引
    idx = np.argmin(np.abs(k_grid - current_k))
    next_k = k_policy[idx]
    k_path.append(next_k)

plt.figure(figsize=(10,6))
plt.plot(k_path)
plt.xlabel('时期')
plt.ylabel('资本存量')
plt.title('资本存量动态路径')
plt.show()

关键修正点说明 #

1. 修正生产函数：将产出计算改为标准柯布道格拉斯形式k**alpha * n**(1-alpha)，同时加入折旧率delta，符合宏观模型的资本积累逻辑。
2. 处理数值边界：将资本和劳动的最小取值设为0.01，避免log(0)或产出为0的情况，减少无效数值。
3. 存储所有最优劳动供给：用数组opt_n_array存储每个资本点对应的最优劳动供给，避免信息丢失。
4. 修正值函数迭代逻辑：遍历所有可能的未来资本点，对每个(k, k_next)组合求解最优劳动供给，确保值函数更新的完整性。
5. 收敛判断时使用数组拷贝：用V = V_new.copy()避免数组引用导致的错误。

内容的提问来源于stack exchange，提问作者Muhammet Ali CANŞI