引言

随着国家不断出台政策调控楼市，房地产行业已然告别高速增长的时代，因此对于房地产上市公司盈利能力的研究尤为重要。通过研究盈利能力，房地产上市公司可以发现自身经营管理中存在的问题，提高自身业绩，而投资者则可以根据盈利能力判断股票价格变化和股息的分配，进行合理的股票投资。

房地产上市公司盈利能力数据是时间序列数据，同时在时间截面上含有多个公司的多个指标，所以应作为纵向数据进行研究。基于纵向数据的特点，目前主要用线性混合效应模型处理纵向数据，但是线性效应模型在很多情况下处理效果不甚理想，本文综合参数混合效应模型与非参数混合效应模型，采用半参数部分线性混合效应模型进行研究。半参数部分线性混合效应模型具有参数模型和非参数模型两者的优点，有效利用数据所含信息的同时考虑了信息不完整的变量，有利于更全面地分析问题。

推断半参数部分线性混合效应模型以往一般采用经典统计方法，主要有截面核法^[1-4]和基于联合最大惩罚似然的光滑样条和惩罚样条方法^[5-7]。相较于这些经典统计方法，贝叶斯统计方法在小样本情况下得到的点估计和区间估计更准确，并且可以通过后验分布将无关的参数通过积分处理掉，有效简化了模型的推断。目前，已经有不少学者利用贝叶斯方法推断统计模型，如王佐仁等^[8]研究了线性混合效应模型的贝叶斯推断，提出了随机效应协方差矩阵的对数的多元正态逼近方法；Li等^[9]研究了纵向数据下半参数混合模型的贝叶斯推断，分析比较了随机效应分别服从正态分布和非参数分布时的估计效果。但是基于共轭先验分布并利用贝叶斯方法推断半参数部分线性混合效应模型方面的研究还未见报道。

本文参考文献[8-9]的方法，选择贝叶斯方法进行模型推断。以我国A股市场的50家房地产上市公司为研究对象，首先根据上市公司盈利指标选取的原则^[10]建立合理的指标体系，将解释变量中与盈利能力的关系比较明确的变量建立参数模型；为了更好地体现时间对于盈利能力的影响，采用非参数函数代替以往文献中常用的线性函数模拟与盈利能力的关系不明确的时间效应，选用半参数部分线性混合效应模型进行分析；选择贝叶斯方法对模型进行推断，利用光滑样条逼近非参数部分，将模型转化为线性混合效应模型，并根据所得后验分布利用Markov Chain Monte Carlo(MCMC)方法^[11]给出了模型中参数的贝叶斯估计；最后通过实证分析证明了所建模型的准确性。

1 变量选取

本文分别从上海证券交易所和深圳证券交易所各选取25家2005年底前上市的房地产公司(不含ST公司)。以这50家具有代表性的来自不同地区的房地产上市公司作为研究对象，截取2006年~2015年共10个年度的纵向数据(所有数据取自国泰安CSMAR金融经济数据库)进行实证分析。

目前，反映上市公司盈利能力的指标主要包括营业利润率、成本费用利润率、盈余现金保障倍数、总资产报酬率、净资产收益率和资本收益率。实际事务中，上市公司一般采用每股收益、每股股利、市盈率、每股净资产等指标评价其获利能力。本文选择净资产收益率作为反映房地产上市公司盈利能力的指标，并从资产结构、负债结构、公司规模、公司风险、发展潜力5个方面选择了流动资产率、资产负债率、总资产、主营业务收入、主营业务收入增长率5个影响盈利能力的因素指标。由于总资产和主营业务收入是数量指标，与其他比率指标相差较大，为了方便分析，取二者的对数值进行研究。所选盈利及其影响因素指标变量如表 1所示。

此外，不同年份的房地产价格走势、市场需求特点也有差别，从而不同的房地产上市公司在不同年份的盈利能力可能存在差异。以往的研究一般直接用线性函数模拟时间效应，但是现实中时间对于盈利能力的影响是很复杂的，一般都是非线性关系，因此本文利用非参数函数刻画时间的影响，进而建立半参数部分线性混合效应模型进行分析。

2 模型描述 2.1 半参数部分线性混合效应模型

