殷亮，王維國

（1.大連民族大學(xué)理學(xué)院與預(yù)科教育學(xué)院，遼寧大連116650；2.東北財經(jīng)大學(xué)經(jīng)濟學(xué)院，遼寧大連116025）

0 引言

在實證分析中，當假設(shè)被解釋變量的分布為連續(xù)型正態(tài)分布的情況時，得到經(jīng)典的線性的計量經(jīng)濟模型，此時模型的預(yù)測精度很高。然而，Escribano[1]表明線性的假定過于嚴格，從大量的經(jīng)濟理論和經(jīng)濟實證中發(fā)現(xiàn)，變量之間的關(guān)系不都是線性的，Gonzalo等[2]指出，使用經(jīng)典計量方法研究變量之間的非線性關(guān)系時，往往得出錯誤的結(jié)論。因此，基于方法論的拓展和實證分析的需要，學(xué)者們將非線性引入到計量模型當中，不僅為計量模型的發(fā)展提供了一個方向，也成為計量經(jīng)濟學(xué)的前沿?zé)狳c問題。

有關(guān)非線性回歸模型的研究，很多成果可以借鑒。從McCullagh和Nelder[3]出版了關(guān)于廣義線性模型（Generalized Linear Model,GLM）的著作以后，學(xué)者們通過各種變換構(gòu)造廣義線性模型對離散的數(shù)據(jù)進行分析，為非線性模型的發(fā)展打下了良好的基礎(chǔ)。之后，GLM推廣到廣義線性混合效應(yīng)模型（Generalized Linear Mixed Model,GLMM），它既可以對來自非正態(tài)分布的數(shù)據(jù)進行擬合，同時對測量數(shù)據(jù)和縱向數(shù)據(jù)模型的內(nèi)在相關(guān)性也可以通過隨機效應(yīng)來刻畫，有關(guān)廣義線性模型的估計方法，McCulloch[4]采用了Monte Carlo Expectation Maximization（MCEM）方法去分析廣義線性混合效應(yīng)模型和廣義線性固定效應(yīng)。

半?yún)?shù)廣義線性（SGLM）模型是非線性回歸的另一個研究方向，半?yún)?shù)模型的參數(shù)部分表示變量間比較確定的線性關(guān)系，非參數(shù)部分表示不確定的部分。同時此模型假設(shè)被解釋變量的條件概率分布屬于指數(shù)分布族，為離散數(shù)據(jù)及屬性數(shù)據(jù)的分析提供了較好的研究方法。

本文結(jié)合Lombardia等[5]的廣義半?yún)?shù)趨勢混合面板模型和Hunsberger[6]的廣義半?yún)?shù)時間模型，建立了廣義半?yún)?shù)趨勢混合效應(yīng)面板模型，其解釋變量包含非線性時間趨勢，為減少其對模型參數(shù)估計的影響，結(jié)合Robinson[7]和Lombardia等[5]的方法，給出模型的加權(quán)極大似然估計量及識別程序。

1 廣義半?yún)?shù)混合效應(yīng)趨勢面板模型

廣義半?yún)?shù)趨勢混合效應(yīng)模型一方面考慮被解釋變量來自指數(shù)分布族，另一方面通過假設(shè)解釋變量包含非線性時間趨勢，是半?yún)?shù)廣義線性模型的特例，也是半?yún)?shù)趨勢面板模型的推廣。其形式如下：

其中，i，t表示橫截面單元和時間維度，為了研究方便，本文的Yit是一維被解釋變量，Xit是已知解釋變量，β(...βd)T是d維未知參數(shù)向量，ft=f(t)，gt=g(t)是未知非線性時間趨勢函數(shù)。嚴格單調(diào)，充分光滑的函數(shù)G稱為聯(lián)系函數(shù)，給定(X，t，αi)時，Y的條件密度所對應(yīng)的σ測度設(shè)為指數(shù)族形式：

