Julia中使用SymPy求解一般均衡问题的索引操作与方程组求解疑问

1. 索引变量赋值/传入表达式的方法 #

SymPy中的IndexedBase是符号索引容器，属于不可变符号对象，不支持[]=直接赋值操作，正确的实现方式是使用符号替换（subs）机制，或是直接按索引生成对应表达式数组：

外生变量校准实现 #

using SymPy

# 设定经济体商品数量N
N = 3
n = symbols("n", integer=true)

# 定义外生符号变量
α = IndexedBase("α")
@syms γ

# 因为N为确定值，直接展开求和比使用SymPy.Sum更便于后续计算
α_sum = sum(α[i] for i in 1:N)
equation_1 = Eq(α_sum + γ, 1)

# 构造校准值替换字典
alpha3 = [1,2,3]
calib_subs = Dict(α[i] => alpha3[i] for i in 1:N)

# 代入校准值得到结果
calibrated_eq1 = subs(equation_1, calib_subs)
# 输出结果为 6 + γ = 1，符合预期逻辑

一阶条件表达式传入 #

对于供需函数这类已推导好的表达式，直接按索引生成表达式数组即可，不需要对IndexedBase对象赋值：

# 定义符号变量
P = IndexedBase("P") # 待求解价格向量
z = IndexedBase("z")
s = IndexedBase("s")
Y = IndexedBase("Y")

# 直接生成各索引对应的供需表达式
S = [s[i] * P[i] for i in 1:N]
D = [z[i] / P[i] for i in 1:N]

2. 带索引内生变量的方程组求解 #

需要先生成所有索引对应的方程，再将所有待求解的价格变量收集为列表传入求解函数即可：

# 生成所有商品的市场出清方程
market_clear_eqs = [Eq(S[i] - D[i], 0) for i in 1:N]
# 生成总收入约束方程
total_constraint = Eq(sum(P[i] * Y[i] for i in 1:N), 0)
# 合并为完整方程组
eq_system = vcat(market_clear_eqs, total_constraint)

# 收集所有待求解的价格变量
solve_vars = [P[i] for i in 1:N]

# 求解方程组（多变量方程组使用solve更适配，solveset主要用于单变量解集求解）
solution = solve(eq_system, solve_vars)

如果后续需要做外生冲击分析，只需要修改校准替换字典中对应外生变量的取值，代入求解结果即可得到新的均衡值。

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