设Y_ij是第i个个体在时间点t_ij的观测值，其中i=1, …, m, j=1, …, n_i，总样本数为n= $\sum\limits_{i = 1}^m {} $ n_i。假设Y_ij服从半参数部分线性混合模型

$ {Y_{ij}} = \mathit{\boldsymbol{X}}_{ij}^{\rm{T}}\mathit{\boldsymbol{\beta }} + f\left( {{t_{ij}}} \right) + \mathit{\boldsymbol{Z}}_{ij}^{\rm{T}}{\mathit{\boldsymbol{b}}_i} + {\varepsilon _{ij}} $

(1)

其中X_ij=(X_ij1, …, X_ijp)^T是与固定效应相关的解释变量构成的p×1维向量；Z_ij=(Z_ij1, …, Z_ijq)^T是与随机效应相关的解释变量构成的q×1维向量；β是p×1维未知参数向量；b_i是第i个个体的q×1维随机效应向量，b_i相互独立且服从协方差D的正态分布N(0, D)；f(t)是关于观测时间t的未知二次可微光滑函数，是模型的非参数部分；ε_ij是随机误差，相互独立，服从正态分布N(0, σ²)且与b_i独立。

2.2 模型参数化

为了对参数β、D和σ²进行估计，首先采用光滑样条方法对未知光滑函数f(·)进行参数化。

利用Li等^[9]的思想，假设f(t)服从Wiener积分过程先验^[12]

$ f\left( t \right) = {\delta _0} + {\delta _1}t + {\lambda ^{ - 1/2}}\int_0^t {\left( {t - u} \right){\rm{d}}W\left( u \right)} ,t \in \left[ {{T_1},{T_2}} \right] $

(2)

其中δ=(δ₀, δ₁)^T是未知的2×1维参数向量，λ是调节参数，T₁和T₂确定t的范围， $\int_0^t {} $ (t－u)dW(u)是标准Wiener过程W(t)的随机积分。由于f(t)已经包含截距项和线性项，本文假设方程(1) 中X_ij不存在与1或者t相关的行。

根据Zhang等^[5]的研究，记T=(1, t⁰)，其中t⁰是由t_ij(i=1, …, m, j=1, …, n_i)中所有离散观测点排序后构成的r×1维向量(r为观测时点的个数)，1是由1构成的r×1维向量。令

$ {\mathit{\boldsymbol{Y}}_i} = {\left( {{T_{i1}}, \cdots ,{Y_{i{n_i}}}} \right)^{\rm{T}}} $

并同样定义X_i、Z_i和ε_i。记N_i为第i个个体的关联矩阵，用来映射t_i=(t_i1, …, t_ini)^T和t⁰。如果t_ij=t_l⁰，则有N_i的第(j, l)个元素为1，否则为0，从而模型(1) 可写成

$ {\mathit{\boldsymbol{Y}}_i} = {\mathit{\boldsymbol{X}}_i}\mathit{\boldsymbol{\beta }} + {\mathit{\boldsymbol{N}}_i}f + {\mathit{\boldsymbol{Z}}_i}{\mathit{\boldsymbol{b}}_i} + {\mathit{\boldsymbol{\varepsilon }}_i} $

(3)

其中f=(f(t₁⁰), …, f(t_r⁰))^T。

记Y=(Y₁^T, …, Y_m^T)^T，并同样定义X、N、b和ε。同时令Z=diag(Z₁, …, Z_m)表示对角元素为Z₁, …, Z_m的块对角矩阵，则有

$ \mathit{\boldsymbol{Y}} = \mathit{\boldsymbol{X\beta }} + \mathit{\boldsymbol{Nf}} + \mathit{\boldsymbol{Zb}} + \mathit{\boldsymbol{\varepsilon }} $

(4)

Li等^[9]指出，f(t)的Wiener积分先验假定等价于对有限维光滑样条f的先验作以下设定：令K表示Wiener积分过程在t⁰处的协方差，若t_j⁰≥t_k⁰，则矩阵K的第(j, k)个和第(k, j)个元素为(t_j⁰)²·(3t_j⁰－t_k⁰)/6。此外令B=L(L^TL)^－1，其中L为满足K=LL^T和L^TT=0的r×(r－2) 维满秩矩阵。