αi，xi是隨機效應(yīng)，這里假定αi～N(0，)，xi～N(0，，vit是平穩(wěn)的殘差序列，vit～N(0，并且限制eit～N(0，表示模型（1）殘差向量。同時假定：

在廣義半?yún)?shù)趨勢混合效應(yīng)面板模型中，Xit通過ηit=β+ft+αi的半?yún)?shù)形式與Yit建立聯(lián)系，且半?yún)?shù)ηit=β+ft+αi通過G與Yit的數(shù)學(xué)期望E(Yit|Xit)相連接，即E(Yit|Xit)=G-1(β+ft+αi)，當G-1為指示函數(shù)時，則模型（1）和模型（2）變?yōu)榘雲(yún)?shù)趨勢隨機效應(yīng)面板模型。

2 廣義半?yún)?shù)趨勢面板模型的加權(quán)似然估計

廣義半?yún)?shù)趨勢混合效應(yīng)面板模型估計的主要問題是識別被解釋變量Yit和解釋變量Xit的關(guān)系，本文的模型估計方法是結(jié)合Robinson[7]的非參數(shù)趨勢面板模型的估計方法和加權(quán)極大似然估計方法[5]，涉及到極大似然函數(shù)的構(gòu)建問題。因為指數(shù)分布族內(nèi)包含多種分布，被解釋變量的不同分布導(dǎo)致了極大似然函數(shù)的不同，所以本文分兩種情況寫出模型的估計方法，第一種是被解釋變量來自一般指數(shù)分布族，第二種是被解釋變量來自指數(shù)分布族的一種分布——指數(shù)分布。

2.1 被解釋變量來自指數(shù)分布族

為了表述方便，對模型中變量表示做如下規(guī)定：Y∈R，X∈Rp，t∈R，考慮αi是隨機影響因素，假設(shè)αi～N(0，σu2)，進一步假設(shè)殘差項=，i=1，...，N，t=1，...，T，令θ=(，δ=(β，θ)，那么對于樣本Yit，Xit，t，(i=1，...,N，t=1，...,T)，Yit在Xit和t下的條件概率密度函數(shù)為：

其中，θit=Xitβ+ft+αi，Xβ是參數(shù)部分，ft=f(t)是非參部分，是充分光滑的函數(shù)，這里，Xit=gt+xi+vit，gt=g(t)也是充分光滑的函數(shù)，E(vit)=0p，E()=。

其中，qt=gtβ+ft是依賴于t的非參數(shù)部分，δ(β，αi)是依賴β，αi的參數(shù)部分，廣義半?yún)?shù)趨勢混合效應(yīng)面板模型（1）和模型（2）的主要結(jié)果是通過估計加權(quán)極大似然函數(shù)得到的β，qt=q(t)和αi。

考慮解釋變量的非線性時間趨勢會對模型的參數(shù)估計產(chǎn)生影響，本文提出兩步驟對廣義半?yún)?shù)趨勢混合效應(yīng)面板模型進行估計。

第一步：采用Robinson[7]提出的非參數(shù)趨勢面板模型的估計方法，對解釋變量Xit進行估計，得到個體效應(yīng)及非線性時間趨勢g=g(t)。即估計模型：

在式（7）中的非線性時間趨勢gt，反映解釋變量受時間的影響而變化的情況，是非線性的不確定的函數(shù)關(guān)系，（因為前文設(shè)定xi～N(0p，，則E(Xit)=gt），將其代替Xit代入θit中，得到：

其中，依賴β，αi的參數(shù)部分是δ(β，αi)，依賴于t的非參數(shù)部分qt=ft，通過估計加權(quán)極大似然函數(shù)得到β和q(t)=f(t)為廣義半?yún)?shù)趨勢混合面板模型（1）和模型（2）的主要結(jié)果。

第二步：采用加權(quán)極大似然估計方法估計模型的參數(shù)和非參數(shù)估計量。

模型的估計問題主要是識別參數(shù)，即參數(shù)部分β和αi，非參數(shù)部分qt=ft是本文要估計的，由后驗分布可知：

根據(jù)極大似然估計的思想，估計的參數(shù)β和αi能使得式（9）值達到最大，只需保證分子達到最大即可。定義：

