AI 驱动的结构估计：把"解模型"变成一次求解

import torch, torch.nn as nn, numpy as np
torch.manual_seed(0); np.random.seed(0)

ALPHA, BETA_DISC, R, PSI = 0.65, 0.96, 0.04, 1.0
NK, NB, NZ = 9, 7, 6          # 候选 k' / 候选 b' / z' 蒙特卡洛抽样数

def profit(k, z):             # 生产利润 π = z * k^α（z 取指数以保证为正）
    return torch.exp(z) * (k ** ALPHA)

# 值函数网络：输入(k,b,z,γ0,δ,ρ,σ) 共7维 → V。参数与价格进状态=扩展状态空间
class ValueNet(nn.Module):
    def __init__(self):
        super().__init__()
        self.net = nn.Sequential(nn.Linear(7,64), nn.SiLU(),
                                 nn.Linear(64,64), nn.SiLU(),
                                 nn.Linear(64,1))
    def forward(self, s): return self.net(s).squeeze(-1)

Vnet, Vtar = ValueNet(), ValueNet()   # Vtar=滞后 target 网络，稳定拟合值迭代
Vtar.load_state_dict(Vnet.state_dict())
optV = torch.optim.Adam(Vnet.parameters(), lr=1e-3)

k_grid = torch.exp(torch.linspace(np.log(0.5), np.log(6.0), NK))
b_grid = torch.linspace(0.0, 3.0, NB)

def sample_state(n):          # 随机采状态 (k,b,z) 与参数 (γ0,δ,ρ,σ)
    return (torch.rand(n)*5+.3, torch.rand(n)*2.5, torch.randn(n)*.5,
            torch.rand(n)*.04+.005, torch.rand(n)*.08+.04,
            torch.rand(n)*.2+.8, torch.rand(n)*.1+.05)

def bellman_residual(n=512):
    k,b,z,g,d,rho,sig = sample_state(n)
    kp, bp = k_grid.view(1,-1), b_grid.view(1,-1)
    inv = kp - (1-d.view(-1,1))*k.view(-1,1)          # 净投资
    c_adj = g.view(-1,1)*(inv.abs()>1e-6).float() + .5*PSI*inv**2  # 固定+二次调整成本
    div_k = profit(k.view(-1,1), z.view(-1,1)) - c_adj - b.view(-1,1)
    div = div_k.unsqueeze(2) + (bp/(1+R)).view(1,1,NB)            # 现金流 (n,NK,NB)
    zp = rho.view(-1,1)*z.view(-1,1) + sig.view(-1,1)*torch.randn(n,NZ)  # z' AR(1)
    sN = torch.zeros(n,NK,NB,NZ,7)
    sN[...,0]=kp.view(1,NK,1,1); sN[...,1]=bp.view(1,1,NB,1)
    sN[...,2]=zp.view(n,1,1,NZ)
    sN[...,3]=g.view(-1,1,1,1); sN[...,4]=d.view(-1,1,1,1)
    sN[...,5]=rho.view(-1,1,1,1); sN[...,6]=sig.view(-1,1,1,1)
    with torch.no_grad():                              # 用滞后网络算 Bellman 目标
        Vnext = Vtar(sN.reshape(-1,7)).reshape(n,NK,NB,NZ).mean(-1)
    V_target = torch.logsumexp((div + BETA_DISC*Vnext).reshape(n,-1), dim=1)  # 平滑 max
    V_cur = Vnet(torch.stack([k,b,z,g,d,rho,sig],1))
    return ((V_cur - V_target)**2).mean()

for ep in range(400):                 # 训练循环：最小化 Bellman 残差
    optV.zero_grad()
    loss = bellman_residual(); loss.backward(); optV.step()
    with torch.no_grad():             # Polyak 缓慢更新 target 网络
        for p,pt in zip(Vnet.parameters(), Vtar.parameters()):
            pt.mul_(0.995).add_(0.005*p)
    if ep % 100 == 0: print(f"epoch {ep}: Bellman residual = {loss.item():.4f}")
# 实跑：残差 35.9 → 0.27，值函数收敛

