Algoritmo GHK

El algoritmo GHK (Geweke, Hajivassiliou y Keane) ^{[ 1 ]} es un método de muestreo de importancia para simular probabilidades de elección en el modelo probit multivariado . Estas probabilidades simuladas pueden utilizarse para recuperar estimaciones de parámetros a partir de la ecuación de máxima verosimilitud mediante cualquiera de los métodos de maximización habituales ( método de Newton , BFGS , etc.). Train ^{[ 2 ]} ofrece pasos bien documentados para implementar este algoritmo en un modelo probit multinomial . Lo que sigue a continuación se aplica al modelo probit multivariado binario.

Consideremos el caso en el que se intenta evaluar la probabilidad de elección de${\displaystyle \Pr(\mathbf {y_ {i}} |\mathbf {X_ {i}\beta } ,\Sigma )}$dónde${\displaystyle \mathbf {y_{i}} =(y_{1},...,y_{J}),\ (i=1,...,N)}$y dónde podemos tomar${\displaystyle j}$como opciones y${\displaystyle i}$como individuos u observaciones, ${\displaystyle \mathbf {X_ {i}\beta}}$es la media y${\displaystyle \Sigma }$es la matriz de covarianza del modelo. La probabilidad de observar la elección${\displaystyle \mathbf {y_ {i}}}$es

${\displaystyle {\begin{aligned}\Pr(\mathbf {y_{i}} |\mathbf {X_{i}\beta } ,\Sigma )=&\int _{A_{J}}\cdots \int _{A_{1}}f_{N}(\mathbf {y} _{i}^{*}|\mathbf {X_{i}\beta } ,\Sigma )dy_{1}^{*}\dots dy_{J}^{*}\\\Pr(\mathbf {y_{i}} |\mathbf {X_{i}\beta } ,\Sigma )=&\int \mathbb {1} _{y^{*}\in A}f_{N}(\mathbf {y} _{i}^{*}|\mathbf {X_{i}\beta } ,\Sigma )d\mathbf {y} _{i}^{*}\end{aligned}}}$

Dónde${\displaystyle A=A_{1}\times \cdots \times A_{J}}$y,

${\displaystyle A_{j}={\begin{cases}(-\infty ,0]&y_{j}=0\\(0,\infty )&y_{j}=1\end{cases}}}$

A menos que${\displaystyle J}$es pequeño (menor o igual que 2) no hay una solución de forma cerrada para las integrales definidas anteriormente (se ha realizado algún trabajo con${\displaystyle J=3}$^{[ 3 ]} ). La alternativa a evaluar estas integrales de forma cerrada o mediante métodos de cuadratura es utilizar la simulación. GHK es un método de simulación para simular la probabilidad anterior utilizando métodos de muestreo de importancia.