令L(Y，α;q，δ)=L1(Y;q，δ)+L2(α，，為了估計非參數(shù)部分qt=ft，本文固定t0點，計算t0點的Y的對數(shù)條件密度函數(shù)：

則加權(quán)極大似然函數(shù)為：

其中，Kh(t0)=K(t-t0h)是帶寬為h的核函數(shù)，滿足在[-1，1]上的lipschitz對稱概率密度函數(shù)。

為表示方便，計算：

此式稱為加權(quán)極大似然函數(shù)，其解為加權(quán)極大似然估計量。

2.2 被解釋變量服從指數(shù)分布

為了使上面的估計方法具體化，選擇指數(shù)分布作為指數(shù)分布族的一類分布，當被解釋變量服從指數(shù)分布時：其概率密度函數(shù)為λe-λy，其中，λ=1 E(Y)，首先給出廣義半?yún)?shù)趨勢混合面板模型的形式，令模型的半?yún)?shù)部分為ηit=Xitβ+ft+αi，被解釋變量的數(shù)學(xué)期望為E(Yit|Xit)=μit，二者通過G建立聯(lián)系，即G(μit)=ηit=Xitβ+ft+αi，其中，嚴格單調(diào)充分光滑的G函數(shù)稱為聯(lián)系函數(shù)。h=G-1，即μit=h(ηit)=1/ηit，或G(μit)=1 μit。此時模型（1）和模型（2）變?yōu)椋?/p>

當被解釋變量服從指數(shù)分布時，借鑒上面給出的模型估計方法的一般情況，得到似然函數(shù)形式如下：

則其加權(quán)似然函數(shù)為：

將采用Robinson的方法對解釋變量的非線性時間趨勢進行分離，將其代替解釋變量代入到加權(quán)似然函數(shù)當中，得到：

其極大似然估計量為?WL(Y;q，δ)/?β=0，求出被解釋變量服從指數(shù)分布時的廣義半?yún)?shù)趨勢混合面板模型的參數(shù)估計量和非參數(shù)估計量。

3 最優(yōu)化算法

半?yún)?shù)趨勢隨機效應(yīng)面板模型的參數(shù)和非參數(shù)估計量是加權(quán)似然函數(shù)的最優(yōu)解，然而它是非顯示解。本文參考Lombardia和Sperlich[5]提出的算法，對其進行估計。

步驟1：對于固定時間的t0，當給定參數(shù)δ（β和αi）時，可以通過下式估計非參數(shù)部分q(t0)：

步驟2：結(jié)合式（20）結(jié)果，可以通過下式估計估計β和αi：

4 蒙特卡羅數(shù)值模擬

本文與Lombardia和Sperlich[5]不同的地方是除了考慮被解釋變量是包含時間趨勢的非平穩(wěn)變量，解釋變量也是如此。如果直接參考Lombardia和Sperlich[5]提出的估計方法估計模型的參數(shù)和非參數(shù)估計量，解釋變量的時間趨勢可能對估計方法產(chǎn)生影響。為了驗證猜想，本文首先考慮解釋變量是不含時間趨勢的平穩(wěn)變量，生成數(shù)據(jù)，對Lombardia和Sperlich[5]的估計方法進行驗證，然后，考慮解釋變量包含非線性時間趨勢，生成數(shù)據(jù)，對本文提出的估計方法和Lombardia和Sperlich[5]提出的估計方法進行對比分析。

4.1 生成數(shù)據(jù)

為了研究方便，在此部分假定Xit為一維變量，則廣義半?yún)?shù)趨勢混合效應(yīng)面板模型形式如下：

其中，Yit來自指數(shù)分布族，因為指數(shù)分布族中的分布的形式過多，同時與上面的估計方法相對應(yīng)，所以本文基于指數(shù)分布族中的一類分布——指數(shù)分布給出廣義半?yún)?shù)趨勢混合效應(yīng)面板模型的加權(quán)極大似然估計方法，在此部分，選用指數(shù)分布作為指數(shù)分布族的特殊分布給出數(shù)據(jù)生成過程。

（1）非線性時間趨勢：本文中非線性時間趨勢是被解釋變量和解釋變量非平穩(wěn)的主要原因，f(u)=2u3+u，g(u)=2sin(πu)，k(u)=2πu為非線性的時間趨勢，反映了政策、法律、法規(guī)等隨機因素對變量的影響。

