引言

证券市场开辟了企业筹资的新渠道，为经济的可持续增长奠定了良好的基础。现金股利是证券市场最常见、最重要的股利支付方式，公司派发现金股利可以带给股东现金收入，增强其投资能力；投资者通过研究上市公司的现金股利支付倾向，可以判断上市公司未来是否有意愿支付现金股利，并根据判断进行合理的投资。

目前对各指标变量影响上市公司支付现金股利概率的研究通常是用logistic回归与线性模型相结合来分析^[1-3]，但是观测数据包含多个上市公司在不同时间点的观测值，而且不同公司间由于规模和效益等条件的不同使支付现金股利的意愿存在着差异；另外时间的影响是非线性的，采用非参数函数刻画更加合理^[4]。广义部分线性混合模型是在logistic模型基础上的进一步推广, 模型以普通的logistic模型为特例，不仅包含参数部分，同时由于加入非参数部分，从而比普通的logistic模型更加灵活，也能更好地拟合数据。作为广义部分线性混合模型的特殊情况，logistic回归模型的统计推断通常采用极大似然、拟似然或广义矩估计等方法^[4-7]，但这些方法在有限样本下相应的估计结果偏差较大，模型参数的估计效率也会降低；而贝叶斯方法引入了研究者的经验信息，可以在样本较少的情况下得到相对可靠的结果。另外，Kleinman等^[8]在随机效应服从非参数先验分布条件下使用Gibbs抽样导出贝叶斯估计；Tang等^[9]利用惩罚样条逼近非参数函数，并结合Gibbs算法^[10]和MH算法^[11]的优势构建混合抽样方法来估计模型中的未知参数。但文献[8-9]的先验设定都比较复杂，导致后验分布复杂度过高，不利于数值计算。

本文以我国12家A股上市公司为研究对象，首先合理选择影响因素指标，建立logistic广义部分线性混合模型进行分析；然后在文献[9]的基础上采用惩罚样条逼近非参数函数，并进一步简化先验设定，推导出logistic部分线性混合模型各参数的后验分布，得到模型参数的贝叶斯估计；最后通过实证分析证明了所建模型在有限样本下的有效性。

1 变量选取

影响现金股利支付倾向的因素主要有公司的成长性、经营状况、股权结构等。本文选择公司规模、成长水平来体现公司成长性，用融资水平、偿债水平、现金流水平来体现公司的经营状况，以股东获利水平、人均GDP和第一大股东持股比例表示公司的股权结构，研究这8个因素对上市公司现金股利分配意愿的影响。考虑到数据的可操作性，对部分数据以取对数或采用比例的形式来进行分析，变量的定义如表 1所示。

2 模型描述 2.1 广义部分线性混合模型

根据文献[8-9]，广义部分线性混合模型定义为：设y_ij是第i个公司在时间点t_ij的观测值，其中i=1, …, m, j=1, …, n_i，随机变量序列y_i1|u_i，y_i2|u_i, …, y_ini|u_i是条件独立的，则y_ij和条件均值、条件方差满足方程(1)

