자기상관

Published

March 27, 2024

1. What is it? What’s the consequence/problem if we do not address it(=if we just use OLS estimates)?

자기상관은 오차 \(u_i\)와 \(u_j\) 간의 공분산이 0이 아닌 것이다. 즉, \[cov(u_i, u_j|x_i, x_j) = E(u_i u_j) \neq 0 \hspace{10pt} i \neq j\] 이다.

말로 표현하자면, 시간 또는 공간 순으로 배열된 관측치 간의 상관관계가 존재하는 것이다.

자기상관이 존재하는 데 불구하고 OLS 추정을 사용하게 된다면 어떻게 되는가?

요약하자면, OLS 추정은 회귀계수의 분산을 과소 또는 과대 추정한다. 추정된 회귀계수는 선형-불편 추정량이지만, 최소분산 추정량은 아니다.

  1. 잔차분산 \(\hat \sigma^2 = \Sigma \frac{\hat u^2_t}{(n-2)}\)가 실제 모분산 \(\sigma^2\)를 과소추정할 수 있다.

  2. 따라서, \(R^2\)를 과대추정한다. (\(R^2 = 1 - \frac{SS_{res}}{SS_{tot}}\))

  3. \(\sigma^2\)가 과소추정되지 않았더라도, \(var(\hat \beta_2)\)는 \(var(\hat \beta_2)_{AR1}\)를 과소 또는 과대추정할 것이다.

  4. 결과적으로, \(t\), \(F\) 테스트가 유효하다고 말할 수 없다. 오히려 결과를 오도할 수 있다.

    \[ t_{score} = \frac{\hat \beta - \beta_0}{SE_{\hat \beta}} \sim \mathcal{T}_{n-2} \] \[ F = \frac{MS_{model}}{MS_{residual}} \]

2. How to detect it?

library(gujarati)
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
as.numeric.df <- function(data){
  df <- data %>% 
    mutate(across(everything(), as.character)) %>% 
    mutate(across(everything(), as.numeric))
  
  return (df)
}

data <- as.numeric.df(Table12_4) 
fit <- lm(log(Y) ~ log(X), data=data)
summary(fit)

Call:
lm(formula = log(Y) ~ log(X), data = data)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.041164 -0.017041  0.001037  0.018077  0.038719 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.60668    0.05471   29.37   <2e-16 ***
log(X)       0.65222    0.01235   52.80   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.02208 on 44 degrees of freedom
Multiple R-squared:  0.9845,    Adjusted R-squared:  0.9841 
F-statistic:  2788 on 1 and 44 DF,  p-value: < 2.2e-16

자기상관을 감지하는 방법에는 Graphical한 방법과 정량적인 검정 방법이 가능하다. 대표적인 시각적인 방법과 정량적 방법인 Durbin-Watson d test에 대해서 알아본다.

plot(data$Year, residuals(fit) * 100, type="l", 
     xlab=NA, ylab=NA, main="Time sequence plot of residuals")
abline(h = 0)

n <- nrow(data)
residuals.sd <- sqrt(sum(residuals(fit)^2)/(n-2))
standardized.residuals <- residuals(fit) / residuals.sd

lines(data$Year, residuals(fit) / residuals.sd, lty="dashed")
legend("topright", c("100*RES", "SDRES"), lty=c("solid", "dashed"))

Graphical 한 방법에는 시간에 따른 잔차의 개형을 보여주는 방법이 있다. 이때, 잔차의 크기는 \(y\)의 값에 의존적이므로 이를 잔차 표준오차로 나누어 표준화 잔차의 값으로 표현할 수 있다.

lag(residuals(fit), 1)
                          1             2             3             4 
           NA -0.0360675258 -0.0307797209 -0.0267240386 -0.0291603896 
            5             6             7             8             9 