（2）隨機效應(yīng)：面板數(shù)據(jù)中，對于特定的個體而言，隨機效應(yīng)是不隨著時間發(fā)生改變的影響因素。αit在N(0，0.25)的正態(tài)分布中隨機抽取，xit在N(0，1)的正態(tài)分布中隨機抽取。

（3）殘差項：eit在N(0，0.25)的正態(tài)分布中提取，序列獨立的序列。vit在N(0，0.6)的正態(tài)分布中提取。

（4）解釋變量與被解釋變量：本文的特別之處就在于考慮了解釋變量是因包含非線性時間趨勢而非平穩(wěn)的變量。在數(shù)據(jù)生成部分，按照兩種情況來生成解釋變量Xit，一種是平穩(wěn)變量，服從正態(tài)分布N(1，0.6)，不包含個體效應(yīng)和時間趨勢，另一種是包含時間趨勢的非平穩(wěn)變量，按照Xit=vit+xi+gt生成。在接下來的蒙特卡羅模擬中，考慮解釋變量的時間趨勢，加入g(u)=2sin(πu)及k(u)=2 πu，被解釋變量Yit按照指數(shù)分布來生成，即：

其中，Yit是以E(Yit)=1/Xitβ+αit+ft為期望的指數(shù)分布生成，ft的形式是按照g(u)及k(u)來構(gòu)造。

4.2 廣義半?yún)?shù)模型蒙特卡羅數(shù)值模擬

按三部分進行模擬實驗：第一，是驗證Lombardia和Sperlich[5]提出的估計方法，第二，從標準差的角度對比分析本文的估計方法和Lombardia和Sperlich[5]的估計方法，第三，從均方誤差的角度對比分析本文的估計方法和Lombardia和Sperlich[5]的估計方法。

（1）首先驗證Lombardia和Sperlich[5]的估計方法，即對廣義半?yún)?shù)混合效應(yīng)面板模型的參數(shù)估計量和非參數(shù)估計量進行估計，分析其有限樣本性質(zhì)。

為了使實驗更一般化，在蒙特卡羅模擬部分，對模型參數(shù)進行兩種假定，即假設(shè)模型參數(shù)分別為β=3和β=5的情況及隨機效應(yīng)的分布方差為=0.0625和=0.09的情況。

考慮到時間維度或截面維度的變化會對廣義半?yún)?shù)混合效應(yīng)面板模型的參數(shù)估計產(chǎn)生影響，首先固定時間維度不變，增加截面維度。觀察模型參數(shù)估計變化情況。然后固定截面維度不變，增加時間維度，觀察模型參數(shù)估計變化情況。

取T=5保持不變，N增加（N=10、N=15、N=20、N=25、N=30），進行500模擬，計算500次模擬參數(shù)估計值的平均值和標準差。為了另取一時間維度不變，增加截面維度，即取T=10保持不變，N增加的情況（N=10、N=15、N=20、N=25、N=30）的情況。進行500次模擬，計算500次模擬參數(shù)估計值的平均值和標準差。結(jié)果如表1和表2所示。

表1 廣義半?yún)?shù)混合效應(yīng)面板模型參數(shù)估計結(jié)果（T=5）

表2 廣義半?yún)?shù)混合效應(yīng)面板模型參數(shù)估計結(jié)果（T=10）

表1和表2列出了廣義半?yún)?shù)混合效應(yīng)面板模型的參數(shù)估計情況，隨著樣本容量的變化，其有限樣本性質(zhì)如下：

當保持時間維度T不變時，增加個體數(shù)N，廣義半?yún)?shù)混合效應(yīng)面板模型的參數(shù)估計結(jié)果有如下特征。第一，模型的參數(shù)估計值的均值變化比較穩(wěn)定。比如，在假定β=3的情況下，保持時間維度T=10不變，N=15時參數(shù)估計值最接近真實值。第二，除個別情況下（T=5，N=25；T=10，N=15），隨著截面?zhèn)€數(shù)N的增加，參數(shù)估計值的標準差有逐漸減小的趨勢。比如，在β=5，T=5的情況下，隨著個體數(shù)N的增加（N=10、N=15、N=20、N=25、N=30），參數(shù)估計值的標準差在逐漸減小。第三，隨機效應(yīng)分布的方差估計的結(jié)果較為理想，在β=3和β=5的假定下，其方差估計的結(jié)果較為準確。