根据文献[13]，非参数分量f可通过一一对应的线性变换表示为

$ \mathit{\boldsymbol{f}} = \mathit{\boldsymbol{T\delta }} + \mathit{\boldsymbol{Ba}} $

(5)

其中δ是一个2×1维模型参数向量。假定服从先验分布δ~N(0, φ I)，φ为已知的特别大的常数；a是一个r×1维随机效应向量，服从正态分布a~N(0, τI)，τ是未知的光滑参数。将f=Tδ+Ba代入式(4)，从而可将式(4) 表示成如式(6) 的线性混合模型

$ \mathit{\boldsymbol{Y}} = {\mathit{\boldsymbol{X}}_ * }{\mathit{\boldsymbol{\beta }}_ * } + {\mathit{\boldsymbol{Z}}^{\left( 1 \right)}}\mathit{\boldsymbol{b}} + {\mathit{\boldsymbol{Z}}^{\left( 1 \right)}}\mathit{\boldsymbol{a}} + \mathit{\boldsymbol{\varepsilon }} $

(6)

其中Z⁽¹⁾=Z，X_*=(X, NT)，β_*=(β^T, δ^T)^T，Z⁽²⁾=NB，b和a为模型的随机效应。

3 贝叶斯分析 3.1 固定效应的条件后验分布

假设固定效应β_*的先验为多元正态分布，即

$ {\rm{ \mathsf{ π} }}\left( {{\mathit{\boldsymbol{\beta }}_ * }} \right) \sim {N_{p + 2}}\left( {{\mathit{\boldsymbol{\beta }}_0},{\sum _0}} \right) $

则β_*的后验分布为

$ \begin{array}{l} {\rm{ \mathsf{ π} }}\left( {{\mathit{\boldsymbol{\beta }}_ * }\left| {\mathit{\boldsymbol{b}},\mathit{\boldsymbol{a}},{\sigma ^2},\mathit{\boldsymbol{Y}}} \right.} \right) \propto \exp \left\{ { - \frac{1}{2}\left( {{{\left( {{\mathit{\boldsymbol{\beta }}_ * },\mathit{\boldsymbol{\hat \beta }}} \right)}^{\rm{T}}}} \right.} \right.\\ \left( {\frac{1}{{{\sigma ^2}}}\mathit{\boldsymbol{X}}_ * ^{\rm{T}}{\mathit{\boldsymbol{X}}_ * }} \right)\left( {{\mathit{\boldsymbol{\beta }}_ * } - \mathit{\boldsymbol{\hat \beta }}} \right) + {\left( {{\mathit{\boldsymbol{\beta }}_ * } - {\mathit{\boldsymbol{\beta }}_0}} \right)^{\rm{T}}}\sum _0^{ - 1}\left( {{\mathit{\boldsymbol{\beta }}_ * } - } \right.\\ \left. {\left. {\left. {{\mathit{\boldsymbol{\beta }}_0}} \right)} \right)} \right\} \end{array} $

其中 ${\mathit{\boldsymbol{\hat \beta }}}$ =(X_*^TX_*)^－1X_*^T(Y－Z⁽¹⁾b－Z⁽²⁾a)。

利用多元配方方法^[14]配方可得

$ \begin{array}{l} {\rm{ \mathsf{ π} }}\left( {{\mathit{\boldsymbol{\beta }}_ * }\left| {\mathit{\boldsymbol{Y}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{a}},{\sigma ^2}} \right.} \right) \propto \exp \left\{ { - \frac{1}{2}\left( {\left( {{\mathit{\boldsymbol{\beta }}_ * } - \left( {{\mathit{\boldsymbol{\beta }}_0} - } \right.} \right.} \right.} \right.\\ \left. {\left. {{{\left. {\left. {{\mathit{\boldsymbol{\beta }}^ * }} \right)} \right)}^{\rm{T}}}\left( {\frac{1}{{{\sigma ^2}}}\mathit{\boldsymbol{X}}_ * ^{\rm{T}}{\mathit{\boldsymbol{X}}_ * } + \sum _0^{ - 1}} \right)\left( {{\mathit{\boldsymbol{\beta }}_ * } - \left( {{\mathit{\boldsymbol{\beta }}_0} - {\mathit{\boldsymbol{\beta }}^ * }} \right)} \right)} \right)} \right\} \end{array} $