$ \left\{ \begin{array}{l} p\left( {{y_i}\left| {{\mathit{\boldsymbol{u}}_i}} \right.,\phi } \right) = \exp \left[ {{\phi ^{ - 1}}\left\{ {{y_{ij}}{\theta _{ij}} - b\left( {{\theta _{ij}}} \right)} \right\} + } \right.\\ \;\;\;\;\;\;\;\left. {c\left( {{y_{ij}},\phi } \right)} \right]\\ {\mu _{ij}} = {\rm{E}}\left( {{y_{ij}}\left| {{\mathit{\boldsymbol{u}}_i}} \right.} \right) = \dot b\left( {{\theta _i}} \right),\dot b\left( \theta \right) = \partial b\left( \theta \right)/\partial \theta \\ {v_{ij}} = \mathit{var}\left( {{y_{ij}}\left| {{\mathit{\boldsymbol{u}}_i}} \right.} \right) = \ddot b\left( {{\theta _i}} \right),\ddot b\left( \theta \right) = {\partial ^2}b\left( \theta \right)/\partial {\theta ^2} \end{array} \right. $

(1)

其中$ \phi $是离差参数，c(·, ·)是仅依赖于y_ij和$ \phi $的函数，θ_ij是典范参数。

存在已知的单调联系函数f(·)，使得式(2)成立

$ f\left( {{\mu _{ij}}} \right) \buildrel \Delta \over = {\eta _{ij}} = \mathit{\boldsymbol{x}}_{ij}^{\rm{T}}\mathit{\boldsymbol{\beta }} + \mathit{\boldsymbol{w}}_{ij}^{\rm{T}}{\mathit{\boldsymbol{u}}_i} + g\left( {{t_{ij}}} \right) $

(2)

其中x_ij是一个由p×1解释变量组成的与固定效应相关的向量，β是未知参数向量，w_ij是由q×1个解释变量组成的与随机效应相关的向量，g(·)是时间效应的未知非参数函数。

当随机效应的先验分布采用非参数DP分布假定时，不但估计过程较复杂，而且随机效应的后验分布的均值不再为零，由此导致非参数函数和回归系数的有偏估计^[12]。因此本文假定随机效应服从正态分布，即u_i服从正态分布N_q(0, D)。

由于每个个体分别有n_i个观测变量，因此可以将模型表示如式(3)

$ \begin{array}{l} \left( {\begin{array}{*{20}{c}} {{\eta _{i1}}}\\ \vdots \\ {{\eta _{i{n_i}}}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {{x_{i11}}}& \cdots &{{x_{i{1_p}}}}\\ \vdots&\ddots&\vdots \\ {{x_{i{n_i}1}}}& \cdots &{{x_{i{n_i}p}}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{\beta _1}}\\ \vdots \\ {{\beta _p}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {g\left( {{t_{i1}}} \right)}\\ \vdots \\ {g\left( {{t_{i{n_i}}}} \right)} \end{array}} \right) + \\ \left( {\begin{array}{*{20}{c}} {{w_{i11}}}& \cdots &{{w_{i{1_q}}}}\\ \vdots&\ddots&\vdots \\ {{w_{i{n_i}1}}}& \cdots &{{w_{i{n_i}q}}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{u_{i1}}}\\ \vdots \\ {{u_{iq}}} \end{array}} \right) \end{array} $

(3)

记η_i=(η_i1, …, η_ini)^T, X_i=(x_i1, …, x_ini)^T, W_i=(w_i1, …, w_ini)^T，则式(3)可简化为

$ {\mathit{\boldsymbol{\eta }}_i} = {\mathit{\boldsymbol{X}}_i}\mathit{\boldsymbol{\beta }} + \mathit{\boldsymbol{g}}\left( {{t_i}} \right) + {\mathit{\boldsymbol{W}}_i}{\mathit{\boldsymbol{u}}_i} $

(4)

进一步有

$ \mathop {\left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{\eta }}_1}}\\ \vdots \\ {{\mathit{\boldsymbol{\eta }}_m}} \end{array}} \right)}\limits_{n \times 1} = \mathop {\left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{X}}_1}}\\ \vdots \\ {{\mathit{\boldsymbol{X}}_m}} \end{array}} \right)}\limits_{n \times p} \mathit{\boldsymbol{\beta + }}\mathop {\left( {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{g}}\left( {{t_1}} \right)}\\ \vdots \\ {\mathit{\boldsymbol{g}}\left( {{t_m}} \right)} \end{array}} \right)}\limits_{n \times 1} + \mathop {\left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{W}}_1}}&0&0\\ \vdots&\ddots&\vdots \\ 0&0&{{\mathit{\boldsymbol{W}}_m}} \end{array}} \right)}\limits_{n \times qm} \mathop {\left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{u}}_1}}\\ \vdots \\ {{\mathit{\boldsymbol{u}}_m}} \end{array}} \right)}\limits_{qm \times 1} $

其中$ n = \sum\limits_{i = 1}^m {{n_i}} $，p为x_ij的维度，q为w_ij的维度。从而将整个模型简写为

$ \mathit{\boldsymbol{\eta }} = \mathit{\boldsymbol{X\beta }} + \mathit{\boldsymbol{g}}\left( t \right) + \mathit{\boldsymbol{Wu}} $

(5)

2.2 非参数部分的样条逼近

首先用惩罚样条来逼近未知函数g(·)^[13]

$ \begin{array}{l} \mathit{\boldsymbol{g}}\left( t \right) = \alpha _0^{\left( 1 \right)} + \alpha _1^{\left( 1 \right)}t + \cdots + \alpha _s^{\left( 1 \right)}{t^s} + \sum\limits_{l = 1}^K {\alpha _l^{\left( 2 \right)}} \\ \left( {t - {\kappa _l}} \right)_ + ^s = {\mathit{\boldsymbol{B}}^{\rm{T}}}\left( t \right)\mathit{\boldsymbol{\alpha }} \end{array} $

且有

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{\alpha }} = \left( {{\mathit{\boldsymbol{\alpha }}^{\left( 1 \right)}},{\mathit{\boldsymbol{\alpha }}^{\left( 2 \right)}}} \right)\\ {\mathit{\boldsymbol{\alpha }}^{\left( 1 \right)}} = \left( {\alpha _0^{\left( 1 \right)}, \cdots ,\alpha _s^{\left( 1 \right)}} \right),{\mathit{\boldsymbol{\alpha }}^{\left( 2 \right)}} = \left( {\alpha _1^{\left( 2 \right)}, \cdots ,\alpha _K^{\left( 2 \right)}} \right)\\ \mathit{\boldsymbol{B}}\left( t \right) = {\left( {1,t, \cdots ,{t^s},\left( {t - {\kappa _l}} \right)_ + ^s, \cdots ,\left( {t - {\kappa _k}} \right)_ + ^s} \right)^{\rm{T}}}\\ a_ + ^s = {\left\{ {\max \left( {a,0} \right)} \right\}^s} \end{array} \right. $

式中，s是多项式的次数，K是节点数量，κ_l代表第l个节点。

对于节点的位置，本文采用样本集T的(l+1)(K+2)次分位数作节点，同时可以通过控制样条节点系数的方差来控制惩罚样条的光滑程度。因此假设

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{\alpha }}_h^{\left( 1 \right)} \sim N\left( {0,\sigma _0^2} \right),h = 1, \cdots ,s\\ \mathit{\boldsymbol{\alpha }}_l^2 \sim N\left( {0,\sigma _\mathit{\boldsymbol{\alpha }}^2\left( {{\kappa _l}} \right)} \right),l = 1, \cdots ,K \end{array} \right. $