@torch.no_grad()
def simulate_moments(g0, d0, r0, s0, T=30, N=16):   # 用训练好的值网络模拟矩
    k = torch.ones(N)*2.0; b = torch.ones(N)*0.5; z = torch.zeros(N)
    inv_hist, lev_hist = [], []
    for t in range(T):
        kp, bp = k_grid.view(1,-1), b_grid.view(1,-1)
        inv = kp - (1-d0)*k.view(-1,1)
        ca = g0*(inv.abs()>1e-6).float() + .5*PSI*inv**2
        div = (profit(k.view(-1,1),z.view(-1,1))-ca-b.view(-1,1)).unsqueeze(-1) \
              + (bp/(1+R)).view(1,1,NB)
        zp = r0*z.view(-1,1) + s0*torch.randn(N,NZ)
        sN = torch.zeros(N,NK,NB,NZ,7)
        sN[...,0]=kp.view(1,NK,1,1); sN[...,1]=bp.view(1,1,NB,1)
        sN[...,2]=zp.view(N,1,1,NZ)
        sN[...,3]=g0; sN[...,4]=d0; sN[...,5]=r0; sN[...,6]=s0
        VN = Vnet(sN.reshape(-1,7)).reshape(N,NK,NB,NZ).mean(-1)
        best = (div + BETA_DISC*VN).reshape(N,-1).argmax(1)   # 最优 (k',b')
        kn, bn = k_grid[best//NB], b_grid[best%NB]
        inv_hist.append(((kn-(1-d0)*k)/k).mean()); lev_hist.append((b/k).mean())
        k,b,z = kn, bn, r0*z + s0*torch.randn(N)
    return torch.tensor([torch.stack(inv_hist[10:]).mean(),
                         torch.stack(lev_hist[10:]).mean(),
                         torch.stack(inv_hist[10:]).std()])

# 伪参数示意：把均衡价格 p 当伪参数，市场出清失衡=0 作为额外矩
def market_clearing_imbalance(pseudo_price, g0,d0,r0,s0):
    lev = simulate_moments(g0,d0,r0,s0)[1]
    return lev - 1.0/pseudo_price          # 均衡要求该失衡 = 0

# --- 采样参数设计集，构造 (参数, 矩) 数据 ---
Nd = 200
g_d,d_d,r_d,s_d = (torch.rand(Nd)*.04+.005, torch.rand(Nd)*.08+.04,
                   torch.rand(Nd)*.2+.8,   torch.rand(Nd)*.1+.05)
M = torch.stack([simulate_moments(g_d[i],d_d[i],r_d[i],s_d[i]) for i in range(Nd)])

class MomentNet(nn.Module):              # 参数→矩 代理网络（输入参数，输出矩）
    def __init__(self):
        super().__init__()
        self.net = nn.Sequential(nn.Linear(4,32), nn.SiLU(),
                                 nn.Linear(32,32), nn.SiLU(), nn.Linear(32,3))
    def forward(self, x): return self.net(x)

mnet = MomentNet(); optM = torch.optim.Adam(mnet.parameters(), lr=1e-3)
X = torch.stack([g_d,d_d,r_d,s_d],1); Xm, Xs = X.mean(0), X.std(0)+1e-8
Mm, Ms = M.mean(0), M.std(0)+1e-8
for ep in range(800):                    # MSE 训练矩映射
    optM.zero_grad()
    loss = ((mnet((X-Xm)/Xs) - (M-Mm)/Ms)**2).mean()
    loss.backward(); optM.step()
with torch.no_grad():                    # 报告样本外/样本内 R^2
    pred = mnet((X-Xm)/Xs)*Ms+Mm
    R2 = 1 - ((M-pred)**2).sum(0)/((M-M.mean(0))**2).sum(0)
print("矩网络 R^2 =", R2.tolist())       # 实跑：[0.995, 0.768, 0.76]

# 真实参数，生成"数据矩"
g_true,d_true,r_true,s_true = 0.02, 0.08, 0.90, 0.12
m_data = simulate_moments(g_true,d_true,r_true,s_true)

# 把待估参数当 leaf variable；通过可微矩网络求梯度（替代每次重新模拟）
th = torch.tensor([0.01,0.10,0.85,0.10], requires_grad=True)
optEst = torch.optim.Adam([th], lr=1e-2)
lo = torch.tensor([.005,.04,.8,.05]); hi = torch.tensor([.045,.12,1.0,.15])
W = torch.diag(1.0/(Ms**2))              # 加权矩阵 = 矩方差逆（对角近似）
for step in range(400):                  # 最小距离估计（SMM），梯度来自矩网络
    optEst.zero_grad()
    m_hat = mnet((th-Xm)/Xs)*Ms+Mm       # 毫秒级前向预测矩
    loss = (m_data-m_hat) @ W @ (m_data-m_hat)   # 加权矩距离（可微）
    loss.backward(); optEst.step()       # 自动微分给出 ∂loss/∂θ
    with torch.no_grad(): th.clamp_(lo,hi)
    if step % 100 == 0:
        print(step, round(loss.item(),4), [round(x,3) for x in th.tolist()])
# 实跑恢复：δ=0.078(真.08), ρ=0.909(真.90)；γ0/σ 因弱识别偏离——见识别诊断