-0.0262463751 -0.0283481204 -0.0175042417 -0.0064187606  0.0070940987 
           10            11            12            13            14 
 0.0184093936  0.0247128352  0.0162890844  0.0253045792  0.0258290878 
           15            16            17            18            19 
 0.0237441145  0.0111308448  0.0183593198  0.0204157821  0.0307809672 
           20            21            22            23            24 
 0.0330231206  0.0316039534  0.0208013425  0.0387191302  0.0144164051 
           25            26            27            28            29 
 0.0017737802  0.0016200723  0.0134709805  0.0137252476  0.0172319138 
           30            31            32            33            34 
-0.0048181019 -0.0062322610 -0.0041175411 -0.0050780893 -0.0106862603 
           35            36            37            38            39 
-0.0235526000 -0.0278739271 -0.0398046993 -0.0411637490 -0.0135764655 
           40            41            42            43            44 
-0.0066740130  0.0108866107  0.0075514987  0.0004531741 -0.0066725611 
           45 
-0.0156503588 
residuals(fit)
            1             2             3             4             5 
-0.0360675258 -0.0307797209 -0.0267240386 -0.0291603896 -0.0262463751 
            6             7             8             9            10 
-0.0283481204 -0.0175042417 -0.0064187606  0.0070940987  0.0184093936 
           11            12            13            14            15 
 0.0247128352  0.0162890844  0.0253045792  0.0258290878  0.0237441145 
           16            17            18            19            20 
 0.0111308448  0.0183593198  0.0204157821  0.0307809672  0.0330231206 
           21            22            23            24            25 
 0.0316039534  0.0208013425  0.0387191302  0.0144164051  0.0017737802 
           26            27            28            29            30 
 0.0016200723  0.0134709805  0.0137252476  0.0172319138 -0.0048181019 
           31            32            33            34            35 
-0.0062322610 -0.0041175411 -0.0050780893 -0.0106862603 -0.0235526000 
           36            37            38            39            40 
-0.0278739271 -0.0398046993 -0.0411637490 -0.0135764655 -0.0066740130 
           41            42            43            44            45 
 0.0108866107  0.0075514987  0.0004531741 -0.0066725611 -0.0156503588 
           46 
-0.0201975369 
plot(lag(residuals(fit) * 100, 1), residuals(fit) * 100,
     xlab="Res1(-1)", ylab="Res1", main="Current residuals versus lagged residuals")

또 다른 방법은 \(t\) 시점의 잔차와 \(t-1\) 시점의 잔차와 비교하여 산점도를 그리는 방법이 있다. 그림과 같이 \(t\) 시점과 \(t-1\) 시점의 산점도 모양이 양의 상관관계를 그리는 것을 직관적으로 확인할 수 있다.

Graphical method는 쉽고 직관적이지만, 정량적인 수치를 제시하지는 않는다.

lmtest::dwtest(fit)

    Durbin-Watson test

data:  fit
DW = 0.21756, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is greater than 0
durbin.watson.test <- function(fit){
  sum((residuals(fit)[-1] - lag(residuals(fit), 1)[-1])^2)/sum(residuals(fit)^2)
}

Durbin-Watson d statistics 는 다음과 같다.

\[d = \frac{\Sigma_{t=2}^{t=n} (\hat u_t - \hat u_{t-1})^2}{\Sigma_{t=1}^{t=n} \hat u_t^2}\] 이는 연속적인 잔차의 차 제곱의 합을 잔차제곱합으로 나눈 것이다. 이 값은 0부터 4의 값을 가질 수 있는데, 0에 가까울수록 양의 자기상관의 증거가 존재하고, 4에 가까울수록 음의 자기상관의 증거가 존재한다고 할 수 있다.

위에서 구한 \(d\) 통계량의 값을 표본의 수 \(n\), 독립변수 \(p\)의 개수에 따라 다른 \(d_L\), \(d_U\)로부터 비교한다.

이외에도, 비모수 검정인 Runs Test, 자기상관 검정의 일반화 방법인 Breusch-Godfrey 방법이 있다.

3. How to fix(address) it?

fit2 <- lm(log(Y) ~ log(X) + Year, data=data)
summary(fit2)

