复现代码，无开源，结果对不上

Question

正在复现这篇文献《信用风险传染效应与外溢冲击研究》的第一部分，模型设定如下图所示：

本文选用产业债信用利差数据构建风险传染网络。具体而言，本文的主要分析变量为煤炭、新能源、房地产、投资平台、非银金融等 27个行业的产业债信用利差指数；分别基于 1-5个交易日、6-20个交易日以及 21-250个交易日三个频段，将产业债间的总体关联分解为短期、中期与长期，从而深入剖析信用风险在不同时期的传染效应。
基于以上，我目前的代码如下，但输出的结果与原文献结果有较大偏差，不知道是哪里出现问题：

import pandas as pd  
import numpy as np  
from statsmodels.tsa.api import VAR  
from scipy.integrate import quad  
  
data = pd.read_excel('行业产业债信用利差.xlsx')   
data.head()

# 设定最大滞后阶数  
maxlags = 5  
  
# 初始化VAR模型并拟合，使用信息准则（如AIC）来选择最优滞后阶数  
model = VAR(data)  
result = model.fit(maxlags=maxlags, ic='aic') 
print('Optimal number of lags:', result.k_ar) 
# 获取移动平均系数矩阵和残差协方差矩阵  
psi_h =result.ma_rep(result.k_ar) #移动平均系数矩阵(3,27,27)
# 获取残差
resid = result.resid
# 计算残差的协方差矩阵
sigma = pd.DataFrame(np.cov(resid.T))  # 获取残差的协方差矩阵(27,27)
psi_h.shape

# 定义函数来计算频率响应函数  
def frequency_response(omega, psi_h):
    n_var = psi_h.shape[1]  # 变量数，即DataFrame中的列数  
    n_lags = result.k_ar # 滞后阶数
    psi_omega = np.zeros((n_var, n_var), dtype=complex) # 初始化频率响应矩阵 
    for h in range(n_lags):
        psi_omega += psi_h[h] * np.exp(-1j * omega * h)  # 计算并累加频率响应
    return psi_omega

# 计算所有频率点的频率响应   
omega_range = np.linspace(0, np.pi / 5, 250)  #短期频段
#omega_range = np.linspace(np.pi / 5, np.pi /20, 1000) 中期频段
#omega_range = np.linspace(np.pi / 20, np.pi /250, 1000) 长期频段
psi_omega_values = np.array([frequency_response(omega, psi_h) for omega in omega_range])  
psi_omega_values.shape

def frequency_response_derivative(omega, psi_h):  
    n_var = psi_h.shape[1]  # 变量数  
    n_lags = psi_h.shape[0]  # 滞后阶数  
    psi_prime_omega = np.zeros((n_var, n_var), dtype=complex)  
    for h in range(n_lags):  
        # 对每一项应用链式法则  
        term_derivative = -1j * h * psi_h[h] * np.exp(1j * omega * h)  # 注意这里 exp 的符号  
        psi_prime_omega += term_derivative  
    return psi_prime_omega

# 定义函数来计算广义因果谱  
def generalized_causality_spectrum(omega, psi_omega, sigma,psi_prime_omega):  
    n_var = psi_omega.shape[0]  
    xi = np.zeros((n_var, n_var))  
    sigma_array = sigma.values
    psi_omega_sigma = np.dot(psi_omega, sigma)
    #psi_omega_sigma_conj_T = np.conjugate(psi_omega).T
    #psi_omega_sigma_psi_omega_T = np.dot(psi_omega_sigma, psi_omega_sigma_conj_T)
    for p in range(n_var):  
        for q in range(n_var):  
            numerator = np.abs(psi_omega_sigma[p, q])**2   
            denominator = np.dot(psi_omega_sigma,psi_prime_omega[p,p])   
            if np.all(denominator != 0):  
                xi = numerator / denominator/sigma_array[p,p]    
    return xi  

def integrand(x, psi_omega_pp):   
    return np.abs(psi_omega_pp)**2  
#定义权重函数  
def calculate_gamma_p(omega_values, psi_h, sigma):  
    n_var = psi_h.shape[1]  
    gamma_p = np.zeros((n_var, len(omega_values)))  
  
    for p in range(n_var):  
        # 初始化 psi_omega 数组以存储所有频率的响应  
        psi_omega_all = np.zeros((len(omega_values), n_var, n_var), dtype=complex)  
        for i, omega in enumerate(omega_values):  
            psi_omega_all[i] = frequency_response(omega, psi_h) 
  
        # 对每个频率计算 psi_omega_pp 的积分  
        integrals = np.zeros(len(omega_values))  
        for i, omega in enumerate(omega_values):  
            psi_omega_pp = psi_omega_all[i, p, p]  
            integrals[i] = quad(lambda x: integrand(x, psi_omega_pp), -np.pi, np.pi)[0] 
  
        # 归一化  
        normalization = np.sum(integrals)  
        if normalization != 0:  
            gamma_p[p, :] = integrals * sigma.diagonal()[p] / normalization  
  
    return gamma_p  

