최근 공공데이터를 활용한 연구를 진행하고 있습니다. 이 데이터의 경우 이상한 보고 규칙 때문에 일부 데이터가 누락된 상태로 공개되었는데 널리 쓰이는 패키지로 적절하게 처리할 수 없어서 직접 코드를 짜서 처리 중입니다. 해당 모델의 경우 Maximum Likelihood Estimation (MLE)로 점추정 (point estimate)을 구하고 있는데 분산추정을 위해서는 데이터의 구조를 적절하게 고려해야만 했습니다. 경제학에서는 지리적으로 인접한 위치에서 추출한 데이터 사이에 공간적 상관관계 (spatial autocorrelation)을 고려하기 위해 추정량의 표준오차 (standard error)을 보정하는 여러 방법을 씁니다. 그 중에서 가장 유명한 것은 Conley (1999)입니다. 동시에 여러 데이터가 한 위치에서 추출됐을 경우 마찬가지의 자기상관 (autocorrelation) 구조가 발생하는데 이를 고려한 표준오차를 주로 cluster robust 분산으로 부릅니다. 가장 유명한 것이 이른바 heteroskedascity and clustering (HAC) robust 추정량입니다. 이 글에서는 이번 프로젝트를 진행하면서 공부한 것을 간략하게 설명할 것입니다. 바로 일반적인 MLE 추정량의 robust 표준오차에 대한 것입니다. 참고로 이 글은 엄밀한 증명을 추구하지 않습니다. 때가 되면 업데이트가 될 수도…?
Standard error estimation in parametric models
모수 $\theta \in \mathbb{R}^p$에 의해 결정되는 일변수 분포 $y \sim f(\cdot|\theta)$를 생각합시다. 확률밀도함수 $f:\mathbb{R} \times \mathbb{R}^p \rightarrow \mathbb{R}$가 두 번 미분가능하고 연속이라는 가정 하에 $\log{f}$의 1계 도함수와 2계 도함수를 정의합니다.
\[g(y;\theta) = \frac{\partial \log{f}}{\partial \theta}\] \[h(y;\theta) = \frac{\partial^2 \log{f}}{\partial \theta^2}\]행렬계산이 쭉 있을텐데 그때 중요한 것 중 하나가 바로 다루는 행렬의 모양입니다. $f$가 $\theta$에 대한 함수이기 때문에 1계 도함수인 $g$는 $1 \times p$ 행렬이고 $h$는 $p \times p$ 행렬입니다. 구체적인 형태를 쓰면
\[g(y;\theta) = \left[ \frac{\partial \log{f}}{\partial \theta_j} \right]_{j} \quad j=1,\ldots ,p\] \[h(y;\theta) = \left[ \frac{\partial^2 \log{f}}{\partial \theta_i \partial \theta_j} \right]_{ij} \quad i=1, \ldots ,p, \, j=1, \ldots ,p\]입니다.
이를 이용해 MLE 추정량 $\hat{\theta}$의 분산추정량을 쓸 수 있는데 (증명은 생략)
\[\mathrm{Var}\hat{\theta} = \left[\sum_{i=1}^n h(y_i;\hat{\theta}) \right]^{-1} \left[\sum_{i=1}^n g(y_i;\hat{\theta})g(y_i;\hat{\theta})^T \right] \left[\sum_{i=1}^n h(y_i;\hat{\theta}) \right]^{-1}\]가 됩니다. 이를 샌드위치 추정량 (sandwich estimator)이라고 부르며 양쪽의 헤시안을 빵, 가운데의 항을 고기로 부르는 관행이 있습니다. 여기서 중요한 가정은 각 $i$에 대한 관찰이 모두 독립이라는 것입니다. 그래서 $f(\mathbf{y};\theta) = \prod_i^n f(y_i;\theta)$로 쪼개서 추정량을 위의 경우처럼 쓸 수 있었던 것입니다. 그런데 이 글을 쓰는 건 실제로 관찰된 대상들이 독립이 아닌 상황을 극복하기 위함이었습니다.
Cluster-robust standard errors
Cluster-robust 표준오차는 각 대상들이 각 cluster $c$ 안에서 서로 독립이 아니지만 $c \neq c’$ 일 때 $c$의 대상과 $c’$의 대상은 독립적으로 발생하는 경우를 다룹니다. 각 cluster $c$의 데이터를 생성하는 밀도함수를 $f_c$로 쓰면 전체 밀도함수는 $f = \prod_{c} f_c$ 가 됩니다. 각 $c$에 대해 $\log f_c$의 1계 도함수와 2계 도함수를 각각 $g_c$ 그리고 $h_c$로 쓰면 분산추정량은 다음과 같이 쓸 수 있습니다.
\[\left[\sum_{c} h_c(y_c;\hat{\theta}) \right]^{-1} \left[\sum_{c} g_c(y_c;\hat{\theta})g_c(y_c;\hat{\theta})^T \right] \left[\sum_{c} h_c(y_c;\hat{\theta}) \right]^{-1}\]Why Generalized Method of Moments (GMM)?
지금까지 익숙한 MLE 추정에 대해 다뤘습니다. MLE 추정은 모든 모수적 추정 방법 중에 점근적으로 가장 효율적인 추정방법임이 알려져 있습니다 (Cramer-Rao Bound). 문제는 MLE는 모델에 대한 제법 강한 가정을 요구하기 때문에 그 가정이 성립하지 않을 때 크고 작은 문제가 생깁니다. Generalized Method of Moments (GMM)은 MLE보다 더 약한 가정을 요구하는 추정 방법으로 대개 1차 혹은 2차 moment 조건만 요구할 뿐 (즉, 평균 혹은 분산에 대한 가정) 구체적인 분포를 가정하지 않습니다. 그 댓가로 추정의 효율은 떨어지지만 훨씬 유연한 추정을 가능케 합니다.
어떤 모수에 의해 결정되는 분포가 있을 때 모수의 참값을 $\theta_0 \in \mathbb{R}^p$라고 하겠습니다. 이때 어떤 함수 $g:\mathbb{R}^q \times \mathbb{R}^p \rightarrow \mathbb{R}^r$가 존재해서
\[m(\theta_0) = \mathbb{E} \left[g(Y, \theta_0) \right] =0\]를 만족할 때를 생각하겠습니다. 대충 말하면 주어진 데이터 $Y \in \mathbb{R}^q$에 대해 실제 모수 $\theta_0$가 어떤 함수적 조건을 만족해야함을 뜻합니다.
이때 GMM 추정량은 다음으로 주어집니다.
\[\hat{\theta} = \mathrm{argmin}_{\theta \in \mathbb{R}^p} \left( \frac{1}{n} \sum_{i=1}^n g(Y_i,\theta) \right)^T W \left( \frac{1}{n} \sum_{i=1}^n g(Y_i,\theta) \right)\]$W$는 임의로 선택된 positive semi-definite matrix 입니다. 저 추정량은 일련의 조건을 만족할 때 consistent하고 점근적으로 정규분포를 띠며 $W$를 잘 고르면 MLE처럼 efficient 함이 알려져 있습니다.
함수 $g$가 실제 모수에서 기댓값이 0이 됨에 착안해서 추정된 평균 ($\frac{1}{n}$)을 낸 $g$의 2-norm이 0에 최대한 가깝게 하는 $\theta$가 실제 모수 $\theta_0$에도 가까울 것이라는 게 대충의 아이디어 입니다. 이 추정량의 분산은 $ G = \mathbb{E} \left[\nabla_{\theta} g(Y,\theta_0) \right]$, $\Sigma = \mathbb{E} \left[ g(Y,\theta_0) g(Y, \theta_0)^T \right]$라고 했을 때
\[\sqrt{n}\left(\hat{\theta} - \theta_0 \right) \rightarrow \mathcal{N}\left[ 0, (G^TWG)^{-1} G^T W \Sigma W^T G (G^TW^TG)^{-1} \right]\]가 됩니다.
그러면 GMM을 이용하여 몇 가지 잘 알려진 추정량을 얻어보겠습니다.
OLS도 GMM이었어 (feat. 배민)
\[y = X^T \boldsymbol{\beta} + \epsilon\]를 생각합시다. Homoskedastic normal linear model은 $\epsilon \sim \mathcal{N}(0, \sigma^2)$과 $\epsilon \perp X$를 가정하지만 OLS를 GMM으로 얻기 위해서는 더 약한 조건인 $\mathbb{E}[X \epsilon] = 0$만 있으면 충분합니다.
\[g(Y,\theta) = g(y, \boldsymbol{\beta}) = X (y-X^T \boldsymbol{\beta})\]를 생각하면 다음이 성립합니다.
\[G = \mathbb{E} \left[\nabla_{\theta}g(Y,\theta) \right] = - \mathbb{E}\left[XX^T\right]\] \[\Sigma = \mathbb{E}\left[ X \epsilon \epsilon^T X^T \right]\] \[\hat{\boldsymbol{\beta}} = \mathrm{argmin}_{\boldsymbol{\beta} \in \mathbb{R}^p} \left( \frac{1}{n} \sum_{i=1}^n (y_i - X_i^T \boldsymbol{\beta})^T X_i^T \right) W \left( \frac{1}{n} \sum_{i=1}^n X_i(y_i - X_i^T \boldsymbol{\beta}) \right)\]. 상단 최적화 문제의 해는 $\nabla_{\boldsymbol{\beta}}$를 해보면 OLS와 동일한 해를 줍니다. 추정에 대한 검정을 위해서 점근분포를 계산할 수 있고 그 결과를 통해 아주 흥미로운 사실 몇 가지를 얻을 수 있습니다. 계산의 편의를 위해 $W = \mathbb{E} \left[XX^T\right]^{-1}$로 두고 점근분포를 계산하면 다음과 같습니다.
그리고 점근분포를 보면
\[\sqrt{n} (\boldsymbol{\hat{\beta}}- \boldsymbol{\beta}) \rightarrow \mathcal{N} \left[0, (\mathbb{E}\left[XX^T\right])^{-1} \Sigma (\mathbb{E}\left[XX^T\right])^{-1} \right]\]가 되고 자연스럽게 heteroskedasity robust variance를 얻을 수 있습니다. $\Sigma$가 기댓값으로 주어진 것에 착안하여
\[\hat{\Sigma} = \frac{1}{n} \sum_{i=1}^n \epsilon_i^2X_iX_i^T\]를 생각하면 되기 때문입니다.
Poisson regresssion과 GMM
저는 count data를 다루고 있어서 Poisson regression을 GMM으로 추정하는 법을 찾고 있습니다. Poisson regression은 대개 다음으로 정의됩니다.
\[y \sim \mathrm{Poisson}(\mu)\] \[\log{\mu} = X^T \beta\]이 모델의 MLE는
\[F(\beta) = \sum_{i=1}^n y_i \log(\mu_i) - \mu_i - \log(y_i !)\]를 최대화하는 $\beta$를 찾는 것인데 이는 $\nabla_{\beta} F =0 $이 되는 $\beta$를 찾는 것입니다. 계산을 통해 이 조건은 다음과 동치임을 알 수 있습니다.
\[\nabla_{\beta} F = \sum_{i=1}^n X_i (y_i -\mu_i) = 0\]만약 데이터가 실제로 위에 적시된 모델을 따를 경우 Score function의 성질에 의해
\[\mathbb{E} \left[ X(y-\mu) \right] =0\]이 성립합니다. Poisson GMM은 이 관찰로부터 $g(Y, \theta) = g(y, \beta) = X(y-\mathrm{exp}(X^T\beta))$를 이용해 추정을 시행합니다. 여기까지 moment condition이 OLS와 동일하게 생겼기 때문에 결과적으로 $\mu$에 대한 식만 다를 뿐 똑같은 점근 추정을 얻을 수 있습니다. 이 경우는 MLE에 기반한 Poisson regression과 다르게 실제 count가 over-dispersed 되어도 올바른 표준오차를 얻을 수 있습니다. 왜냐면 $\mu$에 대한 moment condition만 가정할 뿐 실제로 분포가 Poisson인지는 GMM에 영향을 주지 않기 때문입니다.
Non-spatial autocorrleation
시계열이나 공간데이터를 다룰 경우 표본들이 서로 독립이 아니게 됩니다. 이 경우에는 상단에 기술한 GMM을 바로 적용할 수 없는데 이는 GMM의 점근분포를 얻기 위해서 중심극한정리를 사용하며 중심극한정리가 각 표본이 서로 독립임을 요구하기 때문입니다. 시계열 (time-series)이나 종단 (longitudinal) 데이터의 경우 한 대상에서 복수의 관찰을 얻는데 이 관찰들은 서로 의존적입니다. 반대로 다른 대상에서 온 관찰들은 독립임을 가정할 수 있어서 약간의 변형을 통해 GMM을 적용할 수 있습니다. 총 $n$개의 대상을 $T$번 관찰했을 때 각 $i$에 대하여 길이가 $T$인 확률변수를 다음과 같이 정의합니다.
\[\mathbf{y}_i = \left[ y_{i1} \, \ldots \, y_{iT} \right]\]각 $y_i$는 독립이고 일변수 확률변수였던 이전의 경우와 다르게 벡터 확률변수가 되었음을 알 수 있습니다.
GMM을 적용하기 위해 $g:\mathbb{R}^T \times \mathbb{R}^p \rightarrow \mathbb{R}^{p \times T}$를 다음과 같이 정의하겠습니다.
\[\left[ g(\mathbf{y}_i, \boldsymbol{\beta}) \right]_{t} = X_{it} (y_{it} - \mathrm{exp}(X_{it}^T\boldsymbol{\beta}))\]이제 필요한 항을 하나씩 계산하겠습니다.
\[\left[\nabla_{\boldsymbol{\beta}} g(\mathbf{y}_i, \boldsymbol{\beta})\right] _{t \cdot} = - X_{it} \mathrm{exp}(X_{it}^T \boldsymbol{\beta}) X_{it}^T\]참고로 이 행렬은 $(p \times T) \times p$ 행렬입니다. $t\cdot$ 번째 원소라는 것은 $t$번째 관찰에 해당하는 $g$의 도함수가 $p \times p$ 행렬이 됨을 뜻합니다. $p \times p$ 행렬을 세로로 쭉 쌓으면 $(p \times T) \times p$ 행렬이 됨을 알 수 있습니다.
다음은 $W$를 고를 차례인데 다음과 같은 $(p \times T) \times (p \times T)$ block diagonal matrix로 정하겠습니다.
\[W = \mathrm{diag} \left( \mathbb{E} \left[ \mu_{it}X_{it}X_{it}^T \right]^{-1} \right)_{t=1\ldots T}\]이렇게 정한 이유는
\[G^TWG = \mathbb{E} \left[ \sum_{t=1}^T \mu_{it}X_{it}X_{it}^T \right]\]로 단순화되기 때문입니다.
지금까지 계산 결과를 종합하여 점근분포를 계산하면 다음을 얻습니다.
\[\sqrt{n} \left(\boldsymbol{\hat{\beta}} - \boldsymbol{\beta} \right) \rightarrow \mathcal{N} \left[0, \left(\mathbb{E} \left[ \sum_{t=1}^T \mu_{it}X_{it}X_{it}^T \right]\right)^{-1} \mathbb{E} \left[ \sum_{s} \sum_t \epsilon_{is} \epsilon_{it} X_{is} X_{it}^T \right] \left(\mathbb{E} \left[ \sum_{t=1}^T \mu_{it}X_{it}X_{it}^T \right]\right)^{-1} \right]\]가운데 항을 관찰하면 이전에는 없던 같은 대상 $i$의 관측치들 사이의 autocorrelation이 생겼음을 볼 수 있습니다 ($t \neq s$에 대한 non-zero $\epsilon$ product들). 이는 전형적인 heteroskedasity and cluster robust variance estimator를 줍니다. 구체적으로 적는 건 population mean $\mathbb{E}$를 sample mean $n^{-1} \sum_{i=1}^n$으로 바꾸면 되므로 생략하겠습니다.
(예고편) Spatial autocorrelation
Spatial autocorrelation은 공간적으로 인접한 두 대상 사이의 측정 사이에 autocorrelation이 있는 경우를 뜻합니다. 앞서 언급한 Conley (1999)은 이런 상황에서 GMM을 일반화시켜 robust한 standard error을 얻는 방법을 찾았습니다. 중심극한정리만큼 유명하면서도 많은 사람들이 잘못 알고 있는 수학 정리가 잘 없는데 많은 모델의 검정에 필수적인 역할을 하고 있기 때문에 반드시 한 번 짚고 넘어가야할 것 같습니다.