劉應興的部落格

誤差項相互不獨立的迴歸模型和普通假設誤差項相互獨立的模型一樣地普遍，只要是系統性地取點的資料，如成對時間序列、地理序列，其迴歸模型的誤差項都需要考慮誤差項之間存在相關。因此，迴歸分析的教本固然是以誤差項 i.i.d. 的情況為討論重點，殘差分析也都會談到誤差項之間相關性的檢查，特別是時間序列的資料做迴歸分析時，常會假設誤差項本身的自相關。

假設只有兩變數：反應變數或依變數 Y，解釋變數或自變數 X，在迴歸分析教本上最常考慮的是誤差項屬一階自迴歸模型 (autoregressive model)：

Y(t) = α + β X(t) + ε(t),  ε(t) = ρ ε(t-1) + u(t)

式中 u(t), t=1,...,n 假設是 i.i.d.，稱為白噪音 (white noise)。當然，也可以假設 ε(t)  服從一般的 ARIMA(p,d,q) 模型，那就是時間序（數）列的轉移函數模型 (transfer function model) 了。迴歸分析教本所談的，通常假設 ε(t) 是平穩 (stationary) AR(1) 模型；當然，也可考慮 ε(t) 是一階（可逆 invertible）移動平均模型 MA(1)：

Y(t) = α + β X(t) + ε(t),  ε(t) = u(t) + ψ u(t-1)

「可逆」意味 MA(q) 可表示為 AR(∞)；「平穩」在 ARMA 模型則表示 AR(p) 可表示為 MA(∞)。所以前面的 AR(1) 誤差項可以寫成

ε(t) = u(t) + ρ u(t-1) + ρ^2 u(t-2) + . . . 

至於 ε(t) 是 MA(1) 時則可以寫成

ε(t) = u(t) + ψ ε(t-1) + ψ^2 ε(t-2) + . . .

如果 Var(u(t)) = σ0^2, 則在 AR(1)  得

Var(ε(t)) = ρ^2 Var(ε(t-1)) + σ0^2

但假設 ε(t) 是一平穩時間數列（平穩隨機過程）則必須 |ρ| < 1 且

Var(ε(t)) = σ0^2/(1-ρ^2)

令 σ^2 = Var(ε(t)), 則 σ0^2 = (1-ρ^2)σ^2；而

Cov(ε(t),ε(s)) = σ^2 ρ^{|t-s|}

若 ε(t) 是平穩 MA(1) 數列，則

Cov(ε(t),ε(s)) = σ^2 ψ,  if |t-s| = 1;   = 0,  otherwise

對於誤差項 ε(t) 是 AR(1) 或 MA(1) 的迴歸模型 Y = Xβ + ε（向量式）, 以普通最小平方法 (OLS, ordinary least squares) 估計，則殘差項為

e = (I-M)Y = (I-M)ε,

式中 M = X(X'X)^(-1)X' 是 X 的行空間 c(X) 上的垂直投影矩陣。對殘差平方和取期望值，

E[e'e] = tr((I-M)E[εε']) = tr((I-M)V) = tr((I-M)V(I-M))

式中 V 是 ε 的共變異矩陣，等於 σ^2 乘以 ε 的相關矩陣。以僅一自變數的迴歸模型為例，M 的第 i 列 j 行元素為

(M)ij = (ΣXt^2-XiΣXt-XjΣXt+nXiXj)/[nΣXt^2-(ΣXt)^2]

故殘差平方和期望值為

        σ^2 ΣΣ{δ(i,j)-(ΣXt^2-XiΣXt-XjΣXt+nXiXj)/[nΣXt^2-(ΣXt)^2]}ρ(i,j)

其中 ρ(i,j) = Corr(ε(i),ε(j))。在 AR(1) 或 MA(1) 誤差項模型中，如果可估得 ρ 或 ψ，從上式也就可估得 σ^2，從而得到 V 的估計。並且進一步可改善 β 的估計。取

A = [ I   0 ]    B = [ 0  I ]
A'B = [ 0  I ]
      [ 0  0 ]

矩陣 A, B 中之 I 為 n-1 階單位矩陣，0 為 1×(n-1) 行向量，其元素皆為 0。則 Ae 是 e 去掉 e(n)，Be 是 e 去掉 e(1)，故