Evaluar${\displaystyle \Pr(\mathbf {y_{i}} |\mathbf {X_{i}\beta } ,\Sigma )=\int \mathbb {1} _{y^{*}\in A}f_{N}(\mathbf {y} _{i}^{*}|\mathbf {X_{i}\beta } ,\Sigma )d\mathbf {y} _{i}^{*}}$se simplifica al reconocer que el modelo de datos latentes${\displaystyle \mathbf {y_{i}^{*}} =\mathbf {X_{i}\beta } +\epsilon }$puede reescribirse utilizando una factorización de Cholesky,${\displaystyle \Sigma =CC'}$Esto da${\displaystyle \mathbf {y_{i}^{*}} =\mathbf {X_{i}\beta } +C\eta _{i}}$donde el${\displaystyle \eta _{i}}$Los términos se distribuyen${\displaystyle N(0,\mathbf {I} )}$.

Utilizando esta factorización y el hecho de que la${\displaystyle \eta _{i}}$Se distribuyen independientemente, se pueden simular extracciones de una distribución normal multivariada truncada utilizando extracciones de una distribución normal aleatoria univariada.

Por ejemplo, si la región de truncamiento${\displaystyle \mathbf {A} }$tiene límites inferior y superior iguales a${\displaystyle [a,b]}$(incluyendo a,b =${\displaystyle \pm \infty }$) entonces la tarea se convierte en

${\displaystyle {\begin{array}{lcl}a<&y_{1}^{*}&<b\\a<&y_{2}^{*}&<b\\\vdots &\vdots &\vdots \\a<&y_{J}^{*}&<b\\\end{array}}}$

Nota: ${\displaystyle \mathbf {y_{i}^{*}} =\mathbf {X_{i}\beta } +C\eta _{i}}$, sustituyendo:

${\displaystyle {\begin{array}{lcl}a<&x_{1}\beta _{1}+c_{11}\eta _{1}&<b\\a<&x_{2}\beta _{2}+c_{21}\eta _{1}+c_{22}\eta _{2}&<b\\\vdots &\vdots &\vdots \\a<&x_{J}\beta _{J}+\sum _{k=1}^{J}c_{J,k}\eta _{k}&<b\\\end{array}}}$

Reorganizando lo anterior,

${\displaystyle {\begin{array}{ccc}{\frac {a-x_{1}\beta _{1}}{c_{11}}}&<\eta _{1}<&{\frac {b-x_{1}\beta _{1}}{c_{11}}}\\{\frac {a-(x_{2}\beta _{2}+c_{21}\eta _{1})}{c_{22}}}&<\eta _{2}<&{\frac {b-(x_{2}\beta _{2}+c_{21}\eta _{1})}{c_{22}}}\\\vdots &\vdots &\vdots \\{\frac {a-(x_{J}\beta _{J}+\sum _{k=1}^{J-1}c_{J,k}\eta _{k})}{c_{J,J}}}&<\eta _{J}<&{\frac {b-(x_{J}\beta _{J}+\sum _{k=1}^{J-1}c_{J,k}\eta _{k})}{c_{J,J}}}\\\end{array}}}$

Ahora todo lo que hay que hacer es extraer iterativamente de la distribución normal univariada truncada con los límites dados anteriormente. Esto se puede hacer mediante el método de la función de distribución acumulada inversa y teniendo en cuenta que la distribución normal truncada viene dada por:

${\displaystyle u={\frac {\Phi ({\frac {x-\mu }{\sigma }})-\Phi ({\frac {a-\mu }{\sigma }})}{\Phi ({\frac {b-\mu }{\sigma }})-\Phi ({\frac {a-\mu }{\sigma }})}}}$

Dónde${\displaystyle u}$será un número entre 0 y 1 porque lo anterior es una función de distribución acumulada (CDF). Esto sugiere que para generar extracciones aleatorias de la distribución truncada hay que resolver para${\displaystyle x}$donación,

${\displaystyle x=\sigma F^{-1}(u*(F(\beta )-F(\alpha ))+F(\alpha ))+\mu }$

dónde${\displaystyle \alpha ={\frac {a-\mu }{\sigma }}}$y${\displaystyle \beta ={\frac {b-\mu }{\sigma }}}$y${\displaystyle F}$es la función de distribución acumulada normal estándar. Con tales extracciones se puede reconstruir la${\displaystyle \mathbf {y_{i}^{*}} }$mediante su ecuación simplificada utilizando la factorización de Cholesky. Estas extracciones estarán condicionadas a las extracciones anteriores y, utilizando propiedades de las normales, el producto de las PDF condicionales será la distribución conjunta de la${\displaystyle \mathbf {y_{i}^{*}} }$,

${\displaystyle q(\mathbf {y_{i}^{*}} |\mathbf {X_{1}\beta } ,\Sigma )=q(y_{1}^{*}|\mathbf {X_{1}\beta } ,\Sigma )q(y_{2}^{*}|y_{1}^{*},\mathbf {X_{1}\beta } ,\Sigma )\dots q(y_{J}^{*}|y_{1}^{*},\dots ,y_{J-1}^{*},\mathbf {X_{1}\beta } ,\Sigma )}$

Dónde${\displaystyle q(\cdot )}$es la distribución normal multivariada.

Porque${\displaystyle y_{j}^{*}}$condicionado a${\displaystyle y_{k},\ k<j}$está restringido al conjunto${\displaystyle A}$mediante la configuración utilizando la factorización de Cholesky entonces sabemos que${\displaystyle q(\cdot )}$es una distribución normal multivariada truncada. La función de distribución de una distribución normal truncada es:

${\displaystyle {\frac {\phi ({\frac {x-\mu }{\sigma }})}{\sigma (\Phi ({\frac {b-\mu }{\sigma }})-\Phi ({\frac {a-\mu }{\sigma }}))}}}$

Por lo tanto,${\displaystyle y_{j}^{*}}$tiene distribución,

${\displaystyle {\begin{aligned}q(\mathbf {y_{i}^{*}} |\mathbf {X_{i}\beta } ,\Sigma )&={\frac {{\frac {1}{c_{11}}}\phi _{1}{\Big (}{\frac {y_{j}^{*}-x_{1}\beta }{c_{11}}}{\Big )}}{{\Big (}\Phi _{1}{\Big (}{\frac {b-x_{1}\beta }{c_{11}}}{\Big )}-\Phi _{1}{\Big (}{\frac {a-x_{1}\beta }{c_{11}}}{\Big )}{\Big )}}}\times \dots \times {\frac {{\frac {1}{c_{JJ}}}\phi _{J}{\Big (}{\frac {y_{J}^{*}-(x_{J}\beta +c_{J1}\eta _{1}+c_{J2}\eta _{2}+\dots +c_{JJ-1}\eta _{J-1})}{c_{JJ}}}{\Big )}}{{\Big (}\Phi _{J}{\Big (}{\frac {b-(x_{J}\beta +c_{J1}\eta _{1}+c_{J2}\eta _{2}+\dots +c_{JJ-1}\eta _{J-1})}{c_{JJ}}}{\Big )}-\Phi _{J}{\Big (}{\frac {a-(x_{J}\beta +c_{J1}\eta _{1}+c_{J2}\eta _{2}+\dots +c_{JJ-1}\eta _{J-1}}{c_{JJ}}}{\Big )}{\Big )}}}\\&={\frac {\prod _{j=1}^{J}{\frac {1}{c_{jj}}}\phi _{j}{\Big (}{\frac {y_{j}^{*}-\sum _{k=1}^{k<j}c_{jk}\eta _{k}}{c_{jj}}}{\Big )}}{\prod _{j=1}^{J}{\Big (}\Phi _{j}{\Big (}{\frac {b-\sum _{k=1}^{k<j}c_{jk}\eta _{k}}{c_{jj}}}{\Big )}-\Phi {\Big (}{\frac {a-\sum _{k=1}^{k<j}c_{jk}\eta _{k}}{c_{jj}}}{\Big )}{\Big )}}}\end{aligned}}}$

dónde${\displaystyle \phi _{j}}$es el PDF normal estándar para elegir${\displaystyle j}$.

Porque${\displaystyle y_{j|\{y_{k<j}^{*}\}}^{*}\sim N(\mathbf {X_{i}\beta } +\sum _{k=1}^{k<j}c_{jk}\eta _{k},c_{jj}^{2})}$La estandarización anterior hace que cada término tenga una media de 0 y una varianza de 1.

Sea el denominador${\displaystyle \prod _{j=1}^{J}\Phi _{j}{\Big (}{\frac {b-\sum _{k=1}^{k<j}c_{jk}\eta _{k}}{c_{jj}}}{\Big )}-\Phi {\Big (}{\frac {a-\sum _{k=1}^{k<j}c_{jk}\eta _{k}}{c_{jj}}}{\Big )}=\prod _{j=1}^{J}l_{jj}}$y el numerador${\displaystyle \prod _{j=1}^{J}{\frac {1}{c_{jj}}}\phi _{j}{\Big (}{\frac {y_{j}^{*}-\sum _{k=1}^{k<j}c_{jk}\eta _{k}}{c_{jj}}}{\Big )}=f_{N}(\mathbf {y_{i}^{*}} |\mathbf {X_{i}\beta } ,\Sigma )}$dónde${\displaystyle f_{N}(\cdot )}$es la función de densidad de probabilidad normal multivariada.

Volviendo al objetivo original, evaluar el

${\displaystyle {\begin{aligned}\Pr(\mathbf {y_{i}} |\mathbf {X_{i}\beta } ,\Sigma )=&\int _{A_{j}}f_{N}(\mathbf {y} _{i}^{*}|\mathbf {X_{i}\beta } ,\Sigma )dy_{j}^{*}\\\end{aligned}}}$

Utilizando el muestreo de importancia podemos evaluar esta integral,

${\displaystyle {\begin{aligned}\Pr(\mathbf {y_{i}} |\mathbf {X_{i}\beta } ,\Sigma )=&\int _{A_{j}}f_{N}(\mathbf {y} _{i}^{*}|\mathbf {X_{i}\beta } ,\Sigma )dy_{j}^{*}\\=&\int _{A_{j}}{\frac {f_{N}(\mathbf {y} _{i}^{*}|\mathbf {X_{i}\beta } ,\Sigma )}{q(\mathbf {y_{i}^{*}} |\mathbf {X_{i}\beta } ,\Sigma )}}q(\mathbf {y_{i}^{*}} |\mathbf {X_{i}\beta } ,\Sigma )dy_{j}^{*}\\=&\int _{A_{j}}{\frac {f_{N}(\mathbf {y} _{i}^{*}|\mathbf {X_{i}\beta } ,\Sigma )}{\frac {f_{N}(\mathbf {y} _{i}^{*}|\mathbf {X_{i}\beta } ,\Sigma )}{\prod _{j=1}^{J}l_{jj}}}}q(\mathbf {y_{i}^{*}} |\mathbf {X_{i}\beta } ,\Sigma )dy_{j}^{*}\\=&\mathbb {E} _{\mathbf {q} }{\Big (}\prod _{j=1}^{J}l_{jj}{\Big )}\\\end{aligned}}}$

Esto se aproxima bien mediante${\displaystyle {\frac {1}{S}}\sum _{s=1}^{S}\prod _{j=1}^{J}l_{jj}}$.