在python中实现混频动态因子模型（mixed frequency dynamic factor model）
本文使用的代码引用自 https://github.com/genekindberg/DFM-Nowcaster ，目的是通过对代码的解释和运行加深自己对代码的理解。
不会使用GITHUB的话文末有压缩包。
1.1.1 函数代码
def _InterpolateNaN(Tseries):
not_nan, indices = np.logical_not(np.isnan(Tseries)), np.arange(len(Tseries)).reshape(-1,1)
Tseries_fill = np.interp(indices, indices[not_nan], Tseries[not_nan])
return Tseries_fill1.1.2 命令解释
not_nan 是通过 np.logical_not(np.isnan(Tseries)) 生成的，它首先使用 np.isnan(Tseries) 找到 Tseries 中的NaN值，然后通过 np.logical_not 对其进行逻辑取反，得到一个布尔数组，标记了哪些位置不是NaN。
indices 则是通过 np.arange(len(Tseries)).reshape(-1,1) 生成的，它使用 np.arange(len(Tseries)) 创建了一个包含从0到len(Tseries)-1的整数的一维数组，然后通过 reshape(-1,1) 将其调整为列向量的形状。
indices: 这是插值点的位置，是一个包含整数序列的一维数组。这些整数表示时间序列的索引。
indices[not_nan]: 这是有效插值点的位置，由之前创建的 not_nan 布尔数组确定。只有 not_nan 中对应位置为True的索引才会用于插值。
Tseries[not_nan]: 这是在有效插值点上的时间序列值。由于 not_nan 是一个布尔数组，它会选择那些对应位置为True的时间序列值，即非NaN值。
np.interp 将使用有效插值点的位置 indices[not_nan] 和对应的时间序列值 Tseries[not_nan] 来进行线性插值，生成一个新的时间序列 Tseries_fill。这样，Tseries_fill 中的NaN值就被用相邻的有效值进行了插值填充。
整个操作旨在通过线性插值填充时间序列(Tseries)中的NaN（非数值）值。这段代码是为了保证季度数据和月度数据维度大小一致，但是使用的是线性插值的方法，将用于最终的残差阵的产生。
1.1.3 具体效果
数据集使用了代码作者提供的数据集，一个excel文件，内含四个列表，第一个列表为GDP，第一列为季度数据，第二列为对应GDP数据，后三个列表分别为emp，cons，ipexp，它们的第一列为月度数据，第二第三列为各自对应的数据。
在spyder中读取为dataframe格式，命名为GDP_dat。
单独将第二列GDP数据读取为array格式，命名为GDP。
为进行线性填充，执行下列代码，扩展GDP，命名为GDPnan。
GDPnan = np.asarray(np.kron(GDP, [[np.nan], [np.nan], [1]])
再执行插值代码，线性插值GDPnan，命名为GDP_fill。
GDP_fill = _InterpolateNaN(GDPnan)1.2.1 函数代码
def _NormDframe(df):
df2 = df.copy()
df2 = (df2 - df2.mean())/df2.std()
return df21.2.2 命令解释
df2 = df.copy(): 这行代码创建了数据框 df 的一个副本 df2。这是为了避免修改原始数据框，而是在副本上进行标准化操作。
df2 - df2.mean(): 从每个元素中减去对应列的均值。这样做的效果是将每列的均值变为0。
df2.std(): 除以对应列的标准差。这样做的效果是将每列的标准差变为1。
整个操作实现了标准化，将每个列的数值缩放到均值为0，标准差为1的范围内。这段代码将被用于标准化月度数据，目的是进行主成分分析。
1.2.3 具体效果
执行示例代码。
import pandas as pdimport numpy as np
data = {'A': [1, 2, 3, 4, 5],
'B': [5, 4, 3, 2, 1],
'C': [10, 20, 30, 40, 50]}
df = pd.DataFrame(data)
def _NormDframe(df):
df2 = df.copy()
df2 = (df2 - df2.mean()) / df2.std()
return df2
normalized_df = _NormDframe(df)
print("原始数据框：\n", df)
print("\n标准化后的数据框：\n", normalized_df)此处手动输入的示例数据为dictionary格式，实际运行时没有这个格式转换的过程。
将数据转化为dataframe格式，调用函数标准化，命名为normalized_df。
1.3.1 函数代码
def _PCA_fact_load(data):
pca = PCA(n_components=1)
factor = pca.fit_transform(data)
loadings = pca.components_.T * np.sqrt(pca.explained_variance_)
return factor, loadings1.3.2 命令解释
pca = PCA(n_components=1): 创建了一个PCA对象，设置主成分数量为1。这意味着我们只找到数据的第一主成分。
factor = pca.fit_transform(data): 使用 fit_transform 方法对数据进行主成分分析，找到主成分，并将数据变换为主成分空间。factor 变量现在包含了数据的第一主成分。
loadings = pca.components_.T * np.sqrt(pca.explained_variance_): 计算了主成分的因子载荷。在PCA中，主成分的方向由pca.components_给出，而载荷则表示数据的每个变量对主成分的贡献。np.sqrt(pca.explained_variance_) 用于乘以标准差，确保载荷的缩放是正确的。结果保存在 loadings 变量中。
整个操作目的是对输入的数据进行主成分分析（PCA），并返回主成分分析的结果，包括因子得分（factor）和载荷（loadings）。将用于处理标准化后的月度数据，需要注意的是这里设置了主成分数量为一，所以最终月度数据的数量就等于表的数量，每个表只会被处理出一个主成分因子。
1.3.3 具体效果
此为数据集中的第二个表，月度数据。
第二第三第四个表，读入spyder后为list格式，包含三个dataframe，且不包含日期，命名为MonthlyDat。
执行代码，进行主成分分析。
from sklearn.decomposition import PCA
def _PCA_fact_load(data):
pca = PCA(n_components=1)
factor = pca.fit_transform(data)
loadings = pca.components_.T * np.sqrt(pca.explained_variance_)
return factor, loadings
cutoff = len(GDP_fill) #得到数据的长度
Monthly = MonthlyDat #获取月度数据
MonthlyDat_norm = [] #创建储存表
for df in Monthly: #遍历月度数据逐个进行主成分分析
MonthlyDat_norm.append(_NormDframe(df[0:cutoff]))
loadings = [None]*len(Monthly) #创建储存表
Factors = [None]*len(Monthly) #
for df,num in zip(MonthlyDat_norm, range(len(MonthlyDat_norm))):
Factors[num], loadings[num] = _PCA_fact_load(df[0:cutoff].fillna(method='ffill'))#储存主成分因子和因子载荷首先执行标准化，得到标准化月度数据的list格式，命名为MonthlyDat_norm。
这里对第一个dataframe进行展示，这是标准化后的结果。
最终输出list格式的第一主成分因子，命名为Factors。
对应的因子载荷，命名为loadings。
这里输出的因子载荷是一样的，有些奇怪，但是经过下面的检查测试之后发现确实如此。
import numpy as npfrom sklearn.decomposition import PCA
data = MonthlyDat_norm[1]#或2或3
pca = PCA(n_components=1)
factor = pca.fit_transform(data)
loadings = pca.components_.T * np.sqrt(pca.explained_variance_)
print("因子载荷：", loadings)
输出的因子载荷是相同的，可能是原数据具有完全的相关性。
数据表2的简易散点图。
询问了一位同学，经过她的验证代码本身是没有问题的，数据正交或者维度过低（2）也可能是原因。
1.3.4 主成分分析
主成分分析（Principal component analysis，PCA）又称主分量分析是一种常见的多元统计方法，常用于数据降维。它的如图所示，原始变量显然是线性相关的，黑线表示旋转后的坐标轴。在新的坐标轴下，原始变量可以近似用一个变量替代，从而在不损失较多信息的条件下达到了降低数据维度的目的。基本思想在于利用坐标轴的旋转，将线性相关的变量用少数几个线性无关的变量进行近似表示。
2.1.1 函数代码
def _MakeLags(mat, lags):
cols = mat.shape[1] #获取矩阵 mat 的列数
Matlag = np.zeros((len(mat)-lags,cols*lags)) #创建一个零矩阵，这个矩阵将用于存储创建的滞后矩阵
for ii in range(lags, 0, -1): #迭代从 lags 到 1 的整数
Matlag[:,0+cols*(lags-ii):cols*(lags-ii+1)] = mat[ii-1:(len(mat)-lags+ii-1),:]; #将矩阵 mat 中的每一列进行滞后处理，并存储在 Matlag 的相应位置
mat2 = mat.copy() #复制原始矩阵 mat
mat2 = mat2[lags:,:] #对 mat2 保留从第 lags 行到最后一行的数据，前 lags 行将被丢弃
return mat2, Matlag 2.1.2 命令解释
这个函数 _MakeLags 的目的是将输入的矩阵 mat 进行滞后处理，生成一个包含滞后值的新矩阵 mat2 和一个包含滞后矩阵的矩阵 Matlag。 这段函数在整个代码中将会被多次使用，比如VAR模型中，GIBBS抽样中，以下是一个简单的例子，更具体效果会在函数被使用时再进行解释。滞后为几，生成的矩阵宽度就会乘几倍，长度就会缩短几。
2.1.3 具体效果
import numpy as np
np.random.seed(42)
mat = np.random.rand(5, 3) # 生成一个示例矩阵
lags = 1 #或2或3或4，设置滞后步数
mat2, Matlag = _MakeLags(mat, lags) # 使用 _MakeLags 函数
print("原始矩阵:")# 打印结果
print(mat)
print("\n滞后矩阵:")
print(Matlag)
print("\n移除", lags, "行后的矩阵:")
print(mat2)
原始矩阵:[[0.375 0.951 0.732]
[0.599 0.156 0.156]
[0.058 0.866 0.601]
[0.708 0.021 0.97 ]
[0.832 0.212 0.182]]
滞后矩阵:[[0.375 0.951 0.732]
[0.599 0.156 0.156]
[0.058 0.866 0.601]
[0.708 0.021 0.97 ]]
移除 1 行后的矩阵:[[0.599 0.156 0.156]
[0.058 0.866 0.601]
[0.708 0.021 0.97 ]
[0.832 0.212 0.182]]
滞后矩阵:[[0.599 0.156 0.156 0.375 0.951 0.732]
[0.058 0.866 0.601 0.599 0.156 0.156]
[0.708 0.021 0.97 0.058 0.866 0.601]]
移除 2 行后的矩阵:[[0.058 0.866 0.601]
[0.708 0.021 0.97 ]
[0.832 0.212 0.182]]
滞后矩阵:[[0.058 0.866 0.601 0.599 0.156 0.156 0.375 0.951 0.732]
[0.708 0.021 0.97 0.058 0.866 0.601 0.599 0.156 0.156]]
移除 3 行后的矩阵:[[0.708 0.021 0.97 ]
[0.832 0.212 0.182]]
原始矩阵:[[0.375 0.951 0.732]
[0.599 0.156 0.156]
[0.058 0.866 0.601]
[0.708 0.021 0.97 ]
[0.832 0.212 0.182]]
滞后矩阵:[[0.708 0.021 0.97 0.058 0.866 0.601 0.599 0.156 0.156 0.375 0.951 0.732]]
移除 4 行后的矩阵:[[0.832 0.212 0.182]]2.2.1 VAR模型
VAR 向量自回归模型（Vector Autoregression Model）是一种用于建模多个时间序列变量之间相互关系的统计模型。它是一种常用于时间序列分析和经济学领域的方法。VAR 模型基于向量的形式，能够捕捉变量之间的动态关系和相互影响。VAR 模型的基本思想是将一个时间序列变量的当前值用其过去值的线性组合表示，并考虑多个变量之间的相互作用。具体而言，一个 p 阶的 VAR(p) 模型可以表示。
VAR 模型的阶数 p 决定了模型中考虑的滞后阶数，即过去的观测对当前观测的影响。选择合适的阶数通常需要通过模型拟合的评估指标或信息准则来确定。
VAR 模型的参数估计通常使用最小二乘法或贝叶斯方法。估计得到的系数矩阵可以用于分析变量之间的动态关系，进行预测或冲击响应分析。这里通过最小二乘法估计 VAR 模型的系数矩阵，实际上是在寻找一个最优的线性模型，使得模型的预测值与实际观测值的残差平方和最小。
以下是通过最小二乘法估计 VAR 模型系数的具体步骤。
这样，通过以上步骤，就可以使用最小二乘法估计 VAR 模型的系数矩阵。最终得到的系数矩阵可以用于分析各个变量之间的线性关系，以及进行模型的预测和推断。
2.2.2 函数代码
def _estVAR(Yinput,lags):
Y, X = _MakeLags(Yinput, lags) #创建滞后矩阵，它的输出分别为 Y（处理后的原始数据）和 X（滞后矩阵）
N = Y.shape[0] #获取数据点的数量
X = np.c_[X, np.ones((X.shape[0],1))] #将滞后矩阵 X 增加一列常数项，这是为了估计模型中的截距
B = np.linalg.inv(X.T @ X) @ (X.T @ Y) #计算残差，即实际值与预测值之间的差异
Q = np.divide((e.T @ e),(N-lags)) #计算残差的方差-协方差矩阵的无偏估计， np.divide对两个数组中的元素进行逐元素的除法操作
return B, Q2.2.3 命令解释
这个函数 _estVAR 用于估计VAR（Vector Autoregression）模型的系数矩阵 B 和残差的方差-协方差矩阵 Q。这个函数的输出结果可以用于VAR模型的后续分析和预测。这个函数将被用于且只用于后面的初始化函数中，用于处理已经合并的年度和月度数据，输出两个矩阵将被用于卡尔曼滤波函数。
2.2.4 具体效果
下面是在后文的初始化函数中的使用。GDP_fill在前文中已给出，是经过线性填充后的季度GDP数据。Factors在前文中已给出，是经过主成分分析后的第一主成分因子。
F0 = np.concatenate(Factors,1) #拼接三个主成分因子
F0 = np.hstack((GDP_fill, F0)) #拼接
lags = 6
F0_cut, F0lag = _MakeLags(F0, lags)
B, Q = _estVAR(F0_cut,lags) 拼接后的主成分因子和GDP，格式为array，命名为F0。
滞后设置为6，前五行被切掉后，移除五行后矩阵，格式为array，命名为F0_cut。
滞后矩阵，格式为array，命名为F0lag。
通过最小二乘法估计VAR模型的系数矩阵，格式为array，命名为B。
系数矩阵B是由协方差矩阵的逆矩阵和X、Y之间的关系计算而来。它表示着Y的变化如何由X的变化线性解释。意味着第一个变量受到自己过去值和另外三个过去值的影响。这个矩阵提供了关于变量之间动态关系的信息。正值表示正相关，负值表示负相关。值的绝对值越大，表示影响越显著。
对残差方差-协方差矩阵的无偏估计格式为array，命名为Q。
Q 提供了关于模型残差的信息，即实际观测值与模型预测值之间的差异。方差-协方差矩阵的无偏估计可用于评估模型的拟合质量，了解模型对数据的拟合程度。对角线上的元素表示各个变量的方差，非对角线上的元素表示变量之间的协方差。如果模型拟合得好，那么残差的方差应该足够小，非对角线上的协方差应该接近零，比如前两个变量。如果某个对角线上的值远离零，可能表示模型在该变量上存在较大的残差，即模型未能很好地拟合这个变量，比如后两个变量。最终，Q 将是一个对残差方差-协方差矩阵的无偏估计。这个矩阵在 VAR 模型中非常重要，因为它反映了模型中残差的统计性质，对于后续的模型分析和推断具有重要意义。
2.3.1 BVAR模型
BVAR 模型是贝叶斯向量自回归模型（Bayesian Vector Autoregressive Model）的缩写。它是一种用于描述多个时间序列变量之间动态关系的统计模型。BVAR 模型的主要特点是通过引入先验分布来捕捉未观测的数据结构，从而提供对模型参数后验分布的完整概率描述。
先验分布的设定： 在 BVAR 模型中，通常对系数矩阵的先验使用低秩先验，例如 Minnesota 先验。这种先验对于大规模系统具有稀疏性，促使模型更倾向于简单结构。除了Minnesota先验之外，贝叶斯VAR模型中的其他常见先验方法包括：狄利克雷分布先验（Dirichlet Process Prior）： 该先验用于对VAR模型的系数矩阵引入非参数化先验，允许系数矩阵具有更大的灵活性，适用于不确定性较高的情况。协方差逆矩阵Wishart分布先验： 用于对VAR模型的误差协方差矩阵引入先验信息。通过设定Wishart分布的参数，可以控制误差协方差矩阵的先验强度。脊回归（Ridge Regression）： 将贝叶斯VAR模型与脊回归方法结合，通过引入L2正则化项，可以缓解多重共线性问题。先验因子分析（Bayesian Factor Analysis）： 使用因子分析的思想，对VAR模型的系数矩阵进行建模，引入先验因子以捕捉变量之间的共享信息。
Gibbs 抽样： 为了估计 BVAR 模型，通常使用 Gibbs 抽样等贝叶斯方法，从后验分布中抽取参数值。这使得我们能够获得参数的后验分布，并提供对不确定性的更全面认识。总体而言，BVAR 模型提供了一种灵活的方式来建模多变量时间序列数据，通过引入贝叶斯框架，能够更好地处理参数估计的不确定性，并对模型进行更全面的推断。
2.3.2 函数代码
def _setVARpriors(X,Y,lambda1, lambda3, lambda4, p, alpha0)
n = Y.shape[1] #响应变量的数量
T = X.shape[0] #数据点的数量
k = X.shape[1] #X 的列数
m = k-n*p #截距项的数量
q=n*k; #系数矩阵的总数
arvar = np.empty((n,1)) #生成需要的空矩阵
arvar[:] = np.nan #
ar = np.empty((n,1)) #
ar[:] = np.nan #
for ii in range(0,n,1): #对每个响应变量进行循环
Ytemp=Y[:,ii];
Xtemp=X[:,ii:n*p:n];
B= np.linalg.inv(Xtemp.T @ Xtemp) @ (Xtemp.T @ Ytemp) #估计 AR 模型的系数 B
ar[ii] = B[0] #提取 AR 模型的截距项
eps=Ytemp - (Xtemp @ B) #计算残差
arvar[ii]=np.divide((eps.T @ eps),(T-p)) #计算残差的方差-协方差
beta0=np.zeros((q,1)); #将 AR 模型的截距项设置到 beta0 中
for ii in range(0,n,1): #
beta0[(ii)*k+ii]=ar[ii] #
phi0=np.zeros((k,k)) #创建零矩阵存储模型的自回归（AR）部分的协方差矩阵的先验
for ii in range(0,n,1): #对每个响应变量进行循环
for jj in range(0,p,1): #对滞后阶数中的每个滞后项进行循环
phi0[(jj)*n+ii,(jj)*n+ii]=(1/arvar[ii])*((lambda1/((jj+1)**lambda3))**2)#左侧是位置，右侧是数值
if k>n*p: #检查是否存在外生方差，k 大于 n*p，表示外生变量的数量大于滞后项的数量
m = k-n*p
for ii in range(0,m,1): #对外生变量的数量进行循环，设置截距项的方差-协方差的先验
phi0[k-m+ii,k-m+ii]=(lambda1*lambda4)^2 #左侧是位置，右侧是数值
S0=np.diag(arvar.flatten()); #观测误差的方差-协方差矩阵的先验
return phi0, alpha0, S0, beta02.3.3 命令解释
这个函数的输出结果将被用且仅用于吉布斯抽样中的_EstBVAR函数即设置 VAR 模型的先验参数，注意与_estVAR进行区分。_setVARpriors： 这个函数用于设置 VAR 模型的先验信息，包括自回归部分（系数矩阵）的先验和观测误差的方差-协方差矩阵的先验。它并不执行模型估计，而是为 VAR 模型提供了先验信息，例如对自回归系数和观测误差的先验分布的设定。这个方法是一种设定贝叶斯向量自回归（BVAR）模型的先验（prior）的方法。下面是上述代码中使用的 Minnesota 先验的数学公式。
2.3.4 具体效果
给定先验参数。
lambda1 = 0.5
lambda3 = 1
lambda4 = 1000
p = lags = 6
alpha0 = 1
K = 3
Qs = 1X和Y分别为S_X和S_Y，来源于滞后函数，S[:, 0:K+Qs]：这部分表示从矩阵 S 中选择所有的行 ，以及从第 0 列到第 K+Qs-1 列，形成一个子矩阵。这样，你获得了 S 中的一部分数据，这部分数据包括 K+Qs 列。
S_Y, S_X = _MakeLags(S[:,0:K+Qs], lags)S来源于初始化函数中的卡尔曼滤波返回值--空间状态估计值。
def initializevals(GDP, Monthly, lags, lagsH, K, Qs):
......
S, P = _Kfilter(XY, F, H, Q, R, F0lag)
return XY, H, R, F, Q, S, ints初始化函数将在后面介绍，这里先将代码执行，并观察数据内容，与前文的EstVAR数据类似。
继续执行对应代码得到S_Y, S_X。
至此，我们解释了该函数的输入，接下来我们将给出该部分代码的返回值phi0, alpha0, S0, beta0的一些解释。
n = Y.shape[1] # 时间序列的数量，列数4
T = X.shape[0] # 时间序列的长度，行数297
k = X.shape[1] # 包含所有额外虚拟变量的维度，列数20
m = k-n*p # 系统中未知数的数量
q=n*k; # 需要估计的系数总数alpha0：这个参数用于观测误差的先验。观测误差是模型的输出与真实观测之间的差异。alpha0 表示观测误差的先验信息，由于本段代码并没有实际的使用到这个参数，所以我们暂不进行更详细的展示。
phi0：这是 VAR 模型的自回归（AR）部分的协方差矩阵的先验。这个矩阵表示时间序列变量之间的协方差关系，其中包括滞后项的影响。在贝叶斯 VAR 模型中，先验用于描述我们对这些协方差关系的先验信念。返回的 phi0 就是这些先验信念的表达。
phi0=np.zeros((k,k))
for ii in range(0,n,1):
for jj in range(0,p,1):
phi0[(jj)*n+ii,(jj)*n+ii]=(1/arvar[ii])*((lambda1/((jj+1)**lambda3))**2)
#下面的循环因为m是零，实际上并不会运行
if k>n*p:
m = k-n*p
for ii in range(0,m,1):
phi0[k-m+ii,k-m+ii]=(lambda1*lambda4)^2 
beta0：这是截距项的先验。在 VAR 模型中，包含截距项，beta0 描述了这些截距项的先验分布。
beta0=np.zeros((q,1));
for ii in range(0,n,1):
beta0[(ii)*k+ii]=ar[ii]下表每25行有一个数字。
S0：这是观测误差的方差-协方差矩阵的先验。它描述了观测误差的方差和协方差的先验分布。
S0=np.diag(arvar.flatten())
个别返回值的含义尚不能理解
3.1.1 函数代码
def _EstBVAR(Y,lags, phi0, alpha0, S0, beta0, draws):
Y2, X = _MakeLags(Y, lags) #调用 _MakeLags 函数，创建包含滞后项的矩阵 X。Y2 是移除了滞后期的目标变量矩阵
T = X.shape[0]
n = Y2.shape[1]
k = X.shape[1]
B0=np.reshape(beta0,(k,n), order="F") #将先验系数 beta0 重塑为一个矩阵 B0，其中 k 是所有滞后项的总数，n 是内生变量的数量
invphi0=np.diag(1./np.diag(phi0)) #计算 phi0 的逆矩阵的对角矩阵
invphibar=invphi0+(X.T @ X) #计算先验协方差矩阵的逆矩阵
C=np.linalg.cholesky(invphibar).T #计算 invphibar 的 Cholesky 分解，然后计算其逆矩阵的 Cholesky 分解
invC=np.linalg.inv(C) #
phibar=invC @ invC.T #计算先验协方差矩阵的 Cholesky 分解的逆
Bbar=phibar @ ((invphi0 @ B0) + (X.T @ Y2)) #计算贝叶斯估计的系数均值
betabar=Bbar.flatten().T #将 Bbar 展平为一维数组，得到 betabar
alphabar=T+alpha0 #计算 alphabar，这是先验中参数
Sbar=(Y2.T @ Y2)+S0+(B0.T @ invphi0 @ B0)-(Bbar.T @ invphibar @ Bbar) #算先验中的一个协方差矩阵的函数
sigma_gibbs = np.full((draws,n,n),np.nan) #初始化用于存储抽样结果的数组
beta_gibbs =np.full((draws,k,n),np.nan) #
for ii in range(0,draws):
C=np.linalg.cholesky(Sbar)
Z = np.random.multivariate_normal(np.full(n,0), np.diag(np.full(n,1)), alphabar)#从标准正态分布中抽取随机向量 Z
drawSig=(C @ np.linalg.inv(Z.T @ Z)) @ C.T #计算协方差矩阵 drawSig
sigma_gibbs[ii,:,:]=drawSig #存储 drawSig 到 sigma_gibbs
C=np.linalg.cholesky(drawSig) #重新计算 Cholesky 分解
P=np.linalg.cholesky(phibar)
w=np.random.normal(0, 1, k*n) #生成随机向量 W
W=np.reshape(w,(k,n))
drawBeta=Bbar+((P @ W) @ C.T) #计算系数向量 drawBeta
drawB = np.r_[drawBeta[0:,0:].T, np.c_[np.eye((n)*(lags-1)), np.zeros(((n)*(lags-1),(n)))]]
count = 1
while np.max(np.abs(np.linalg.eigvals(drawB)))>0.999: #将系数矩阵 drawB 合并
W = np.reshape(np.random.normal(0, 1, k*n),(k,n))
drawBeta=Bbar+((P @ W) @ C.T)
drawB = np.r_[drawBeta[0:,0:].T, np.c_[np.eye((n)*(lags-1)), np.zeros(((n)*(lags-1),(n)))]]
count = count+1
if count>10000: #检查 drawB 是否满足平稳性条件，若不满足则重新生成，最多尝试 10000 次
raise Exception('VAR not stationary!')
beta_gibbs[ii,:,:]=drawBeta #存储 drawBeta 到 beta_gibbs
return sigma_gibbs, beta_gibbs #返回抽样结果 sigma_gibbs 和 beta_gibb3.1.2 命令解释
这是一个贝叶斯向量自回归（BVAR）模型的蒙特卡罗马尔科夫链蒙特卡罗（MCMC）抽样函数。在这个函数中，BVAR模型的参数（包括系数矩阵和协方差矩阵）通过MCMC方法从后验分布中抽样得到。以下是主要步骤的简要解释：
数据准备： 使用_MakeLags函数为目标变量创建包含滞后项的矩阵。然后，将目标变量矩阵移除滞后期，得到Y2。这些数据将用于贝叶斯估计。
先验设置： 使用_setVARpriors函数设置VAR模型的先验。该函数计算AR模型的系数、截距项和协方差矩阵的先验分布。
参数初始化： 将先验系数beta0重塑为矩阵B0。计算先验协方差矩阵的逆矩阵invphi0、invphibar和其Cholesky分解。
MCMC抽样： 对于每次MCMC迭代，进行以下步骤：
计算协方差矩阵Sbar，其中包括数据的协方差、先验截距项的协方差和先验系数的协方差。从标准正态分布中抽取随机向量Z，然后计算协方差矩阵drawSig。使用Cholesky分解重新计算drawSig的Cholesky分解C。生成随机向量W，计算系数向量drawBeta。将系数矩阵drawB合并，确保其特征值的绝对值小于0.999，最多尝试10000次。
存储抽样结果。
返回结果： 返回抽样得到的协方差矩阵序列sigma_gibbs和系数矩阵序列beta_gibbs。
这个函数的目的是获得贝叶斯估计的系数和协方差矩阵，通过MCMC方法从后验分布中抽样得到。
3.2.1 函数代码
def _Kfilter(Y, F, H, Q, R, F0):
n = Y.shape[1] # 内生变量的数量
k = F0.shape[1] # 状态向量的维度
T = Y.shape[0] # 观测序列的长度
P = np.full((T, k**2), np.nan) # 用于存储系统方程误差协方差矩阵的数组
S = np.full((T, k), np.nan) # 用于存储隐藏状态的数组
Sp = F0[0, :] # 初始状态向量
Pp = np.eye(k) # 初始系统方程误差协方差矩阵
for i in range(0, T, 1):
y = Y[i, :].T.copy() # 当前时刻的观测向量
if np.isnan(y).any():
replace = H @ Sp # 如果观测值中存在缺失值，用系统方程的预测值替代
y[np.where(np.isnan(y))] = replace[np.where(np.isnan(y))]
nu = y - (H @ Sp) # 观测残差
f = ((H @ Pp) @ H.T) + R # 观测方程误差协方差矩阵
finv = np.linalg.inv(f) # 观测方程误差协方差矩阵的逆矩阵
Stt = Sp + (Pp @ H.T) @ (finv @ nu) # 隐藏状态的更新
Ptt = Pp - ((Pp @ H.T) @ (finv @ H) @ Pp)# 系统方程误差协方差矩阵的更新
if i < T-1:
Sp = F @ Stt # 下一时刻的预测状态
Pp = ((F @ Ptt) @ F.T) + Q # 下一时刻的系统方程误差协方差矩阵
S[i, :] = Stt.T # 存储当前时刻的隐藏状态
P[i, :] = np.reshape(Ptt, (k**2)) # 存储当前时刻的系统方程误差协方差矩阵
return S, P3.2.2 命令解释
该函数执行卡尔曼滤波（Kalman Filter）操作，用于估计状态空间模型的隐藏状态。
3.3.1 函数代码
def _Ksmoother(F, H, Q, R, S, P, lags, n):
T = S.shape[0] # 观测序列的长度
k = n * lags # 状态向量的维度
Qstar = Q[0:n, 0:n]
Fstar = F[0:n, :]
Pdraw = np.full((T, k**2), np.nan) # 用于存储系统方程误差协方差矩阵的数组
Sdraw = np.full((T, k), np.nan) # 用于存储平滑后的隐藏状态的数组
for ii in range(T-1, -1, -1):
Sf = S[ii, 0:n].T # 当前时刻的卡尔曼滤波估计值
Stt = S[ii-1, :] # 上一时刻的隐藏状态
Ptt = np.reshape(P[ii-1, :], (k, k)) # 上一时刻的系统方程误差协方差矩阵
f = (Fstar @ Ptt @ Fstar.T) + Qstar # 系统方程误差协方差矩阵
finv = np.linalg.inv(f) # 系统方程误差协方差矩阵的逆矩阵
nu = Sf - (Fstar @ Stt) # 残差
Smean = Stt + (Ptt @ Fstar.T @ finv @ nu) # 平滑后的隐藏状态的均值
Svar = Ptt - (Ptt @ Fstar.T @ finv @ Fstar @ Ptt) # 平滑后的隐藏状态的方差
Sdraw[ii, :] = Smean.T # 存储平滑后的隐藏状态
Pdraw[ii, :] = np.reshape(Svar, (k**2)) # 存储平滑后的系统方程误差协方差矩阵
return Sdraw, Pdraw, Svar3.3.2 命令解释
总体来说，该函数实现了卡尔曼平滑算法，通过递归地从后向前处理估计状态向量，得到每个时间点的平滑处理结果。函数的输入包括状态转移矩阵 F、观测矩阵 H、状态方差 Q、观测方差 R、以及卡尔曼滤波的输出 S、P。函数的输出是平滑处理后的估计状态向量、状态协方差矩阵和最终的协方差矩阵。
前向滤波 (Kalman Filtering)。
目的： 用于递归地估计系统的状态，即给定当前的观测值和之前的状态估计，计算下一个时刻的状态估计。
过程： 通过观测值和状态方程，计算当前时刻的预测（先验）状态和预测（先验）协方差矩阵。然后，通过与观测值的比较，计算卡尔曼增益，最终得到当前时刻的状态估计和协方差矩阵。
后向平滑 (Kalman Smoothing)。
目的： 在得到全部观测信息后，通过后向迭代来提高系统状态的估计精度。
过程： 从最后一个时刻开始，使用后向迭代递归地计算每个时刻的平滑（后验）状态和平滑（后验）协方差矩阵。这一过程考虑了整个时间序列的信息，从而提供更准确的状态估计。
这两个步骤共同构成了卡尔曼滤波与平滑算法，允许在动态系统中进行实时状态估计，并在观测完整序列后对状态进行更精确的估计。
3.3.4 具体效果
卡尔曼相关概念不易理解，这里举一个简单的例子。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42) # 生成模拟数据
F = np.array([[1, 1],
[0, 1]]) # 状态转移矩阵 F
Q = np.array([[0.01, 0.],
[0., 0.01]]) # 状态转移协方差矩阵 Q
H = np.array([[1, 0]]) # 观测矩阵 H
R = np.array([[0.1]]) # 观测协方差矩阵 R
initial_state = np.array([0, 0]) # 初始状态
T = 100 # 生成状态和观测序列
true_states = np.zeros((T, 2))
observations = np.zeros((T, 1))
state = initial_state
for t in range(T):
true_states[t] = state
observation = H @ state + np.random.normal(0, np.sqrt(R[0, 0]))
observations[t] = observation
state = F @ state + np.random.multivariate_normal([0, 0], Q)
def kalman_filter(Y, F, H, Q, R, initial_state): # 卡尔曼滤波
T = len(Y)
k = len(initial_state)
S = np.zeros((T, k))
P = np.zeros((T, k, k))
S[0] = initial_state
P[0] = Q
for t in range(1, T):
S[t] = F @ S[t-1] # 预测
P[t] = F @ P[t-1] @ F.T + Q
K = P[t] @ H.T @ np.linalg.inv(H @ P[t] @ H.T + R) # 更新
S[t] = S[t] + K @ (Y[t] - H @ S[t])
P[t] = P[t] - K @ H @ P[t]
return S, P
filtered_states, _ = kalman_filter(observations, F, H, Q, R, initial_state)# 使用卡尔曼滤波估计状态
def kalman_smoother(F, P, filtered_states): # 卡尔曼平滑
T, k = filtered_states.shape
S_smooth = np.zeros((T, k))
S_smooth[-1] = filtered_states[-1]
for t in range(T-2, -1, -1):
J = P[t] @ F.T @ np.linalg.inv(F @ P[t] @ F.T + Q)
S_smooth[t] = filtered_states[t] + J @ (S_smooth[t+1] - F @ filtered_states[t])
return S_smooth
smoothed_states = kalman_smoother(F, filtered_states, filtered_states) # 使用卡尔曼平滑估计状态
plt.figure(figsize=(10, 6)) # 绘制结果
plt.plot(true_states[:, 0], label='True Position', linestyle='--')
plt.plot(observations, label='Observations', marker='o', linestyle='None')
plt.plot(filtered_states[:, 0], label='Filtered Position', linestyle='-', alpha=0.7)
plt.plot(smoothed_states[:, 0], label='Smoothed Position', linestyle='-', linewidth=2)
plt.title('Kalman Filter and Smoother Example')
plt.legend()
plt.show()
4.1.1 函数代码
def _RemoveOutliers(Data, SDs, Qs):
for col in range(0,Data.shape[1]): #对数据矩阵中的每一列进行循环
Data[(Data[:,col] - np.nanmean(Data[:,col])) > SDs * np.nanstd(Data[:,col]), col] = np.nanmean(Data[:,col])+SDs * np.nanstd(Data[:,col])
Data[(Data[:,col] - np.nanmean(Data[:,col])) < -(SDs * np.nanstd(Data[:,col])), col] = np.nanmean(Data[:,col])-SDs * np.nanstd(Data[:,col])
return Data4.1.2 命令解释
Data[(Data[:,col] - np.nanmean(Data[:,col])) > SDs * np.nanstd(Data[:,col]), col] = np.nanmean(Data[:,col])+SDs * np.nanstd(Data[:,col])：对当前列的数据进行处理：计算数据列的均值 np.nanmean(Data[:,col]) 和标准差 np.nanstd(Data[:,col])。将大于指定标准差倍数 SDs 的值替换为均值加上标准差倍数。这一步的目的是将超过指定标准差倍数的异常值截断到一个阈值。
Data[(Data[:,col] - np.nanmean(Data[:,col])) < -(SDs * np.nanstd(Data[:,col])), col] = np.nanmean(Data[:,col])-SDs * np.nanstd(Data[:,col])：对当前列的数据进行处理：将小于指定标准差倍数的值替换为均值减去标准差倍数。这一步的目的是将低于指定标准差倍数的异常值截断到一个阈值。
总体来说，该函数使用指定的标准差倍数 SDs 来检测和处理每一列数据中的异常值。将大于指定倍数的异常值替换为均值加上倍数的标准差，将小于指定倍数的异常值替换为均值减去倍数的标准差。这样的处理方式有助于减小异常值对整体分析的影响。
4.1.3 具体效果
执行一个随机的例子展示函数效果。
import numpy as np # 生成一些包含异常值的示例数据
np.random.seed(42)
data = np.random.normal(loc=0, scale=1, size=(100, 3))
data[10, 0] = 10 # 添加一个异常值
print("Original Data Statistics:") # 查看原始数据的统计信息
print("Mean:", np.nanmean(data, axis=0))
print("Standard Deviation:", np.nanstd(data, axis=0))
Original Data Statistics:
Mean: [ 0.19778304 -0.18323331 0.07482166]
Standard Deviation: [1.28030189 0.97419224 1.10733139]
SDs = 2.0
Qs = 0
modified_data = _RemoveOutliers(data, SDs, Qs) # 应用 _RemoveOutliers 函数
print("\nModified Data Statistics:") # 查看修改后的数据的统计信息
print("Mean:", np.nanmean(modified_data, axis=0))
print("Standard Deviation:", np.nanstd(modified_data, axis=0))
Modified Data Statistics:
Mean: [ 0.12895238 -0.17944938 0.05886791]
Standard Deviation: [0.85149081 0.92719664 1.03718499]生成数据中的异常值10已被替换。
4.2.1 函数代码
def _SampleFactorsSimple(S,Svar, K, Qs):
keep = np.setdiff1d(np.arange(0,S.shape[1]), np.arange(0,S.shape[1],(K+Qs)))
Sdraw = S.copy() #复制原始的状态向量矩阵 S，以便对其进行修改
for ii in range(S.shape[0]):
Sdraw[ii,keep] = np.random.multivariate_normal(S[ii,keep].T, np.diag(np.diag(Svar)[keep]), 1)
return Sdraw4.2.2 命令解释
这个函数 _SampleFactorsSimple 用于在给定的卡尔曼平滑结果 S 和状态协方差矩阵 Svar 的基础上，从每个时间点的状态向量中抽取因子。下面是对一些命令的详细解释。
keep = np.setdiff1d(np.arange(0,S.shape[1]), np.arange(0,S.shape[1],(K+Qs)))：计算一个索引数组 keep，其中包含了将被保留的因子的索引。np.arange(0,S.shape[1],(K+Qs)) 生成从0开始，以 (K+Qs) 为步长的索引，表示要保留的因子的位置。np.setdiff1d 用于获取不在这个列表中的索引，即剩余的因子的索引。
Sdraw[ii,keep] = np.random.multivariate_normal(S[ii,keep].T, np.diag(np.diag(Svar)[keep]), 1)：对于每个时间点，从剩余因子的正态分布中抽取新的因子。S[ii,keep].T 是当前时间点保留因子的均值向量。np.diag(np.diag(Svar)[keep]) 是协方差矩阵中仅包含剩余因子的对角元素。np.random.multivariate_normal 用于生成一个多元正态分布的随机样本，以便替换原始矩阵中的对应位置。
return Sdraw：返回更新后的状态向量矩阵 Sdraw。
总体来说，该函数通过从剩余的因子中重新抽取值，对给定的卡尔曼平滑结果进行了简单的因子抽样。这有助于在一定程度上考虑不确定性，并生成具有相似统计特性的新状态向量矩阵。
具体而言，它从多元正态分布中抽取新的状态向量，其中均值和协方差由平滑后的状态向量和协方差矩阵确定。在这个例子中，我们可以考虑一个具体的情景。
import numpy as np
# 假设有5个状态变量，每隔2个变量为一组，且有2组（K=2，Qs=0）
# S的形状为(时间点数, 5)
S = np.array([
[1.0, 2.0, 3.0, 4.0, 5.0],
[2.0, 3.0, 4.0, 5.0, 6.0],
[3.0, 4.0, 5.0, 6.0, 7.0],
])
Svar = np.eye(5) # 假设协方差矩阵Svar是单位矩阵
K = 2 # K表示每隔多少个变量一组，这里为2
Qs = 0 # Qs表示有多少组变量，这里为0
keep = np.setdiff1d(np.arange(0, S.shape[1]), np.arange(0, S.shape[1], (K+Qs)))# 选取需要保留的状态变量的索引，这里保留的是每一组的第一个变量
Sdraw = S.copy() # 对状态向量进行抽样
for ii in range(S.shape[0]):
Sdraw[ii, keep] = np.random.multivariate_normal(S[ii, keep].T, np.diag(np.diag(Svar)[keep]), 1)
# 打印结果
print("原始状态向量：")
print(S)
print("\n抽样后的状态向量：")
print(Sdraw)
原始状态向量：
[[1. 2. 3. 4. 5.]
[2. 3. 4. 5. 6.]
[3. 4. 5. 6. 7.]]
抽样后的状态向量：
[[1. 1.171 3. 3.44 5. ]
[2. 3.747 4. 5.61 6. ]
[3. 3.979 5. 6.117 7. ]]在卡尔曼滤波和平滑算法中，我们通常希望对隐藏状态进行估计，并且希望得到这些估计的不确定性信息。卡尔曼平滑算法可以提供在给定观测序列的情况下，对隐藏状态的后验均值和协方差矩阵的估计。在实际应用中，我们可能对隐藏状态的分布感兴趣，而不仅仅是均值和方差。通过抽样，我们可以得到隐藏状态的不同样本，这有助于我们更全面地理解隐藏状态的不确定性和分布。抽样也可以用于蒙特卡洛方法，例如蒙特卡洛积分，以更好地估计某些函数的期望值。
在上述代码中，对状态向量进行抽样的目的是得到隐藏状态的不同可能值，以更好地理解隐藏状态的变异性和分布。这在分析卡尔曼平滑结果时可能是有用的。
4.2.3 具体效果
具体使用场景在GIBBS抽样中的此处。
S, P = _Kfilter(XY, F, H, Q, R, S)
S, P, Svar = _Ksmoother(F, H, Q, R, S, P, lags, n)
S = _SampleFactorsSimple(S,Svar, K, Qs)4.3.1 函数代码
def _SampleLoadingsLags(S, XY, s0, alpha0, L_var_prior, lags, interv, Qs, K, lagsH):
K = int(S.shape[1] / lags - Qs) # 计算滞后阶数下的状态向量维度
T = XY.shape[0] # 观测序列的长度
N = XY.shape[1] # 观测变量的数量
XYnan = XY[~np.isnan(XY[:, Qs:]).any(axis=1)].copy() # 删除观测数据中包含 NaN 值的行
Beta = [] # 存储回归系数的列表
cumsum = 0
for ii in range(0, len(interv)): # 计算回归系数 Beta
Beta.append(np.linalg.inv(S[0:XYnan.shape[0], Qs+ii:Qs+ii+(lagsH-1)*(Qs+K)+1:Qs+K].T @ S[0:XYnan.shape[0], Qs+ii:Qs+ii+(lagsH-1)*(Qs+K)+1:Qs+K]) @
(S[0:XYnan.shape[0], Qs+ii:Qs+ii+(lagsH-1)*(Qs+K)+1:Qs+K].T @ XYnan[:, Qs+cumsum:Qs+cumsum+interv[ii]]))
cumsum = cumsum + Beta[ii - Qs].shape[1]
Lf = np.zeros((N - Qs, K * lagsH + Qs * lagsH)) # 初始化负载荷矩阵
for jj in range(0, lagsH):
cumsize = 0
for ii in range(len(interv) + Qs):
if ii == 0 or ii % (K + Qs) == 0:
Lf[:, ii:ii+1] = np.zeros((N - Qs, 1))
else:
Lf[cumsize:cumsize + Beta[ii - Qs].shape[1], jj * (K + Qs) + ii] = Beta[ii - Qs][jj, :].reshape(-1, 1).T
cumsize = cumsize + Beta[ii - Qs].shape[1]
Lf_con = np.c_[np.tile(np.c_[1/3, np.zeros((1, K))], (1, 3)), np.tile(np.zeros((1, K + Qs)), (1, lags - 3))]
Lf = np.r_[Lf_con, np.c_[Lf, np.zeros((N - Qs, (lags - lagsH) * (K + Qs)))]]
StoreE = np.empty([XY.shape[0], XY.shape[1]]) # 存储观测残差的数组
R = np.zeros((XY.shape[1], XY.shape[1])) # 初始化观测误差协方差矩阵
for n in range(0, XY.shape[1]): # 对每个观测变量进行处理
ed = XY[:, n] - S[:, :] @ Lf[n, :].T # 计算观测残差
StoreE[:, n] = ed
ed = ed[~np.isnan(ed)] # 去除残差中的 NaN 值
if n == 0: # 对观测误差协方差矩阵元素进行采样
s0 = 0.01
else:
s0 = 5
R_bar = s0 + (ed.T @ ed) + Lf[n, :] @ np.linalg.inv(L_var_prior + np.linalg.inv(S.T @ S)) @ Lf[n, :].T
Rd = np.random.normal(0, 1, T + alpha0)
Rd = Rd.T @ Rd
Rd = np.divide(R_bar, Rd)
R[n, n] = Rd
H = Lf # 负载荷矩阵
return H, R, StoreE4.3.2 命令解释
这段代码的目的是对负载荷矩阵（观测变量与状态变量之间的映射关系）和观测误差协方差矩阵进行采样。这个步骤涉及到贝叶斯统计的概念，主要有以下几个原因。
不确定性建模： 通过采样得到负载荷矩阵和观测误差协方差矩阵的后验分布，可以更全面地了解这些参数的不确定性。在贝叶斯统计中，参数不是固定的值，而是随机变量，其分布描述了对这些参数的信念或不确定性。
先验信息： 通过引入先验信息，例如观测误差协方差矩阵的先验信息，可以在采样过程中更好地结合先验信念和观测数据。这有助于提高对模型参数的估计精度，特别是当观测数据较少或不够准确时。
避免过拟合： 贝叶斯方法的采样过程中引入了随机性，这有助于避免过拟合，即在模型参数估计中考虑过多的噪声或特定样本的影响。通过引入随机性，模型更能适应不同的数据集。
模型诊断： 采样过程中得到的后验分布可以用于进行模型诊断。通过观察后验分布的性质，例如均值、置信区间等，可以发现模型中可能存在的问题，或者确定模型参数的合理范围。
总体而言，这一步是贝叶斯建模中的一个重要环节，通过融合观测数据和先验信息，得到模型参数的后验分布，提供了更为全面和灵活的建模手段。
4.3.3 具体效果
使用在吉布斯抽样中。
H, R, StoreE = _SampleLoadingsLags(S, XY, s0, alpha0, L_var_prior, lags, Ints, Qs, K, lagsH)
if ii >= burn:
Hdraw[:,:,ii-burn] = H.copy()
Qdraw[:,:,ii-burn] = Q.copy()
Fdraw[:,:,ii-burn] = F.copy()
Pdraw[:,:,ii-burn] = P.copy()
Rdraw[:,:,ii-burn] = R.copy()
Sdraw[:,:,ii-burn] = S.copy()
return Hdraw, Qdraw, Fdraw, Pdraw, Rdraw, Sdraw 这些返回值，在后面经过简单处理后，输入卡尔曼滤波函数，将会得到最终输出的预测值。
5.1.1 函数代码
def _SampleLoadingsLagsPhi(S, XY, s0, alpha0, L_var_prior, lags, interv, Qs, K, lagsH, Phi):
K = int(S.shape[1]/lags-Qs) #计算滞后阶数下的状态向量维度。
T = XY.shape[0] #获取观测序列的长度。
N = XY.shape[1] #获取观测变量的数量。
Y_phi = _quasidiff(XY,Phi, Qs) #计算差分操作
Beta = [] #初始化存储回归系数的列表。
cumsum = 0 #初始化累计变量。
cumsize = 0
for ii in range(0,len(interv)):#对每个滞后期的间隔进行循环。
Betatemp = np.zeros((lagsH,interv[ii]))
for jj in range(0,interv[ii]):
notnan = ~np.isnan(Y_phi[:,Qs:]).any(axis=1)
Sreg = _quasidiff(S[notnan,Qs+ii:Qs+ii+(lagsH)*(Qs+K):Qs+K],np.tile(Phi[cumsize+jj,:],(lagsH,1)),0)
notnan[0] = False
Best = np.linalg.inv(Sreg.T @ Sreg) @ (Sreg.T @ Y_phi[notnan,Qs+cumsum+jj:Qs+cumsum+jj+1])
Betatemp[:,jj] = Best.reshape((1,-1))
Beta.append(Betatemp) #对每个观测变量执行回归操作，计算回归系数，并将其存储在Beta列表中。
cumsize = cumsize+interv[ii]
Lf = np.zeros((N-Qs,K*(lagsH+1)+Qs*(lagsH+1))) #初始化负载荷矩阵。
Lfraw = np.zeros((N-Qs,K*(lagsH)+Qs*(lagsH))) #初始化未经处理的负载荷矩阵。
for jj in range(0,lagsH): #对每个时间滞后期进行循环。
cumsize = 0
for ii in range(K+Qs):
if ii==0 or ii+1 % (K+Qs)==0:
Lf[:,ii:ii+1] = np.zeros((N-Qs,1))
else:
Lf[cumsize:cumsize+Beta[ii-Qs].shape[1],jj*(K+Qs)+ii] = Lf[cumsize:cumsize+Beta[ii-Qs].shape[1],jj*(K+Qs)+ii]+Beta[ii-Qs][jj,:].reshape(-1,1).T
Lfraw[cumsize:cumsize+Beta[ii-Qs].shape[1],jj*(K+Qs)+ii] = Beta[ii-Qs][jj,:].reshape(-1,1).T
Lf[cumsize:cumsize+Beta[ii-Qs].shape[1],(jj+1)*(K+Qs)+ii] = (-Phi[cumsize:cumsize+Beta[ii-Qs].shape[1],:].T)*Beta[ii-Qs][jj,:].reshape(1,-1)
cumsize = cumsize+Beta[ii-Qs].shape[1] #对每个观测变量和每个滞后期执行操作，将回归系数添加到负载荷矩阵中。
Lf_con = np.c_[np.tile(np.c_[1/3, np.zeros((1,K))],(1,3)), np.tile(np.zeros((1,K+Qs)),(1,lags-3))]#构建负载荷矩阵的前置部分。
Lf = np.r_[Lf_con, np.c_[Lf, np.zeros((N-Qs,(lags-lagsH-1)*(K+Qs)))]] #将负载荷矩阵的前置部分和计算得到的部分合并。
Lfraw = np.r_[Lf_con, np.c_[Lfraw, np.zeros((N-Qs,(lags-lagsH)*(K+Qs)))]] #将未经处理的负载荷矩阵的前置部分和计算得到的部分合并。
StoreE = np.empty([Y_phi.shape[0],Y_phi.shape[1]]) #初始化存储观测残差的数组。
R = np.zeros((XY.shape[1],XY.shape[1])) #初始化观测误差协方差矩阵。
for n in range(0,Y_phi.shape[1]): #对每个观测变量进行循环。
ed=Y_phi[:,n] - S[:,:] @ Lf[n,:].T #计算观测残差，即观测值减去预测值。
edraw = XY[1:,n]-S[:,:] @ Lfraw[n,:].T #将观测残差存储在StoreE数组中。
StoreE[:,n] = edraw
ed = ed[~np.isnan(ed)]
if n==0: #如果是第一个观测变量，为观测误差协方差矩阵的元素进行采样。
s02=s0/10
else:
s02=s0*5
R_bar = s02 +(ed.T @ ed)+Lf[n,:] @ np.linalg.inv(L_var_prior+np.linalg.inv(S.T @ S)) @ Lf[n,:].T
Rd = np.random.normal(0,1,T+alpha0) #通过采样得到的元素填充观测误差协方差矩阵。
Rd = Rd.T @ Rd
Rd = np.divide(R_bar,Rd)
R[n,n]=Rd
H = Lf #将负载荷矩阵作为输出。
return H, R, Stor5.1.2 命令解释
这个函数 _SampleLoadingsLagsPhi 的目的是通过贝叶斯方法抽样得到因子载荷、测量误差协方差矩阵和测量残差。与上一个函数类似，不过添加了MA项。
Phi: 季度观测数据的 MA 系数。
生成准差分数据 Y_phi：使用 Phi 对原始数据进行准差分。
5.2.1 函数代码
def _SamplePhi(StoreE,R, Qs):
Phi = np.zeros((StoreE.shape[1]-Qs,1))
for ii in range(Qs,StoreE[:,Qs:].shape[1]):
eps = StoreE[~np.isnan(StoreE[:,ii]),ii].reshape(-1,1)
Eps, EpsL = _MakeLags(eps, 1)
Phi[ii-1,0] = np.linalg.inv(EpsL.T @ EpsL) @ (EpsL.T @ Eps)
return Phi 5.2.2 命令解释
这个函数 _SamplePhi 用于通过贝叶斯方法抽样得到季度数据的 MA 模型的系数 Phi，将被用于带MA项的吉布斯抽样中，以下是每个命令的详细解释。
初始化 Phi：Phi 是一个存储 MA 模型系数的向量。
遍历残差列：对于每个残差列（除去季度可观测 GDP 系列），通过建立自回归模型进行贝叶斯抽样。
计算 MA 参数 Phi：对于每个残差列，取出非缺失值的残差，并将其构建为一个列向量。使用 _MakeLags 函数，将残差序列构建成具有一阶滞后的矩阵 EpsL。通过线性回归计算 MA 参数 Phi。
返回结果：返回抽样得到的 MA 参数 Phi。
总体来说，该函数执行了季度数据的 MA 模型参数的贝叶斯抽样。
5.3.1 函数命令
def _quasidiff(Y,Phi,Skip):
Yphi = Y[1:,Skip:]-Y[0:-1,Skip:] * Phi.T
if Skip>0:
Yphi = np.c_[Y[1:,0:Skip],Yphi]
return Yphi5.3.2 命令解释
这个函数 _quasidiff 是一个实现准差分的工具函数，用于在时间序列上应用季节性调整。仅在结合MA方法时，才能使用。以下是每个命令的详细解释：
输入参数：Y：输入的时间序列数据。Phi：MA 模型的系数。Skip：需要跳过的列数。
应用准差分：函数首先从第二行开始计算 Y 的准差分。对于每一列（时间序列），它减去上一行的值乘以对应的 Phi 系数。如果存在需要跳过的列数（Skip > 0），则在准差分结果的左侧添加相应的列（原始数据中需要跳过的列）。
返回结果：返回应用准差分后的时间序列。
总体来说，该函数实现了一个简单的准差分操作，通过减去前一期的值乘以 MA 模型的系数来调整时间序列。
5.3.3 具体效果
def _SamplePhi(StoreE,R, Qs):
Phi = np.zeros((StoreE.shape[1]-Qs,1))
for ii in range(Qs,StoreE[:,Qs:].shape[1]):
eps = StoreE[~np.isnan(StoreE[:,ii]),ii].reshape(-1,1)
Eps, EpsL = _MakeLags(eps, 1)
Phi[ii-1,0] = np.linalg.inv(EpsL.T @ EpsL) @ (EpsL.T @ Eps)
return Phi
Phi = _SamplePhi(StoreE,R, Qs)
Phidraw[:,:,ii-burn] = Phi.copy()
if self.MAterm==1:
Phidraw = self.Phidraw
Phifor = np.median(Phidraw[:,:,:],axis=2)
XYf = _quasidiff(XYf,Phifor,Qs)6.1.1 函数代码
def initializevals(GDP, Monthly, lags, lagsH, K, Qs):
GDPnan = np.asarray(np.kron(GDP, [[np.nan], [np.nan], [1]])) #将季度 GDP 数据扩展为包含 NaN 值的矩阵，使其具有与月度数据相同的形状
GDP_fill = _InterpolateNaN(GDPnan) #通过 _InterpolateNaN 函数填充 NaN 值，线性插值
cutoff = len(GDP_fill) #取处理后的数据的长度
MonthlyDat_norm = [] #生成空表，储存标准化后的数据表
for df in Monthly: #遍历每个列表并进行标准化
MonthlyDat_norm.append(_NormDframe(df[0:cutoff])) #
loadings = [None]*len(Monthly) #生成空表，用来储存因子载荷
Factors = [None]*len(Monthly) #生成空表，用来储存主成分因子
for df,num in zip(MonthlyDat_norm, range(len(MonthlyDat_norm))): #遍历标准化后的月度数据列表，进行主成分分析
Factors[num], loadings[num] = _PCA_fact_load(df[0:cutoff].fillna(method='ffill')) #储存因子载荷和因子
F0 = np.concatenate(Factors,1) #将主成分因子进行拼接，得到一个矩阵
N = len(np.concatenate(loadings,0)) #获取所有加载矩阵垂直拼接后的总行数
Lf = np.zeros((N,len(loadings))) #创建一个零矩阵，用于存储载荷矩阵，这个矩阵的每一列对应一个月度数据的载荷矩阵
cumsize = 0 #初始化一个变量 cumsize 用于累积每个月度因子的因子载荷数量
for ii in range(len(loadings)):
Lf[cumsize:cumsize+len(loadings[ii]),ii] = loadings[ii].reshape(-1,1).T
cumsize = cumsize+len(loadings[ii]) #Lf 处理
F0 = np.hstack((GDP_fill, F0)) #将主成分因子与处理后的 GDP 组合
F0_cut, F0lag = _MakeLags(F0, lags) #滞后处理
Lf_con = np.c_[np.tile(np.c_[1/3, np.zeros((1,K))],(1,3)), np.tile(np.zeros((1,K+Qs)),(1,lags-3))]
H = np.r_[Lf_con, np.c_[np.zeros((N,1)), Lf, np.zeros((N,(lags-1)*(K+1)))]]
Monthly_dat = np.concatenate(MonthlyDat_norm,axis=1) #将标准化后的月度数据进行拼接
e = np.c_[GDP_fill[lags:,0:], Monthly_dat[lags:,0:]] - F0lag[:,] @ H.T #计算残差矩阵，通过矩阵运算计算模型的预测值，见示例四，从实际值中减去预测值得到残差
e = e[~np.isnan(e).any(axis=1)] #移除包含 NaN 值的行，确保残差矩阵 e 不包含任何缺失值
R=np.divide((e.T @ e),(Monthly_dat.shape[0]-lags)) #e.T @ e是残差的平方和，除以自由度（Monthly_dat.shape[0]-lags计算均值
R = np.diag(np.diag(R)) #将协方差矩阵的非对角元素置零，得到一个对角协方差矩阵
B, Q = _estVAR(F0_cut,lags) #通过时间序列数据 F0_cut 估计向量自回归模型的系数和协方差
Q=np.c_[Q, np.zeros((K+1,(K+1)*(lags-1)))] #调整Q的维度
Q = np.r_[Q, np.zeros(((K+1)*(lags-1),(K+1)*(lags)))] #
F=np.r_[B[0:-1,0:].T, np.c_[np.eye((K+Qs)*(lags-1)), np.zeros(((K+Qs)*(lags-1),(K+Qs)))]]
XY = np.c_[GDPnan, Monthly_dat] #未线性插值的 GDP 与月度数据拼接
S, P = _Kfilter(XY, F, H, Q, R, F0lag) #通过 Kalman 滤波法，返回状态的估计和协方差。
ints = [listel.shape[0] for listel in loadings] #loadings 包含因子载荷矩阵的列表，ints 中的每个元素就是对应因子载荷矩阵的行数。
return XY, H, R, F, Q, S, ints6.1.2 命令解释
GDP读取文件输入原始季度GDP数据，Monthly读取文件输入原始月度数据，示例输入了三个表，每个表两列数据，共六列，K月度因子数量，此处为表的数量，注意不是数据的数量，Qs季度因子数量，此处为一，lags转移矩阵的滞后阶数，lagsH可观测滞后阶数。
这个 initializevals 函数的目的是初始化模型的一些参数，以便进行后续的贝叶斯VAR估计。
GDP填充将用于残差计算，GDP扩展将用于状态估计。展示数据处理过程。
构建输出参数：返回初始化后的数据 XY，载荷矩阵 H，观测残差协方差矩阵 R，通过时间序列数据估计向量自回归模型的系数和协方差F，Q，拼接未线性插值的 GDP 与月度数据，进行 Kalman 滤波，返回状态的估计和协方差S，P。ints 是一个列表，获取每个因子载荷矩阵的行数。
6.2.1 函数代码
def Gibbs_loop(XY,F, H, Q, R, S, lags, lagsH, K, Qs, s0, alpha0, L_var_prior, Ints, burn, save, GDPnorm):
Hdraw = np.empty([H.shape[0],Q.shape[1],save])
Qdraw = np.empty([Q.shape[0],Q.shape[1],save])
Fdraw = np.empty([F.shape[0],F.shape[1],save])
Pdraw = np.empty([S.shape[0],S.shape[1]*S.shape[1],save])
Rdraw = np.empty([R.shape[0],R.shape[1],save])
Sdraw = np.empty([S.shape[0],S.shape[1],save])
iter = burn+save
n = K+Qs
lambda1 = 0.5
lambda3 = 1
lambda4 = 1000
n = (K+Qs)
alpha0=1
S_Y, S_X = _MakeLags(S[:,0:K+Qs], lags)
phi0, alpha0, S0, beta0 = _setVARpriors(S_X,S_Y,lambda1, lambda3, lambda4, lags, alpha0)
keep = np.setdiff1d(np.arange(0,S.shape[1]), np.arange(0,S.shape[1],n))
for ii in range(iter):
if ii % 10==0:
print('Completed ', str(ii/iter*100), "% of the loop")
S, P = _Kfilter(XY, F, H, Q, R, S)
S, P, Svar = _Ksmoother(F, H, Q, R, S, P, lags, n)
S = _SampleFactorsSimple(S,Svar, K, Qs)
if GDPnorm:
S[:,:] = (S[:,:]-np.mean(S[:,:], axis=0))/np.std(S[:,:],axis=0)
else:
S[:,keep] = (S[:,keep]-np.mean(S[:,keep], axis=0))/np.std(S[:,keep],axis=0)
draws = 1
sigma_gibbs, beta_gibbs = _EstBVAR(S[:,0:K+Qs],lags, phi0, alpha0, S0, beta0, draws)
sigma = sigma_gibbs.reshape((sigma_gibbs.shape[1],sigma_gibbs.shape[2]))
beta = beta_gibbs.reshape((beta_gibbs.shape[1],beta_gibbs.shape[2]))
Q = np.c_[sigma, np.zeros((K+1,(K+1)*(lags-1)))]
Q = np.r_[Q, np.zeros(((K+1)*(lags-1),(K+1)*(lags)))]
F=np.r_[beta[0:,0:].T, np.c_[np.eye((K+Qs)*(lags-1)), np.zeros(((K+Qs)*(lags-1),(K+Qs)))]]
H, R, StoreE = _SampleLoadingsLags(S, XY, s0, alpha0, L_var_prior, lags, Ints, Qs, K, lagsH)
if ii >= burn:
Hdraw[:,:,ii-burn] = H.copy()
Qdraw[:,:,ii-burn] = Q.copy()
Fdraw[:,:,ii-burn] = F.copy()
Pdraw[:,:,ii-burn] = P.copy()
Rdraw[:,:,ii-burn] = R.copy()
Sdraw[:,:,ii-burn] = S.copy()
return Hdraw, Qdraw, Fdraw, Pdraw, Rdraw, Sdraw 6.2.2 命令解释
这是一个 Gibbs 采样的循环函数，用于估计模型的参数。以下是每个命令的详细解释：
输入参数：XY：包含月度数据的矩阵。F：状态转移矩阵的初始值。H：加载矩阵的初始值。Q：状态转移方差项的初始值。R：加载方差项的初始值。S：状态矩阵的初始值。lags：状态转移矩阵的滞后阶数。lagsH：可观测数据的滞后阶数。K：月度因子数量。Qs：季度因子数量。s0：因子加载方差的先验。alpha0：因子加载方差的缩放参数。L_var_prior：因子加载系数的先验权重。Ints：因子加载的限制 - 每个因子上加载的数量。burn：燃烧期，抽样过程中的迭代次数。save：保存期，保存的迭代次数。GDPnorm：是否对 GDP 进行标准化。
初始化存储结果的矩阵：Hdraw：存储加载矩阵的采样结果。Qdraw：存储状态转移方差项的采样结果。Fdraw：存储状态转移矩阵的采样结果。Pdraw：存储状态的平滑方差项的采样结果。Rdraw：存储加载方差项的采样结果。Sdraw：存储状态矩阵的采样结果。
Gibbs 采样循环：对每次迭代执行以下步骤：Kalman 滤波和平滑，计算状态矩阵和平滑方差项。通过随机冲击重新采样状态矩阵。标准化状态矩阵（除了 GDP）。随机生成 BVAR 模型，得到状态转移方差项 Q 和系数矩阵 F。为因子加载进行采样。将采样结果存储起来。
返回结果：返回加载矩阵、状态转移方差项、状态转移矩阵、平滑方差项、加载方差项和状态矩阵的采样结果。
吉布斯抽样（Gibbs Sampling）是一种马尔可夫链蒙特卡洛（MCMC）的方法，用于从多维分布中抽取样本。它属于马尔可夫链的一种特殊类型，通过在给定条件下对每个变量进行逐个抽样，从而形成一个完整的马尔可夫链。吉布斯抽样的过程如下。
初始化参数： 随机初始化模型参数。
逐个变量抽样： 对于每个变量，依据给定条件分布进行抽样。在吉布斯抽样中，每个变量的抽样都是在其他所有变量已知的条件下进行的。
更新参数： 将抽样得到的新值更新到模型参数中。
重复步骤 2 和 3： 重复逐个变量抽样和更新参数的过程，直到收敛为止。
Gibbs_loop 函数就是吉布斯抽样的一个实现。在每次循环中，该函数对隐藏状态进行 Kalman 滤波和平滑，然后采样因子和更新模型参数（包括负载荷矩阵、观测误差协方差矩阵等）。这一过程通过马尔可夫链的方式迭代进行，最终得到的样本就是从联合分布中抽取的。吉布斯抽样的优势之一是它对高维参数空间的建模很有效，尤其在无法直接从联合分布中抽样的情况下。通过在给定条件下逐个变量进行抽样，吉布斯抽样能够在高维空间中灵活地移动，并提供对参数的后验分布的估计。
6.3.1 函数代码
def Gibbs_loopMA(XY,F, H, Q, R, S, lags, lagsH, K, Qs, s0, alpha0, L_var_prior, Ints, burn, save, GDPnorm):
Hdraw = np.empty([H.shape[0],Q.shape[1],save])
Qdraw = np.empty([Q.shape[0],Q.shape[1],save])
Fdraw = np.empty([F.shape[0],F.shape[1],save])
Pdraw = np.empty([S.shape[0]-1,S.shape[1]*S.shape[1],save])
Rdraw = np.empty([R.shape[0],R.shape[1],save])
Sdraw = np.empty([S.shape[0]-1,S.shape[1],save])
Phidraw = np.empty([XY.shape[1]-Qs,1,save])
iter = burn+save
n = K+Qs
lambda1 = 0.5
lambda3 = 1
lambda4 = 1000
n = (K+Qs)
alpha0=1
S_Y, S_X = _MakeLags(S[:,0:K+Qs], lags)
phi0, alpha0, S0, beta0 = _setVARpriors(S_X,S_Y,lambda1, lambda3, lambda4, lags, alpha0)
Phi = np.zeros((XY.shape[1]-Qs,1))
keep = np.setdiff1d(np.arange(0,S.shape[1]), np.arange(0,S.shape[1],n))
for ii in range(iter):
if ii % 10==0:
print('Completed ', str(ii/iter*100), "% of the loop")
XYquas = _quasidiff(XY,Phi, Qs)
S, P = _Kfilter(XYquas, F, H, Q, R, S)
S, P, Svar = _Ksmoother(F, H, Q, R, S, P, lags, n)
S = _SampleFactorsSimple(S,Svar, K, Qs)
if GDPnorm:
S[:,:] = (S[:,:]-np.mean(S[:,:], axis=0))/np.std(S[:,:],axis=0)
else:
S[:,keep] = (S[:,keep]-np.mean(S[:,keep], axis=0))/np.std(S[:,keep],axis=0)
draws = 1
sigma_gibbs, beta_gibbs = _EstBVAR(S[:,0:K+Qs],lags, phi0, alpha0, S0, beta0, draws)
sigma = sigma_gibbs.reshape((sigma_gibbs.shape[1],sigma_gibbs.shape[2]))
beta = beta_gibbs.reshape((beta_gibbs.shape[1],beta_gibbs.shape[2]))
Q = np.c_[sigma, np.zeros((K+1,(K+1)*(lags-1)))]
Q = np.r_[Q, np.zeros(((K+1)*(lags-1),(K+1)*(lags)))]
F=np.r_[beta[0:,0:].T, np.c_[np.eye((K+Qs)*(lags-1)), np.zeros(((K+Qs)*(lags-1),(K+Qs)))]]
H, R, StoreE = _SampleLoadingsLagsPhi(S, XY, s0, alpha0, L_var_prior, lags, Ints, Qs, K, lagsH, Phi)
Phi = _SamplePhi(StoreE,R, Qs)
if ii >= burn:
Hdraw[:,:,ii-burn] = H.copy()
Qdraw[:,:,ii-burn] = Q.copy()
Fdraw[:,:,ii-burn] = F.copy()
Pdraw[:,:,ii-burn] = P.copy()
Rdraw[:,:,ii-burn] = R.copy()
Sdraw[:,:,ii-burn] = S.copy()
Phidraw[:,:,ii-burn] = Phi.copy()
return Hdraw, Qdraw, Fdraw, Pdraw, Rdraw, Sdraw, Phidraw 6.3.2 命令解释
这是一个与之前相似的 Gibbs 采样循环函数，用于估计包含季度可观测数据和MA（移动平均）项的贝叶斯VAR（向量自回归）模型的参数。Phidraw：存储MA系数的采样结果，具体步骤与前一函数相同，仅多了一步对MA系数进行采样。
7.1.1 函数代码
class DynamicFactorModel():
def __init__(self, Q_GDP, Monthly, Dates, K, Qs, lags, lagsH, MAterm, normGDP):
self.GDP = Q_GDP #读取文件输入原始季度GDP数据
self.Monthly = Monthly #读取文件输入原始月度数据，示例输入了三个表，每个表两列数据，共六列
self.Dates = Dates #读取文件输入原始季度GDP的日期数据
self.K = K #月度因子数量，此处为表的数量，注意不是数据的数量
self.Qs = Qs #季度因子数量，此处为一
self.lags = lags #转移矩阵的滞后阶数，如何设置合适？
self.lagsH = lagsH #可观测滞后阶数，如何设置合适？
self.MAterm = MAterm #是否使用移动平均项，0表示不使用，1表示使用，使用后预测将会变得平滑，但准确度会下降
self.normGDP = normGDP #是否对GDP进行标准化，0表示不使用，1表示使用，差别不大
def estimateGibbs(self, burn, save, s0 = 0.1, alpha0 = 1, L_var_prior = None):
#burn一次性初始抽取的次数，如何设置合适？
#save吉布斯抽签储存，如何设置合适？
self.s0 = s0 #s0: 因子加载方差的先验，如何设置合适？
self.alpha0 = alpha0 #alpha0: 因子加载方差的缩放参数，如何设置合适？
self.L_var_prior = L_var_prior #L_var_prior: 因子加载系数的先验权重，默认为单位矩阵。
if L_var_prior is None:
L_var_prior = np.identity((self.K+self.Qs)*self.lags) #生成一个月数加季数*转移滞后的单位矩阵
self.XY, H, R, F, Q, S, ints = initializevals(self.GDP, self.Monthly, self.lags, self.lagsH, self.K, self.Qs) #将原始数据输入初始化函数，具体见本系列文（6）
if self.normGDP: #判断是否对GDP进行标准化
self.GDPmean = np.nanmean(self.GDP) #
self.GDPstd = np.nanstd(self.GDP) #
self.XY[0:,0] = (self.XY[0:,0] - self.GDPmean)/self.GDPstd #
XYcomp = self.XY[0:self.GDP.shape[0]*3,:].copy() #创建一个季度数据的三倍（月度）的表，切片处理后的年度数据
if self.MAterm == 0: #选择是否加入MA项，进而调用不同的基布斯抽样
self.Hdraw, self.Qdraw, self.Fdraw, self.Pdraw, self.Rdraw, self.Sdraw = Gibbs_loop(self.XYcomp,F, H, Q, R, S, self.lags, self.lagsH, self.K \
, self.Qs, s0, alpha0, L_var_prior, ints, burn, save, self.normGDP) #
return self.Hdraw, self.Qdraw, self.Fdraw, self.Pdraw, self.Rdraw, self.Sdraw #
elif self.MAterm == 1: #
self.Hdraw, self.Qdraw, self.Fdraw, self.Pdraw, self.Rdraw, self.Sdraw, self.Phidraw = Gibbs_loopMA(self.XYcomp,F, H, Q, R, S, self.lags, self.lagsH, self.K \
, self.Qs, s0, alpha0, L_var_prior, ints, burn, save, self.normGDP) #注意缩进格式
def Nowcast(self, start, Horz):
#start: 开始日期
#Horz: 预测的时间跨度（季度）
XY = self.XY #处理后的年度数据
Dates = self.Dates #原始季度日期数据
Hdraw = self.Hdraw #基布斯抽样的返回（下同），分别是？？
Qdraw = self.Qdraw #
Fdraw = self.Fdraw #
Rdraw = self.Rdraw #
Sdraw = self.Sdraw #
Qs = self.Qs #季度因子
K = self.K #月度因子
Hfor = np.median(Hdraw[:,:,:],axis=2) #从三维数组中取一个二维的中位数矩阵（下同）
Qfor = np.median(Qdraw[:,:,:],axis=2) #
Ffor = np.median(Fdraw[:,:,:],axis=2) #
Rfor = np.median(Rdraw[:,:,:],axis=2) #
Sfor = np.median(Sdraw[:,:,:],axis=2) #
Total = Dates.shape[0]*3+Horz*3 #生成一个原始季度数据的三倍加预测季度三倍大小的总矩阵
if XY.shape[0]<Total: #用空数据将原始数据扩展以容纳预测数据
NewNans = np.empty((Total-XY.shape[0],XY.shape[1]))#
NewNans[:] = np.nan #
XYf = np.r_[XY, NewNans] #
if self.MAterm==1: #添加移动平均项
Phidraw = self.Phidraw #有MA的基布斯多返回的一个值
Phifor = np.median(Phidraw[:,:,:],axis=2) #
XYf = _quasidiff(XYf,Phifor,Qs) #使用之前编写的差分函数进行差分处理
Findex = int(np.where(Dates==start)[0]) #找到指定的开始预测的日期的位置
Quartersback = (Dates.shape[0]-Findex) #得到从开始预测日期到最后有多少季度
self.Quartersback = Quartersback
Fcast_current = np.empty((Quartersback+1, 3)) #创建四个空数组，形状为最终输出的季度数加一再乘三
Fcast_next = np.empty((Quartersback+1, 3)) #
Outturn_current = np.empty((Quartersback+1,3)) #
Outturn_next = np.empty((Quartersback+1,3)) #
DateFcast = np.arange(start,start+Quartersback/4,0.25) #输出预测日期，注意缩进格式
for ii in range(Quartersback+1):
if ii ==Quartersback: #检查当前的迭代是否达到了最后一个季度。如果是的话，就会为变量 newvar 赋值为 1
newvar = 1 #
PriorS = Sfor[-(Quartersback)*3-1+ii*3,:].reshape(1,-1) #从处理后的数据中摘取一部分？？
XYsnap = XYf[-(Quartersback+Horz)*3+ii*3:-(Quartersback+Horz)*3+(ii*3+3),:] #
for jj in range(3):
Outturn_current[ii,jj] = XYsnap[2,0] #从 XYsnap 数组中提取的第三行和第一列的元素值
if ii<Quartersback-1:
Outturn_next[ii,jj] = XYf[-(Quartersback+Horz)*3+(ii*3+5):-(Quartersback+Horz)*3+(ii*3+6),0]
Feed = XYsnap.copy()
if jj<2:
Feed[-2+jj:,:] = np.nan
Feed[2,0] = np.nan
S_f, P_f = _Kfilter(Feed, Ffor, Hfor, Qfor, Rfor, PriorS) #通过卡尔曼滤波法得到预测值
Fcast_current[ii,jj] = np.mean(S_f[0:3,0]) #从状态矩阵S_f中取出第1至3行、第1列的元素，即取出状态矩阵的前三个元素，计算平均值
FeedNext = np.r_[Feed,np.nan*np.empty((3,Feed.shape[1]))] #FeedNext 是在 Feed 的末尾添加三行 NaN 值
S_fn, P_fn = _Kfilter(FeedNext, Ffor, Hfor, Qfor, Rfor, PriorS) #通过卡尔曼滤波法得到预测值
Fcast_next[ii,jj] = np.mean(S_fn[3:,0]) #从状态矩阵中取出第1至3行、第1列的元素，即取出状态矩阵的前三个元素，计算平均值
if self.normGDP: #将输出值进行标准化处理
Fcast_current = Fcast_current*self.GDPstd+self.GDPmean #
Fcast_next = Fcast_next*self.GDPstd+self.GDPmean #
Outturn_current = Outturn_current*self.GDPstd+self.GDPmean #
Outturn_next = Outturn_next*self.GDPstd+self.GDPmean #
RMSE = np.empty((3,2)) #创建一个用于计算误差的空数组
for ii in range(3): #
CurrentErr = np.square(Fcast_current[:,ii] - Outturn_current[:,ii]) #输出当前预测的误差
NextErr = np.square(Fcast_next[:,ii] - Outturn_next[:,ii]) #输出下一季预测的误差
RMSE[ii,:] = np.c_[np.sqrt(np.nanmean(CurrentErr)), np.sqrt(np.nanmean(NextErr))] #第一列是当前季度的均方根误差，第二列是下一季度的均方根误差
Datesaug = np.r_[Dates, Dates[-1:]+0.25] #为日期输出增加一个季度
self.S_f = S_f
self.P_f = P_f
self.S_fn = S_fn
self.S_fn = S_fn
self.Fcast_current = Fcast_current
self.Fcast_next = Fcast_next
self.Outturn_current = Outturn_current
self.Outturn_next = Outturn_next
self.RMSE = RMSE
self.Datesaug = Datesaug
def PlotFcast(self, Month):
titlelist = ['Month 1', 'Month 2', 'Month 3']
fig1 = plt.figure(figsize=(15,15)) #创建一个图形对象，设置图形大小为 (15,15)
ax1 = fig1.add_subplot(111) #在一个 1x1 的网格中创建第一个（唯一的）子图
ax1.plot(self.Datesaug[-self.Quartersback-1:], self.Fcast_current[:,Month-1],marker='o', color='olive', linewidth=2) #选择数组中第 Month 列的所有行
ax1.plot(self.Datesaug[-self.Quartersback-1:], self.Outturn_current[:,Month-1], marker='o', color='blue', linewidth=2)
ax1.legend(['Nowcast', 'Actual']) #添加图例，显示 'Nowcast' 表示预测值，'Actual' 表示实际值
ax1.set_title(titlelist[Month-1]) #设置子图的标题，标题内容从 titlelist 中选择，根据传入的 Month 参数来确定
plt.show()7.1.2 命令解释
这是一个用于实现动态因子模型的 Python 类 DynamicFactorModel。以下是该类的主要方法和属性：
__init__ 方法： 初始化类的实例。接受季度GDP数据（Q_GDP）、月度数据（Monthly）、日期数据（Dates）、月度因子数量（K）、季度因子数量（Qs）、转移矩阵滞后阶数（lags）、可观测数据滞后阶数（lagsH）、是否使用移动平均项（MAterm）、是否对GDP进行标准化（normGDP）等参数。
estimateGibbs 方法： 使用Gibbs采样法估计模型参数。接受燃烧期（burn）、保存期（save）、因子加载方差的先验（s0）、因子加载方差的缩放参数（alpha0）、因子加载系数的先验权重（L_var_prior）等参数。
lags（转移矩阵的滞后阶数）： 这个参数指定了动态因子模型中状态转移矩阵的滞后阶数。滞后阶数决定了过去时刻的状态对当前时刻的状态的影响程度。选择适当的滞后阶数通常涉及平衡模型的复杂性和拟合数据的能力。较高的滞后阶数可能更灵活，但也可能导致过拟合。通常，可以通过交叉验证或信息准则（例如，AIC，BIC）来选择一个合适的滞后阶数。
lagsH（可观测滞后阶数）： 这个参数指定了可观测部分（例如，月度数据）的滞后阶数。类似于状态转移矩阵的滞后阶数，可观测滞后阶数决定了过去时刻的观测数据对当前时刻的观测的影响。选择适当的滞后阶数同样可以通过交叉验证或信息准则来进行，以确保在模型中捕获重要的时间动态特征。
burn（燃烧期）： 这个参数决定了在采样过程中丢弃的初始迭代次数（燃烧期）。在燃烧期内，采样器可能还没有收敛到目标分布。通常的做法是将burn设置为一个合理的值（例如，50），以便在收集样本之前使采样器稳定下来。
save（保存期）： 这个参数定义了在燃烧期之后保存的迭代次数。保存的样本然后用于推断。较大的save值将提供更多用于分析的样本，但也会增加计算时间。通常有益的做法是检查跟踪图或其他诊断图，以确定马尔可夫链是否已经收敛，然后相应地设置save。
s0： 因子加载方差的先验。选择s0取决于您对因子加载变异性的先验信念。较小的s0表示较小的变异性，而较大的s0则允许更多的变异性。通过尝试不同的s0值并观察其对结果的影响，可以进行敏感性分析。
alpha0： 因子加载方差的缩放参数。与s0类似，选择alpha0影响对因子加载方差的先验强度。较小的alpha0值表示较小的方差，而较大的值表示更大的方差。同样，敏感性分析可以帮助选择适当的值。
L_var_prior： 这是因子加载的先验权重矩阵。默认的是单位矩阵。该矩阵影响因子加载的先验结构。您可能需要根据先验知识或模型需求调整这个矩阵。
返回加载矩阵（Hdraw）、状态转移方差项（Qdraw）、状态转移矩阵（Fdraw）、平滑方差项（Pdraw）、加载方差项（Rdraw）、状态矩阵（Sdraw）等采样结果。
Nowcast 方法：生成模型的实时预测结果。接受开始日期（start）和预测的时间跨度（Horz）等参数。返回当前月和下月的预测值（Fcast_current 和 Fcast_next）、当月和下月的真实值（Outturn_current 和 Outturn_next）、RMSE（RMSE）等信息。
PlotFcast 方法：根据输入的月份，绘制当月和下月的实时预测值与真实值的对比图。接受月份参数（Month）。
这个类封装了动态因子模型的估计和预测过程，提供了一种方便的方式来执行这些任务。
向量自回归（VAR）： 这是一个时间序列建模方法，用于描述多个变量之间的动态关系。在VAR中，变量的当前值被其过去的值线性回归。你可以使用VAR模型来捕捉观测数据的动态特征，例如，将经济指标、收入等作为变量。卡尔曼滤波： 卡尔曼滤波是一种递归估计技术，用于估计系统的状态变量，尤其是当系统的状态是不完全观测时。在时间序列预测中，你可以使用卡尔曼滤波来动态地估计系统状态，以便更准确地进行预测。吉布斯抽样： 吉布斯抽样是一种贝叶斯统计的方法，用于从参数的后验分布中抽取样本。在这个框架中，你可以使用吉布斯抽样来对模型的未知参数进行贝叶斯推断，获取参数的后验分布。
综合这三个方法，其中VAR模型用于建模动态关系，卡尔曼滤波用于动态状态估计，吉布斯抽样用于更新参数。
7.2.1 函数代码
梳理一下整个模型的预测执行过程，先载入相关包，其中DFM是代码集合。
import pandas as pd
import os
import matplotlib.pyplot as plt
import numpy as np
import DFM as DF读取代码文件目录中的数据集，只需要将数据文件放在代码同一个文件夹里就好。
dirname = os.path.dirname(__file__) #读取文件路径
excel_file = os.path.join(dirname, r'ImportFile.xlsx') #读取数据文件
GDP_dat = pd.read_excel(excel_file, sheet_name='GDP') #读取数据文件中的每一个表
Emp_dat = pd.read_excel(excel_file, sheet_name='Emp') #
Cons_dat = pd.read_excel(excel_file, sheet_name='Cons') #
IPExp_dat = pd.read_excel(excel_file, sheet_name='IPExp') #
GDP = GDP_dat['GDP'][0:].values.reshape(-1,1) #读取除日期外的GDP数据
Dates = GDP_dat.Date #读取季度日期
MonthlyDat = [Emp_dat.drop(columns=['Date']), #储存月度数据
Cons_dat.drop(columns='Date'), IPExp_dat.drop(columns='Date')]
MAterm = 0 #是否使用移动平均项，0表示不使用，1表示使用
lags = 6 #转移矩阵的滞后阶数
lagsH = 1 #可观测滞后阶数
K = 3 #月度因子数量，与月度表格数量一致
Qs = 1 #季度因子数量
normGDP = 1 #是否对GDP进行标准化
burn = 50 #一次性初始抽取的次数
save = 50 #吉布斯抽签储存首先训练模型将上述值输入，这些值将会作为程序内的一些初始值。
DynamicFac = DF.DynamicFactorModel(GDP, MonthlyDat, Dates, K, Qs, lags, lagsH, MAterm, normGDP)
#输入基布斯抽样需要的初始值:
DynamicFac.estimateGibbs(burn, save)
#输入预测的起始年份和向后预测的季度数:
DynamicFac.Nowcast(2008, 2)
print('Plot quasi-out of sample month 2 nowcasts against data (these use full estimation sample to estimate parameters and previous quarter state)')
#根据数据绘制准样本外第2个月的nowcasts(这些使用完整的估计样本来估计参数和上一季度状态)
DynamicFac.PlotFcast(3)
np.set_printoptions(precision=3) #浮点数精确度
print('View RMSEs, as current quarter, next quarter (columns), and months (1,2,3) of the quarter (rows).')
#查看RMSEs，如当前季度、下一季度(列)和该季度的月份(1、2、3)(行)。
print(DynamicFac.RMSE)
print('Fcast_current')
#当前季度的系列即时预测(季度行，列，每个季度的月份)
print(DynamicFac.Fcast_current)
print('Fcast_next')
#下一季度系列即时预测(季度行，列，每个季度的月份)
print(DynamicFac.Fcast_next)
print('Outturn_current')
print(DynamicFac.Outturn_current)
print('Dates series')
print(DynamicFac.Datesaug)