Call:
lm(formula = log(Y) ~ log(X) + Year, data = data)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.034870 -0.010276  0.000406  0.011741  0.049155 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 14.868718   2.712260   5.482 2.06e-06 ***
log(X)       1.028295   0.077552  13.259  < 2e-16 ***
Year        -0.007528   0.001539  -4.890 1.45e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.0179 on 43 degrees of freedom
Multiple R-squared:   0.99, Adjusted R-squared:  0.9896 
F-statistic:  2132 on 2 and 43 DF,  p-value: < 2.2e-16
durbin.watson.test(fit2)
[1] 0.4497413
fit3 <- lm(Y ~ X + I(X^2), data=data)
summary(fit3)

Call:
lm(formula = Y ~ X + I(X^2), data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.3873 -1.1568  0.5073  1.3592  2.7745 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 10.9680426  3.6757881   2.984  0.00468 ** 
X            1.1850370  0.0842133  14.072  < 2e-16 ***
I(X^2)      -0.0028466  0.0004614  -6.169 2.07e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.757 on 43 degrees of freedom
Multiple R-squared:  0.9875,    Adjusted R-squared:  0.987 
F-statistic:  1705 on 2 and 43 DF,  p-value: < 2.2e-16
durbin.watson.test(fit3)
[1] 0.3038069
  1. pure autocorrlation임을 확인한다.

중요한 변수가 생략되었거나 함수의 형태가 잘못되어서 잔차 양상이 관찰될 수 있다. 따라서, 모형에 대해 추가적인 작업을 수행해본다.

  1. pure autocorrlation이라면, 기존 모델의 적절한 조작을 통해 autocorrlation이 없는 모델로 만든다.

이는 자기상관 계수 \(\rho\)를 알고 있을 때와 \(\rho\)를 모를 때로 나뉜다.

  1. \(\rho\)를 알 때

\(\rho\)를 알 떄에는, 수정된 OLS 회귀를 수행할 수 있다.

\[Y_t = \beta_1 + \beta_2 X_t + u_t, \hspace{10pt} u_t = \rho u_{t-1} + \epsilon_t \hspace{5pt} -1 \lt \rho \lt 1\] \[Y_{t-1} = \beta_1 + \beta_2 X_{t-1} + u_{t-1}\] \[\rho Y_{t-1} = \rho \beta_1 + \rho \beta_2 X_{t-1} + \rho u_{t-1}\]

\[(Y_t - \rho Y_{t-1}) = \beta_1 (1 - \rho) + \beta_2 (X_t - \rho X_{t-1}) + \epsilon_t\]

\[Y_t^* = \beta_1^* + \beta_2^* X_t^* + \epsilon_t\] \[Y_t^* = Y_t - \rho Y_{t-1}, \hspace{2pt} \beta_1^* = \beta_1 (1 - \rho),\hspace{2pt} \beta_2^* = \beta_2 \hspace{2pt}, X_t^* = X_t - \rho X_{t-1} \]

  1. \(\rho\)를 모를 때

\(\rho\)를 모를 때는 \(\rho\)가 극단적으로 -1이나 1의 값을 가진다고 생각하고, 1차 차분에 대한 회귀를 수행할 수 있다. 또는, Durbin-Watson d 통계량과 \(\rho\)의 아래 관계를 사용하여 \(\rho\)를 추정한다.

\[\hat \rho \approx 1 - \frac{d}{2}\]

u.hat <- residuals(fit)[-1]
u.lag.hat <- lag(residuals(fit),1)[-1]

fit6 <- lm(u.hat ~ u.lag.hat + 0)

summary(fit6)

Call:
lm(formula = u.hat ~ u.lag.hat + 0)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0197735 -0.0066147  0.0000806  0.0070989  0.0221492 

Coefficients:
          Estimate Std. Error t value Pr(>|t|)    
u.lag.hat  0.86789    0.06814   12.74 2.34e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.009884 on 44 degrees of freedom
Multiple R-squared:  0.7866,    Adjusted R-squared:  0.7818 
F-statistic: 162.2 on 1 and 44 DF,  p-value: 2.337e-16

또 다른 방법으로는 잔차 \(\hat u_{t}\) 를 \(\hat u_{t-1}\)에 대해 회귀하는 것이다. 즉, 다음과 같은 회귀를 수행한다.

\[\hat u_{t} = \rho \hat u_{t-1} + v_t\] 이외에도 반복적인 \(\rho\)의 추정이 가능하다.

fit4 <- data %>% 
  mutate(lag_Y = lag(Y, 1), 
         lag_X = lag(X, 1),
         dY = log(Y) - log(lag_Y), dX = log(X) - log(lag_X)) %>% 
  na.omit() %>% 
  lm(dY ~ dX + 0, data=.)
summary(fit4)

Call:
lm(formula = dY ~ dX + 0, data = .)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.024361 -0.005615 -0.000190  0.006854  0.027542 

Coefficients:
   Estimate Std. Error t value Pr(>|t|)    
dX  0.65386    0.05734    11.4 9.97e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.0103 on 44 degrees of freedom
Multiple R-squared:  0.7472,    Adjusted R-squared:  0.7415 
F-statistic: 130.1 on 1 and 44 DF,  p-value: 9.972e-15
consumption <- as.numeric.df(Table10_7)

fit6 <- lm(log(C) ~ log(Yd) + log(W) + I, data=consumption)

summary(fit6)

Call:
lm(formula = log(C) ~ log(Yd) + log(W) + I, data = consumption)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.018441 -0.010001  0.000337  0.007039  0.032578 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.4677111  0.0427780 -10.933 7.33e-15 ***
log(Yd)      0.8048729  0.0174979  45.998  < 2e-16 ***
log(W)       0.2012700  0.0175926  11.441 1.43e-15 ***
I           -0.0026891  0.0007619  -3.529 0.000905 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.01193 on 50 degrees of freedom
Multiple R-squared:  0.9996,    Adjusted R-squared:  0.9995 
F-statistic: 3.783e+04 on 3 and 50 DF,  p-value: < 2.2e-16
nrow(consumption)
[1] 54
lmtest::dwtest(fit6)

    Durbin-Watson test

data:  fit6
DW = 1.2892, p-value = 0.0009444
alternative hypothesis: true autocorrelation is greater than 0
durbin.watson.test(fit6)
[1] 1.289219
u.hat <- residuals(fit6)[-1]
u.lag.hat <- lag(residuals(fit6),1)[-1]

fit7 <- lm(u.hat ~ u.lag.hat + 0)

summary(fit7)

Call:
lm(formula = u.hat ~ u.lag.hat + 0)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0215502 -0.0078554  0.0002098  0.0066877  0.0305022 

Coefficients:
          Estimate Std. Error t value Pr(>|t|)  
u.lag.hat   0.3247     0.1351   2.402   0.0199 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.01092 on 52 degrees of freedom
Multiple R-squared:  0.09991,   Adjusted R-squared:  0.0826 
F-statistic: 5.772 on 1 and 52 DF,  p-value: 0.01989
rho.hat <- 1 - durbin.watson.test(fit6)/2; rho.hat
[1] 0.3553907
fit5 <- lm(Y ~ X, data=data)
summary(fit5)

Call:
lm(formula = Y ~ X, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.7247 -1.9665  0.3394  2.0095  3.9586 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 32.74190    1.39402   23.49   <2e-16 ***
X            0.67041    0.01567   42.78   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.384 on 44 degrees of freedom
Multiple R-squared:  0.9765,    Adjusted R-squared:  0.976 
F-statistic:  1830 on 1 and 44 DF,  p-value: < 2.2e-16
library(sandwich)
NeweyWest(fit5)
            (Intercept)          X
(Intercept)  3377.60155 -42.325262
X             -42.32526   0.530507
sqrt(diag(vcovHAC(fit5)))
(Intercept)           X 
 3.34556545  0.03700348
  1. 큰 샘플 사이즈라면 Newey-West 방법을 통해 자기상관이 조정된 OLS 추정치의 표준 오차를 구할 수 있다.