# 示例：计算一系列频率下的 gamma_p  
sigma_array = sigma.values
gamma_p = calculate_gamma_p(omega_range, psi_h, sigma_array)  
print(gamma_p.shape)
gamma_p

omega = np.pi / 5 
psi_omega = frequency_response(omega, psi_h) 
psi_prime_omega = frequency_response_derivative(omega, psi_h)
xi_pq_values = np.array([generalized_causality_spectrum(omega, psi_omega, sigma,psi_prime_omega) for omega in omega_range])
print(xi_pq_values.shape)
xi_pq_values 

n_var = psi_omega.shape[0]
theta_d_pqs = np.zeros((n_var, n_var))
# 遍历 xi_pq_values 的后两个维度（即变量对）  
for i in range(n_var):  # 遍历变量 p  
    for j in range(n_var):  # 遍历变量 q  
        # 提取对应的 Gamma_p_values（假设第一个维度对应不同的 p）  
        Gamma_p_for_this_pq = gamma_p[i, :]  
  
        # 提取 xi_pq_values 中的对应 (p, q) 元素  
        xi_pq_for_this_pq = xi_pq_values[:, i, j]  
  
        # 数值积分计算 (theta_d)_{p,q}   
        theta_d_pq = (1 / (2 * np.pi)) * np.trapz(Gamma_p_for_this_pq * xi_pq_for_this_pq, dx=omega_range[1] - omega_range[0])  
  
        # 存储结果  
        theta_d_pqs[i, j] = theta_d_pq  
print(f"The spectral representation of prediction error variance matrix is:\n{theta_d_pqs}")  

# 标准化传染效应矩阵  
theta_d_normalized = theta_d_pqs / theta_d_pqs.sum()  
theta_d_normalized

Answer 1

以下回复参考：皆我百晓生、券券喵儿等免费微信小程序作答:

根据您提供的代码，我发现可能的问题是在计算广义因果谱时，分母可能出现零值。为了解决这个问题，您可以在计算广义因果谱之前检查分母是否为零，并在必要时进行修正。以下是修改后的代码：

def generalized_causality_spectrum(omega, psi_omega, sigma, psi_prime_omega):  
    n_var = psi_omega.shape[0]  
    xi = np.zeros((n_var, n_var))  
    sigma_array = sigma.values
    psi_omega_sigma = np.dot(psi_omega, sigma)
    for p in range(n_var):  
        for q in range(n_var):  
            numerator = np.abs(psi_omega_sigma[p, q])**2   
            denominator = np.dot(psi_omega_sigma, psi_prime_omega[p, p])   
            if np.all(denominator != 0):  
                xi[p, q] = numerator / denominator / sigma_array[p, p]    
    return xi  

此外，我还注意到在计算权重函数 calculate_gamma_p 时，您使用了 sigma.diagonal()[p] 作为归一化因子。这可能会导致问题，因为 sigma 是一个协方差矩阵，其对角线元素表示各个变量的方差。在这种情况下，您应该使用 np.sqrt(sigma.diagonal()) 作为归一化因子。以下是修改后的代码：

def calculate_gamma_p(omega_values, psi_h, sigma):  
    n_var = psi_h.shape[1]  
    gamma_p = np.zeros((n_var, len(omega_values)))  
  
    for p in range(n_var):  
        # 初始化 psi_omega 数组以存储所有频率的响应  
        psi_omega_all = np.zeros((len(omega_values), n_var, n_var), dtype=complex)  
        for i, omega in enumerate(omega_values):  
            psi_omega_all[i] = frequency_response(omega, psi_h) 
  
        # 对每个频率计算 psi_omega_pp 的积分  
        integrals = np.zeros(len(omega_values))  
        for i, omega in enumerate(omega_values):  
            psi_omega_pp = psi_omega_all[i, p, p]  
            integrals[i] = quad(lambda x: integrand(x, psi_omega_pp), -np.pi, np.pi)[0] 
  
        # 归一化  
        normalization = np.sum(integrals)  
        if normalization != 0:  
            gamma_p[p, :] = integrals * np.sqrt(sigma.diagonal()) / normalization  
  
    return gamma_p  

希望这些修改能解决您的问题。如果还有其他问题，请随时告诉我。