其中β^*=(X_*^TX_*+σ²Σ₀^－1)^－1X_*^TX_*(β₀－ ${\mathit{\boldsymbol{\hat \beta }}}$ )。通过核函数可得β_*的后验分布为

$ \begin{array}{l} \;\;\;\;\;{\mathit{\boldsymbol{\beta }}_ * }\left| {\mathit{\boldsymbol{Y}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{a}},{\sigma ^2}} \right. \sim {N_{p + 2}}\left( {\left( {{\mathit{\boldsymbol{\beta }}_0} - {\mathit{\boldsymbol{\beta }}^ * }} \right),\left( {{\sigma ^{ - 2}}\mathit{\boldsymbol{X}}_ * ^{\rm{T}}{\mathit{\boldsymbol{X}}_ * } + } \right.} \right.\\ \left. {{{\left. {\sum _0^{ - 1}} \right)}^{ - 1}}} \right) \end{array} $

(7)

3.2 随机效应的条件后验分布 3.2.1 随机效应b_i的条件分布

已知b_i服从多元正态分布N_q(0, D)，不妨令

$ {\mathit{\boldsymbol{w}}_i} = {\mathit{\boldsymbol{Y}}_i} - {\mathit{\boldsymbol{X}}_i}{\mathit{\boldsymbol{\beta }}_ * } - {\mathit{\boldsymbol{Z}}^{\left( 2 \right)}}\mathit{\boldsymbol{a}} $

则第i个个体似然方程为

$ \begin{array}{l} \;\;\;\;\;\;\;l\left( {{\mathit{\boldsymbol{\beta }}_ * },{\mathit{\boldsymbol{b}}_i},\mathit{\boldsymbol{a}},{\sigma ^2}\left| {{\mathit{\boldsymbol{Y}}_i}} \right.} \right) \propto \exp \left\{ { - \frac{1}{{2{\sigma ^2}}}{{\left( {{\mathit{\boldsymbol{b}}_i} - {{\mathit{\boldsymbol{\hat b}}}_i}} \right)}^{\rm{T}}}} \right.\\ \left. {{{\left( {\mathit{\boldsymbol{Z}}_i^{\left( 1 \right)}} \right)}^{\rm{T}}}\mathit{\boldsymbol{Z}}_i^{\left( 1 \right)}\left( {{\mathit{\boldsymbol{b}}_i} - {{\mathit{\boldsymbol{\hat b}}}_i}} \right)} \right\} \end{array} $

其中 ${{\mathit{\boldsymbol{\hat b}}}_i}$ =((Z_i⁽¹⁾)^TZ_i⁽¹⁾)^－1(Z_i⁽¹⁾)^Tw_i。

同样利用多元配方法^[14]配方可得第i个个体的随机效应b_i的后验分布为

$ \begin{array}{l} {\rm{ \mathsf{ π} }}\left( {{\mathit{\boldsymbol{b}}_i}\left| {{\mathit{\boldsymbol{Y}}_i},\mathit{\boldsymbol{\beta }},\mathit{\boldsymbol{D}},{\sigma ^2}} \right.} \right) \propto \exp \left\{ { - \frac{1}{{2{\sigma ^2}}}{{\left( {{\mathit{\boldsymbol{b}}_i} - {{\mathit{\boldsymbol{\hat b}}}_i}} \right)}^{\rm{T}}}} \right.\\ \left. {{{\left( {{\mathit{\boldsymbol{Z}}_1}} \right)}^{\rm{T}}}\mathit{\boldsymbol{Z}}_i^{\left( 1 \right)}\left( {{\mathit{\boldsymbol{b}}_i} - {{\mathit{\boldsymbol{\hat b}}}_i}} \right)} \right\}\exp \left\{ { - \frac{1}{2}\mathit{\boldsymbol{b}}_i^{\rm{T}}{\mathit{\boldsymbol{D}}^{ - 1}}{\mathit{\boldsymbol{b}}_i}} \right\} \propto \exp \\ \left\{ { - \frac{1}{2}{{\left( {{\mathit{\boldsymbol{b}}_i} - \mathit{\boldsymbol{b}}_i^ * } \right)}^{\rm{T}}}\left( {\frac{1}{{{\sigma ^2}}}{{\left( {\mathit{\boldsymbol{Z}}_i^{\left( 1 \right)}} \right)}^{\rm{T}}}\mathit{\boldsymbol{Z}}_i^{\left( 1 \right)} + {\mathit{\boldsymbol{D}}^{ - 1}}} \right)\left( {{\mathit{\boldsymbol{b}}_i} - } \right.} \right.\\ \left. {\left. {\mathit{\boldsymbol{b}}_i^ * } \right)} \right\} \end{array} $