所以有

$ \mathit{\boldsymbol{\alpha }} = {\left( {{\mathit{\boldsymbol{\alpha }}^{\left( 1 \right)}},{\mathit{\boldsymbol{B}}^{\left( 2 \right)}}} \right)^{\rm{T}}} \sim {N_{s + K + 1}}\left( {{\bf{0}},\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}} \right) $

其中，Λ=diag(Λ₁, Λ₂)，Λ₁=diag(σ₀², …, σ₀²)，Λ₂=diag(σ_α²(κ₁), …, σ_α²(κ_K))，σ₀²是已知常数，σ_α²(κ_l)是光滑参数。

为了估计光滑参数，对σ_α²(κ_l)进行建模

$ \left\{ \begin{array}{l} \lg \left\{ {\sigma _\mathit{\boldsymbol{\alpha }}^2\left( x \right)} \right\} = \delta _0^{\left( 1 \right)} + \delta _1^{\left( 1 \right)}x + \cdots + \delta _q^{\left( 1 \right)}{x^q} + \\ \;\;\;\sum\limits_{k = 1}^{{K_\alpha }} {\delta _k^{\left( 2 \right)}\left( {x - \kappa _k^\mathit{\boldsymbol{\alpha }}} \right)_ + ^q} + v\left( x \right) \buildrel \Delta \over = {\mathit{\boldsymbol{C}}^{\rm{T}}}\left( x \right)\mathit{\boldsymbol{\delta }} + v\left( x \right)\\ \mathit{\boldsymbol{\delta }} \sim {N_q}\left( {0,\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}} \right),v\left( x \right) \sim N\left( {0,\sigma _v^2} \right) \end{array} \right. $

式中δ=(δ⁽¹⁾, δ⁽²⁾)^T，δ⁽¹⁾=(δ₀⁽¹⁾, …, δ_q⁽¹⁾)，δ⁽²⁾=(δ₀⁽²⁾, …, δ⁽²⁾_Kα)，q是多项式样条的次数，C(x)=(1, x, …, x^q, (x－κ₁^α)₊^q, …, (x－κ_Kα^α)₊^q)^T，κ₁^α < … < κ_Kα^α是样条节点，Ω=diag(σ_0δ²I_p+1, γI_Kα)。γ是未知参数，σ_0δ², σ_v²是已知常数，通常取值很小(如0.01)，作为误差项存在。

3 模型的Bayes分析 3.1 模型的全条件概率分布

由2.1节和2.2节总结可得，已知的先验分布分别为

$ p\left( \mathit{\boldsymbol{D}} \right) \sim I{W_r}\left( {{\rho _0},{R_0}} \right),p\left( \phi \right) \sim \Gamma \left( {{a_\phi },{b_\phi }} \right) $

(6)

$ p\left( {\mathit{\boldsymbol{\beta }}\left| \phi \right.} \right) \sim {N_p}\left( {{\mathit{\boldsymbol{\beta }}_0},{\phi ^{ - 1}}{\mathit{\boldsymbol{H}}_0}} \right),p\left( \gamma \right) \sim \Gamma \left( {{a_\gamma },{b_\gamma }} \right) $

(7)

其中，ρ₀, R₀, a$ _\phi $, b$ _\phi $, a_γ, b_γ, β₀, H₀都是已知的超参数。从而全条件概率分布为

$ \begin{array}{l} p\left( {\mathit{\boldsymbol{\beta }},\phi ,\mathit{\boldsymbol{D}},\mathit{\boldsymbol{u}},\mathit{\boldsymbol{\alpha }},{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_2},\mathit{\boldsymbol{\delta }},\gamma \left| {\mathit{\boldsymbol{X}},\mathit{\boldsymbol{Y}},\mathit{\boldsymbol{W}},\mathit{\boldsymbol{T}},\sigma _v^2} \right.} \right) = p\\ \left( \mathit{\boldsymbol{u}} \right)p\left( \mathit{\boldsymbol{\delta }} \right)p\left( \mathit{\boldsymbol{ \boldsymbol{\varLambda} }} \right)p\left( \mathit{\boldsymbol{\alpha }} \right)p\left( \mathit{\boldsymbol{Y}} \right)p\left( {\mathit{\boldsymbol{\beta }},\phi ,\mathit{\boldsymbol{D}},\gamma } \right) \propto {\left| \mathit{\boldsymbol{D}} \right|^{ - \frac{1}{2}}}\\ {\left| \mathit{\boldsymbol{ \boldsymbol{\varOmega} }} \right|^{ - \frac{1}{2}}}{\left| \mathit{\boldsymbol{ \boldsymbol{\varLambda} }} \right|^{ - \frac{1}{2}}}\exp \left( {{\mathit{\boldsymbol{u}}^{\rm{T}}}{\mathit{\boldsymbol{D}}^{ - 1}}\mathit{\boldsymbol{u}} + {\mathit{\boldsymbol{\delta }}^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1}}\mathit{\boldsymbol{\delta }} + {\mathit{\boldsymbol{\alpha }}^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}^{ - 1}}\mathit{\boldsymbol{\alpha }} + } \right.\\ \sigma _v^{ - 2}\sum\limits_{l = 1}^K {{{\left[ {\ln \sigma _v^2\left( {{\kappa _l}} \right) - {\mathit{\boldsymbol{C}}^{\rm{T}}}\left( {{\kappa _l}} \right)\mathit{\boldsymbol{\delta }}} \right]}^2}} + \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^{{n_i}} {\left( {{\phi ^{ - 1}}} \right.} } \\ \left. {\left( {{y_{ij}}{\theta _{ij}} - b\left( {{\theta _{ij}}} \right)} \right)} \right)p\left( {\mathit{\boldsymbol{\beta }},\phi ,\mathit{\boldsymbol{D}},\gamma } \right) \end{array} $