當保持截面維度N不變，增加時期數(shù)T（從T=5到T=10）時，廣義半?yún)?shù)混合效應(yīng)面板模型的參數(shù)估計結(jié)果有如下特征。第一，模型的參數(shù)估計值的均值變化的穩(wěn)定性提高，即在T=5下的情況下，β=3的均值變化幅度較大，而在T=10時，其均值變化比較穩(wěn)定。第二，模型參數(shù)估計值的標準差在逐漸減小，即在T=5時，β=3的參數(shù)估計值的標準差在0.9～1.9之間，而在T=10時，β=3的參數(shù)估計值的標準差在0.7～1.1之間。第三，估計的隨機效應(yīng)分布的方差結(jié)果較為理想，在β=3和β=5的假定情況下的假定情況下，模型的隨機效應(yīng)方差估計的結(jié)果都比較準確。

（2）表1和表2是在假設(shè)解釋變量為平穩(wěn)的變量的情況下得到的結(jié)論，驗證了廣義半?yún)?shù)混合面板模型的有限樣本性質(zhì)，下面考慮解釋變量為包含時間趨勢的非平穩(wěn)變量時，對比分析Lombardia和Sperlich[5]提出的廣義半?yún)?shù)混合效應(yīng)面板模型估計方法和本文提出的兩步驟模型估計方法。表3和表4表明Lombardia和Sperlich[5]提出的方法的估計結(jié)果。

表3 廣義半?yún)?shù)趨勢混合效應(yīng)面板模型參數(shù)估計結(jié)果（T=5）

表4 廣義半?yún)?shù)趨勢混合效應(yīng)面板模型參數(shù)估計結(jié)果（T=10）

從表3和表4的參數(shù)估計結(jié)果發(fā)現(xiàn)，當解釋變量是含有非線性時間趨勢的非平穩(wěn)變量時，采用Lombardia和Sperlich[5]提出的極大似然估計方法估計廣義半?yún)?shù)趨勢混合面板模型的參數(shù)，參數(shù)估計結(jié)果較好。第一，隨著樣本容量N、T的變化，參數(shù)估計值變化比較穩(wěn)定。第二，除個別情況下，隨著樣本容量N、T的變化，參數(shù)估計值的標準差有逐漸減小的趨勢。第三，隨機效應(yīng)分布的方差估計結(jié)果較為理想。

盡管Lombardia和Sperlich[5]提出的極大似然估計方法依然可以對廣義半?yún)?shù)趨勢混合面板模型參數(shù)進行估計，然而，由于解釋變量從平穩(wěn)變量變成包含時間趨勢的非平穩(wěn)變量，Lombardia和Sperlich[5]提出的極大似然估計方法是否是最佳的？是本文需要研究的問題。接下來，采用本文提出的兩步驟估計方法對廣義半?yún)?shù)趨勢混合面板模型進行估計，即首先采用Robinson(2012)[7]提出的非參數(shù)趨勢面板模型估計方法對解釋變量的時間趨勢進行分離，然后用其代替解釋變量代入到廣義半?yún)?shù)混合效應(yīng)面板模型當中進行加權(quán)極大似然估計，觀察估計結(jié)果。如表5和下頁表6所示。

表5 廣義半?yún)?shù)趨勢混合效應(yīng)面板模型參數(shù)估計結(jié)果（T=5調(diào)整）

表6 廣義半?yún)?shù)趨勢混合效應(yīng)面板模型參數(shù)估計結(jié)果（T=10調(diào)整）

