随机吃蛋糕问题MATLAB实现：值函数与策略函数绘图异常求助

问题描述 #

在MATLAB中实现动态随机吃蛋糕问题，目标是计算值函数与策略函数并绘制对应图像。代码无语法错误，但绘图结果为恒定值的水平线，而非预期曲线。调试发现内循环运行正常，但外循环中V(i_s,:)和pol_fun(i_s,:)的值始终恒定，导致绘图异常。

原始代码 #

% taste schock values
y_L = 0.75;
y_H = 1.25;

% Parameters
sigma = 2;
beta = 0.95;
n_s = 100;
n_c = 120;
n_y = 140;

% CRRA utility function
u = @(c) (c.^(1-sigma))/(1-sigma);

% Transition matrix
Pi = [0.95 0.05; 0.1 0.9];

% Grid for cake size
X = linspace(0, 1, n_s);

% initializing values for consumption
X_L = 0.5;

% Grid for taste shocks
Y = linspace(y_L, y_H, n_y);
pol_fun=zeros(n_s, n_y);

% Initializing value function
V = zeros(n_s, n_y);

% Value function iteration
for i_s = 1:n_s

    aux = zeros(1, n_c);

    c_L = X_L; %min value for consumption
    c_H = X(i_s); %max value for consumption
    c_grid = linspace(c_L, c_H, n_c);
    for i_c = 1:n_c
         EnextV = zeros(1, n_y);
         for i_y = 1:n_y
             nextX = beta * (X(i_s) - c_grid(i_c)) + Y(i_y);
             if nextX >= 0 && nextX <= 1
                nextV = interp1(X, V(:,i_y), nextX, 'spline');
             else
                nextV = -Inf;
            end
            EnextV(i_y) = nextV;
        end
        aux(i_c) = u(c_grid(i_c)) + beta * sum(max(EnextV) .* Pi(:,2)');
    end
    [V(i_s,:), pol_fun(i_s,:)] = max(aux); % current state
end

% Plotting value function
figure;
hold on;
for i_s = 1:n_s
    plot(Y, V(i_s,:), 'LineWidth', 2);
end
xlabel('Taste shock')
ylabel('Value function')
title('Value function')
hold off;

% Plotting policy function
figure;
hold on;
for i_s = 1:n_s
    plot(Y, c_L + (c_H - c_L)/(n_c-1) * (pol_fun(i_s,:) - 1), 'LineWidth', 2);
end
xlabel('Taste shock')
ylabel('Consumption')
title('Policy function')
hold off;

核心问题分析 #

1. 仅执行一轮值函数迭代：值函数迭代需要多轮循环直到收敛，原始代码只遍历一次蛋糕网格X。初始V全为0，第一轮计算时所有未来值均为0，导致aux的计算结果不随冲击状态i_y变化，最终max(aux)返回单一值，使得V(i_s,:)全为常数。
2. 期望计算逻辑错误：sum(max(EnextV) .* Pi(:,2)')完全不符合随机动态规划的期望计算规则——应该对当前冲击状态下的所有未来冲击状态值函数，用转移矩阵加权求和，而非取EnextV的最大值后再乘转移矩阵列向量。
3. 冲击设定矛盾：Y是连续的140个点的网格，但Pi是2×2的离散转移矩阵，两者维度不匹配，无法正确计算跨期期望。
4. 消费网格无效：c_L固定为0.5，当X(i_s) < 0.5时，linspace(c_L, c_H, n_c)会生成反向的消费网格，导致后续计算出现异常值。
5. 状态维度匹配错误：值函数V的维度为n_s × n_y，对应蛋糕状态和冲击状态，但原始代码计算未来值时未正确匹配冲击状态的转移关系。

修正后的代码实现 #

% 冲击状态设定（离散两状态，匹配转移矩阵）
y_vals = [0.75, 1.25];
n_y = length(y_vals);

% 参数定义
sigma = 2;
beta = 0.95;       % 贴现因子
delta = 0.95;      % 蛋糕折旧率（与贴现因子分开，更符合模型逻辑）
n_s = 100;         % 蛋糕存量网格数
n_c = 120;         % 消费网格数
tol = 1e-6;        % 收敛阈值
max_iter = 1000;   % 最大迭代次数

% CRRA效用函数（消费为0时设为负无穷）
u = @(c) (c.^(1-sigma))/(1-sigma);
u(0) = -Inf;

% 转移矩阵：Pi(i,j)表示从冲击状态i转移到j的概率
Pi = [0.95 0.05; 0.1 0.9];

% 蛋糕存量网格
X = linspace(0, 1, n_s);

% 初始化值函数与策略函数
V = zeros(n_s, n_y);
pol_fun = zeros(n_s, n_y);

% 值函数迭代
for iter = 1:max_iter
    V_new = zeros(n_s, n_y);
    pol_fun_new = zeros(n_s, n_y);
    
    for i_s = 1:n_s
        x = X(i_s);
        c_grid = linspace(0, x, n_c);  % 消费网格从0到当前蛋糕存量
        
        for i_y = 1:n_y
            aux = zeros(1, n_c);
            
            for i_c = 1:n_c
                c = c_grid(i_c);
                if c > x
                    aux(i_c) = -Inf;
                    continue;
                end
                
                % 计算未来蛋糕存量的基础值（折旧后的剩余蛋糕）
                next_x_base = delta * (x - c);
                % 计算所有未来冲击状态下的值函数
                next_v = zeros(1, n_y);
                for i_y_next = 1:n_y
                    next_x = next_x_base + y_vals(i_y_next);
                    if next_x < 0 || next_x > 1
                        next_v(i_y_next) = -Inf;
                    else
                        % 插值得到未来值函数，超出范围返回负无穷
                        next_v(i_y_next) = interp1(X, V(:, i_y_next), next_x, 'spline', -Inf);
                    end
                end
                
                % 计算期望未来值（当前冲击状态下的加权平均）
                E_next_v = Pi(i_y, :) * next_v';
                % 当前效用+贴现后的期望未来值
                aux(i_c) = u(c) + beta * E_next_v;
            end
            
            % 记录最优值与对应消费索引
            [V_new(i_s, i_y), pol_fun_new(i_s, i_y)] = max(aux);
        end
    end
    
    % 检查收敛
    if max(max(abs(V_new - V))) < tol
        fprintf('收敛于第%d次迭代\n', iter);
        break;
    end
    
    V = V_new;
    pol_fun = pol_fun_new;
    
    if iter == max_iter
        fprintf('达到最大迭代次数未收敛\n');
    end
end

% 计算最优消费值
opt_c = zeros(n_s, n_y);
for i_s = 1:n_s
    x = X(i_s);
    c_grid = linspace(0, x, n_c);
    for i_y = 1:n_y
        opt_c(i_s, i_y) = c_grid(pol_fun(i_s, i_y));
    end
end

% 绘制值函数
figure;
hold on;
plot(X, V(:, 1), 'r-', 'LineWidth', 2, 'DisplayName', '低冲击状态');
plot(X, V(:, 2), 'b-', 'LineWidth', 2, 'DisplayName', '高冲击状态');
xlabel('蛋糕存量');
ylabel('值函数');
title('值函数随蛋糕存量的变化');
legend;
hold off;

% 绘制策略函数
figure;
hold on;
plot(X, opt_c(:, 1), 'r-', 'LineWidth', 2, 'DisplayName', '低冲击状态');
plot(X, opt_c(:, 2), 'b-', 'LineWidth', 2, 'DisplayName', '高冲击状态');
xlabel('蛋糕存量');
ylabel('最优消费');
title('最优消费策略随蛋糕存量的变化');
legend;
hold off;

关键修正说明 #

• 统一冲击设定：将连续冲击网格改为离散两状态，与转移矩阵维度匹配，确保期望计算逻辑正确。若需连续冲击，需改用AR(1)过程的转移密度进行积分计算期望。
• 多轮迭代收敛：加入收敛判断循环，让值函数迭代直到相邻两轮的差值小于阈值，保证计算结果收敛到真实值函数。
• 修正期望计算：对每个当前冲击状态，用转移矩阵加权所有未来冲击状态的值函数，得到符合动态规划要求的期望未来值。
• 修复消费网格：消费下限设为0，避免出现无效的反向网格，确保消费值始终在合理范围内。
• 匹配状态维度：明确值函数的维度对应关系（蛋糕状态×冲击状态），未来值的插值和期望计算均严格匹配状态转移逻辑。

内容的提问来源于stack exchange，提问作者LioM