(8)

3.2 固定效应的条件后验分布

由p(β| $ \phi $)~N_p(β₀, $ \phi $^－1H₀)、公式(8)和贝叶斯公式可得固定效应的条件后验分布为

$ \begin{array}{l} p\left( {\mathit{\boldsymbol{\beta }}\left| {\mathit{\boldsymbol{Y}},\mathit{\boldsymbol{X}},\mathit{\boldsymbol{W}},\mathit{\boldsymbol{T}},\mathit{\boldsymbol{u}},\phi } \right.} \right) \propto \exp \left\{ {{\phi ^{ - 1}}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^{{n_i}} {\left( {{y_{ij}}{\theta _{ij}} - } \right.} } } \right.\\ \left. {\left. {b\left( {{\theta _{ij}}} \right)} \right) - 0.5\phi {{\left( {\mathit{\boldsymbol{\beta }} - {\mathit{\boldsymbol{\beta }}_0}} \right)}^{\rm{T}}}\mathit{\boldsymbol{H}}_0^{ - 1}\left( {\mathit{\boldsymbol{\beta }} - {\mathit{\boldsymbol{\beta }}_0}} \right)} \right\} \end{array} $

(9)

式(9)的条件分布不是标准分布，因此根据文献[11]的方法，采用MH算法从后验分布p(β|Y, X, W, T, u, $\phi $)中抽取样本观察值，给定当前状态β^(m)，从正态分布N(β^(m), σ_β²Ω_β)中抽取一个潜在的β，其接受概率为

$ \min \left\{ {1,\frac{{p\left( {\mathit{\boldsymbol{\beta }}\left| {\mathit{\boldsymbol{Y}},\mathit{\boldsymbol{X}},\mathit{\boldsymbol{W}},\mathit{\boldsymbol{T}},\mathit{\boldsymbol{u}},\phi } \right.} \right)}}{{p\left( {{\mathit{\boldsymbol{\beta }}^{\left( m \right)}}\left| {\mathit{\boldsymbol{Y}},\mathit{\boldsymbol{X}},\mathit{\boldsymbol{W}},\mathit{\boldsymbol{T}},\mathit{\boldsymbol{u}},\phi } \right.} \right)}}} \right\} $

正态分布中，${\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_\mathit{\boldsymbol{\beta }}} = {\phi ^{ - 1}}{\left( {\frac{1}{2}\sum\limits_{i = 1}^n {\mathit{\boldsymbol{X}}_i^{\rm{T}}{{\mathit{\boldsymbol{\bar V}}}_i}{\mathit{\boldsymbol{X}}_i} + \mathit{\boldsymbol{H}}_0^{ - 1}} } \right)^{ - 1}}, {{\mathit{\boldsymbol{\bar V}}}_i} = {E_{{y_i}}}\left( {{\partial ^2}b\left( {{\mathit{\boldsymbol{\theta }}_i}} \right)/\partial {\mathit{\boldsymbol{\theta }}_i}\partial \mathit{\boldsymbol{\theta }}_i^{\rm{T}}} \right)|\mathit{\boldsymbol{\beta }} = {\mathit{\boldsymbol{\beta }}^{(m)}} $。同时选择适当的方差调节系数σ_β²，使得潜在转移点的平均接受概率保持在0.25左右。

3.3 随机效应的条件后验分布

对于p(u|Y, X, W, T, D, $ \phi $, β)，由全条件概率分布可得其条件后验分布为

$ \begin{array}{l} p\left( {\mathit{\boldsymbol{\beta }}\left| {\mathit{\boldsymbol{Y}},\mathit{\boldsymbol{X}},\mathit{\boldsymbol{W}},\mathit{\boldsymbol{T}},\mathit{\boldsymbol{D}},\phi ,\mathit{\boldsymbol{\beta }}} \right.} \right) \propto \prod\limits_{i = 1}^n {\exp \left\{ { - \frac{1}{2}} \right.} \\ \left. {\left[ {\sum\limits_{j = 1}^{{n_i}} {{\phi ^{ - 1}}\left( {{y_{ij}}{\theta _{ij}} - b\left( {{\theta _{ij}}} \right)} \right)} + \mathit{\boldsymbol{u}}_i^{\rm{T}}{\mathit{\boldsymbol{D}}^{ - 1}}{\mathit{\boldsymbol{u}}_i}} \right]} \right\} \propto \exp \\ \left\{ { - \frac{1}{2}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^{{n_i}} {{\phi ^{ - 1}}\left( {{y_{ij}}{\theta _{ij}} - b\left( {{\theta _{ij}}} \right)} \right)} } + \mathit{\boldsymbol{u}}_i^{\rm{T}}{\mathit{\boldsymbol{D}}^{ - 1}}{u_i}} \right\} \end{array} $

