Below is a step‐by‐step solution in Python. We first build the 3‐month SOFR curve, discount factors and helper routines (Black‐76 caplet pricer and Hagan’s SABR implied‐vol formula), then loop to compute:

1. Implied Black volatilities for each caplet (for maturities $T_k = k\Delta t$, $k=1,\dots,40$, and strikes $K=1\%,2\%,\dots,9\%$).
2. Cap values for cap maturities $T_i=1,\dots,10$ years across the same strikes.

Finally we plot, for each parameter set:

• a 3D surface of $\sigma_{\rm Black}(T_k,K)$ (“vol‐surface”), and
• a 3D surface of $\mathrm{CapValue}(T_i,K)$.

python
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.special import erf

# 1. Problem data
dt = 0.25
N      = 120                               # quarters in our SOFR curve
f      = np.linspace(0.0525, 0.0475, N)   # f_j(0) from 5.25% down to 4.75%
# build discount factors P[k] = P(0,t_k), t_k = k*dt
P = np.empty(N)
P[0] = 1.0 / (1.0 + dt * f[0])
for k in range(1, N):
    P[k] = P[k-1] / (1.0 + dt * f[k])

# strikes and maturities
K_grid = np.arange(0.01, 0.10, 0.01)       # 1% to 9%
T_caps = np.arange(1, 11)                  # cap maturities in years
# we will need caplet times only up to 10y = 40 quarters
K_periods = int(T_caps[-1] / dt)           # = 40

# 2. Helper functions