從表5和表6與表3和表4的估計結(jié)果發(fā)現(xiàn)：當解釋變量是包含時間趨勢的非平穩(wěn)變量時，本文提出的兩步驟估計方法與Lombardia和Sperlich[5]提出的加權(quán)極大似然估計方法對比結(jié)果如下：從共同點來說，除個別情況下，兩種參數(shù)估計結(jié)果都保持了參數(shù)估計值變化比較穩(wěn)定，估計值的標準差隨著樣本容量的增大而減小的特征。從不同點來說，若直接采用廣義半?yún)?shù)模型進行分析采用加權(quán)極大似然法估計，參數(shù)估計值均值估計較穩(wěn)定的情況下，參數(shù)估計值的標準差結(jié)果好于先對解釋變量進行分離然后再進行加權(quán)極大似然估計的估計結(jié)果，后者均值的標準差在整體水平上有明顯的增大趨勢。

（3）從以上兩種方法的估計結(jié)果來看，參數(shù)估計值均值變化都較穩(wěn)定，本文提出的估計方法的參數(shù)估計值的標準差雖然在整體水平上略大于Lombardia和Sperlich[5]提出的方法，但是相差不大。僅從估計值的標準差大小來判斷兩種估計方法哪個更可靠，不夠客觀，下面計算兩種估計方法的估計值的均方誤差，觀察計算結(jié)果，進一步比較兩種估計方法的優(yōu)良。

對于參數(shù)β的估計值的均方誤差計算公式如下：

均方誤差刻畫的是估計值的損失函數(shù)，均方誤差越小，估計效果越好。為了清楚地對比兩種方法均方誤差的大小，下面列出β=3時，兩種估計方法的參數(shù)估計值的均方誤差，用L(2008)表示采用Lombardia和Sperlich[5]提出的加權(quán)極大似然估計方法，用L&R(2012)表示本文提出的估計方法，即先采用Roinson[7]對解釋變量X的時間趨勢進行分離，然后采用加權(quán)極大似然法進行估計。計算結(jié)果見表7。從表7的結(jié)果可以看出，本文提出的估計方法的參數(shù)估計值的均方誤差略小于Lombardia和Sperlich[5]提出的加權(quán)極大似然估計方法。

表7 參數(shù)估計值的均方誤差

5 結(jié)論

本文在Lombardia和Sperlich[5]和Hunsberger[6]的基礎(chǔ)上，建立了解釋變量為包含時間趨勢的廣義半?yún)?shù)趨勢混合效應(yīng)面板模型，并給出加權(quán)極大似然估計的有效估計量及識別程序。在仿真實驗部分中，當被解釋變量服從指數(shù)分布時，對于解釋變量為非平穩(wěn)的時間序列，本部分采用了Lombardia和Sperlich[5]提出的加權(quán)極大似然估計方法和本文提出的估計方法進行對比分析了兩種估計方法的估計結(jié)果。

兩種估計方法得到的參數(shù)估計結(jié)果較為理想，都保持了參數(shù)估計值變化穩(wěn)定，估計值的標準差隨著樣本容量增加而減小的特征。從參數(shù)估計值標準差上來看，本文提出的估計方法的參數(shù)估計值的標準差略大于Lombardia和Sperlich[5]提出的加權(quán)極大似然估計方法的標準差，從均方誤差來看，本文提出的估計方法的均方誤差略小于加權(quán)極大似然估計方法的均方誤差。

[1] Escribano A.Nonlinear Error Correction:The Case of Money Demand in the UK(1878-2000)[J].Macroeconomic Dynamics，2004，(8).

[2] Gonzalo J,Pitarakis J Y.Threshold Effects in Cointegrating Relationships[J].Oxford Bulletin of Economics and Statistics,2006，(68).

[3] McCullagh P,Nelder J A.Generalized Linear Models(second edition)[M].London:Chapman and Hall,1989.

[4] McCulloch C E,Natarajan R.A Note on the Existence of the Posterior Distribution for a Class of Mixed Models for Binomial Responses[J].Bimometrika，1995，(82).

[5] Lombardia M J,Sperlich S.Semi-parametric Inference in Generalized Mixed Effects Models Working Paper[J].Journal of the Royal Statisti?cal Society,2008,70(5).

[6] Hunsberger S.Semi-parametric Regression in Likelihood-Bsed Mod?els[J].Journal of the American Stastical Association,1994,89,(428).

[7] Robinson P M.Non-parametric Trending Regression With Cross-sec?tional Dependence[J].Journal of Econometrics,2012,(169).