(10)

因公式(10)的条件分布不是标准分布，所以同样根据后验分布p(u|Y, X, W, T, D, $\phi $, β)采用MH算法抽取样本观察值，在给定当前状态u^(m)下从正态分布N(u^(m), σ_u²D)中抽取一个潜在的u，其接受概率为

$ \min \left\{ {1,\frac{{p\left( {\mathit{\boldsymbol{u}}\left| {\mathit{\boldsymbol{Y}},\mathit{\boldsymbol{X}},\mathit{\boldsymbol{W}},\mathit{\boldsymbol{T}},\mathit{\boldsymbol{D}},\phi ,\mathit{\boldsymbol{\beta }}} \right.} \right)}}{{p\left( {{\mathit{\boldsymbol{u}}^{\left( m \right)}}\left| {\mathit{\boldsymbol{Y}},\mathit{\boldsymbol{X}},\mathit{\boldsymbol{W}},\mathit{\boldsymbol{T}},\mathit{\boldsymbol{D}},\phi ,\mathit{\boldsymbol{\beta }}} \right.} \right)}}} \right\} $

同时，选择适当的方差调节系数σ_u²使得潜在转移点的平均接受概率保持在0.25左右。

3.4 方差项的条件后验分布

首先对于p(D|u)，由D的先验分布、公式(8)和贝叶斯公式可得其后验分布为

$ \begin{array}{l} p\left( {\mathit{\boldsymbol{D}}\left| \mathit{\boldsymbol{u}} \right.} \right) = p\left( {\mathit{\boldsymbol{u}}\left| \mathit{\boldsymbol{D}} \right.} \right)p\left( \mathit{\boldsymbol{D}} \right) \propto {\left| \mathit{\boldsymbol{D}} \right|^{ - \frac{{n + {R_0} + r + 1}}{2}}}\exp \\ \left\{ { - \frac{1}{2}tr\left[ {{\mathit{\boldsymbol{D}}^{ - 1}}\left( {\sum\limits_{i = 1}^n {{\mathit{\boldsymbol{u}}_i}\mathit{\boldsymbol{u}}_i^{\rm{T}}} + {\rho _0}} \right)} \right]} \right\} \end{array} $

(11)

所以有

$ p\left( {\mathit{\boldsymbol{D}}\left| \mathit{\boldsymbol{u}} \right.} \right) \sim I{W_r}\left( {n + {R_0},{\rho _0} + \sum\limits_{i = 1}^n {{\mathit{\boldsymbol{u}}_i}\mathit{\boldsymbol{u}}_i^{\rm{T}}} } \right) $

且有

$ \begin{array}{l} tr\left( {\sum\limits_{i = 1}^n {\mathit{\boldsymbol{u}}_i^{\rm{T}}{\mathit{\boldsymbol{D}}^{ - 1}}{\mathit{\boldsymbol{u}}_i}} } \right) = tr\left( {\sum\limits_{i = 1}^n {{\mathit{\boldsymbol{D}}^{ - 1}}{\mathit{\boldsymbol{u}}_i}\mathit{\boldsymbol{B}}_i^{\rm{T}}} } \right) = \\ tr\left( {{\mathit{\boldsymbol{D}}^{ - 1}}\sum\limits_{i = 1}^n {{\mathit{\boldsymbol{u}}_i}\mathit{\boldsymbol{u}}_i^{\rm{T}}} } \right) \end{array} $

其次对于p(Λ₂|T, α, δ, σ_v²)，由公式(8)可得其后验分布为

$ \begin{array}{l} p\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_2}\left| {\mathit{\boldsymbol{\alpha }},\mathit{\boldsymbol{\delta }},\sigma _v^2} \right.} \right) \propto p\left( {{\mathit{\boldsymbol{\alpha }}^{\left( 2 \right)}}\left| {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_2}} \right.} \right)p\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_2}} \right) \propto \\ \prod\limits_{l = 1}^K {p\left( {\alpha _l^{\left( 2 \right)}\left| {\alpha _l^{\left( 2 \right)}\left( {{\kappa _l}} \right)} \right.} \right)p\left( {\sigma _\mathit{\boldsymbol{\alpha }}^2\left( {{\kappa _l}} \right)\left| {\mathit{\boldsymbol{\delta }},\mathit{\boldsymbol{T}},\sigma _v^2} \right.} \right)} \end{array} $

于是有

$ \begin{array}{l} p\left( {\sigma _\mathit{\boldsymbol{\alpha }}^2\left( {{\kappa _l}} \right)\left| {\mathit{\boldsymbol{\delta }},\mathit{\boldsymbol{T}},\sigma _v^2} \right.} \right) \propto {\left( {\sigma _\mathit{\boldsymbol{\alpha }}^2\left( {{\kappa _l}} \right)} \right)^{ - 0.5}}\exp \\ \left\{ { - \frac{{\alpha _l^{{{\left( 2 \right)}^2}}}}{{2\sigma _\mathit{\boldsymbol{\alpha }}^2\left( {{\kappa _l}} \right)}} - \frac{{{{\left( {\ln \sigma _\alpha ^2\left( {{\kappa _l}} \right) - {\mathit{\boldsymbol{C}}^{\rm{T}}}\left( {{\kappa _l}} \right)\mathit{\boldsymbol{\delta }}} \right)}^2}}}{{2\sigma _v^2}}} \right\} \end{array} $