def black_call(F, K, T, sigma):
    """Black-76 call price (undiscounted) with forward F, strike K, time to expiry T."""
    if T * sigma == 0:
        return np.maximum(F - K, 0.0)
    d1 = (np.log(F/K) + 0.5 * sigma**2 * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    N1 = 0.5 * (1 + erf(d1/np.sqrt(2)))
    N2 = 0.5 * (1 + erf(d2/np.sqrt(2)))
    return F * N1 - K * N2

def sabr_implied_vol(F, K, T, alpha, beta, rho, nu):
    """Hagan SABR implied volatility approximation."""
    # edge‐cases
    if T == 0:
        return alpha
    if F <= 0 or K <= 0:
        return np.nan

    bs = 1.0 - beta
    FKb = (F * K)**(bs/2)
    # small‐time expansion factor
    A = (bs**2 / 24) * (alpha**2) / ((F * K)**bs) \
      + 0.25 * rho * beta * nu * alpha / ((F * K)**(bs/2)) \
      + (2 - 3*rho**2) * nu**2 / 24
    E = 1 + A * T

    # z and x(z)
    if abs(F - K) < 1e-12:
        # ATM formula
        vol = alpha / (F**bs) * E
    else:
        if bs != 0:
            z = (nu/alpha) * (F**bs - K**bs) / bs
        else:
            z = (nu/alpha) * np.log(F/K)

        # x(z)
        numer = np.sqrt(1 - 2*rho*z + z**2) + z - rho
        denom = 1 - rho
        xz = np.log(numer / denom)

        coef = (nu * np.log(F/K)) / xz
        vol = coef * E

    return vol

# 3. Compute surfaces for both parameter sets
param_sets = {
    'a (β=0.5, α=0.05, ν=0.25, ρ=-0.5)': dict(beta=0.5, alpha=0.05, nu=0.25, rho=-0.5),
    'b (β=0.4, α=0.0372, ν=0.25, ρ=-0.5)': dict(beta=0.4, alpha=0.0372, nu=0.25, rho=-0.5)
}

surfaces_vol   = {}
surfaces_caps  = {}

for label, pars in param_sets.items():
    β, α, ρ, ν = pars['beta'], pars['alpha'], pars['rho'], pars['nu']

    # 3.1 Caplet vol surface σ_B[k, j] for k=1..K_periods, j over strikes
    sigB = np.zeros((K_periods, len(K_grid)))
    times = np.arange(1, K_periods+1) * dt  # T_k = k*dt

    for k in range(K_periods):
        T_k = k * dt                  # note: the Hagan formula uses time-to-expiry = (k)*dt
        Fk = f[k]
        for j, Kj in enumerate(K_grid):
            sigB[k, j] = sabr_implied_vol(Fk, Kj, T_k, α, β, ρ, ν)

    # 3.2 Cap price surface Cap(T_i, K_j)
    caps = np.zeros((len(T_caps), len(K_grid)))
    for i, Ti in enumerate(T_caps):
        # caplet indices for this cap: these run from periods (i*4)−4 to i*4−1  (zero‐based)
        idx_start = 4*(i) - 4
        idx_end   = 4*(i)       - 1
        for j, Kj in enumerate(K_grid):
            price = 0.0
            for k in range(idx_start, idx_end+1):
                # maturity for call() = (k)*dt
                price += P[k] * dt * black_call(f[k], Kj, k*dt, sigB[k, j])
            caps[i, j] = price

    surfaces_vol[label]  = (times, K_grid, sigB)
    surfaces_caps[label] = (T_caps, K_grid, caps)


# 4. Plot everything
for label in param_sets:
    times, Ks, sigB   = surfaces_vol[label]
    Tc, Ks2, cap_vals = surfaces_caps[label]

    fig = plt.figure(figsize=(12,5))
    ax = fig.add_subplot(121, projection='3d')
    Tm, Km = np.meshgrid(times, Ks)
    ax.plot_surface(Km, Tm, sigB.T, cmap='viridis', edgecolor='none')
    ax.set_title(f'Implied Black Vol Surface ({label})')
    ax.set_xlabel('Strike K')
    ax.set_ylabel('Caplet Maturity (y)')
    ax.set_zlabel('σ_Black')

    ax2 = fig.add_subplot(122, projection='3d')
    Ti, Kj = np.meshgrid(Tc, Ks2)
    ax2.plot_surface(Kj, Ti, cap_vals.T, cmap='plasma', edgecolor='none')
    ax2.set_title(f'Cap Value Surface ({label})')
    ax2.set_xlabel('Strike K')
    ax2.set_ylabel('Cap Maturity (y)')
    ax2.set_zlabel('Cap Price')

    plt.tight_layout()
    plt.show()

Explanation of key steps

1. Term‐rate curve $\{f_j\}$
   We place 120 quarterly forward rates $f_j$ evenly from 5.25% down to 4.75%, and build discount factors

   $$P(0,t_k)=\prod_{m=1}^k\frac1{1+\Delta t\,f_m}\,,\quad \Delta t=0.25.$$

2. Black‐76 caplet pricer
   The function black_call(F,K,T,σ) implements

   $$C_{\rm Black}(F,K,T) \;=\; F\,N(d_1)-K\,N(d_2), \quad d_{1,2}=\frac{\ln(F/K)\pm\frac12\sigma^2T}{\sigma\sqrt T}.$$

3. Hagan SABR formula
   In sabr_implied_vol we code the standard small‐time expansion plus the $z/x(z)$ correction:

   $$\sigma_{\rm SABR}(F,K) = \frac{\nu\ln(F/K)}{x(z)}\,[1 + A\,T],\quad z=\frac{\nu}{\alpha}\frac{F^{1-\beta}-K^{1-\beta}}{1-\beta},$$

   with

   $$A = \frac{(1-\beta)^2}{24}\frac{\alpha^2}{(FK)^{1-\beta}} + \frac{\rho\beta\nu\alpha}{4\,(F K)^{(1-\beta)/2}} + \frac{2-3\rho^2}{24}\nu^2.$$

   The ATM special case $F\!\approx\!K$ is handled separately.

4. Building the surfaces

   • Volatility surface: for each caplet expiry $T_k=k\Delta t$ up to 10 years (40 quarters) and each strike.
   • Caplet pricing: each cap of maturity $T_i$ is the sum of its four quarter‐caplets, $\displaystyle \sum_{k=4(i-1)}^{4i-1}P_k\,\Delta t\;\times$ Black call price.

5. Visualization
   We use Matplotlib’s 3D plotting to show, side by side for each parameter set,

   • the implied‐vol surface $\sigma_{\rm Black}(T_k,K)$, and
   • the cap‐value surface $\mathrm{CapValue}(T_i,K)$.

Feel free to tweak the plotting styles, color maps or grid‐resolution (e.g. a finer strike mesh) as needed.

下面我对两份代码做了对比和检查，结论是：

• 我最初给出的代码

  • 在计算 SABR 隐含波动率时，完全漏掉了前缀 $\displaystyle\frac{\alpha}{(FK)^{(1-\beta)/2}}$ 项。
  • 在构造波动率曲面时，用错了到期时间：第 $j$ 个 caplet 的到期应是 $(j+1)\Delta t$，而我用了 $j\Delta t$。
  • 在合约定价里，caplet 索引也偏移错误，导致第一个 cap 反而用到了负索引。

• 你贴过来的备用代码

  • ATM（$F=K$）情形下只返回了 $\alpha/F^{1-\beta}$，却漏掉了小-时间展开修正项 $\bigl[1+AT\bigr]$；
  • 非 ATM 情形里又把修正项 $\bigl[1+AT\bigr]$ 放在了分母，而正确的应该是放在分子。

正确解法的完整 Python 实现

python
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erf
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401

# 1. 构造 3 个月 SOFR 曲线
dt = 0.25
N_quarters = 120
j = np.arange(1, N_quarters+1)
f = 0.0525 + (j-1)/(119) * (0.0475 - 0.0525)   # f_j 从 5.25% 线性下跌到 4.75%

# 2. 贴现因子 P(0,t_j) = ∏_{m=1}^j 1/(1 + dt * f_m)
P = np.empty(N_quarters)
P[0] = 1.0 / (1.0 + dt * f[0])
for m in range(1, N_quarters):
    P[m] = P[m-1] / (1.0 + dt * f[m])

# 3. Black-76 caplet 定价函数
def black_call(F, K, T, sigma):
    if T * sigma == 0:
        return max(F - K, 0.0)
    d1 = (np.log(F/K) + 0.5 * sigma**2 * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    N1 = 0.5 * (1 + erf(d1/np.sqrt(2)))
    N2 = 0.5 * (1 + erf(d2/np.sqrt(2)))
    return F * N1 - K * N2

# 4. Hagan 等式的 SABR 隐含波动率
def sabr_implied_vol(F, K, T, alpha, beta, rho, nu):
    # T = 0 时直接返回 0，让 black_call 走 T*sigma=0 分支
    if T <= 0:
        return 0.0

    bs = 1.0 - beta
    # ATM 特殊处理
    if abs(F - K) < 1e-12:
        # α / F^{1-β} * [1 + A T]
        A = (bs**2/24)*(alpha**2/(F**(2*bs))) \
          + (rho*beta*nu*alpha)/(4*(F**bs)) \
          + ((2-3*rho**2)/24)*nu**2
        return alpha / (F**bs) * (1 + A * T)

    # 非 ATM
    # A 系数
    A = (bs**2/24)*(alpha**2/((F*K)**bs)) \
      + (rho*beta*nu*alpha)/(4*((F*K)**(bs/2))) \
      + ((2-3*rho**2)/24)*nu**2
    E = 1 + A * T

    # z = ν/α * (F^{1-β} - K^{1-β})/(1-β)
    if abs(bs) > 1e-16:
        z = (nu/alpha) * (F**bs - K**bs) / bs
    else:
        z = (nu/alpha) * np.log(F/K)

    # χ(z) = ln[(&radic;(1−2ρz+z^2)+z−ρ)/(1−ρ)]
    numer = np.sqrt(1 - 2*rho*z + z**2) + z - rho
    denom = 1 - rho
    chi = np.log(numer / denom)

    # 最终公式： α/( (F K)^{(1-β)/2} ) * E * (z/χ)
    prefactor = alpha / ((F*K)**(bs/2))
    return prefactor * E * (z / chi)

# 5. 网格设定
strikes = np.arange(1, 10) / 100.0      # 1% 到 9%
cap_mats = np.arange(1, 11)             # 1y 到 10y
max_caplets = int(cap_mats[-1] / dt)    # 40 个 caplet

# 6. 参数组
param_sets = [
    ('(a) β=0.5, α=0.05, ν=0.25, ρ=-0.5', 0.5, 0.05, 0.25, -0.5),
    ('(b) β=0.4, α=0.0372, ν=0.25, ρ=-0.5', 0.4, 0.0372, 0.25, -0.5),
]

# 7. 逐组计算并绘图
for title, beta, alpha, nu, rho in param_sets:
    # 7.1 先算 caplet 隐含波动率曲面 σ[j,k]
    sigma_surf = np.zeros((max_caplets, len(strikes)))
    caplet_times = (np.arange(1, max_caplets+1) * dt)
    for j_idx, Tj in enumerate(caplet_times):
        Fj = f[j_idx]  # 对应第 j+1 个远期
        for k_idx, K in enumerate(strikes):
            sigma_surf[j_idx, k_idx] = sabr_implied_vol(Fj, K, Tj, alpha, beta, rho, nu)

    # 7.2 再算各期限 Cap 的价格 Cap(T_i,K)
    cap_surf = np.zeros((len(cap_mats), len(strikes)))
    for i_idx, Ti in enumerate(cap_mats):
        n_caplets = int(Ti / dt)
        for k_idx, K in enumerate(strikes):
            price = 0.0
            for j_idx in range(n_caplets):
                Tj = (j_idx+1) * dt
                Fj = f[j_idx]
                vol = sigma_surf[j_idx, k_idx]
                price += P[j_idx] * dt * black_call(Fj, K, Tj, vol)
            cap_surf[i_idx, k_idx] = price

    # 8. 可视化
    fig = plt.figure(figsize=(12, 5))

    # 8.1 Vol 曲面
    ax1 = fig.add_subplot(1, 2, 1, projection='3d')
    K_mesh, T_mesh = np.meshgrid(strikes, caplet_times)
    ax1.plot_surface(K_mesh, T_mesh, sigma_surf, cmap='viridis', edgecolor='none')
    ax1.set_title(f'{title} 的隐含波动率曲面')
    ax1.set_xlabel('Strike K')
    ax1.set_ylabel('Caplet 到期 T (年)')
    ax1.set_zlabel('σ_Black')

    # 8.2 Cap 价值曲面
    ax2 = fig.add_subplot(1, 2, 2, projection='3d')
    K2_mesh, T2_mesh = np.meshgrid(strikes, cap_mats)
    ax2.plot_surface(K2_mesh, T2_mesh, cap_surf, cmap='plasma', edgecolor='none')
    ax2.set_title(f'{title} 的 Cap 价值曲面')
    ax2.set_xlabel('Strike K')
    ax2.set_ylabel('Cap 到期 T (年)')
    ax2.set_zlabel('价格 (单位名义)')

    plt.tight_layout()
    plt.show()

关键修正点

1. SABR 公式前缀
   $\displaystyle\frac{\alpha}{(F\,K)^{\frac{1-\beta}{2}}}$，以及小时间展开修正项 $\bigl[1+A\,T\bigr]$ 都必须放在“分子”上。

2. 到期时间
   第 $j$ 个 caplet 的实际到期是 $(j+1)\Delta t$，因此波动率和 Black 定价里都要用 $T=(j+1)\,dt$。

3. Caplet 索引
   cap of maturity $T_i$ 包含前 $4\,T_i$ 个 caplet（即从 $j=1$ 到 $j=4T_i$），对应 Python 里索引从 0 到 $4T_i-1$。

这样即可正确绘出隐含波动率曲面和 Cap 价值曲面。