简单CGE模型实施冲击后结果无变化的技术求助

问题描述 #

我是Python新手，正在求解一个非常简单的CGE模型。但在对经济施加技术冲击（shock_A=1.1）后，求解结果与基准情况几乎一致，想知道是否需要调整代码写法、更换优化器或采取其他措施。

原模型类代码 #

import numpy as np
from scipy import optimize
from scipy.optimize import fsolve

class SimpleEquilibriumModel:
    def __init__(self):
        # Define objects and dimensions
        self.Y = 0  # Output (Endogenous)
        self.L = 0  # Labor (Endogenous)
        self.C = 0  # Consumption (Endogenous)
        self.W = 0  # Wage (Endogenous)
        self.P = 0  # Price (Endogenous)
        self.A = 0  # Technology (Exogenous)
        self.N = 0  # Population (Exogenous)

        # Define equations
        def E_L():
            return self.L == self.Y / self.A

        def E_NP():
            return self.P * self.Y == self.W * self.L

        def E_W():
            return self.L == self.N

        def E_C():
            return self.C == self.W / self.P * self.N

        def E_P():
            return self.P == 1
        
        def E_A():
            return self.A == self.W / self.P

        self.equations = [E_L, E_NP, E_W, E_C, E_P, E_A]

        # Define national accounts parameter
        self.NR = 1000

    def solve(self):
        # Initialization and calibration
        self.W = 1
        self.P = 1
        self.L = self.NR
        self.Y = self.NR
        self.N = self.L
        self.A = self.W / self.P


        # Solve the model
        def model_equations(x):
            self.W, self.P, self.L, self.Y, self.N, self.A = x
            return [E() for E in self.equations]

        self.W, self.P, self.L, self.Y, self.N, self.A = fsolve(model_equations, [1, 1, self.L, self.Y, self.N, self.A])

        self.C = self.W / self.P * self.N

        return self.Y, self.L, self.C, self.W, self.P
    
    def solve_shock(self, shock_A = None):
        # Initialization and calibration
        self.W == 2
        self.P == 2
        self.L = self.NR
        self.Y = self.NR
        self.N = self.L
        self.A = (self.W / self.P)*shock_A


        # Solve the model
        def model_equations(x):
            self.W, self.P, self.L, self.Y, self.N, self.A = x
            return [E() for E in self.equations]

        self.W, self.P, self.L, self.Y, self.N, self.A = fsolve(model_equations, [1, 1, self.L, self.Y, self.N, self.A])

        self.C = self.W / self.P * self.N

        return self.Y, self.L, self.C, self.W, self.P

原笔记本调用代码 #

import numpy as np
from scipy import optimize
from . import SimpleEquilibriumModel

output, employment, consumption, wage, price = model.solve()

print("Equilibrium output:", output)
print("Equilibrium employment:", employment)
print("Equilibrium consumption:", consumption)
print("Equilibrium wage:", wage)
print("Equilibrium price:", price)

output, employment, consumption, wage, price = model.solve_shock(shock_A = 1.1)

print("New equilibrium output:", output)
print("New equilibrium employment:", employment)
print("New equilibrium consumption:", consumption)
print("New equilibrium wage:", wage)
print("New equilibrium price:", price)

原求解输出 #

Equilibrium output: 1000.0000018014192
Equilibrium employment: 1000.0000056491614
Equilibrium consumption: 999.9999869002581
Equilibrium wage: 0.9999999925494194
Equilibrium price: 1.0000000149011612
New equilibrium output: 1000.0
New equilibrium employment: 1000.0000149011612
New equilibrium consumption: 999.9999701976781
New equilibrium wage: 0.9999999850988388
New equilibrium price: 1.0000000149011612

问题根源分析 #

1. 赋值语法错误：solve_shock方法中self.W == 2和self.P == 2是比较运算符，而非赋值（应使用=），导致初始值未被正确修改，仍沿用基准场景的数值。
2. 外生变量被内生化：原模型通过E_A()方程将外生技术变量A与内生变量W/P绑定，且把A纳入求解向量x中，导致技术冲击的初始设置被求解过程覆盖，无法体现冲击效果。
3. 方程返回类型错误：原方程返回布尔值（==判断），但fsolve需要返回数值残差（等式两边的差值，目标是残差为0），布尔值转成0/1会干扰求解精度。
4. 冗余约束：6个方程对应6个变量，但N是固定外生变量，E_P()强制P=1，实际存在冗余约束，直接抵消了技术冲击的影响。

修正后的模型代码 #

import numpy as np
from scipy.optimize import fsolve

class SimpleEquilibriumModel:
    def __init__(self):
        # 内生变量
        self.Y = 0  # 产出
        self.L = 0  # 劳动
        self.C = 0  # 消费
        self.W = 0  # 工资
        self.P = 0  # 价格
        # 外生变量
        self.A = 1  # 初始技术水平
        self.N = 1000  # 固定人口

        # 定义均衡方程（返回数值残差，而非布尔值）
        def E_L():
            return self.L - self.Y / self.A  # 生产函数残差：L = Y/A → L - Y/A = 0

        def E_NP():
            return self.P * self.Y - self.W * self.L  # 利润为零残差：PY = WL → PY - WL = 0

        def E_W():
            return self.L - self.N  # 劳动市场均衡残差：L = N → L - N = 0

        def E_C():
            return self.C - self.W / self.P * self.N  # 消费均衡残差：C = (W/P)*N → C - (W/P)N = 0

        def E_P():
            return self.P - 1  # 价格标准化残差：P = 1 → P - 1 = 0

        self.equations = [E_L, E_NP, E_W, E_C, E_P]

    def solve(self):
        # 内生变量初始猜测：W, P, L, Y, C
        init_guess = [1, 1, self.N, self.A * self.N, self.N]

        def model_equations(x):
            self.W, self.P, self.L, self.Y, self.C = x
            return [eq() for eq in self.equations]

        # 求解均衡
        self.W, self.P, self.L, self.Y, self.C = fsolve(model_equations, init_guess)
        return self.Y, self.L, self.C, self.W, self.P
    
    def solve_shock(self, shock_A=1.1):
        # 直接修改外生技术变量，施加冲击
        self.A = shock_A
        # 复用基准求解逻辑
        return self.solve()

修正后的调用代码 #

from SimpleEquilibriumModel import SimpleEquilibriumModel

# 初始化模型
model = SimpleEquilibriumModel()

# 求解基准场景
output, employment, consumption, wage, price = model.solve()
print("基准均衡产出:", output)
print("基准均衡就业:", employment)
print("基准均衡消费:", consumption)
print("基准均衡工资:", wage)
print("基准均衡价格:", price)

# 求解冲击场景
output_shock, employment_shock, consumption_shock, wage_shock, price_shock = model.solve_shock(shock_A=1.1)
print("\n冲击后均衡产出:", output_shock)
print("冲击后均衡就业:", employment_shock)
print("冲击后均衡消费:", consumption_shock)
print("冲击后均衡工资:", wage_shock)
print("冲击后均衡价格:", price_shock)

预期修正后输出 #

技术冲击（A提升10%）会带来符合理论预期的结果：

基准均衡产出: 1000.0
基准均衡就业: 1000.0
基准均衡消费: 1000.0
基准均衡工资: 1.0
基准均衡价格: 1.0

冲击后均衡产出: 1100.0
冲击后均衡就业: 1000.0
冲击后均衡消费: 1100.0
冲击后均衡工资: 1.1
冲击后均衡价格: 1.0

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