(12)

因公式(12)的条件分布不是标准分布，所以同样采用MH算法从后验分布p(σ_α²(κ_l)|δ, T, σ_v²)中抽取样本观察值，在给定当前状态为σ_α^2(m)(κ_l)时，从正态分布N(σ_α^2(m)(κ_l), σ²_κlI)中抽取一个潜在的σ_α²(κ_l)，其接受概率为

$ \min \left\{ {1,\frac{{p\left( {\sigma _\mathit{\boldsymbol{\alpha }}^2\left( {{\kappa _l}} \right)\left| {\mathit{\boldsymbol{\delta }},\mathit{\boldsymbol{T}},\sigma _v^2} \right.} \right)}}{{p\left( {\sigma _\mathit{\boldsymbol{\alpha }}^{2\left( m \right)}\left( {{\kappa _l}} \right)\left| {\mathit{\boldsymbol{\delta }},\mathit{\boldsymbol{T}},\sigma _v^2} \right.} \right)}}} \right\} $

同时，选择适当的方差调节系数σ²_κl使得潜在转移点的平均接受概率保持在0.25左右。

最后对于p(γ|δ)，考虑关于γ^－1的条件概率，有

$ p\left( {{\gamma ^{ - 1}}\left| \mathit{\boldsymbol{\delta }} \right.} \right) \propto p\left( {\mathit{\boldsymbol{\delta }}\left| \gamma \right.} \right)p\left( {{\gamma ^{ - 1}}} \right) $

由此可得

$ p\left( {{\gamma ^{ - 1}}\left| \mathit{\boldsymbol{\delta }} \right.} \right) \sim \mathit{\Gamma }\left( {{a_\gamma } + {K_\alpha }/2,{\tau _\gamma } + {\mathit{\boldsymbol{\delta }}^{{{\left( 2 \right)}^{\rm{T}}}}}{\mathit{\boldsymbol{\delta }}^{\left( 2 \right)}}/2} \right) $

(13)

3.5 离差参数的后验条件分布

由全条件概率分布和p(β| $\phi $)~N_p(β₀, $\phi $^－1H₀)，p($\phi $)~Γ(a $_\phi $, b$_\phi $)，可得离差参数$\phi $的后验条件分布为

$ \begin{array}{l} p\left( {\phi \left| {\mathit{\boldsymbol{Y}},\mathit{\boldsymbol{X}},\mathit{\boldsymbol{W}},\mathit{\boldsymbol{T}},\mathit{\boldsymbol{u}},\mathit{\boldsymbol{\beta }},\mathit{\boldsymbol{D}}} \right.} \right) \propto {\phi ^{0.5p + {\alpha _\varphi } - 1}}\exp \left\{ { - \phi } \right.\\ \left[ {{b_\phi } + 0.5{{\left( {\mathit{\boldsymbol{\beta }} - {\mathit{\boldsymbol{\beta }}_0}} \right)}^{\rm{T}}}\mathit{\boldsymbol{H}}_0^{ - 1}\left( {\mathit{\boldsymbol{\beta }} - {\mathit{\boldsymbol{\beta }}_0}} \right)} \right] + {\phi ^{ - 1}}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^{{n_i}} {} } \\ \left. {\left( {{y_{ij}}{\theta _{ij}} - b\left( {{\theta _{ij}}} \right)} \right) + \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^{{n_i}} {c\left( {{y_{ij}},\phi } \right)} } } \right\} \end{array} $

(14)

公式(14)的条件分布不是标准分布，所以采用MH算法从后验分布p($\phi $|Y, X, W, T, u, β, D)中抽取样本观察值，在给定当前状态为$\phi $^(m)时，从正态分布N($\phi $^(m), σ_{$\phi $}²I)中抽取一个潜在的$\phi $，其接受概率为

$ \min \left\{ {1,\frac{{p\left( {\phi \left| {\mathit{\boldsymbol{Y}},\mathit{\boldsymbol{X}},\mathit{\boldsymbol{W}},\mathit{\boldsymbol{T}},\mathit{\boldsymbol{u}},\mathit{\boldsymbol{\beta }},\mathit{\boldsymbol{D}}} \right.} \right)}}{{p\left( {{\phi ^{\left( m \right)}}\left| {\mathit{\boldsymbol{Y}},\mathit{\boldsymbol{X}},\mathit{\boldsymbol{W}},\mathit{\boldsymbol{T}},\mathit{\boldsymbol{u}},\mathit{\boldsymbol{\beta }},\mathit{\boldsymbol{D}}} \right.} \right)}}} \right\} $

同时，选择适当的方差调节系数σ_{$\phi $}²使得潜在转移点的平均接受概率保持在0.25左右。

3.6 惩罚样条参数项α、δ后验条件分布

首先，基于条件分布的定义可得

