第 6 章 修正
自從 Petersen (2009) 一文，Panel data 具異質變異和序列相關的穩健標準差計算的修正式， 就變得很重要。尤其在財務與會計學門使用panel data的研究很多，但是， 多半沒有處理穩健共變異問題。基本上的問題，就是有非同質變異時，卻沒有修正的話，會低估變異數，導致係數會相當顯著。相關經典論文，如 Stock and Watson (2008) 和 Thompson (2011) 有都提出一些處理方式。本章分兩節介紹這個問題：第1節是修正係數估計式，第2節是指修正變異數估計式。
6.1 具序列相關修正之模型估計
前章第4節介紹了序列相關的檢定，有的序列相關存在與否，和設定的推論有關。如同一般時間序列迴歸問題，當迴歸殘差有序列相關時，會導致估計的參數不是最佳的(Non-Optimality)。如同前面的實證範例，我們發現，不管是單維或是雙維，固定效果和隨機效果對x3係數的估計結果，相差甚遠，除了一正一負，顯著性也極端不同。隨機效果就很顯著，固定效果就不顯著。這問題該怎麼辦？
根據我們在Hausman Test的說明，如果只是內生性問題，兩者的估計參數結果會很近，差別是變異數的效率問題。因此，這很有可能是殘差有序列相關，但是，未受修正。因此，我們必須使用類似標準時間序列使用的Cochrane-Orcutt 演算來修正這個問題。在Panel Data，Baltagi and Li (1991) 提出兩步驟修正方法：
第1步先用LS估計原式，在計算殘差序列相關係數；若其為一階，就是AR(1)係數；
第2步則利用了Prais-Winsten5轉換矩陣(簡稱PW)來修正估計式。在Panel Data分析中，PW轉換比 Cochrane and Orcutt (1949) 要有效率，因為PW轉換利用了每一個時間點t所有的觀察值，但是，Cochrane-Orcutt的修正方法，需要扣掉一筆資料。在T相對很少的追蹤資料結構，Cochrane-Orcutt的方法顯的昂貴。承前面的符號，panel data 迴歸殘差的AR(1)迴歸如下：
\[
\tag{6.1}
{{\varepsilon }_{it}}=\rho {{\varepsilon }_{it-1}}+{{\eta }_{it}}
\]
上式中：
\(\left| \rho \right|<1, \;\; {{\eta }_{it}}\sim iid(0,\sigma _{\eta }^{2})\)
且
\({{\varepsilon }_{it}}\sim iid(0,\frac{\sigma _{\eta }^{2}}{1-{{\rho }^{2}}})\)
我們利用 \(i=1\) ，也就是時間序列的例子，說明此變異數的代數結構如下：
\[
\begin{align}
& \sigma _{\varepsilon }^{2}=E{{(\rho {{\varepsilon }_{t-1}}+{{\eta }_{t}})}^{2}}-{{(E(\rho {{\varepsilon }_{t-1}}+{{\eta }_{t}}))}^{2}} \\
& \because E(\rho {{\varepsilon }_{t-1}}+{{\eta }_{t}})=0 \\
& \therefore \sigma _{\varepsilon }^{2}=E{{(\rho {{\varepsilon }_{t-1}}+{{\eta }_{t}})}^{2}}={{\rho }^{2}}E\varepsilon _{t-1}^{2}+2\rho E({{\varepsilon }_{t-1}}{{\eta }_{t}})+E\eta _{t}^{2} \\
& \Rightarrow \sigma _{\varepsilon }^{2}={{\rho }^{2}}\sigma _{\varepsilon }^{2}+\sigma _{\eta }^{2} \\
\end{align}
\]
故 \(\sigma _{\varepsilon }^{2}=\frac{\sigma _{\eta }^{2}}{1-{{\rho }^{2}}}\)
也就是說，如果殘差有序列相關 \(\rho \ne 0\) ，但是卻假設它沒有序列相關 \(\rho= 0\) ，會低估標準差。被低估的標準差，會產生估計係數的 \(t\) 檢定相對顯著的假象。嚴重時，就會做出錯誤的統計推論，以為資料解釋理論解釋的相當好。過去數十年，許多研究偏好使用pooled LS的顯著結果，來支持理論假設，其實是一種計量上的錯誤。所以，檢定後的修正是必然的。
因此，利用iterative procedure，一個一般化的自我共變異函數(Auto-Covariance Function) 為 \(\operatorname{cov}({{\varepsilon }_{t}},{{\varepsilon }_{t+h}})=\frac{{{\rho }^{h}}}{1-{{\rho }^{2}}}\) 故這個AR(1)的variance-covariance 矩陣 \(\Omega\) 可以寫成
\[
\Omega =\left[ \begin{matrix}
\frac{1}{1-{{\rho }^{2}}} & \frac{\rho }{1-{{\rho }^{2}}} & \frac{{{\rho }^{2}}}{1-{{\rho }^{2}}} & \cdots & \frac{{{\rho }^{T-1}}}{1-{{\rho }^{2}}} \\
\frac{\rho }{1-{{\rho }^{2}}} & \frac{1}{1-{{\rho }^{2}}} & \frac{\rho }{1-{{\rho }^{2}}} & \cdots & \frac{{{\rho }^{T-1}}}{1-{{\rho }^{2}}} \\
\frac{{{\rho }^{2}}}{1-{{\rho }^{2}}} & \frac{\rho }{1-{{\rho }^{2}}} & \frac{1}{1-{{\rho }^{2}}} & \cdots & \frac{{{\rho }^{T-3}}}{1-{{\rho }^{2}}} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
\frac{{{\rho }^{T-1}}}{1-{{\rho }^{2}}} & \frac{{{\rho }^{T-2}}}{1-{{\rho }^{2}}} & \frac{{{\rho }^{T-3}}}{1-{{\rho }^{2}}} & \cdots & \frac{1}{1-{{\rho }^{2}}} \\
\end{matrix} \right]
\]
在LS迴歸之中，將此共變異數矩陣的逆矩陣，內插入LS公式就得到修正AR(1)之PW的 GLS估計式 \({{({X}'{{\Omega }^{-1}}X)}^{-1}}{X}'{{\Omega }^{-1}}Y\)
承前，Baltagi and Li (1991) 兩步驟修正方法接著：
第1步。因為 \(\Omega\) 是對稱正定矩陣，故其逆矩陣可分割為\(\Omega^{-1}=C'C\)
利用PW提出的轉換矩陣 C，來轉換原殘差，將之轉換成無序列相關殘差。C矩陣如下
\[
\mathbf{C}=\left[ \begin{matrix}
\sqrt{1-{{\rho }^{2}}} & 0 & 0 & \cdots & \cdots & 0 & 0 & 0 \\
-\rho & 1 & \cdots & \cdots & \cdots & 0 & 0 & 0 \\
0 & -\rho & 1 & \cdots & \cdots & \cdots & \cdots & \vdots \\
\vdots & 0 & -\rho & 1 & \cdots & \cdots & \cdots & \vdots \\
\vdots & \cdots & 0 & -\rho & 1 & \cdots & \cdots & \vdots \\
\vdots & \cdots & \cdots & \cdots & \ddots & \ddots & 0 & 0 \\
\vdots & \cdots & \cdots & \cdots & 0 & -\rho & 1 & 0 \\
0 & \cdots & \cdots & \cdots & 0 & 0 & -\rho & 1 \\
\end{matrix} \right]
\]
轉換後的新殘差為
\[
\tag{6.2}
{{\xi }_{it}}=({{\mathbf{I}}_{\mathbf{N}}}\otimes \mathbf{C})u=({{\mathbf{I}}_{\mathbf{N}}}\otimes \mathbf{C})\mu +({{\mathbf{I}}_{\mathbf{N}}}\otimes \mathbf{C})\varepsilon
\]
令 \({{C}_{\tau T}}=(1-\rho )\tau _{T}^{\alpha }\) ，且 \(\tau _{T}^{\alpha }=(\alpha ,{{{\tau }'}_{T-1}}),\alpha =\sqrt{\frac{1+\rho }{1-\rho }}\) ，則(6.2)可以寫成
\[
\tag{6.3}
{{\xi }_{it}}=(1-\rho )({{\mathbf{I}}_{\mathbf{N}}}\otimes \tau _{T}^{\alpha })\mu +({{\mathbf{I}}_{\mathbf{N}}}\otimes \mathbf{C})\varepsilon
\]
因此，可以得到一個新的共變異數矩陣
\[
\tag{6.4}
\mathbf{\tilde{\Omega }}=E(\mathbf{\xi {\xi }'})=\sigma _{\mu}^{2}{{(1-\rho )}^{2}}({{\mathbf{I}}_{\mathbf{N}}}\otimes \tau_{T}^{\alpha}{\tau'}_{T}^{\alpha })
({{\mathbf{I}}_{\mathbf{N}}}\otimes {{\mathbf{I}}_{\mathbf{T}}})\sigma _{\varepsilon }^{2}
\]
第2步，將原始資料前乘矩陣 \({{\mathbf{\sigma }}_{\mathbf{\eta }}}\sqrt{{\mathbf{\tilde{\Omega }}}}\left( {{\mathbf{I}}_{\mathbf{N}}}\otimes \mathbf{C} \right)\) 。例如，被解釋變數向量y，轉換後變成 \({{\mathbf{\sigma }}_{\mathbf{\eta }}}\sqrt{{\mathbf{\tilde{\Omega }}}}\left( {{\mathbf{I}}_{\mathbf{N}}}\otimes \mathbf{C} \right)\mathbf{y}\) ，用同樣的方法轉換解釋變數。之後，再利用固定效果或隨機效果的估計式，去重新估計迴歸方程式。
這樣轉換後的迴歸殘差，就不具一階序列相關。需要進一步理論細節的讀者，或序列相關的結構是其他形式，例如，AR(2), MA(1), 或季節性結構AR(4)等等，請參考 Wooldridge (2010) 第10章或 Baltagi (2013) 第5章。
模組plm內建的估計方法，沒有對估計式做兩階段序列相關修正，這問題在R裡面很簡單。我們利用R內建的gls()函數即可，關鍵就是裡面的 corAR1()函數了。
固定效果 AR(1)修正的R 程式如下：
library(nlme)
gsp_feAR1=gls(as.formula(Eq),data=temp1, correlation = corAR1(0,form=~year|state))
summary(gsp_feAR1)$tTable## Value Std.Error t-value p-value
## (Intercept) 2.805534534 0.212847679 13.180950 4.424391e-36
## x1_hwy 0.070161815 0.041038954 1.709639 8.771603e-02
## x2_water 0.045945405 0.021187335 2.168532 3.040902e-02
## x3_other -0.007173820 0.028522088 -0.251518 8.014775e-01
## x4_private 0.066371301 0.021567665 3.077352 2.158893e-03
## x5_emp 0.879342712 0.031705093 27.735062 1.615147e-119
## x6_unemp -0.005256745 0.000739034 -7.112995 2.500988e-12correlation=corAR1(0,form=~year|state)這個函數的說明有兩點，瞭解後就更知道如何使用它。
第1.
corAR1(0,form=~year|state)省略了一個內定宣告，就是
corAR1(0,form=~year|state, fixed=FALSE)
前面的數值0，是說序列相關數值演算的起始值(starting values)，關鍵是後面的
fixed=FALSE
FALSE是說這個值不固定，須要估計。如果宣告TRUE就是不修正，也就是無序列相關。內定是需要估計，所以一般是省略了。
第2. 參數form=~year|state是說，這個殘差的序列相關，是如何修正，也就是由時間(year)去排列殘差，然後，grouping 是state，也就是panel data 的 \(N\) ：在每一個state內，year都沒有重複。
第3. 更廣泛的高階序列相關結構，是可以利用
corARMA(p=, q=,form=~year|state)
p和q分別是AR和MA的階次。所以，如果需要更高階的序列相關修正，例如，AR(2)，季節結構，就可以使用。
估計結果摘要如下：
## Generalized least squares fit by REML
## Model: as.formula(Eq)
## Data: temp1
## AIC BIC logLik
## -3694.67 -3652.408 1856.335
##
## Correlation Structure: AR(1)
## Formula: ~year | state
## Parameter estimate(s):
## Phi
## 0.9877442
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 2.8055345 0.21284768 13.180950 0.0000
## x1_hwy 0.0701618 0.04103895 1.709639 0.0877
## x2_water 0.0459454 0.02118734 2.168532 0.0304
## x3_other -0.0071738 0.02852209 -0.251518 0.8015
## x4_private 0.0663713 0.02156767 3.077352 0.0022
## x5_emp 0.8793427 0.03170509 27.735062 0.0000
## x6_unemp -0.0052567 0.00073903 -7.112995 0.0000
##
## Correlation:
## (Intr) x1_hwy x2_wtr x3_thr x4_prv x5_emp
## x1_hwy -0.736
## x2_water 0.105 -0.094
## x3_other 0.140 -0.404 -0.223
## x4_private -0.309 -0.103 -0.153 -0.128
## x5_emp 0.339 -0.309 -0.265 -0.291 -0.301
## x6_unemp 0.237 -0.165 -0.175 -0.257 -0.222 0.668
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.0756893 -0.6420634 -0.1824753 0.2992425 4.8399730
##
## Residual standard error: 0.1393526
## Degrees of freedom: 816 total; 809 residual
以上還有一些outputs，於此省略。檢視這個結果，phi=0.988就是殘差序列相關係數，可見得序列相關程度嚴重。如果不修正的話，結果會不可靠。當然，這問題也接近了非定態問題。後面我們可以再討論。
隨機效果 AR(1)修正的R 程式如下：
隨機效果修正AR(1)序列相關的結果，和固定效果就相當一致了。
## Linear mixed-effects model fit by REML
## Data: temp1
## AIC BIC logLik
## -3692.67 -3645.712 1856.335
##
## Random effects:
## Formula: ~1 | state
## (Intercept) Residual
## StdDev: 1.695642e-05 0.1393526
##
## Correlation Structure: AR(1)
## Formula: ~year | state
## Parameter estimate(s):
## Phi
## 0.9877442
## Fixed effects: as.formula(Eq)
## Value Std.Error DF t-value p-value
## (Intercept) 2.8055345 0.21284768 762 13.180950 0.0000
## x1_hwy 0.0701618 0.04103895 762 1.709639 0.0877
## x2_water 0.0459454 0.02118734 762 2.168532 0.0304
## x3_other -0.0071738 0.02852209 762 -0.251518 0.8015
## x4_private 0.0663713 0.02156767 762 3.077352 0.0022
## x5_emp 0.8793427 0.03170509 762 27.735062 0.0000
## x6_unemp -0.0052567 0.00073903 762 -7.112995 0.0000
## Correlation:
## (Intr) x1_hwy x2_wtr x3_thr x4_prv x5_emp
## x1_hwy -0.736
## x2_water 0.105 -0.094
## x3_other 0.140 -0.404 -0.223
## x4_private -0.309 -0.103 -0.153 -0.128
## x5_emp 0.339 -0.309 -0.265 -0.291 -0.301
## x6_unemp 0.237 -0.165 -0.175 -0.257 -0.222 0.668
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.0756893 -0.6420634 -0.1824753 0.2992425 4.8399730
##
## Number of Observations: 816
## Number of Groups: 48## Value Std.Error DF t-value p-value
## (Intercept) 2.805534532 0.212847678 762 13.180950 6.904951e-36
## x1_hwy 0.070161816 0.041038953 762 1.709639 8.773968e-02
## x2_water 0.045945406 0.021187335 762 2.168532 3.042698e-02
## x3_other -0.007173820 0.028522088 762 -0.251518 8.014814e-01
## x4_private 0.066371301 0.021567665 762 3.077352 2.163278e-03
## x5_emp 0.879342711 0.031705093 762 27.735062 1.360803e-117
## x6_unemp -0.005256745 0.000739034 762 -7.112995 2.620496e-12
我們也可以執行一個ANOVA分析的比較
## Model df AIC BIC logLik Test L.Ratio p-value
## gsp_reAR1 1 10 -3692.67 -3645.712 1856.335
## gsp_feAR1 2 9 -3694.67 -3652.408 1856.335 1 vs 2 6.6349e-08 0.9998
ANOVA 的p-value指出兩個幾乎沒有差別。
最後，我們再介紹如何比較LSDV之下的固定效果，和修正序列相關的結果。因為需要ANOVA比較，所以，plm的架構不能使用。故都一致使用gls()函數估計一個LSDV的固定效果，和前面的AR(1)修正結果比較。
固定效果 \(i.i.d.\) vs. 固定效果 AR(1)修正的R程式碼：
先寫新公式myFormula1。將共同截距去除，改以用state產生虛擬變數矩陣as.factor(state)。如果不習慣這樣處理公式，用人工key 也可以。
myFormula1=as.formula(paste(paste(Eq,"-1"),"as.factor(state)",sep="+"))
gsp_feIID=gls(myFormula1,data=temp1)
k=length(names(myData1)[4:9])
\(k\) 計算解釋變數個數。取出部份係數摘要。需要這樣做，是因為估計結果中有許多截距式的固定效果
## Value Std.Error t-value p-value
## x1_hwy 0.076752199 0.031243700 2.456566 1.424921e-02
## x2_water 0.078697075 0.015003077 5.245396 2.022587e-07
## x3_other -0.114778244 0.018147120 -6.324874 4.317272e-10
## x4_private 0.234986423 0.026214245 8.964074 2.382829e-18
## x5_emp 0.801172149 0.029756978 26.923841 1.005712e-112
## x6_unemp -0.005179978 0.000979674 -5.287450 1.621685e-07
這個做法使用了比較進階的公式建立方法，我也將state轉成因子，這樣才能產生within group的效果。我們看看ANOVA比較的結果
## Warning in nlme::anova.lme(object = gsp_feIID, gsp_feAR1): fitted objects with
## different fixed effects. REML comparisons are not meaningful.## Model df AIC BIC logLik Test L.Ratio p-value
## gsp_feIID 1 55 -2610.425 -2355.448 1360.212
## gsp_feAR1 2 9 -3694.670 -3652.408 1856.335 1 vs 2 992.2455 <.0001
看最後面的值，可以知道序列相關修正是必要的。同時也知道，修正後，固定效果和隨機效果的結果會趨於一致。前章發現的問題，至此就沒有了。所以，序列相關的檢定和修正，就時十分重要。但是，本書用這個高序列相關的極端例子，只是為了凸顯修正的重要，以及非iid的問題。一般的資料，如果序列相關不是很嚴重，例如，0.2，造成的偏誤不會像本例。
當然，我們也發現，在有序列相關的情形，係數值問題最嚴重的是隨機效果，固定效果似乎沒有太大影響，固定效果的問題在於顯著性。這個結論，用代數方式說明如下。
假設一個 DGP為 \(Y=X \beta+\varepsilon\)
實證上，\(Y=Xb+e\) 是 \(Y=X \beta+\varepsilon\) 的樣本估計式
\[
\begin{align}
b&={{({X}'X)}^{-1}}{X}'Y={{({X}'X)}^{-1}}{X}'(X\beta +\varepsilon) \\
& ={{({X}'X)}^{-1}}{X}'X\beta +{{({X}'X)}^{-1}}{X}'\varepsilon \\
\end{align}
\]
故 \(b=\beta +{{({X}'X)}^{-1}}{X}' \varepsilon\)
如果無內生性， \(EX'\varepsilon=0\) ，則 \(Eb=\beta\) 的結果，和殘差序列相關就無關。 如果有內生性， \(EX'\varepsilon \ne0\) ， \(b\) 就不再是 \(\beta\) 的不偏估計式了。根據內生性有無和固定效果與隨機效果設定的關係，有關序列相關的檢定與修正的重要性，就十分清楚了：一旦有內生性，序列相關就會影響到參數估計，此時使用隨機效果問題最大，因為隨機效果必須在無內生性時，才具有最佳性質。影響的過程，只需要把序列相關遞迴結構寫出來就知道，如下：
\[
\begin{align}
e_{t} &= \rho e_{t-1}+ \varepsilon_t \\ \\
e_{1} &= \rho e_{0}+ \varepsilon_1 \\
e_{2} &= \rho e_{1}+ \varepsilon_2 \\
...
\end{align}
\]
展開後帶入，計算變異數，就可以證明序列相關程度的影響。此時的參數估計式，包含了一項嚴重的干擾參數，除非序列相關係數很小，不然在Panel Data 短短的時間， 這個影響還蠻大的。
6.2 異質殘差的穩健修正 Robust Coefficient Covariance
我們先瞭解，不修正殘差異質性，會怎樣？假設一個普通線性迴歸
\[
y_{i}=a+bx_{i}+e_{i}
\]
\(e_{i} \sim IID(0,\sigma^2_e)\) 為獨立同分布，也就是同質變異。此時，我們計算y的變異數
\[
Var(y_i)=Ey_i^2-(Ey_i)^2=E(a+bx_{i}+e_{i})^2-(a+bx_{i})^2=Ee_{i}^2
\]
如果殘差真的是\(i.i.d.\)，則 \(Ee_{i}^2 = \sigma^2_e\)
若殘差和解釋變數有關，也就是
\(e_i=\lambda x_i+\varepsilon_i \\ \varepsilon_i \sim IID(0,\sigma^2_{\varepsilon})\)
則，
\[
\begin{align}
Ee_{i}^{2} &=E{{(\lambda {{x}_{i}}+{{\varepsilon }_{i}})}^{2}} \\
& ={{\lambda }^{2}}Ex_{i}^{2}+2\lambda E({{x}_{i}}{{\varepsilon }_{i}})+E\varepsilon _{i}^{2} \\
& ={{\lambda }^{2}}Ex_{i}^{2}+\sigma _{\varepsilon }^{2} \\
& \because E({{x}_{i}}{{\varepsilon }_{i}})=0 \\
\end{align}
\]
如果同質變異 \(\lambda=0\) ，則\(\sigma^2_{e}=\sigma^2_{\varepsilon}\) 。
如果\(\lambda \ne 0\) ，卻假設\(\lambda=0\) ，則變異數就不會是一個常數，會有低估變異數的情形，這低估程度隨 \(x^2\) 的期望值變動；因此，傳統的 \(t \; test\) 檢定會偏高，使得統計顯著性大增。
異質變異的穩健共變異數修正方式，一般迴歸的處理方式，第1章第5節已經有提到，這一節我們介紹在panel data的處理方式。在panel data，異質變異的基本來源就是和資料維度有關， 一則是來自 \(N\) ，一則是來自 \(T\) 。依照文獻，我們將之分成三類：
第1類：White型
這是基於已故計量大師Halbert White之兩篇論文 White (1984) 和 White (1980) ， 以及 Arellano (1987) 的論文所做的延伸。White的修正方式有兩種：第一種是從 \(N\) 的方面去修正一般化的異質變異，但假設殘差無序列相關；第二種則是從 \(T\) 的方面去修正異質變異，但假設沒有橫段面相關。White的修正方式看起來有其侷限，必須假設序列相關或橫段面相依不存在。因此，Arellano方法則是一般化的情況：假設橫段面相依和殘差的序列相關皆有異質變異。
在R內的修正異質變異函數是vcovHC()。
第2類：PCSE型
這種形式是以估計一個無條件 (unconditional) 的穩健共變異數矩陣為基礎，文獻上由 Beck and Katz (1995) 延伸到 Panel data模型。Beck-Katz方法使用了Panel Corrected Standard Errors (PCSE)。PCSE方法從 \(N\) 或 \(T\) 的方向，計算了殘差共變異數的無條件估計式。若從 \(N\) 的方向，藉以調整每個樣本 \(i\) 的殘差序列相關和timewise異質性；若從 \(T\) 的方向，藉以調整在每一個時間點，橫段面相關性和groupwise異質性。也就是可以計算timewise或 groupwise HC 變異數。
使用上必須知道，Beck-Katz公式是建立在 \(N \to \infty\) 或 \(T \to \infty\) 的漸近性質，不是每一種資料結構都可以用。
Beck-Katz在R內的函數是vcovBK()。
第3類：無母數方法
Driscoll and Kraay (1998) 提出在序列相關和橫段面相依之下，無母數方法的異質變異修正。Driscoll-Kraay方法的漸近性質是從 \(T \to \infty\) 的，與 \(N\) 的維度無關。
Driscoll-Kraay方法在R內的函數是vcovSCC()。
對於考慮殘差異質性和序列相關的穩健共變異數，R有2類模組。sandwich，另有pcse，pcse處理的情形是pooled OLS。我們接下來的程式以 sandwich 為主， sandwich全部都可以做。
R做穩健共變數估計選項，非常簡易，只需要將plm的迴歸物件用vcovHC處理即可重新計算共變異數矩陣。請看以下程式碼：
RLab: 估計隨機效果panel data迴歸
coeftest(gsp_re)標準無HC修正的變異數，如下：
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.1677401 0.1431346 15.1448 < 2.2e-16 ***
## x1_hwy 0.0620968 0.0222839 2.7866 0.005451 **
## x2_water 0.0755775 0.0139882 5.4029 8.628e-08 ***
## x3_other -0.0983994 0.0170746 -5.7629 1.175e-08 ***
## x4_private 0.2732151 0.0202804 13.4718 < 2.2e-16 ***
## x5_emp 0.7491027 0.0253493 29.5512 < 2.2e-16 ***
## x6_unemp -0.0058945 0.0008935 -6.5971 7.572e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
接下來修正異質變異數：
R 的各種方法如下三類穩健共變異的內定HC型，也就是: HC, BK和SCC。透過 “vcov=” 宣告。
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.1677401 0.1913600 11.3281 < 2.2e-16 ***
## x1_hwy 0.0620968 0.0503246 1.2339 0.21759
## x2_water 0.0755775 0.0309106 2.4450 0.01470 *
## x3_other -0.0983994 0.0531029 -1.8530 0.06425 .
## x4_private 0.2732151 0.0423981 6.4440 1.997e-10 ***
## x5_emp 0.7491027 0.0683882 10.9537 < 2.2e-16 ***
## x6_unemp -0.0058945 0.0022746 -2.5914 0.00973 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.1677401 0.2254488 9.6152 < 2.2e-16 ***
## x1_hwy 0.0620968 0.0450704 1.3778 0.1686547
## x2_water 0.0755775 0.0291932 2.5889 0.0098021 **
## x3_other -0.0983994 0.0398472 -2.4694 0.0137387 *
## x4_private 0.2732151 0.0397729 6.8694 1.286e-11 ***
## x5_emp 0.7491027 0.0504770 14.8405 < 2.2e-16 ***
## x6_unemp -0.0058945 0.0017016 -3.4640 0.0005601 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.1677401 0.3619031 5.9898 3.158e-09 ***
## x1_hwy 0.0620968 0.0499687 1.2427 0.2143333
## x2_water 0.0755775 0.0249281 3.0318 0.0025084 **
## x3_other -0.0983994 0.0273614 -3.5963 0.0003424 ***
## x4_private 0.2732151 0.0512384 5.3322 1.260e-07 ***
## x5_emp 0.7491027 0.0980089 7.6432 6.000e-14 ***
## x6_unemp -0.0058945 0.0012683 -4.6475 3.922e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
以下是上述三種修正中的進一步處理方式。必須注意，除了vcovHC有三種method可以選，vcovBK和vcovSCC沒有，然後vcovSCC也沒有cluster選項，因為它同時考慮兩種clusters情況的異質變異修正。
vcovHC()的一般形式，有很多選項。解釋如下。
vcovHC(FH_2vlm, method="white1", type="HC0", cluster="time")
method: 有3種處理方法= c(“white1”, “white2”, “arellano”)。
“white1”處理一般異質性，沒有序列相關。“white2”則是以white1為基礎，多考慮了每一個橫斷面 \(i\) 內的共同變異數。“arellano” 則考慮了相關於異質性和序列相關的充分一般可能性。
type: 有5種處理方法= c(“HC0”, “HC1”, “HC2”, “HC3”, “HC4”)，細節如前。
cluster: 有2種處理方法= c(“time”, “group”)。
“time” 則考慮序列相關(serial correlation)；
“group” 則考慮了橫斷面相關。
如果是選method=“white1”，則和 “time”, “group” 這兩個選項無關。
下面列出實際比較，發現不修正異質變異時，所有係數都很顯著，考慮穩健共變異數後的標準差會比較大，使的顯著性檢定的結果趨於保守，不會太輕易就拒絕虛無假設，也使的迴歸結果較為可靠。當然，所謂「較可靠」試一種實證研究的哲學，不是一個為真(TRUE)的說法。好比電影上常常在演，如果有賭客一直贏，就會被懷疑有詐賭；計量經濟學全班期末考都考90分以上，老師雖嘴裡說信任同學，事實上卻不可鄉愿。一個太好的結果，不是不對，只是提醒研究者必須審慎。
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.16774014 0.14990930 14.4603 < 2.2e-16 ***
## x1_hwy 0.06209679 0.02373783 2.6159 0.009064 **
## x2_water 0.07557748 0.01480763 5.1040 4.149e-07 ***
## x3_other -0.09839937 0.01908805 -5.1550 3.190e-07 ***
## x4_private 0.27321506 0.02158193 12.6594 < 2.2e-16 ***
## x5_emp 0.74910272 0.02910093 25.7415 < 2.2e-16 ***
## x6_unemp -0.00589453 0.00088998 -6.6232 6.406e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.1677401 0.2403878 9.0177 < 2.2e-16 ***
## x1_hwy 0.0620968 0.0233249 2.6623 0.0079160 **
## x2_water 0.0755775 0.0176527 4.2814 2.080e-05 ***
## x3_other -0.0983994 0.0229575 -4.2861 2.037e-05 ***
## x4_private 0.2732151 0.0436371 6.2611 6.201e-10 ***
## x5_emp 0.7491027 0.0511126 14.6559 < 2.2e-16 ***
## x6_unemp -0.0058945 0.0016981 -3.4713 0.0005454 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.1677401 0.3649469 5.9399 4.234e-09 ***
## x1_hwy 0.0620968 0.0504183 1.2316 0.2184445
## x2_water 0.0755775 0.0250672 3.0150 0.0026502 **
## x3_other -0.0983994 0.0275648 -3.5697 0.0003784 ***
## x4_private 0.2732151 0.0516901 5.2856 1.613e-07 ***
## x5_emp 0.7491027 0.0990148 7.5656 1.051e-13 ***
## x6_unemp -0.0058945 0.0012820 -4.5979 4.950e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1最後，異質變異修正的方法這麼多，應該怎麼選擇？每一種作法，都有特定控制的層面，例如，假設無序列相關，或無橫段面相依。實證上建議取用 Driscoll and Kraay (1998) 的vcovSCC()，理由是它考慮了橫段面相依和序列相關，再去修正異質變異。這樣的一般性當然比較好。但是，修正程度不一定會很完整。當然，一個實證的處理方式，應該是都做一做，綜合比較後，再取出適合最大交集的結果。
最後，我們介紹 Petersen (2009) 的方法。
使用時須注意：估計的模型，必須以lm()執行模型的線性迴歸所產生的物件，不能用plm()產生的物件，也就是說， Petersen (2009) 異質修正的模型，必需是pooled model。如下，我們先檢視原始無修正的估計結果：
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9260003 0.0525031 36.6836 < 2.2e-16 ***
## x1_hwy 0.0588882 0.0154114 3.8211 0.0001431 ***
## x2_water 0.1185816 0.0123566 9.5966 < 2.2e-16 ***
## x3_other 0.0085530 0.0123541 0.6923 0.4889360
## x4_private 0.3120204 0.0110875 28.1417 < 2.2e-16 ***
## x5_emp 0.5496945 0.0155369 35.3800 < 2.2e-16 ***
## x6_unemp -0.0072715 0.0013836 -5.2553 1.892e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1本範例須要讀取 petersen.src。這個程式是 Petersen 網站上提供的， 本書將之編成 source code 載入使用
Petersen 共變異數修正：Clustered by \(N\)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9260003 0.2142841 8.9881 < 2.2e-16 ***
## x1_hwy 0.0588882 0.0507781 1.1597 0.2465064
## x2_water 0.1185816 0.0345015 3.4370 0.0006182 ***
## x3_other 0.0085530 0.0406234 0.2105 0.8332975
## x4_private 0.3120204 0.0467833 6.6695 4.755e-11 ***
## x5_emp 0.5496945 0.0676971 8.1199 1.735e-15 ***
## x6_unemp -0.0072715 0.0029460 -2.4682 0.0137850 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Petersen 共變異數修正：Clustered by \(T\)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9260003 0.0686890 28.0394 < 2.2e-16 ***
## x1_hwy 0.0588882 0.0142107 4.1439 3.773e-05 ***
## x2_water 0.1185816 0.0040413 29.3422 < 2.2e-16 ***
## x3_other 0.0085530 0.0083878 1.0197 0.3082
## x4_private 0.3120204 0.0063710 48.9753 < 2.2e-16 ***
## x5_emp 0.5496945 0.0210498 26.1140 < 2.2e-16 ***
## x6_unemp -0.0072715 0.0018526 -3.9250 9.411e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Petersen 共變異數修正：Clustered by both \(N\) and \(T\)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9260003 0.2170679 8.8728 < 2.2e-16 ***
## x1_hwy 0.0588882 0.0504639 1.1669 0.2435796
## x2_water 0.1185816 0.0328660 3.6080 0.0003275 ***
## x3_other 0.0085530 0.0395645 0.2162 0.8289038
## x4_private 0.3120204 0.0453753 6.8764 1.227e-11 ***
## x5_emp 0.5496945 0.0681980 8.0603 2.729e-15 ***
## x6_unemp -0.0072715 0.0032344 -2.2482 0.0248330 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
其實看完前面關於異質修正的介紹，cl()和mcl()兩個修正函數應該很容易理解。有興趣的讀者不妨和 Driscoll and Kraay (1998) 的vcovSCC()比較， 其實， 大同小異。
References
S. J. Paris and C. B. Winstein “Trend Estimators and Serial Correlation,” Unpublished Cowles Commission, Discussion Paper, Chicago, 1954.↩︎