其中b_i^*=((Z_i⁽¹⁾)^TZ_i⁽¹⁾+σ²D^－1)^－1(Z_i⁽¹⁾)^T(Y_i－X_iβ－Z⁽²⁾a)。

根据核函数，可得随机效应b_i的后验分布为

$ \begin{array}{l} \;\;\;\;\;\;\;{\mathit{\boldsymbol{b}}_i}\left| {{\mathit{\boldsymbol{Y}}_i},{\mathit{\boldsymbol{\beta }}_ * },\mathit{\boldsymbol{D}},{\sigma ^2}} \right. \sim {N_q}\left( {\mathit{\boldsymbol{b}}_i^ * ,\left( {{\sigma ^{ - 2}}{{\left( {\mathit{\boldsymbol{Z}}_i^{\left( 1 \right)}} \right)}^{\rm{T}}}\mathit{\boldsymbol{Z}}_i^{\left( 1 \right)} + } \right.} \right.\\ \left. {{{\left. {{\mathit{\boldsymbol{D}}^{ - 1}}} \right)}^{ - 1}}} \right) \end{array} $

(8)

3.2.2 随机效应a的条件分布

令v=Y－X_*β_*－Z⁽¹⁾b，则似然方程可写为

$ \begin{array}{l} \;\;\;\;\;\;\;l\left( {{\mathit{\boldsymbol{\beta }}_ * },\mathit{\boldsymbol{b}},\mathit{\boldsymbol{a}},{\sigma ^2}\left| \mathit{\boldsymbol{Y}} \right.} \right) \propto \exp \left\{ { - \frac{1}{{2{\sigma ^2}}}{{\left( {\mathit{\boldsymbol{a}} - \mathit{\boldsymbol{\hat a}}} \right)}^{\rm{T}}}} \right.\\ \left. {{{\left( {{\mathit{\boldsymbol{Z}}^{\left( 2 \right)}}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{Z}}^{\left( 2 \right)}}\left( {\mathit{\boldsymbol{a}} - \mathit{\boldsymbol{\hat a}}} \right)} \right\} \end{array} $

其中 ${\mathit{\boldsymbol{\hat a}}}$ =((Z⁽²⁾)^TZ⁽²⁾)^－1(Z⁽²⁾)^Tv。

同时a~N_r－2(0, τI)，则a的条件后验分布为

$ \begin{array}{l} {\rm{ \mathsf{ π} }}\left( {\mathit{\boldsymbol{a}}\left| {\mathit{\boldsymbol{Y}},{\mathit{\boldsymbol{\beta }}_ * },\mathit{\boldsymbol{b}},\mathit{\boldsymbol{D}},\mathit{\boldsymbol{\tau }},{\sigma ^2}} \right.} \right) \propto \exp \left\{ { - \frac{1}{2}{{\left( {\mathit{\boldsymbol{a}} - {\mathit{\boldsymbol{a}}^ * }} \right)}^{\rm{T}}}} \right.\\ \left. {\left( {\frac{1}{{{\sigma ^2}}}{{\left( {{\mathit{\boldsymbol{Z}}^{\left( 2 \right)}}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{Z}}^{\left( 2 \right)}} + {{\left( {\tau \mathit{\boldsymbol{I}}} \right)}^{ - 1}}} \right)\left( {\mathit{\boldsymbol{a}} - {\mathit{\boldsymbol{a}}^ * }} \right)} \right\} \end{array} $

其中a^*=((Z⁽²⁾)^TZ⁽²⁾+σ²(τI)^－1)^－1(Z⁽²⁾)^T(Y－ X_*β_*－Z⁽¹⁾b)。

根据核函数，可得a的后验分布为

$ \begin{array}{l} \;\;\;\;\;\;\mathit{\boldsymbol{a}}\left| {\mathit{\boldsymbol{Y}},{\mathit{\boldsymbol{\beta }}_ * },\mathit{\boldsymbol{b}},\mathit{\boldsymbol{D}},\mathit{\boldsymbol{\tau }},{\sigma ^2}} \right. \sim {N_{r - 2}}\left( {{\mathit{\boldsymbol{a}}^ * },\left( {{\sigma ^{ - 2}}{{\left( {{\mathit{\boldsymbol{Z}}^{\left( 2 \right)}}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{Z}}^{\left( 2 \right)}} + } \right.} \right.\\ \left. {{{\left. {{\tau ^{ - 1}}\mathit{\boldsymbol{I}}} \right)}^{ - 1}}} \right) \end{array} $

(9)

3.3 随机效应方差项的条件后验分布 3.3.1 D的条件分布

本文考虑逆Wishart分布(公式中用IW表示)作为D的先验分布，即D~IW_q(ρ₀, R₀)(ρ₀和R₀为超参数)，则

$ p\left( \mathit{\boldsymbol{D}} \right) \propto {\left| \mathit{\boldsymbol{D}} \right|^{ - \frac{{{p_0} + q + 1}}{2}}}\exp \left\{ { - \frac{1}{2}{\rm{tr}}\left( {{\mathit{\boldsymbol{D}}^{ - 1}}{\mathit{\boldsymbol{R}}_0}} \right)} \right\} $

又由b_i|D~N_q(0, D)可得

$ p\left( {{\mathit{\boldsymbol{b}}_i}\left| \mathit{\boldsymbol{D}} \right.} \right) \propto {\left| \mathit{\boldsymbol{D}} \right|^{ - \frac{m}{2}}}\exp \left\{ { - \frac{1}{2}\sum\limits_{i = 1}^m {\mathit{\boldsymbol{b}}_i^{\rm{T}}{\mathit{\boldsymbol{D}}^{ - 1}}{\mathit{\boldsymbol{b}}_i}} } \right\} $

则D的后验分布为

$ \begin{array}{l} \;\;\;\;{\rm{ \mathsf{ π} }}\left( {\mathit{\boldsymbol{D}}\left| \mathit{\boldsymbol{b}} \right.} \right) \propto {\left| \mathit{\boldsymbol{D}} \right|^{ - \frac{m}{2}}}\exp \left\{ { - \frac{1}{2}\sum\limits_{i = 1}^m {\mathit{\boldsymbol{b}}_i^{\rm{T}}{\mathit{\boldsymbol{D}}^{ - 1}}{\mathit{\boldsymbol{b}}_i}} } \right\}\\ {\left| \mathit{\boldsymbol{D}} \right|^{ - \frac{{{p_0} + q + 1}}{2}}}\exp \left\{ { - \frac{1}{2}{\rm{tr}}\left( {{\mathit{\boldsymbol{D}}^{ - 1}}{\mathit{\boldsymbol{R}}_0}} \right)} \right\} \propto {\left| \mathit{\boldsymbol{D}} \right|^{ - \frac{{m + {p_0} + q + 1}}{2}}}\exp \\ \left\{ { - \frac{1}{2}{\rm{tr}}\left( {{\mathit{\boldsymbol{D}}^{ - 1}}\left( {{\mathit{\boldsymbol{R}}_0} + \sum\limits_{i = 1}^m {{\mathit{\boldsymbol{b}}_i}\mathit{\boldsymbol{b}}_i^{\rm{T}}} } \right)} \right)} \right\} \end{array} $

所以D的条件后验分布为

$ \mathit{\boldsymbol{D}}{\left| \mathit{\boldsymbol{b}} \right._i} \sim I{W_q}\left( {m + {\rho _0},{\mathit{\boldsymbol{R}}_0} + \sum\limits_{i = 1}^m {{\mathit{\boldsymbol{b}}_i}\mathit{\boldsymbol{b}}_i^{\rm{T}}} } \right) $

(10)

3.3.2 τ的条件分布

对于τ，同样考虑服从先验分布IW(ρ₁, R₁)(ρ₁和R₁为已知的超参数)，则

$ p\left( \tau \right) \propto {\left| \tau \right|^{ - \frac{{{\rho _1} + 1 + 1}}{2}}}\exp \left\{ { - \frac{1}{2}{\rm{tr}}\left( {{\tau ^{ - 1}}{\mathit{\boldsymbol{R}}_1}} \right)} \right\} $

又根据a|τ~N_r－2(0, τI)，则

$ p\left( {\mathit{\boldsymbol{a}}\left| \tau \right.} \right) \propto {\left| \tau \right|^{ - \frac{{r - 2}}{2}}}\exp \left\{ { - \frac{1}{2}{\mathit{\boldsymbol{a}}^{\rm{T}}}{{\left( {\tau \mathit{\boldsymbol{I}}} \right)}^{ - 1}}\mathit{\boldsymbol{a}}} \right\} $

通过推导可得τ的后验分布为

$ \begin{array}{l} \;\;\;\;\;\;\;{\rm{ \mathsf{ π} }}\left( {\tau \left| \mathit{\boldsymbol{a}} \right.} \right) \propto {\left| \tau \right|^{ - \frac{{r - 2 + {\rho _1} + 1 + 1}}{2}}}\exp \left\{ { - \frac{1}{2}{\rm{tr}}\left( {{\tau ^{ - 1}}\left( {{\mathit{\boldsymbol{R}}_1} + } \right.} \right.} \right.\\ \left. {\left. {\left. {{\mathit{\boldsymbol{a}}^{\rm{T}}}\mathit{\boldsymbol{a}}} \right)} \right)} \right\} \end{array} $

所以τ的条件后验分布为

$ \tau \left| \mathit{\boldsymbol{a}} \right. \sim IW\left( {r - 2 + {\rho _1},{\mathit{\boldsymbol{R}}_1} + {\mathit{\boldsymbol{a}}^{\rm{T}}}\mathit{\boldsymbol{a}}} \right) $

(11)

3.4 随机误差方差项的条件后验分布

对于随机误差的方差σ²，采用扩散先验分布，即π(σ²)~1/σ²，则σ²的条件后验分布为

$ \begin{array}{l} {\rm{ \mathsf{ π} }}\left( {{\sigma ^2}\left| {{\mathit{\boldsymbol{\beta }}_ * },\mathit{\boldsymbol{b}},\mathit{\boldsymbol{a}},\mathit{\boldsymbol{Y}}} \right.} \right) \propto {\left( {{\sigma ^2}} \right)^{ - \left( {\frac{n}{2} + 1} \right)}}\exp \left\{ { - \frac{1}{{{\sigma ^2}}}} \right.\\ \left. {\left( {\frac{1}{2}{\mathit{\boldsymbol{u}}^{\rm{T}}}\mathit{\boldsymbol{u}}} \right)} \right\} \end{array} $

其中u=Y－X_*β_*－Z⁽¹⁾b－Z⁽²⁾a。

不难发现上述结果是逆Gamma分布(简记为IG)的核函数，所以σ²的后验分布为

$ {\sigma ^2}\left| {\mathit{\boldsymbol{Y}},{\mathit{\boldsymbol{\beta }}_ * },\mathit{\boldsymbol{b}},\mathit{\boldsymbol{a}}} \right. \sim IG\left( {\frac{n}{2},\frac{1}{2}{\mathit{\boldsymbol{u}}^{\rm{T}}}\mathit{\boldsymbol{u}}} \right) $

(12)

注意，逆Gamma分布是一个易于处理的后验分布，可以有效简化计算步骤。

3.5 贝叶斯推断

记{(β^(φ), D^(φ), τ^(φ), σ^2(φ)):φ=l, …, M}是利用MCMC方法生成的{β, σ², D, τ}的观测值。同时，令L=M－l，l表示迭代过程中的迭代次数，则β, D, τ, σ²的贝叶斯估计分别为

$ \begin{array}{l} \mathit{\boldsymbol{\hat \beta }} = \frac{1}{L}\sum\limits_{t = l}^M {{\mathit{\boldsymbol{\beta }}^{\left( \mathit{\boldsymbol{\varphi }} \right)}}} ,\mathit{\boldsymbol{\hat D}} = \frac{1}{L}\sum\limits_{t = l}^M {{\mathit{\boldsymbol{D}}^{\left( \mathit{\boldsymbol{\varphi }} \right)}}} \\ \hat \tau = \frac{1}{L}\sum\limits_{t = l}^M {{\tau ^{\left( \mathit{\boldsymbol{\varphi }} \right)}}} ,{{\hat \sigma }^2} = \frac{1}{L}\sum\limits_{t = l}^M {{{\left( {{\sigma ^2}} \right)}^{\left( \mathit{\boldsymbol{\varphi }} \right)}}} \end{array} $

此外，通过{β^(φ):φ=l, …, M}也可以得到后验协方差Var(β|Y, X, Z)的估计

$ \widehat {Var}\left( {\mathit{\boldsymbol{\beta }}\left| {\mathit{\boldsymbol{Y}},\mathit{\boldsymbol{X}},\mathit{\boldsymbol{Z}}} \right.} \right) = \frac{1}{{L - 1}}\sum\limits_{t = m}^M {\left( {{\mathit{\boldsymbol{\beta }}^{\left( \mathit{\boldsymbol{\varphi }} \right)}} - \mathit{\boldsymbol{\hat \beta }}} \right){{\left( {{\mathit{\boldsymbol{\beta }}^{\left( \mathit{\boldsymbol{\varphi }} \right)}} - \mathit{\boldsymbol{\hat \beta }}} \right)}^{\rm{T}}}} 。$

4 实证分析

将50家房地产上市公司每年的净资产收益率作为响应变量Y，以流动资产率(X₁)、资产负债率(X₂)、总资产的对数(X₃)、主营业务收入的对数(X₄)、主营业务收入增长率(X₅)作为解释变量建立半参数部分线性混合效应模型：

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;{Y_{ij}} = {\beta _1}{X_{1ij}} + {\beta _2}{X_{2ij}} + {\beta _3}{X_{3ij}} + {\beta _4}{X_{4ij}} + {\beta _5}{X_{5ij}} + \\ f\left( {{t_{ij}}} \right) + {b_i} + {\varepsilon _{ij}}\left( {i = 1, \cdots ,50,j = 1, \cdots ,10} \right) \end{array} $

(13)

其中，f(t)是一个光滑函数；b_i是相互独立的随机截距，服从正态分布N(0, θ)。

将数据代入模型(13)，通过贝叶斯分析和MCMC方法进行模拟可得模型(13) 的参数估计，具体结果见表 2、3。

由表 2可以看出，代表公司资产结构的流动资产率的系数0.1466为正值，说明流动资产率正面影响房地产上市公司的盈利能力，即公司的流动资产率越高，其盈利能力越好，这与房地产行业盈利能力强的公司都保持着良好的流动资产率的情况是相符的；资产负债率的系数为0.0246，表明资产负债率可以促进房地产公司的盈利能力，即资产负债率越高公司盈利能力越强，这与现状也是一致的。总资产对数的系数为0.0412，说明公司总资产越大，净资产收益率越高，即房地产公司规模越大，盈利能力越强，这与实际相符；主营业务收入的系数为0.0003，说明主营业务收入的与数对房地产上市公司的盈利能力呈正相关，较高的主营业务收入体现着房地产公司较好的盈利能力；代表公司发展潜力的主营业务收入增长率与净资产收益率存在正相关关系，说明具有发展潜力的公司受市场和投资者的认可，市场前景好，因此利润率相对较高。

5 结束语

本文根据贝叶斯理论对半参数部分线性混合效应模型进行了统计推断，可以发现模型固定效应、随机效应及其方差以及随机误差在先验分布下的后验分布均为标准分布，有效简便了模型估计的统计计算。同时，纵向数据下利用半参数部分线性混合效应模型对房地产上市公司盈利能力相关影响因素的实证结果表明这些因素对于房地产公司盈利能力有着不同程度的影响，实证分析结论也与房地产行业的实际情况相符，说明所提模型在实际应用中是可行的。但是在实际情况中，房地产上市公司的盈利能力不仅受到文中所述的解释变量影响, 还受到国家政策、宏观经济状况等外部因素的影响，从而造成模型的精确度降低，如果可以量化这些因素并纳入模型中将会有更加精准的结果，这也是本文后续进一步研究的内容。