$ \begin{array}{l} p\left( {\mathit{\boldsymbol{\alpha }}\left| {\mathit{\boldsymbol{Y}},\mathit{\boldsymbol{X}},\mathit{\boldsymbol{W}},\mathit{\boldsymbol{T}},\mathit{\boldsymbol{u}},\phi ,\mathit{\boldsymbol{\beta }},{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_2}} \right.} \right) \propto \exp \\ \left\{ {{\phi ^{ - 1}}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^{{n_i}} {\left( {{y_{ij}}{\theta _{ij}} - b\left( {{\theta _{ij}}} \right)} \right)} } - 0.5{\mathit{\boldsymbol{\alpha }}^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}^{ - 1}}\mathit{\boldsymbol{\alpha }}} \right\} \end{array} $

(15)

公式(15)的条件分布不是标准分布，所以采用MH算法从p(α|Y, X, W, T, u, $\phi $, β, Λ₂)中抽取样本观察值，在给定当前状态为α^(m)时，从正态分布N(α^(m), σ_α²I)中抽取一个潜在的α，其接受概率为

$ \min \left\{ {1,\frac{{p\left( {\mathit{\boldsymbol{\alpha }}\left| {\mathit{\boldsymbol{Y}},\mathit{\boldsymbol{X}},\mathit{\boldsymbol{W}},\mathit{\boldsymbol{T}},\mathit{\boldsymbol{u}},\phi ,\mathit{\boldsymbol{\beta }},{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_2}} \right.} \right)}}{{p\left( {{\mathit{\boldsymbol{\alpha }}^{\left( m \right)}}\left| {\mathit{\boldsymbol{Y}},\mathit{\boldsymbol{X}},\mathit{\boldsymbol{W}},\mathit{\boldsymbol{T}},\mathit{\boldsymbol{u}},\phi ,\mathit{\boldsymbol{\beta }},{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_2}} \right.} \right)}}} \right\} $

同时，选择适当的方差调节系数σ_α²使得潜在转移点的平均接受概率保持在0.25左右。

其次，对于δ直接计算可得其后验分布为

$ \begin{array}{l} p\left( {\mathit{\boldsymbol{\delta }}\left| {\mathit{\boldsymbol{T}},{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_2},\gamma ,\sigma _v^2} \right.} \right) \propto \exp \left\{ { - \frac{1}{{2\sigma _v^2}}\sum\limits_{l = 1}^K {} } \right.\\ \left. {{{\left( {\ln \sigma _\alpha ^2\left( {{\kappa _l}} \right) - {\mathit{\boldsymbol{C}}^{\rm{T}}}\left( {{\kappa _l}} \right)\mathit{\boldsymbol{\delta }}} \right)}^2} - 0.5{\mathit{\boldsymbol{\delta }}^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1}}\mathit{\boldsymbol{\delta }}} \right\} \end{array} $

于是有

$ p\left( {\mathit{\boldsymbol{\delta }}\left| {\mathit{\boldsymbol{T}},{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_2},\gamma ,\sigma _v^2} \right.} \right) \sim {N_{q + {K_\mathit{\boldsymbol{\alpha }}} + 1}}\left( {{\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{\delta }}},{\mathit{\boldsymbol{M}}_\mathit{\boldsymbol{\delta }}}} \right) $

(16)

其中${\mathit{\boldsymbol{M}}_\mathit{\boldsymbol{\delta }}} = {\left( {{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^{ - 1}} + \sum\limits_{l = 1}^K {\mathit{\boldsymbol{C}}\left( {{\mathit{\boldsymbol{\kappa }}_l}} \right){\mathit{\boldsymbol{C}}^{\rm{T}}}\left( {{\mathit{\boldsymbol{\kappa }}_l}} \right)/\sigma _v^2} } \right)^{ - 1}}, {\rm{ }}{\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{\delta }}} = {\mathit{\boldsymbol{M}}_\mathit{\boldsymbol{\delta }}}\sum\limits_{l = 1}^K {\mathit{\boldsymbol{C}}\left( {{\mathit{\boldsymbol{\kappa }}_l}} \right){\rm{ln}}\sigma _\mathit{\boldsymbol{\alpha }}^2\left( {{\mathit{\boldsymbol{\kappa }}_l}} \right)/\sigma _v^2} $。

3.7 Bayes推断

利用Gibbs与MH混合抽样算法获得抽样样本后，可以得到参数和非参数以及随机效应的标准差和贝叶斯估计。

$\{ {\phi ^{(m)}},{\mathit{\boldsymbol{\beta }}^{(m)}},{\mathit{\boldsymbol{\delta }}^{(m)}},{\mathit{\boldsymbol{\alpha }}^{(m)}},{\mathit{\boldsymbol{u}}^{(m)}},{\gamma ^{(m)}},{\mathit{\boldsymbol{D}}^{(m)}},\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_2^{(m)}\} , $ m=1, …, M是由混合算法产生的{ $ \phi $, β, δ, α, u, γ, D, Λ₂}的随机样本，所以$ \phi $, β, δ, α, u, γ, D, Λ₂的贝叶斯估计可以表示为

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{\hat \beta }} = \frac{1}{M}\sum\limits_{m = 1}^M {{\mathit{\boldsymbol{\beta }}^{\left( m \right)}}} ,\hat \phi = \frac{1}{M}\sum\limits_{m = 1}^M {{\phi ^{\left( m \right)}}} ,\hat \gamma = \frac{1}{M}\sum\limits_{m = 1}^M {{\gamma ^{\left( m \right)}}} \\ \mathit{\boldsymbol{\hat \alpha }} = \frac{1}{M}\sum\limits_{m = 1}^M {{\mathit{\boldsymbol{\alpha }}^{\left( m \right)}}} ,\mathit{\boldsymbol{\hat \delta }} = \frac{1}{M}\sum\limits_{m = 1}^M {{\mathit{\boldsymbol{\delta }}^{\left( m \right)}}} ,\mathit{\boldsymbol{\hat u}} = \frac{1}{M}\sum\limits_{m = 1}^M {{\mathit{\boldsymbol{u}}^{\left( m \right)}}} \\ {{\mathit{\boldsymbol{ \boldsymbol{\hat \varLambda} }}}_2} = \frac{1}{M}\sum\limits_{m = 1}^M {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_2^{\left( m \right)}} ,\mathit{\boldsymbol{D}} = \frac{1}{M}\sum\limits_{m = 1}^M {{\mathit{\boldsymbol{D}}^{\left( m \right)}}} \end{array} \right. $

