Kernel Density Estimation
引言
在讨论核密度估计前，我们先回顾一下直方图。
直方图
直方图是一种估计概率密度函数的简单方法。直方图通过将数据划分为若干个离散的区间来统计频率。假设样本数据\(X_1,...,X_n\)独立同分布(i.i.d),我们可以对区间\([x_0,x_0+h]\)中样本计数，并通过每个区间频率来估计\(x\)在\([x_0,x_0+h)\)上的密度\(f(x)\)如下： \[\begin{split}
f(x_0)=&F'(x_0)\\
=&\lim_{h\rightarrow 0^+}\frac{F(x_0+h)-F(x_0)}{h}\\
=&\lim_{h\rightarrow 0^+}\frac{\mathbb{P}[x_0<X<x_0+h]}{h}\\
\end{split}\] 具体来说，给定初始点\(t_0\)和带宽\(h>0\)，直方图通过计算每个区间内的采样点的数量，在区间\({B_k:=[t_k,t_{k+1}):t_k=t_0+hk,k\in\mathbb{Z}}\)内构建一个分段常数函数。这些分割区间长度固定，计数因此已经标准化，以便在固定间隔中计算样本的相对频率。点x处的直方图定义为：
\[\hat{f}_H(x;t_0,h):=\frac{1}{nh}\sum_{i=1}^n
\mathbb{1}_{\{X_i\in B_i:x_i\in B_k\}}\]
由于位于同一组内所有点的直方图密度估计均相等，直方图所对应的分布函数\(\hat{F}_h(x)\)是单调递增的阶梯函数，这与经验分布函数形状类似。
假设样本\(X_1,...,X_n\)在\(B_k\)区间的观测值个数为\(v_k\),则上述直方图也可以表示为： \[
\hat{f}_H(x;t_0,h)=
\begin{cases}
\frac{v_k}{nh},& if\ x\in B_k\ for\ a\ certain\ k\in \mathbb{Z}\\
0,& else
\end{cases}
\]
直方图的高度表示给定数据点位于该区间内的概率贡献。通常，如果区间过小或过多，则仅需要少量的数据点即可产生”噪声”，反之，如果区间过大，则可能会损失分布的显著特征。因此，直方图估计通常不能很好地捕捉数据分布的细节，因为它们的形状依赖于选定的区间宽度和位置。
KDE根据已知数据集计算一个连续的概率密度函数，反映随机变量的概率分布。在KDE中，每个数据点都被认为服从一个核函数分布，这些分布的和组成了估计概率密度函数。KDE提供了一个连续的、光滑的曲线来描述数据的密度分布，它能更好地反映数据集中的细节和规律。
因此，核密度估计通常比直方图估计更加精确和灵活。
KDE
从直方图到核密度估计
上节介绍的方法求密度时需要依赖于某一数据点\(x_0\)的选取。不妨将上述方法稍做改进，将相应区间选取为\((x-h,x+h)\)并通过计算样本\(X_1,...,X_n\)在其中频率来估计\(x\)处的密度如下：
\[\begin{split}
f(x_0)=&F'(x)\\
=&\lim_{h\rightarrow 0^+}\frac{F(x+h)-F(x-h)}{2h}\\
=&\lim_{h\rightarrow 0^+}\frac{\mathbb{P}[x-h<X<x+h]}{2h}\\
\end{split}\] 与直方图不同，这种方法下选择的区间只与\(x\)有关且关于\(x\)对称，这使得密度估计仅依赖于\(F\)在\(x\)处对称的导数而不依赖于\(x_0\)处导数。
同样地，给定带宽\(h>0\)，可以通过考虑\(X_1,…,X_n\)在\((x−h,x +h)\)
内的相对频率来构建一个分段常数函数并做如下等价变换： \[\begin{split}\hat{f}_N(x;h):=&\frac{1}{2nh}\sum_{i=1}^n
\mathbb{1}_{\{x-h<X_i<x+h\}}\\
=&\frac{1}{nh}\sum_{i=1}^n
\frac{1}{2}\mathbb{1}_{\{x-h<X_i<x+h\}}\\
=&\frac{1}{nh}\sum_{i=1}^n K(\frac{x-X_i}{h})
\end{split}\] 其中\(K(z)=\frac{1}{2}\mathbb{1}_{\{-1<z<1\}}\)。
由于\(\mathbb{P}[x-h<X<x+h]=\mathbb{P}[-1<\frac{x-X}{h}<1]\)，可以把上述\(K\)看作在区间(-1,1)上的均匀密度。对上述密度函数进行积分：
\[\begin{split}\int_{-\infty}^{+\infty}\hat{f}(x)dx=&\frac{1}{nh}\sum_{i=1}^n\int_{-\infty}^{+\infty}
K(\frac{x-X_i}{h})dx\\
=&\frac{1}{n}\sum_{i=1}^n\int_{-\infty}^{+\infty} K(t)dt\\
=&\int_{-\infty}^{+\infty} K(t)dt
\end{split}\] 显然，只要函数\(K\)的积分等于1，就能保证估计出来的密度函数积分为1。我们很自然地对上式进行推广：用任意分布密度函数替换K。比如，我们用标准正态分布的密度函数作为\(K\)，估计就变成\(\hat{f}(x)=\frac{1}{nh}\sum_{i=1}^n\phi(\frac{x-x_i}{h})\)
除此之外，若\(K\)还满足以下条件则为核函数。
* ：对于每个\(x\)都有\(K(x)\geq0\)；
* ：对于每个\(x\)都有\(K(x)=K(-x)\);
* ：对于每个\(x>0\)都有\(K'(x)\leq0\)。
下面几幅图将简单展示核密度估计的基本思想，其中红色点为样本数据。简单而言，核函数就是在每个数据点上构造一个核函数再加总归一化使得最后函数与坐标轴围成面积为1。
大样本性质
估计偏误
通常，增加模型复杂度可以降低偏差但会同时导致更大的估计方差。在构建非参数模型时需要同时衡量估计偏差和估计方差(也就是variance-bais
tradeoff)，进而使得以均方误差(MSE)为模型总体拟合程度的指标达到最小值。
估计偏误可以写做： \[Bias(\hat{f}(x))=E(\hat{f}(x))-f(x)\]
其中第一项核函数估计值的期望为： \[E\frac{1}{h}k(\frac{X_i-x}{h})=\int_{-\infty}^{\infty}\frac{1}{h}k(\frac{z-x}{h})f(z)dz\]
对核函数括号内变量进行替换\(u=(z-x)/h\),上式等价于 \[\int_{-\infty}^{\infty}k(u)f(x+hu)du\]
上述积分无法求出数值解，又由于\(h\rightarrow
0\)，可以使用泰勒展开对\(f(x)\)在\(x\)处进行二阶展开。可以得到： \[f(x+hu)=f(x)+f^{(1)}(x)hu+\frac{1}{2}f^{(2)}(x)h^2u^2+o(h^2)\]
余项为\(o(h^2)\)。 \[\begin{split}\int_{-\infty}^{+\infty}k(u)f(x+hu)du=&f(x)+f^{(1)}(x)h\kappa_{1}(k)+\frac{1}{2}f^{(2)}(x)h^2\kappa_{2}(k)+o(h^2)\\
=&f(x)+\frac{1}{2 !}f^{(2)}(x)h^{2}\kappa_{2}(k)+o(h^{2})
\end{split}\] 其中\(\kappa_1=\int_{-\infty}^{\infty}k(u)u=0\)，\(\kappa_2=\int_{-\infty}^{\infty}k(u)u^2\)
\[Bias(\hat{f}(x))=E(\hat{f}(x))-f(x)=\frac{1}{2
!}f^{(2)}(x)h^{2}\kappa_{2}(k)+o(h^2)\]
从偏误表达式可以看出估计偏差随着带宽平方增加而增加，因此选择更小的带宽可以减少偏差。同时，偏差也与密度函数的二阶导数\(f^{(2)}(x)\)有关。
估计方差
由于核估计是线性估计量以及\(k(\frac{X_i-x}{h})\)是独立同分布的，可以得到
\[\begin{split}Var\left(\hat{f}(x)\right)=&Var\left(\frac{1}{nh}\sum_{i=1}^n
k\left(\frac{X_i-x}{h}\right)\right)\\
=&\frac{1}{n^2h^2}Var\left(\sum_{i=1}^nk\left(\frac{X_i-x}{h}\right)\right)\\
=&\frac{1}{n^2h^2}\times n\times
Var\left(k\left(\frac{X_i-x}{h}\right)\right)\\
=&\frac{1}{nh^2}Var\left(k\left(\frac{X_i-x}{h}\right)\right)\\
=&\frac{1}{nh^2}E
k\left(\frac{X_i-x}{h}\right)^2-\frac{1}{n}\left(\frac{1}{h} E
k\left(\frac{X_i-x}{h}\right)\right)^2\end{split}\]
根据上节证明过的核函数估计的渐进无偏性，可以得到\(\frac{1}{n}\left(\frac{1}{h} E
k\left(\frac{X_i-x}{h}\right)\right)^2=\frac{1}{n}\left(f(x)+\frac{1}{2!}f^{(2)}\kappa_2(k)+o(h^2)\right)^2\),因此第二项为\(\frac{1}{h}E
k\left(\frac{X_i-x}{h}\right)=f(x)+o(1)=O(\frac{1}{n})\)。
对于第一项，可以将期望用积分形式表达再同样利用变量变换和一阶泰勒展开得到：
\[\begin{split} \frac{1}{n^2h}
Ek\left(\frac{X_i-x}{h}\right)^2=&\frac{1}{n^2h}\int_{-\infty}^{+\infty}k\left(\frac{z-x}{h}\right)^2f(z)dz\\
=&\frac{1}{nh}\int_{-\infty}^{+\infty}k(u)^2f(x+hu)du\\
=&\frac{1}{nh}\int_{-\infty}^{+\infty}k(u)^2(f(x)+O(h))du\\
=&\frac{\left(f(x)+O(h)\right)R(k)}{nh}
\end{split}\] 其中，\(R(k)=\int_{-\infty}^{+\infty}k(u)^2du\)。
由此，我们可以得到 \[Var\left(\hat{f}(x)\right)=\frac{f(x)R(k)}{nh}+O(\frac{1}{n})\]
由于\(h\rightarrow
0,h^{-1}\rightarrow\infty\)，余项\(O(\frac{1}{n})\)是比\(O(\frac{1}{nh})\)还小的级数。
均方误差（MSE）和渐进积分平方误差（AMISE）
通常使用均方误差对模型的准确度进行刻画，其定义为 \[\begin{split}
MSE(\hat{f}(x))=&E(\hat{f}(x)-f(x))^2\\
=&Bias\left(\hat{f}(x)\right)^2+Var\left(\hat{f}(x)\right)\\
= &\left(\frac{1}{2
!}f^{2}(x)h^{2}\kappa_{2}(k)\right)^2+\frac{f(x)R(k)}{nh}\\
=& \frac{\kappa_{2}^2(k)}{4}f^{(2)}(x)^2h^{4}+\frac{f(x)R(k)}{nh}\\
=&AMSE\left(\hat{f}(x)\right)
\end{split}\]
由于以上近似是基于渐进展开结果，因此被称为渐进均方误差(AMSE)。可以注意到渐进均方误差与样本量n，带宽h，核函数有关，并且当\(x\)改变时随着\(f^{(2)}(x)\)和\(f(x)\)改变。
另一个衡量全局拟合准确度的标准为渐进积分平方误差（AMISE）： \[\begin{split}AMISE=&\int_{-\infty}^{+\infty}AMSE\left(\hat{f}(x)\right)dx\\
=&\frac{\kappa_{2}^2(k)}{4}R(f^{(2})h^{4}+\frac{R(k)}{nh}\\
\end{split}\] 其中，\(R(f^{(2)})=\int_{-\infty}^{+\infty}\left(f^{(2)}(x)\right)^2dx\)，\(\kappa_2=\int_{-\infty}^{\infty}k(u)u^2\)。
收敛速度
通常情况下，参数方法按\(\sqrt{n}\)的速度收敛.非参数方法收敛速度小于这个速度。通过前文分析可知核密度估计的收敛速度为\(\sqrt{nh}\),由于\(h\rightarrow 0\)，核密度估计的收敛速度远慢于普通参数估计的收敛速度。因此，相对于参数方法，非参数估计需要更大样本量n。
参数选择
核密度估计的优点之一是可以直观地可视化局部最大值和最小值的位置。然而，核密度估计在近似低频的尾部数据时可能会导致误差增加，并且在具有高峰的多峰分布时也会出现误差。为了解决这个问题，需要使用不同的核函数和带宽估计核密度函数。
带宽的选择
为了求渐进最优带宽，很自然想到对渐进积分平方误差（AMISE）关于带宽求导：
\[\begin{split}\frac{d}{dh}AMISE=&\frac{d}{dh}\left(\frac{\kappa_{2}^2(k)}{(2!)^2}R(f^{(2)})h^{4}+\frac{R(k)}{nh}\right)\\
=&h^{3}\kappa^2_{2}(k)R\left(f^{(2)}\right)-\frac{R(k)}{nh^2}\\
=&0
\end{split}\] 得到解： \[h=\left(\frac{R(k)}{nR(f^{(2)})\kappa_2^2(k)}\right)^{\frac{1}{5}}\]
其中\(\kappa_2=\int_{-\infty}^{\infty}k(u)u^2\)。
最优带宽与\(n^{-1/5}\)正相关，数量级为\(O(n^{-1/5})\)。当带宽取到最优带宽\(h_0\)时，渐进积分平方误差可以写做： \[AMISE_0(k)=\frac{5}{4}\left(\kappa_2^2(k)R(k)^4R(f^{(2)})\right)^{1/5}n^{-4/5}\]
大的带宽可能会导致过平滑（over-smooth），掩盖数据的真实结构；小的带宽会产生尖刻多峰的密度估计，不具有良好数学性质。可以将带宽理解为一个控制核在样本点分布的宽度的缩放因子。在选取带宽时可以采取一些经验性规则，如：
：如果数据分布比较集中，可以选择较小的带宽；如果数据分布比较分散，则需要选择较大的带宽。
：如果问题中存在噪声或者确信密度变化非常显著的区域，则需要使用较小的带宽来避免过拟合；如果数据分布比较光滑，则可以使用较大的带宽。
最优带宽依赖于未知的\(\int_{-\infty}^{+\infty}\left(f^{(2)}(x)\right)^2dx\)，其衡量了核函数平滑度。Plug-in方法使用插值的非参数方法对其进行估计。这种方法不依赖于对分布的先验性假设，具有较强稳健性，但较为复杂。
可以用已知样本标准差\(\hat{\sigma}\)的参考密度分布估计值\(g_{\sigma}\)来代替未知的\(f^{(\nu)}\)。最标准的做法是令\(g_{\sigma}=\phi_{\hat{\sigma}}\)，假设g服从正态分布\(N(0,\hat{\sigma}^2)\)。这个近似方法的思路是：如果真实密度函数是正态分布密度函数，那么计算出的带宽就是最优带宽；如果真实密度函数接近正态分布密度函数，那么计算出的带宽也会接近最优带宽。
所谓经验法则是指选取带宽\(h=\hat{\sigma}C(k)n^{-1/(5)}\)，其中\(\hat{\sigma}\)是样本标准差。
\[\begin{split}C(k)=&R\left(\phi^{(2)}\right)^{-1/5}A(k)\\
=&2\left(\frac{\pi^{1/2}(2!)^3R(k)}{4(4)!\kappa^2_{2}(k)}\right)^{1/5}
\end{split}\] 其中，\(\kappa_2=\int_{-\infty}^{\infty}k(u)u^2\)。下面列举常见核函数对应的\(C_{\nu}(k)\)择以做参考:
| 核函数 | \(C(k)\) |
|---|---|
| Epanechnikov | 2.34 |
| Biweight | 2.78 |
| Triweight | 3.15 |
| Gaussian | 1.06 |
特别地，如果选用高斯函数来近似单变量函数，则最小化积分平方误差的\(h\)最佳选择是 \[h=(\frac{4\hat{\sigma}^5}{3n})^{\frac{1}{5}}\approx 1.06 \hat{\sigma}n^{-\frac{1}{5}}\] 其中\(\hat{\sigma}\)是样本的标准差，这种近似被称为正态分布近似。 除此之外，Silverman(1986)还提供了另一种对数据中异常值更加稳健的经验性法则：将\(\hat{\sigma}\)替换为\(A=min(standard\quad deviation,interquantile\quad range/1.34)\)
交叉验证法与前面所提到的方法有所不同，其直接考虑最小化MISE,是一种靠数据驱动的带宽选择方法。其思想是使用样本两次:一次用于计算KDE，另一次用于评估其在估计f时的性能，两次使用样本均不相同。采用交叉验证的方式可以有效避免结果对样本选取的依赖性。
首先将\(MISE\)的平方项展开 \[\begin{split}
MISE[\hat{f(\cdot;h)}]=&E\left[\int_{-\infty}^{+\infty}\left(\hat{f}(x;h)-f(x)\right)^2dx\right]\\
=&E\left[\int_{-\infty}^{+\infty}
\hat{f}(x;h)^2dx\right]-2E\left[\int_{-\infty}^{+\infty}\hat{f}(x;h)f(x)dx\right]+\int_{-\infty}^{+\infty}
f(x)^2dx
\end{split}\] 其中，\(\kappa_2=\int_{-\infty}^{\infty}k(u)u^2\)。
由于最后一项与带宽\(h\)无关，则最小化\(MISE\left[\hat{f}(\cdot;h)\right]\)问题等价于最小化\(E\left[\int
\hat{f}(x;h)^2dx\right]-2E\left[\int\hat{f}(x;h)f(x)\right]\)。
虽然这两项大小未知，但根据前面证明过的核密度估计的渐进无偏性,可以证明LSCV是其无偏估计量(此处证明略)：
\[LSCV(h):=\int\hat{f}(x;h)^2dx-2n^{-1}\sum_{i=1}^n\hat{f}_{-i}(X_i;h)\]
其中，\(\hat{f}_{-i}(\cdot;h)\)是按照”去一法”(leave-one-out
method)去掉样本点\(X_i\)计算出来的KDE
\[\hat{f}_{-i}(x;h)=\frac{1}{n-1}\sum_{j=1\\j\neq
i}^n K_h(x-X_j)\]
对于LSCV的构造很直观：第一项在构造时就是无偏估计，第二项可以通过样本点\(X_1,...,X_n～f\)的蒙特卡洛模拟来近似\(\int\hat{f}(x;h)f(x)dx\)即，可以通过\(f(x)dx=dF(x)\)来替换\(dF_n(x)\)。那么有， \[\int \hat{f}(x;h)f(x)dx\approx
\frac{1}{n}\sum_{i=1}^n\hat{f}(X_i;h)\]
为了减少对样本的依赖性，采用去一法，进一步用\(\hat{f}_{-i}(X_i;h)\)替代\(\hat{f}(X_i;h)\)。与之前的plug-in方法相反，最小二乘交叉验证法(Least
Square
Cross-Validation)对目标函数的形状几乎没有要求。最小二乘交叉验证方法也被称为无偏交叉验证法，其旨在于寻找使得LSCV最小的带宽h：
\[\hat{h}_{LSCV}:=\arg\min\limits_{h>0}
LSCV\]
然而，LSCV函数的数值优化具有挑战性。在实际应用中，可能存在多个局部极小值，目标函数的粗糙度随n和f的不同而有显著变化。因此，优化例程可能陷入伪解。为了安全起见，总是建议(在可能的情况下)通过绘制范围为h的LSCV(h)来检查解，或者在带宽网格中执行搜索。
似然交叉验证法是Duin(1976)中提出的另一种由数据驱动的选择平滑参数h的方法。与最小二乘交叉验证法不同的是，这种方法旨在选出使得去一法后的对数似然函数最大化的带宽h：
\[\mathcal{L}=lnL=\sum_{i=1}^nln\hat{f}_{-i}(X_i)\]
其中，\(\hat{f}_{-i}(X_i)\)是去掉样本数据\(X_i\)后估计的核密度估计量。
似然交叉验证的主要问题是它受到\(f(x)\)的尾部行为的严重影响，并且在使用常用核函数时可能导致厚尾分布的结果不一致(见Hall
(1987a,1987b))。因此，似然交叉验证方法适用于一系列标准分布(即薄尾分布)，这种方法在统计文献中并不常用。
有偏交叉验证法（Biased
Cross-Validation）结合了plug-in和交叉验证思想。为了计算简单不妨首先考虑二阶核函数(\(\nu=2\))的AMISE表达式: \[AMISE\left[\hat{f}(\cdot;h)\right]=\frac{1}{4}\kappa_2^2(k)R(f'')h^4+\frac{R(k)}{nh}\]
接着使用plug-in方法基于对\(R\left(\hat{f}(\cdot;h)\right)\)的修正来估计\(R(f'')\)。这个修正为： \[\tilde{R}(f''):=R(\hat{f}(\cdot;h))-\frac{R(k'')}{nh^5}=
\frac{1}{n^2}\sum_{i=1}^n\sum_{j=1\\ j\neq i}^n ({k''_h} *
{k''_h})(X_i-X_j)\] 由于\(E\left[R\left(\hat{f}''(\cdot;h)\right)\right]=R(f'')+\frac{R(k'')}{nh^5}+O(h^2)\),对\(R\left(\hat{f}(\cdot;h)\right)\)的修正是为了减少减少偏差。将上式带入AMISE表达式，可以得到有偏交叉验证的目标函数以及最优带宽：
\[\begin{split}
BCV(h):=&\frac{1}{4}\kappa_2^2(k)\tilde{R}(f'')h^4+\frac{R(k)}{nh}\\
\hat{h}_{BCV}:=&\arg loc\min\limits_{h>0} BCV(h)
\end{split}\] 其中，\(\arg
loc\min\limits_{h>0}
BCV(h)\)代表使得BCV在局部取得最小值的带宽h。
对于固定n，当\(h\rightarrow\infty\)时，\(BCV\rightarrow 0\),这意味着\(R(f'')\)趋于0的速度大于\(O(h^{-4})\)。因此，当最小化BCV(h)时，需要注意：因为我们实际上感兴趣的是获得最小的局部最小化。因此，需要避免带宽网格过大的的上极端值，因为它们的存在会错过局部最小值，而会在\(h\rightarrow\infty\)处获得全局最小值。
\(\hat{h}_{BCV}\)相对于\(\hat{h}_{LSCV}\)而言，它的方差比\(\hat{h}_{LSCV}\)小得多。但方差的减少是以偏差增加为代价的，这往往使\(\hat{h}_{BCV}\)大于\(\hat{h}_{MISE}\)。
最后再简单讨论一下自适应带宽的核密度估计方法(Adaptive Kernel
Density)。自适应带宽的核密度估计方法是在固定带宽核密度函数的基础上，通过修正带宽参数而得到的，核函数在每个数据点处的带宽都可以不相同。其思想为：对于样本数据密度较小的点可以选取较大带宽，对于数据点较为紧密的区域选取较小带宽。其形式为：
\[\hat{f}(x)=\frac{1}{n}\sum_{i=1}^n\frac{1}{h_i}K\left(\frac{x-X_i}{h_i}\right)\]
具体带宽估计方法为：首先在某一固定带宽\(h\)下计算出密度函数的初始值\(f_h(x)\)，则可以利用\(F:=\left(\prod_{i=1}^nf_h(x_i)\right)^{1/n}\)和\(h_i=h\sqrt{\frac{F}{f_h(x_i)}}\)计算出上式。
可以参考相关论文：
Bootstrap
choice of bandwidth for density estimation(JASA,1990)
Bootstrap
bandwidth selection
核函数的选择
核的(moment)可以定义为\(\kappa_j(k)=\int_{-\infty}^{+\infty}u^jk(u)du\)，由此可知,核的所有奇数矩均为0。事实上，大多数非参数估计都使用对称的核函数。核的(order),\(\nu\)定义为第一个非零矩的阶数。比如说，如果\(\kappa_1(k)=0\)，但\(\kappa_2(k)>0\)，则\(k\)是一个二阶核；如果，\(\kappa_1(k)=\kappa_2(k)=\kappa_3(k)=0\)，但\(\kappa_4(k)>0\)，那么\(k\)是一个四阶核。因此，对称核函数的阶数永远都是偶数。二阶核为对称非负核，阶数高于二阶的核函数被统一称作高阶核函数，这些核函数存在负数部分且不是密度分布函数，它们也被称为偏误缩减核(bias-reducing
kernels)。
定义核函数的(roughness)为\(R(g)=\int_{-\infty}^{\infty}g(u)^2du\)，即核函数的形状和幅度如何影响平滑和峰度。一个平滑的核函数会产生光滑的密度估计，而一个不光滑的核函数会产生峰状的密度估计。在密度估计中，核函数的roughness对于选择可靠的模型和模型的正确性至关重要。通常，越光滑的核函数具有更好的bias-variance
trade-off，但它们可能会失去一些分布的细节。较粗糙的核函数具有更高的方差，并更容易产生过拟合的结果，但它们可以更好地捕捉分布的局部细节。
假定已经根据样本数据选定了核函数阶数，接下来如何选择具体的核函数呢？通过上节中给出的\(AMISE_0\)表达式可以发现在给定\(\nu\)时，核函数通过\(\kappa_{\nu}(k)R(k)^{\nu}\)来影响渐进估计准确度。在其他影响因素固定时，可以通过选择最优核函数使得AMISE达到最小。对于核函数而言，只有shape重要而scale不重要，因此不妨假设\(\kappa_{\nu}=1\),因此问题变为了如何最小化核函数平滑度\(R(k)=\int_{-\infty}^{+\infty}k(u)^2du\)，约束条件为\(\int_{-\infty}^{+\infty}k(u)du=1\)和\(\int_{-\infty}^{\infty}u^{\nu}k(u)du=1\)。Muller(1984)中证明了这个问题的解也就是最优核函数是一种经过缩放调整的\(k_{\nu,1}\)。由于是否缩放对于核函数而言不重要，这意味着高阶Epanechikov核函数\(k_{\nu,1}\)可能在合适带宽下有最小AMISE。由于这个原因，Epanechikov核函数也被称为最优核函数（optimal
kernel）。
待选择的带宽宽度为最佳宽度的情况下，核函数的形式对MISE的影响甚小，因此在选择核函数上，需要更加注重于核函数的数学性质，比如其求解的难度和导数的阶数。下面列举一些常见的二阶核函数及部分核函数图：
| 二阶核函数 | 表达式 |
|---|---|
| Epeanechnikov | \(k_{2,1}=\frac{3}{4}(1-u^2)_+\) |
| Biweight | \(k_{2,2}=\frac{3}{4}(1-u^2)^2_+\) |
| Triweight | \(k_{2,3}=\frac{35}{32}(1-u^2)^3_+\) |
| Gaussian | \(k_{2,\phi}=\frac{1}{\sqrt{2\pi}}exp(-\frac{u^2}{2})\) |
| Uniform | \(k_{2,5}=\frac{1+u}{2}\) |
| Triangle | \(k_{2,6}=(1-|u|)_+\)) |
事实上，核密度估计继承了其所选用核的全部性质。朴素Rosenblatt核是不连续的，ASH三角核的导数不连续，而Cauchy核函数不存在矩。因此，在选择核函数时需要注意所选择核函数的特殊性质。最保守的选择是一个平滑的、具有明确函数形式的关于原点对称的二阶核函数。
高阶核函数通常是由二阶核函数与多项式相乘得到。通过前文对核函数偏误表达式的推导过程，我们可以回答为什么高阶核函数可以缩小偏误的问题。高阶核函数意味着密度函数有更多阶数的导数，估计偏误级数小于前面讨论的核函数偏误级数\(O_p(h^2)\)。因此，使用高阶核函数的估计误差小于使用二阶核函数时的估计误差，因此高阶核函数也被称为”bias-reducing
kernels”。高阶核函数具有更高的复杂性和自由度，可以更好地拟合非线性模型。常用的四阶核和六阶核列出如下：
| 名称 | 四阶核函数 | 六阶核函数 |
|---|---|---|
| Epeanechnikov | \(k_{4,1}=\frac{15}{8}(1-\frac{7}{3}u^2)k_{2,1}(u)\) | \(k_{6,1}=\frac{175}{64}(1-6u^2+\frac{33}{5}u^4)k_{2,1}(u)\) |
| Biweight | \(k_{4,2}=\frac{7}{4}(1-3u^2)k_{2,2}(u)\) | \(k_{6,2}=\frac{315}{128}(1-\frac{22}{3}u^2+\frac{143}{15}u^4)k_{2,2}(u)\) |
| Triweight | \(k_{4,3}=\frac{27}{16}(1-\frac{11}{3}u^2)k_{2,3}(u)\) | \(k_{6,3}=\frac{297}{128}(1-\frac{26}{3}u^2+13u^4)k_{2,3}(u)\) |
| Gaussian | \(k_{4,\phi}=\frac{27}{16}(1-\frac{11}{3}u^2)k_{2,\phi}(u)\) | \(k_{6,3}=\frac{297}{128}(1-\frac{26}{3}u^2+13u^4)k_{2,\phi}(u)\) |
一方面，使用高阶核函数可以减少估计偏误；而另一方面，使用高阶核函数可能会出现局部负函数值的情况。随着核函数阶数增加，所需要的样本量也随之增加，收敛速率增加最终接近参数速率\(n^{-1}\)，一定程度说明了使用高阶核函数能够减轻非参数估计造成的收敛速度缓慢问题。收敛速度的提高要求密度函数非常平滑，导数可以达到\((\nu+1)\)阶。当密度曲线变得越来越平滑时，估计更接近于参数估计。高阶核函数在局部区域内产生的扰动很小，因此可能会错过数据的细节特征，而把噪声或者本不应该存在的模式也检测出来了。此外，高阶核函数在边缘处会产生奇异点，这可能会影响结果的准确性。因此，计算高阶核密度估计时需要特别注意处理这些问题，例如考虑使用带有带宽参数的正则化方法来控制过拟合和奇异点的影响。
至于如何均衡利弊，Bruce
Hansen对于核函数阶数选择的建议是：在小样本中，二阶内核将是最佳选择，在中等样本中，四阶内核，在较大样本中，可以使用六阶内核。
多元KDE
多元KDE
基于与一元情形相同的原理，核密度估计可以扩展到估计\(\mathbb{R}^P\)中的多元密度函数\(f\)。对于\(\mathbb{R}^P\)中的一个样本\(X_1\),…,\(X_n\)，在\(x\in\mathbb{R}^P\)处的KDE定义为: \[\hat{f}(\mathbf{x})=\frac{1}{n|\mathbf{H}^{1/2}|}\sum_{i=1}^n K(\mathbf{H}^{-1/2}(\mathbf{x}-\mathbf{X}_i))\] 其中，\(K\)是一个依赖于\(p\times p\)维对称正定带宽矩阵\(H\)的多元核函数，带宽矩阵非对角线元素正负号决定多元核密度函数的偏斜方向，对角线元素则控制每个维度上的平滑度。
积核(Product Kernel)
积核是指将多变量密度函数分解为多个一元密度函数的乘积形式。它是多元密度函数的一种推广形式，通过将多元密度函数分解为较小的单元，可以有效简化参数估计以及带宽选择。此时，带宽矩阵取\(H=diag(h_1^2,...,h_p^2)\)，核函数分解为若干个一元核函数的乘积\(K(u)=k(u_1)...k(u_p)\)。在积核方法中，可以选择不同的核函数来构建多元密度函数。积核估计在每个维度下所使用的核函数与一元情形下相同，但是具有不同的平滑参数。
以二元情形为例，积核可以定义为： \[\hat{f}(\mathbf{x})=\frac{1}{nh_1h_2}\sum_{i=1}^n\left[k(\frac{x_1-x_{1i}}{h_1})k(\frac{x_2-x_{2i}}{h_2})\right]\]
\(x_{1i}\)和\(x_{2i}\)分别是，\(\mathbf{x}=(x_1,...,x_d)^T\)的密度函数的估计是逐点的。下图展示了几个常见二元积核。
估计误差和方差和渐进平均积分平方误差（AMISE）
当多元核函数\(K(u)\)取积核形式时，可以计算出估计量的估计偏误、方差及渐进平均积分平方误差分别为：
\[\begin{split}Bias\left(\hat{f}(x)\right)=&\frac{\kappa^2_{\nu}(k)}{\nu!}\sum_{j=1}^q\frac{\partial^\nu}{\partial
x_j^{\nu}}f(x)h_j^{\nu}+o(h_1^{\nu}+...+h_q^{\nu})\\
Var\left(\hat{f}(x)\right)=&\frac{f(x)R(k)}{n|\mathbf{H}^{1/2}|}+O(\frac{1}{n})\\
=&\frac{f(x)R(k)^q}{nh_1h_2...h_q}+O(\frac{1}{n})\\
AMISE\left(\hat{f}(x)\right)=&\frac{\kappa_2^2(k)}{(\nu!)^2}\int\left(\sum_{j=1}^q
\frac{\partial^\nu}{\partial x_j^{\nu}
}f(x)h_j^{\nu}\right)^2\end{split}\]
由于最小化AMISE问题没有数值解，无法通过其求解最优带宽。可以发现，与一元形式相同，多元核函数的AMISE只与\(R(k)\)和\(\kappa^2_\nu(k)\)有关，因此，最优核可以最小化\(R(k)\)。
参数选择
最优带宽选择
假设：密度函数\(f\)平方可积，二阶连续可微以及所有二阶偏导数都是平方可积的；核函数K是球形对称的并且密度函数二阶矩有界、平方可积；\(\mathbf{H}=\mathbf{H_n}\)是一个确定性对称正定矩阵，使得\(n\rightarrow \infty\)时，\(vecH\rightarrow0\)以及\(n|\mathbf{H}|^{1/2}\rightarrow\infty\)。
与一元情形类似，\(MISE[\hat{f}(\cdot;H)]=E[\int
\left(\hat{f}(x;H)-f(x)\right)^2dx]\)，定义使得MISE取得最小值的最优带宽矩阵为
\[\mathbf{H}_{MISE}=\arg\min\limits_{\mathbf{H} \in
SPD_p}MISE\left[\hat{f}(\cdot;H) \right]\] 其中\(SPD_p\)时p阶对称正定矩阵的集合。
事实上，MISE的数值解较难求得，通常用AMISE代替。在以上假设下， \[AMISE[\hat{f}(\cdot;\mathbf{H})]=\frac{1}{4}\kappa_2^2(k)R\left(tr\left(\mathbf{H}Hf(\cdot)\right)\right)+\frac{R(k)}{n|\mathbf{H}|^{1/2}}\]
在\(n\rightarrow \infty\)时，\(AMISE[\hat{f}(\cdot;H)]\rightarrow
0\)。通过最小化AMISE可以求出最优窗宽矩阵。进一步假设一种特殊情形，即当\(\mathbf{H}=h^2I_p\)时，得到最优带宽为 \[h_{AMISE}=\left[\frac{pR(k)}{\kappa_2^2(k)R(tr(Hf(\cdot)))n}\right]^{1/(p+4)}\]
由上述表达式可以看出最优带宽大小为\(O(n^{-1/(p+4)})\)；对于带宽矩阵而言，其大小为\(O(n^{-2/(p+4)})\)。随着维度p的增加，最优带宽也相应增加。
Rule of Thumb法则
与一元情形类似，可以用缩放的正态核函数的平滑度代替未知密度函数的平滑度，得到此时经验性的带宽： \[\mathbf{H}_{NS}=(4/(p+2))^{2/(p+4)}n^{-2/(p+4)}\Sigma\] 在实际估计时可以用正态核的估计样本方差矩阵\(S\)代替\(\Sigma\)来估计\(\hat{\mathbf{H}}_{NS}\)。
最小二乘交叉验证法
回顾一元核密度函数估计的最小二乘交叉估计的出发点：最小化MISE \[MISE(x;h)=\int_{-\infty}^{\infty}\left(\hat{f}(x)-f(x)\right)^2dx=\int_{-\infty}^{\infty}\hat{f}(x)^2(dx)-2\int_{-\infty}^{\infty}\hat{f}(x)f(x)(dx)+\int_{-\infty}^{\infty}f(x)^2(dx)\] 最后一项与h无关，则只需最小化前两项即可。在多元情形下，分别使用其无偏估计进行替代，得到最小化目标函数为 \[LSCV(h_1,...,h_q)=\frac{1}{n^2|\mathbf{H}^{1/2}|}\sum_{i=1}^n\sum_{j=1}^n\bar{K}\left(\mathbf{H}^{-1/2}(X_i-X_j)\right)-\frac{2}{n(n-1)|\mathbf{H}^{1/2}|}\sum_{i=1}^n\sum_{i\neq j}^nK\left(\mathbf{H}^{-1/2}(X_j-X_i)\right)\]
极大似然交叉验证法
与一元情形类似。
结合参数Copula的多元密度估计方法
Copula函数可以对多个随机变量的边缘分布函数建模以考察相关性，Copula也被称为“连接函数”。Sklar定理指出，当\(F\)为多个边缘分布函数\(F_1(\cdot),...,F_p(\cdot)\)的联合分布函数，存在一个Copula函数C使得\(F(x_1,...,x_p)=C(F_1(x_1),...,F_p(x_p))\)。对应的联合密度函数表示为可以反应多个随机变量之间相依性结构的Copula函数C的密度函数与边缘密度函数的乘积形式:\[f(x_1,...,x_p)=c\left(F_1(x_1),...,F_p(x_p)\right)\prod_{i=1}^pf_i(x_i)\]
其中，\(c\left(F_1(x_1),...,F_p(x_p)\right)=\frac{\partial
C\left(F_1(x_1),...,F_p(x_p)\right)}{\partial F_1(x_1)...\partial
F_p(x_p)}\)。
Sklar定理将多元变量联合分布分解为多个边缘分布来对变量间相依性进行建模，同时也为多元密度函数估计提供了另一种半参数视角：将Copula函数的参数估计与边缘密度函数的非参数估计进行结合。这种建模思路有效克服了高维情形下参数过多估计困难的问题，并一定程度上改善了非参数方法收敛慢而导致所需样本容量过多的缺点。
通常，参数Copula常用极大似然（ML）方法进行估计。考虑最简单情形，即各个变量间独立同分布(i.i.d)时，可以将上述联合密度函数取对数表示为：
\[L=\sum_{k=1}^n\log
f(x_1^{(k)},...,x_p^{(k)})\] 其中\(x^{(k)}=(x_1^{(k)},...,x_p^{(k)}),k=1,2,...,n\)。
那么取对数后的联合概率密度分布被分为了两部分的加和\(L=L_C+\sum_{n=1}^N L_n\)。第一部分为\(L_C=\sum_{k=1}^n\log
c(F_1(x_1^{(k)}),...,F_p(x_p^{(k)}))\)，代表由Copula函数C所表示的相依性结构的对数似然贡献；第二部分为\(L_n=\sum_{k=1}^n\log
f_i(x_i^{(k)}),i=1,2,...,p\)，代表每个边缘密度分布的对数似然贡献。在各个变量不相关的情形下，第二部分就是各个密度函数的对数似然函数之和。
但在边缘密度函数不能通过参数方法求得时，上述极大似然方法第二项无法估计，控制变量间相依关系的参数\(\alpha\)必须与各个边缘函数无关。Genest等人(1995)提出了一种通过最大化伪对数似然函数(pseudo
log-likelihood)的半参数方法实现了对Copula参数的估计，并证明了提出估计量的一致性和渐进正态性。具体而言，实现Copula参数\(\alpha\)一致性估计要求最大化如下伪似然函数：
\[L(\alpha)=\sum_{k=1}^n\log\left[c_{\alpha}\{F_{1n}(x_{1k}),...,F_{pn}(x_{pk})\}\right]\]
其中，\(F_{in}\)代表第i个边缘经验分布函数乘\(n/(n+1)\)。
这种缩放调整有效避免了由于某个边缘分布函数接近于1而导致\(\log
c_{\alpha}\{F_{1n}(x_{1k}),...,F_{pn}(x_{pk})\}\)存在的潜在无边界问题。
Copula的分类以及其他估计方法可以参考【Book】Dependence
Modeling with Copulas，本节不多赘述。
KDE在经济学中的应用
应用领域
核密度估计是一种经济学中常用的非参数统计方法，它广泛应用于估计概率密度函数、模拟随机变量和处理高维数据等问题。以下是核密度估计在经济学中的一些应用及例子：
：金融市场变动对企业的财务风险具有显著影响，核密度估计可以帮助企业进行风险管理。例如，可以利用核密度方法估计市场风险，通过对概率密度函数进行模拟，评估市场风险的变化状态。
：在消费者索赔中，核密度估计可用于估计产品价格的概率密度函数，可为消费者提供有关产品定价的信息，以确定索赔金额。通过分析数据样本，可以建立价格分布的概率密度函数，并根据索赔的情况，计算出相应的索赔金额。
：核密度估计可以帮助研究者对经济增长问题进行分析。例如，可以利用核密度方法估计经济增长的概率密度函数，从而找到经济增长的概率分布规律，并据此制定相应的经济政策。
：区域经济差异是研究者关心的重要问题，核密度估计可用于研究区域变量在不同时间间隔内的变化。例如，可以利用核密度方法估计一个城市的人口密度或企业分布，或分析不同区域的收入分布情况。
：在战略分析中，核密度估计常用于研究市场规律，评估竞争格局。例如，核密度估计可用于估计市场中销售额、市场份额和价格等变量的概率密度函数，从而帮助企业制定相应的市场营销策略。
论文实例
本文着重展示KDE在金融中的应用。
Time-varying general dynamic factor models and the measurement of financial connectedness【JOE 2021】
这篇论文针对高维局部平稳时间序列(locally
stationary)，提出了一种新的时变广义动态因子模型。估计基于动态主成分分析结合奇异VAR估计，并将Forni
et al.(2017)
针对平稳数据提出的单边估计方法推广到局部平稳情况。文章证明了在时间序列的样本量T和维数n都增长到无穷大情形下时变脉冲响应函数的估计量的一致性。文章最后将该方法用于实证应用，使用329家标普500成分股的调整对数价差数据研究不同时段、不同频率、不同产业各个公司的时变关联性。
文章重点在于对时变脉冲相应函数的估计，其中估计谱密度时使用了核密度估计方法。总体来说，估计分为主要四个步骤：
1）谱密度估计
对于给定\(\tau\in(0,1)\)和某一带宽\(b_T\in[0,\frac{1}{2}]\)，定义\(T_1(\tau)=\lfloor T\tau\rfloor-\lfloor
Tb_T\rfloor+1\)和\(T_2(\tau)=\lfloor
T\tau \rfloor+\lfloor Tb_T\rfloor\)。令\(M_T=2\lfloor Tb_T\rfloor\)，定义\(X_n(\tau)\)的间隔为\(l\)的自协方差矩阵的局部估计量为： \[\hat{\Gamma}_{n,T}(\tau;l)=
\begin{cases}
\frac{1}{M_T}\sum_{s=T_1(\tau)+l}^{T_2 (\tau)} J(\frac{s-\lfloor \tau
T\rfloor}{M_T}) X_{n,s-l}X'_{ns}, &\tau\in(0,1), 0\leq l\leq
(M_T-1)\\
\hat{\Gamma}_{n,T}^X (\tau;-l), &(-M_T+1)\leq l\leq -1
\end{cases}\] 接着，对于一些\(m_T\in
\mathbb{N}_0\)，使得\(m_T<M_T\),定义\(X_n(\tau)\)的谱密度函数矩阵局部估计量为：
\[\hat{\Sigma}_{n,T}^X(\tau;\theta)=\frac{1}{2\pi}\sum_{l=-m_T}^{m_T}
K(\frac{|l|}{m_T})\hat{\Gamma}_{n,T}(\tau;l)e^{-\iota l
\theta},\tau\in(0,1),\theta\in[-\pi,\pi]\] 其中，\(\iota=\sqrt{-1}\)为虚数单位；本文选取的核函数分别为：
\[J(u)=
\begin{cases}
1, &|u|\leq \frac{1}{2}\\
0,& else
\end{cases}\] \[K(u)=
\begin{cases}
1-|u| &|u|\leq 1\\
0,& else
\end{cases}\] 2）动态主成分分析
\(\hat{\Sigma}_{n,T}^X(\frac{t}{T};\theta_j)\)的第j个特征值和相应的n维特征向量分别为\(\hat{\lambda}_{j;n,T}^X(\frac{t}{T};\theta_j)\)和\(\hat{P}_{j;n,T}^X(\frac{t}{T};\theta_j)\),则对于给定的q个因子，公共部分谱密度估计量如下
\[\hat{\chi}_{j;n,T}^X(\frac{t}{T};\theta_j)=\sum_{j=1}^q
\hat{\lambda}_{j;n,T}^X(\frac{t}{T};\theta_j)\hat{P}_{j;n,T}^X(\frac{t}{T};\theta_j)\hat{P}_{j;n,T}^{X{\dagger}}(\frac{t}{T};\theta_j),\frac{M_T}{2}\leq
t\leq (T-\frac{M_T}{2}),|j|\leq m_T\]
最后通过傅立叶逆变换，可以得到公共部分的局部自协方差矩阵估计量为：
\[\hat{\Gamma}_{n,T;l}^{\chi}(\frac{t}{T};l)=\frac{2\pi}{m_T}\sum_{j=-m_T}^{m_T}\hat{\Sigma}_{n,T}^{\chi}(\frac{t}{T};\theta_j)e^{-\iota
l \theta_j},\frac{M_T}{2}\leq t\leq
(T-\frac{M_T}{2}),l\in\mathbb{Z}\] 3）VAR滤波
假设存在某一整数\(m\)使得\(n=m(q+1)\)，则\(m\)个自回归模型为： \[A_n^{(k)}(\frac{t}{T};L){\chi}_{nt}^{(k)}(\frac{t}{T})=H_{nt}^{(k)}(\frac{t}{T})u_t\]
基于自协方差矩阵估计量并使用AIC准则决定VAR阶数，计算自回归滤波\(A_n^{(k)}(\frac{t}{T};L)\)的Yule-Walker估计值\(\hat{A}_n^{(k)}(\frac{t}{T};L)\)。然后，构建含有m个对角区块的\(n\times n\)分块对角滤波\(\hat{A}_{n,T}(\frac{t}{T};L)\)为\(\hat{A}_{n,T}^{(1)}(\frac{t}{T};L),...,\hat{A}_{n,T}^{(m)}(\frac{t}{T};L)\)；且滤波过程
\[\hat{Z}_{nt}(\frac{t}{T})=\hat{A}_{n,T}(\frac{t}{T};L)X_{nt},\frac{M_T}{2}\leq
t\leq (T-\frac{M_T}{2})\] 4）主成分分析
\(\hat{Z}_{nt}(\frac{t}{T})\)的协方差估计量为
\[\hat{\Gamma}_{n,T}^{\hat{Z}}(\frac{t}{T})=\frac{1}{M_T}\sum_{s=T_1(t/T)}^{T_2(t/T)}J(\frac{s-t}{M_T})\hat{Z}_{ns}(\frac{t}{T})\hat{Z}_{ns}(\frac{t}{T})',\frac{M_T}{2}\leq
t\leq (T-\frac{M_T}{2})\] 令\(\hat{\mu}_{j;n,T}^{\hat{Z}}(\frac{t}{T})\)和\(\hat{P}_{j;n,T}^{\hat{Z}}(\frac{t}{T})\)分别是\(\hat{\Gamma}_{n,T}^{\hat{Z}}(\frac{t}{T})\)的第\(j\)个特征值和标准化的\(n\)维特征向量,且定义 \[\hat{R}_{j;n,T}(\frac{t}{T})=(\hat{R}_{1;n,T}(\frac{t}{T}),...,\hat{R}_{q;n,T}(\frac{t}{T}))=\hat{P}_{j;n,T}^{\hat{Z}}(\frac{t}{T})\sqrt{\hat{\mu}_{j;n,T}^{\hat{Z}}(\frac{t}{T})},\frac{M_T}{2}\leq
t\leq (T-\frac{M_T}{2}),1\leq j \leq q\]
由此可以得到时变脉冲响应函数的估计量为： \[\hat{C}_{n,T}(t;L)=[\hat{A}_{n,T}(\frac{t}{T};L)]^{-1}\hat{R}_{n,T}(\frac{t}{T}),\frac{M_T}{2}\leq
t\leq (T-\frac{M_T}{2})\] 其\((i,j)\)元素\(\hat{c}_{ij;n,T}(t;L)=\sum_{k=0}^{\infty}c_{ijk;n,T}^*(t)L^k\)
在实证中，根据股市交易机制，作者选取带宽参数分别为\(m_T=5\)和\(M_T=22\)分别对应一周交易天数和一个月的交易天数。
Granger causality in risk and detection of extreme risk spillover between financial markets【JOE 2009】
这篇文章将Granger因果检验拓展到风险和极端风险溢出监测领域，提出了一种基于核方法的检测市场间极端风险溢出Granger因果统计量，克服了长期滞后和微弱溢出所带来的传统统计量检验势下降的不足，并贴合与金融市场相对于受远期过去事件影响小于近期事件影响的特征事实。该检验统计量的蒙特卡洛模拟结果显示：在有限样本下，其有较正确的检验尺度，并对各类备择假设具有良好的检验功效。该检验对于研究不同市场不同国家之间的风险共振及风险传染现象具有重大意义。最后，作者使用提出统计量对欧元-美元汇率和日本-美元汇率的极端风险溢出做了一个简单实证，与领域内文献结论一致。
文章重点在于对风险Granger因果检验统计量的构造。首先作者定义一个风险示性变量：
\[Z_{it}:=\mathbf{1}(Y_{it}<-V_{it})\]
其中\(\mathbf{1}(\cdot)\)是示性函数。
当实际损失超过对应分位数点的VaR时取1，否则取0。定义\(t-1\)时期的信息集\(I_{t-1}:=(I_{1(t-1)},I_{2(t-1)})\),其中，\(I_{1(t-1)}=\{Y_{1(t-1)},...,Y_{11}\},I_{2(t-1)}=\{Y_{2(t-1)},...,Y_{21}\}\)。同样可以依照Granger因果检验的传统表达形式写出原假设和备择假设：
\[\mathbb{H}_1^O:E(Z_{1t}|I_{1(t-1)})=E(Z_{1t}|I_{t-1})\quad
versus\quad \mathbb{H}_1^A:E(Z_{1t}|I_{1(t-1)})\neq
E(Z_{1t}|I_{t-1})\] 因此\(\{Y_{1t}\}\)和\(\{Y_{1t}\}\)之间的Granger因果关系可以看作是示性变量\(\{Z_{1t}\}\)和\(\{Z_{1t}\}\)之间的Granger因果关系。与Granger(1969)引入Granger因果关系定义类似，我们可以考虑\(\{Z_{1t}\}\)和\(\{Z_{1t}\}\)的交叉谱。对于双变量方差平稳序列\(\{Z_{1t},Z_{2t}\}\)，正交化的谱密度估计为
\[f(\omega):=\frac{1}{2\pi}\sum_{j=-\infty}^{\infty}\rho(j)e^{-ij\omega},\quad
\omega\in[-\pi,\pi],i=\sqrt{-1}\] 其中，\(\rho(j):=corr(Z_{1t},Z_{2(t-j)})\)。由于，\(\rho(j)\neq \rho(-j)\)，\(f(\omega)\)的值通常是复数。由于\(\rho(\omega)\)和\(f(\omega)\)互为傅立叶变换的关系，它们包含了关于\(\{Z_{1t}\}\)与\(\{Z_{2t}\}\)之间关系的相同信息。因此，无论使用\(\rho(\omega)\)还是\(f(\omega)\)去检验原假设都可行。由于\(f(\omega)\)具有更好的数学性质，在这篇文章中作者使用\(f(\omega)\)带入检验。
在原假设下，对于所有\(j\),有\(\rho(j)=0\)，\(f(\omega)\)可以写作： \[f_1^0(\omega)=\frac{1}{2\pi}\sum_{j=-\infty}^0\rho(j)e^{-ij\omega},\quad
\omega\in[-\pi,\pi]\] 因此，可以根据\(f(\omega)\)和\(f_1^0(\omega)\)之间的差异构造检验统计量。接下来，作者只需通过非参数方法对未知\(f(\omega)\)和\(f_1^0(\omega)\)进行一致估计（此处省去原文证明）即可。
对此，作者首先假设了一个参数VaR模型\(V_{lt}\theta_{l}:=V_1(I_{l(t-1)},\theta_l),\quad
l=1,2,\)并通过模型估计出\(\hat{\theta}_l\)。接着，作者定义\(\hat{Z}_{lt}:=\hat{Z}_{lt}(\hat{\theta}_l),l=1,2\)，其中\(\hat{Z}_{lt}(\hat{\theta}_l):=\mathbf{1}[Y_{lt}<-V_{lt}(\theta_l)]\)。\(\{Z_{1t}\}\)和\(\{Z_{1t}\}\)之间的样本交叉协方差函数可以写作：
\[\hat{C}(j):= \begin{cases}
T^{-1}\sum_{t=1+j}^T(\hat{Z}_{1t}-\hat{\alpha}_1)(\hat{Z}_{2(t-j)}-\hat{\alpha}_2),
& 0\leq j\leq T-1,\\
T^{-1}\sum_{t=1+j}^T(\hat{Z}_{1(t-j)}-\hat{\alpha}_1)(\hat{Z}_{2t}-\hat{\alpha}_2),
& {1-T} \leq j\leq 0,
\end{cases}\] 其中，\(\hat{\alpha}_l:=T^{-1}\sum_{t=1}^T
\hat{Z}_{it}\)。
\(\{Z_{1t}\}\)和\(\{Z_{1t}\}\)之间的样本交叉相关性函数为
\[\hat{\rho}(j):=\hat{C}(j)/\hat{S}_1\hat{S}_2,\quad
j=0,j=0,\pm1,...,\pm(T-1)\] 其中\(\hat{S}_l:=\hat{\alpha}_l(1-\hat{\alpha_l})\)是\(\{\hat{Z}_{lt}\}\)的样本方差。
交叉谱密度\(f(\omega)\)和\(f_1^0(\omega)\)的核估计量可以用\(\hat{\rho}(l)\)表示为： \[\hat{f}(\omega):=\frac{1}{2\pi}\sum_{j=1-T}^{T-1}k(j/M)\hat{\rho}(j)e^{-ij\omega}\]
\[\hat{f}_1^0
(\omega):=\frac{1}{2\pi}\sum_{j=1-T}^{0}k(j/M)\hat{\rho}(j)e^{-ij\omega}\]
作者使用平方距离衡量\(f(\omega)\)和\(f_1^0(\omega)\)间距离： \[L^2(\hat{f},\hat{f}_1^0):=2\pi\int_{-\pi}^\pi
|f(\omega)-f_1^0(\omega)|d\omega=\sum_{j=1}^{T-1}k^2(j/M)\hat{\rho}^2(j)\]
最后，作者构造统计量如下： \[Q_1(M):=[T\sum_{j=1}^{T-1}k^2(j/M)\hat{\rho}^2(j)-C_{1T}(M)]/D_{1T}(M)^{\frac{1}{2}}\]
其中，两个常数为平方距离的渐进均值和方差，分别为： \[C_{1t}(M):=\sum_{j=1}^{T-1}(1-j/k)k^2(j/M)\]
\[D_{1t}(M):=2\sum_{j=1}^{T-1}(1-j/k)(1-(j+1)/T)k^4(j/M)\]
Hong(1996)证明有一大类核函数只要满足以下条件就可以保证统计量的检验势较高：
\[\mathbb{K}(\tau)=\{k(\cdot):k(0)=1,\int_{-\infty}^{\infty}k(z)dz<\infty,k_2:=\lim\limits_{z
\to \infty}\frac{1-k(z)}{z^2}\in(0,\infty) \}\]
这些核函数包括了Parzen核、二次谱核（但不包括Barlett核）。这其中Daniel核\(k_D(z)=sin(\pi z)/\pi
z\)最大化了检验统计量的检验势。但是，无论采用除了truncated（如uniform）核以外的何种核函数对结果影响都不大，再次验证了上节结论。
使用核函数方法天然具有可以对于时间距离更远减少权重对于时间距离更近增加权重的优势，相对于一般方法更符合金融市场的事实情况，一定程度提高了估计精度。
KDE的R实现
补充：R中实现直方图
本例采用R中自带数据集faithful。该数据集中包括了美国黄石公园Old
Faithful温泉喷发持续时间和间隔时间两个变量。输出数据集前五行数据并画出对应散点图如下：
## eruptions waiting
## 1 3.600 79
## 2 1.800 54
## 3 3.333 74
## 4 2.283 62
## 5 4.533 85
## 6 2.883 55
通过散点图可以看出数据集中分布在两块区域，表明原始数据分布存在“双峰”特征。接着画出喷发时间的直方图和密度直方图。
#密度直方图
histograph2=hist(faithful$eruptions,main="histograph of eruptions",
xlab="eruptions",probability = TRUE)
最后，根据喷发时间的最小值和最大值，设置组距bins，进一步画密度直方图。
Bk=seq(min(faithful$eruptions),max(faithful$eruptions),by=0.25)
histograph3=hist(faithful$eruptions,main="histograph of eruptions",
xlab="eruptions",probability = TRUE,breaks = Bk)
#拟合出样本点
rug(faithful$eruptions)
R中实现KDE
在R语言中，画核密度估计函数的包较多，列出但不限于如下几个：
1）density:
这是R内部自带的包，可以使用density()函数绘制核密度估计曲线。
2）KernSmooth:
这个包提供了不同的平滑核函数和带宽选择方法，可以使用bkde()函数绘制核密度估计曲线。
3）MASS:
这个包是用于统计学习和数据分析的包，可以使用kde2d()函数绘制二维核密度估计曲线。
4）ggplot2:
这个包是非常受欢迎的数据可视化包，可以使用geom_density()函数绘制核密度估计曲线。
5）ks:
这个包括了单变量和多变量数据的核平滑估计，包括密度函数、密度导数、累积分布、模态聚类、判别分析、显著模态区域和两样本假设检验等。
本文只列举了一些常见的包，实际上还有许多其他的包也可以用于核密度估计绘图。下面就部分常用包继续使用faithful数据集使用进行R语言演练。
R自带density函数
简单的绘制核密度图可以用 \(\textbf{plot(density(x))}\)命令实现。R默认的核密度估计函数是高斯密度函数；默认的带宽采用Silverman经验法则（Silverman’s rule of thumb）。核函数也可以从”gaussian”, “epanechnikov”, “rectangular”, “triangular”, “biweight”,“cosine”,“optcosine”中进行选择。一行命令即可画出喷发时间变量以biweight为核函数的核密度图如下：
不妨通过作图对比不同窗宽选择对密度函数估计的影响。
par(mfrow=c(2,2))
plot(density(faithful$eruptions,kernel="biweight",bw=8),main="bandwidth=8")
plot(density(faithful$eruptions,kernel="biweight",bw=0.8),main="bandwidth=0.8")
plot(density(faithful$eruptions,kernel="biweight",bw=0.08),main="bandwidth=0.08")
plot(density(faithful$eruptions,kernel="biweight",bw=0.008),main="bandwidth=0.008")
在选取较大窗宽时，核密度估计不能很好刻画原始数据的“双峰”特征。随着窗宽减小，数据特征逐步趋于明显，局部呈现崎岖不平的尖锐特征，最终表现为过拟合。
KernSmooth包
# 若没有安装KernSmooth包需要先安装载入包
#install.packages('KernSmooth', repos='http://mirrors.tuna.tsinghua.edu.cn/CRAN/');
library(KernSmooth);## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009# 计算核密度估计值
density <- bkde(faithful$eruptions)
# 绘制核密度估计函数
plot(density,type = 'l', col = 'blue', lwd = 2,
main = 'Kernel Density Estimate',
xlab = 'Sepal Length', ylab = 'Density',
xlim = range(faithful$eruptions),
ylim = c(0, max(density$y) * 1.2))
KernSmooth 包中的 kb
函数提供了多种核函数，这些核函数如下所示:
* Gaussian kernel (gaussian)
* Epanechnikov kernel (epanechnikov)
* Uniform kernel (rectangular)
* Triangular kernel (triangular)
* Biweight kernel (biweight)
* Cosine kernel (cosine)
* Optcosine kernel (optcosine)
* Logistic kernel (logistic)
* Silverman’s “rule of thumb” optimal kernel (bkde)
不妨在控制其他参数不变的情形下，改变估计所使用核函数类型，分别设置为Epanechnikov核函数、正态核函数、Biweight核函数以及Rectangular核函数。
par(mfrow=c(2,2))
plot(density(faithful$eruptions,kernel="epanechnikov"), main = "Epanechnikov Kernel")
plot(density(faithful$eruptions,kernel="gaussian"), main = "Normal Kernel")
plot(density(faithful$eruptions,kernel="biweight"), main = "Biweight Kernel")
plot(density(faithful$eruptions,kernel="rectangular"), main = "Rectangular Kernel") 
上图可以很好验证前文所提到的核函数选择并不是影响估计结果的最重要的因素，除了Rectangular Kernel估计结果明显有不光滑的凸起外，其他三种核函数估计出来的结果几乎一致。
ks包
ks包功能十分强大，不仅提供了核密度估计的内置函数，包含了前文提到的Silverman经验法则交叉验证以及plug-in插值法等带宽选择方法函数，还涵盖了与核密度估计有关的其他统计方法，如核判别分析、核Copula密度估计等，具体说明详见RDocument说明文档。下面以ks包为例，本节将展示不同带宽选择方法下的核密度估计结果。
##
## The downloaded binary packages are in
## /var/folders/08/fdkvz7q525bcg7dx6zw92p6h0000gn/T//Rtmp0kXr6G/downloaded_packageslibrary(ks) # 加载ks包
par(mfrow=c(2,2))
bw=hpi(faithful$eruptions,1)# 使用plug-in方法选择窗宽迭代阶数为1
plot(density(faithful$eruptions,bw=bw),main="plug_in");
bww=hscv(faithful$eruptions,1);
plot(density(faithful$eruptions,bw=bww),main="Smoothed CV");
bwww=hlscv(faithful$eruptions);
plot(density(faithful$eruptions,bw=bwww),main="LSCV");
bwwww=hns(faithful$eruptions,deriv.order = 2)
plot(density(faithful$eruptions,bw=bwwww),main="Normal Scale");
上述结果表明不同的窗宽选择方法也会对估计结果造成影响，其本质在于不同窗宽选择方法选择了不同窗宽。比较上述四种窗宽选择方法的估计结果可以看出plug-in方法估计结果最接近数据真实密度分布。
编程实现核密度函数
在R中，我们可以使用上节提到的包和函数来进行核密度估计。但是，如果不使用任何包和函数，可以用以下代码实现核密度估计：
根据数据的范围和带宽长度创建一组均匀网格。
数据的范围是指数据的最小值和最大值之间的范围。
关于带宽长度的选择，最常用的是标准得到的带宽长度，常常是Silverman规则；例如，设
\(h = 1.06\hat{\sigma} n^{-1/5}\)，其中
\(\hat{\sigma}\)
是数据的标准偏差，\(n\)
是数据的样本大小。
对于每个网格点，计算该点左右范围内的数据点的密度。
密度可以使用高斯函数（正态分布函数）来实现。高斯函数的公式为：\[\frac{1}{\sqrt{2\pi}h}\exp\left(-\frac{(x -
x_i)^2}{2h^2}\right)\] 其中，\(x_i\)表示当前网格点的值，\(x\)表示数据中的每个点，\(h\)表示带宽长度。
计算每个网格点的密度之和，并将其除以数据点的数量和带宽长度的乘积。
将这些估计值组合在一起，形成一个估计密度曲线。
# 生成一组数据
data <- rnorm(100, mean = 10, sd = 2)
# 选择带宽长度
h <- 1.06 * sd(data) * length(data)^(-1/5)
# 创建均匀网格
x <- seq(min(data), max(data), length.out = 1000)
# 对每个网格点计算密度
dens <- rep(0, length(x))
for (i in seq_along(dens)) {
for (j in seq_along(data)) {
dens[i] <- dens[i] + dnorm(x[i], mean = data[j], sd = h)
}
}
# 将密度除以带宽长度和数据点数量
density <- dens / (length(data) * h)
# 绘制密度估计图
plot(x, density, type = "l")
虽然此代码结果与density()函数的结果相似，但其优化性能远比现成函数差，而且灵活性较差。
参考文献
[1] Scott D W.
Multivariate density estimation: theory, practice, and visualization[M].
John Wiley & Sons, 2015.
[2] Devroye
L. Nonparametric density estimation[J]. The L_1 View, 1985.
[3] Joe
H. Dependence modeling with copulas[M]. CRC press, 2014.
[4] Silverman
B W. Density estimation for statistics and data analysis[M]. CRC press,
1986.
[5] Wand
M P, Jones M C. Kernel smoothing[M]. CRC press, 1994.
[6] Li
Q, Racine J S. Nonparametric econometrics: theory and practice[M].
Princeton University Press, 2007.
[7] Chacón
J E, Duong T. Multivariate kernel smoothing and its applications[M]. CRC
Press, 2018.
[8] Genest
C, Ghoudi K, Rivest L P. A semiparametric estimation procedure of
dependence parameters in multivariate families of distributions[J].
Biometrika, 1995, 82(3): 543-552.
[9] Marron
J S, Wand M P. Exact mean integrated squared error[J]. The Annals of
Statistics, 1992, 20(2): 712-736.
[10] Scott
D W, Terrell G R. Biased and unbiased cross-validation in density
estimation[J]. Journal of the american Statistical association, 1987,
82(400): 1131-1146.
[11] Sheather
S J, Jones M C. A reliable data‐based bandwidth selection method for
kernel density estimation[J]. Journal of the Royal Statistical Society:
Series B (Methodological), 1991, 53(3): 683-690.
[12] Notes On
Nonparametric Density Estimation-Berkeley
[13] Physics 509:
Kernel Density Estimation
[14] Copula
Estimation - Harvard University
[15] Lecture
Notes on Nonparametrics-Bruce Hansen
[16] Chapter
1 Density estimation
[17] Lecture
7: Density Estimation-Washington University
[18] Math
6070 Elements of Density Estimation-University of Utah
[19] 统计计算课件-李东风
[20] Rdocument-ks
[21] Rdocument-KernSmooth