E[Σe(t)e(t-1)] = E[e'A'Be] = tr(A'B(I-M)V(I-M)）

矩陣 A'B 的作用是把後面的矩陣第 2 至 n 列往上提（把第一列去掉），再補上一列 0。所以：上列期望值計算結果是取

V_ = (I-M)V(I-M)

主對角線之下一元素的和，而前面殘差平方和的期望值是取上列 V_ 主對角線元素的和。如果沒有前後乘以 I-M，我們有

ρ = [tr(A'BV)/(n-1)/[tr(V)/n]；但現在 V 換成 V_，是否可用

r' = [nΣe(t)e(t-1)]/[(n-1)Σe(t)^2]

來估計 ρ？如果我們拿 e(t) 建立自迴歸模型，則得 ρ 的 OLS 估計：

r = Σe(t)e(t-1)/Σ_{t>1}e(t)^2

這相當於在 ρ = [tr(A'BV)/(n-1)/[tr(V)/n] 式中把分母的 tr(V) 去掉第一項，因而不需除以其項數。

事實上整個程序要複雜些：在 AR(1) 誤差模型

Y(t) = α + X(t)β + ε(t),   ε(t) = ρ ε(t-1) + u(t)

式中 X(t) 為列向量（單一自變數則 X(t) 是純量）。Cochrane–Orcutt (1949) 首先用 OLS 配適原迴歸模型，並由殘差 e(t) 如上述計算出 r 做為 ρ 的初估值後，考慮

    Y'(t) = α(1-ρ) + X'(t)β + u(t),    Y'(t) = Y(t)-ρY(t-1), X'(t) = X(t)-ρX(t-1)

重新估計 α, β；檢定新模型的殘差是否還存在序列相關，如果證實 u(t) 已是白燥音，則用上列楔型之參數估計轉計算原模型之參數估計並進行後績統計推論；否則就依新的迴歸係數估計計算新的殘差，並得到新的 ρ 的估計，再轉換新變數 Y', X', 計算新的迴歸參數估計。在 u(t) 被「接受」是純燥音時，意味我們也接受了原先的 r 是 ρ 的合理估計；但如發現對應 u(t) 項的殘差值證實 u(t) 是白噪音的假說不成立，也就是原來的 r 不是 ρ 的適當估計，如何由變數轉換後的迴歸模型係數轉回原模型的係數？要計算 ρ 的新估計需要原模型的新殘差，而兩模型的關係：

Y(t) = α + X(t)β + ε(t);   Y'(t) = α(1-ρ) + X'(t)β + u(t)

由 α(1-ρ) 的估計轉成原變數模型的 α 似乎需要 ρ，這不成死結？其實沒事，在原變數模型，可取

α^ = Ybar - (Xbar)β^

也就是只需要變數轉換後模型的 β 估計即可。

檢定模型的殘差是否存在序列相關？有的作者提出的是疊代收斂準則，看疊代過程中的一些估計結果是否穩定；Durbin–Watson 則提出一個假說檢定方法，考慮統計量

d = Σ(e(t)-e(t-1))^2/Σ(e(t))^2

以矩陣形式表示，則是

d= e'(B-A)'(B-A)e/e'e = ε'(I-M)F(I-M)ε/ε'(I-M)ε

式中 F 的第 (i,j) 元素為

Fij = 2  if i = j;  = -1  if |i-j| = 1;  = 0  otherwise

所以其分子期望值是上列 F 作用在 V_  = (I-M)V(I-M) 後取主對角線元素和，分母期望值則如先前說的，是對 V_ 取主對角線元素。如果沒有 I-M 的作用，分子期望值是 2nσ^2(1-ρ), 分母期望值是 nσ^2，也就是說 d ≒ 2(1-ρ)。因 -1≦ρ≦1, 故 d 「應」介於 0 和 4 之間。做檢定時，如果要檢定 Ha: ρ > 0, 則直接 d 與臨界值比較，若 d 值小於較小臨界值，表示 ρ > 0 應接受；若 d 值大於較高臨界值，表示可以認為 Ha 不成立；若 d 值介於兩臨界值之間，則此檢定法無法判定。若對立假說是 Ha: ρ < 0, 則以 4-d 取代 d 與臨界值比較。

總結一下 Cochrane–Orcutt 程序：

(1) 用 OLS 估計 Y(t) = α + X(t)β + ε(t),  t = 1,...,n
(2) e(t) = Y(t) - (α^) - X(t)(β^)
(3) ρ 的估計, 用
  　　  r = Σe(t)e(t-1)/Σ_{t>1}e(t)^2, 或
　　    r = [n/(n-1)]Σe(t)e(t-1)/Σe(t)^2
(4) 計算
　　    Y'(t) = Y(t) - r Y(t-1)
　　    X'(t) = X(t) - r X(t-1)
(5) 用 OLS 配適 Y'(t) = α' + X'(t)β + u(t)
(6) 做 Durbin-Watson 檢定:
   　　 d = Σ(u(t)^ -u(t-1)^)^2/Σ(u(t)^)^2
    d 與 d_L, d_U 比較以判斷是否接受 Ha: ρ_u > 0
    4-d 與 d_L, d_U 比較判斷是否接受 Ha: ρ_u < 0
(7) 若判定 ρ_u = 0 則 Y(t) = α + X(t)β + ε(t) 之
　 　   β^ 為 (5) 所得
  　　  α^ = Ybar - Xbar．β^
　　　       或 α'^/(1-ρ^)
　　    ρ^ = r
　　    σ0^2 = Var(u(t)) 之估計用 (5) 之 SSEu/(df)
　　    σ^2 = Var(ε(t)) 之估計用 σ0^2/(1-ρ^2) 之估計。
(8) 若在 (6) 判定 ρ>0 或 ρ<0 則用 (5) 之 β^ 及
　　    α^ = Ybar - Xbar．β^
    代入 (2) 回到 (2)。

上述 Cochrane–Orcutt 程序除了 ρ>0 時用 r 估計有低估問題以外，也被認為失去了一些效率，因為在變數轉換模型只有 n-1 個資料點，Prais–Winsten (1954) 程序則補上

        √(1-ρ^2) Y(1) = α√(1-ρ^2) + √(1-ρ^2) X(1) + √(1-ρ^2) ε(1)

因為 ε(1) 與 u(2),...,u(n) 是相互獨立的。但是，此資料項的常數項與其他資料項不等：

√(1-ρ^2) α ≠ (1-ρ) α

如果單獨列一個效應項，同樣也是效率喪失；如果強把它們視為相等，結果製造出偏誤。不過，Prais–Winsten 程序似乎不是用轉換變數變成適用 OLS 的模型，而是估計 ρ 後做 GLS？沒找到原始論文，有找到相關的卻看不懂，網路上介紹的也看不懂 . . .

如果 ρ 已知，則可計算誤差項共變異矩陣 V 的反矩陣，其 (i,j) 元素是 c(i,j)/[σ^2(1-ρ^2)], 其中

c(i,j) = 1,        if i = j = 1 or n;
       = 1+ρ^2,    if 1 < i = j < n;
       = -ρ,       if |i-j| = 1;
       = 0,        ortherwise.   

當我們有 ρ 的良好估計 r 時，把 ρ = r 代入求得 V 的反矩陣的估計，而後計算模型的參數估計

Y = Xβ + ε,  β^ = (X'V*^(-1)X)^(-1)X'V*^(-1)Y

以 P* 表示投影矩陣 X(X'V*^(-1)X)^(-1)X'V*^(-1)，則殘差

e = (I-P*)Y = (I-P*)ε

其平方和期望值為

E[e'e] = tr((I-P*)V(I-P*)')

但在 GLS 我們考慮 e'V^(-1)e, 現在以 V* 代 V，得

E[e'V*^(-1)e]  = tr((I-P*)'V*^(-1)(I-P*)V) = tr((I-P*)'V*^(-1)V)

當 ρ 用 r 估計而 σ^2 未知時，V* 用的是估計的相關矩陣，若 r ≒ ρ 則 V ≒ σ^2 V*，所以

MSE = e'V*^(-1)e/tr(I-P*) = e'V*^(-1)e/rank(I-P*)

將是 σ^2 的適當估計。而 Σe(t)e(t-1) = e'A'Be 之期望值為

    E[Σe(t)e(t-1)] = tr(A'B(I-P*)V(I-P*)') ≒ σ^2 tr(A'BV*(I-P*)')

可以用 r = Σe(t)e(t-1)/MSE 取代 r = Σe(t)e(t-1)/Σ_{t>1}e(t)^2 做為 ρ 的新估計。這樣的程序可以疊代至疊代收斂準則符合，或用 Durbin-Watson 檢定判定

u(t)^ = e(t) - ρ e(t-1),  t = 2,...,n

沒有序列相關。注意是採用原先配適模型的 ρ（估計或設定值）而非根據殘差重新計算的 r，因為要檢定的是原先設定的模型適當與否。如果檢定結果 u(t) 之間存在序列相關，即表示原先用以估計 α, β, ρ 及 σ^2 的模型（假設或估計的 ρ 值）不正確，需要用新估計的 ρ 值重新配適模型。

總結用 GLS 程序疊代如下：

(1) 設定一個 ρ 的初值, 誤差項的相關矩陣 Ω = [ρ^|i-j|],
    或設定 Ω = I, P = X(X'Ω^(-1)X)^(-1)X'
(2) 用 GLS 配適模型
         Y(t) = α + X(t)β + ε(t),  t = 1,...,n
(3) 用 MSE = e'Ω^(-1)e/tr(I-P) = e'Ω^(-1)e/rank(I-P)
    估計 σ^2; 用
        r = Σe(t)e(t-1)/MSE 或
        r = Σe(t)e(t-1)/Σ_{t>1}e(t)^2 或
        r = [n/(n-1)]Σe(t)e(t-1)/Σe(t)^2
    做為 ρ 的新估計.
(4) 依疊代停止準則或用 Durbin-Watson 檢定於
        u(t)^ = e(t) - ρ^ e(t-1)  (用舊的 ρ 估計)
    檢定 u(t) 是否存在序列相關。若未達疊代停止準則或 u(t)
    尚有序列相關則回到 (2)（用新的 ρ 估計值）。
(5) 若 r ≒ 1 則改配適差分模型:
      Y(t) - Y(t-1) = α* + (X(t)-X(t-1))β + u(t)

如「線性模型：當誤差非等幅變異時」文中說的，在 V 未知時用 OLS 或用 V* 取代 V 做 GLS 結果 E[Y] = Xβ 的估計仍是不偏的，只是可能不是 BLUE；而最大的問題只是誤差的估計錯亂了，因此對可估函數 λ'β 或 Λ'β 估計式的評估或假說檢定，也都不適當。而另一方面，Xβ 及其他可估函數的估計相對地還是較穩健的，因此 Newey–West 把重點放在誤差的估算。

在談 Cochrane–Orcutt 程序時我們提到該程序可能低估自相關幅度；Prais–Winsten 嫌該法捨棄第一個資料點導致效率損失而對該點提出不同的變換，但個人不瞭解其實際操作；Hildreth–Lu 則提出另一套程序：選擇 ρ 的估計以極小化 u(t)^ 的平方和，也就是極小化

SSEu(ρ) = Σ[Y(t)-ρY(t-1)-(1-ρ)α-(Xt-ρX(t-1))β]^2

具體做法是假設不同 ρ 值，用 OLS 計算對應的 SSEu(ρ) 而後取其最小的。

如果估計的 ρ 接近 1，表示 ρ = 1 可能是合理或正確的。在 ρ = 1 時，原迴歸式誤差項

ε(t) = ε(t-1) + u(t)

是不平穩的時間數列，Var(ε(t)) 是常數的假設不成立。但此時

Y(t) - Y(t-1) = α* + (X(t)-X(t-1))β + u(t)

符合 OLS 的假設，估計了 β 和 σ0^2 = Var(u(t)) 後也可以代回原迴歸式：

Y(t) = α + X(t)β + ε(t)

其中常數項利用迴歸函數通過平均點的特性可以算出來。

我們沒有考慮 ε(t) = -ε(t-1) + u(t) 或 ε(t) = kε(t-1)+u(t) 其中 |k| > 1 的情況，當然如果碰到這樣的情形，或其他或許更複雜的模型，就需要針對該模型進行不同於以上方法的處理了。