(17)

同时，非参数部分g(t)可以通过$\hat g\left( t \right) = {\mathit{\boldsymbol{B}}^{\rm{T}}}\left( t \right)\mathit{\boldsymbol{\hat \alpha }}$来估计。类似地，参数的后验协方差矩阵的相合估计也可以通过随机样本的样本协方差矩阵得到

$ \begin{array}{l} \mathop {\mathit{var}}\limits^ \wedge \left( {\mathit{\boldsymbol{\beta }}\left| {\mathit{\boldsymbol{Y}},\mathit{\boldsymbol{X}},\mathit{\boldsymbol{W}},\mathit{\boldsymbol{T}}} \right.} \right) = \frac{1}{{M - 1}}\sum\limits_{m = 1}^M {\left( {{\mathit{\boldsymbol{\beta }}^{\left( m \right)}} - \mathit{\boldsymbol{\hat \beta }}} \right)} \\ {\left( {{\mathit{\boldsymbol{\beta }}^{\left( m \right)}} - \mathit{\boldsymbol{\hat \beta }}} \right)^{\rm{T}}} \end{array} $

4 实证分析

本文从中国A股上市公司中选取12家有代表性的上市公司，根据这12家上市公司2004~2011年的数据进行实证分析。

根据文献[3]的方法，以现金股利支付意愿(Y)作为响应变量，以第一大股东持股比例(X₁)、公司规模(X₂)、公司成长性(X₃)、融资能力(X₄)、偿债能力(X₅)、现金流能力(X₆)、股东获利能力(X₇)和人均GDP(X₈)作为解释变量，建立如式(18)的logistic部分线性混合模型来刻画解释变量对支付倾向的影响

$ \left\{ \begin{array}{l} P\left( {{Y_{ij}} = 1} \right) = \frac{1}{{1 + \exp \left( { - {\eta _{ij}}} \right)}}\\ \begin{array}{*{20}{c}} {{\eta _{ij}} = {\beta _1}{X_{1ij}} + {\beta _2}{X_{2ij}} + {\beta _3}{X_{3ij}} + {\beta _4}{X_{4ij}} + {\beta _5}{X_{5ij}} + }\\ {{\beta _6}{X_{6ij}} + {\beta _7}{X_{7ij}} + {\beta _8}{X_{8ij}} + g\left( {{t_{ij}}} \right) + {\mathit{\boldsymbol{u}}_i}} \end{array} \end{array} \right. $

(18)

式中，g(t_ij)是关于时间的光滑函数，作为模型的非参数部分；u_i为相互独立的随机截距项，其先验信息服从正态分布N(0, D)。

将12家上市公司2004~2011年的数据代入公式(16)，通过贝叶斯分析和混合算法抽样，得到参数和超参数的估计结果如表 2和表 3所示。

由表 3可以看出，各个变量都对是否进行现金股利分配产生了影响，说明所建模型是有效的。同时代表第一大股东持股比例、公司规模、公司成长性、偿债能力和股东获利能力的参数系数都为正值，说明这5个因素正面影响上市公司的现金股利分配意愿。原因是第一大股东持股比例越大，上市公司现金股利支付水平越高，公司拥有更大的发展机会、能力和收入，上市公司的现金股利支付率也就越高。另外，现金流能力、融资能力和人均GDP 3个因素的系数都为负值，说明这3个因素与现金股利分配意愿成反比例，这与已有结论相一致，即公司偿债能力、资本市场发展水平和经济发展水平越高，公司现金股利支付率越低。

时间t定义为t=(t_year－2003)/5，其中t_year为年份。然后根据非参数估计结果画图，得到g(t)的样条逼近如图 1。由图可知，时间函数的影响效应曲线是不规则的连续函数，因此需要利用非参数函数刻画。这也印证了本文所建模型的正确性。

5 结束语

本文简化了随机效应部分的先验分布，并采用Gibbs与MH抽样相结合的混合抽样算法，更高效地获得了参数的贝叶斯估计。所选择的影响因素指标对于上市公司的现金股利分配意愿有着不同程度的影响，且影响效果与显示情况基本相符，说明所构建的logistic广义部分线性混合模型在实际应用中是可行的。

另外本文方法也可以针对其他6种不同的股利支付方式进行进一步研究。