Algoritmo GHK

Método de muestreo de importancia

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 se pueden utilizar para recuperar estimaciones de parámetros de la ecuación de verosimilitud maximizada utilizando cualquiera de los métodos de maximización conocidos ( método de Newton , BFGS , etc.). Train ^[2] tiene pasos bien documentados para implementar este algoritmo para un modelo probit multinomial. Lo que sigue aquí se aplicará al modelo probit binario multivariado.

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

${\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 y, ${\displaystyle A=A_{1}\times \cdots \times A_{J}}$

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

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

La evaluación se simplifica al reconocer que el modelo de datos latentes se puede reescribir utilizando una factorización de Cholesky, . Esto da como resultado dónde se distribuyen los términos . ${\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}^{*}}$ ${\displaystyle \mathbf {y_{i}^{*}} =\mathbf {X_{i}\beta } +\epsilon }$ ${\displaystyle \Sigma =CC'}$ ${\displaystyle \mathbf {y_{i}^{*}} =\mathbf {X_{i}\beta } +C\eta _{i}}$ ${\displaystyle \eta _{i}}$ ${\displaystyle N(0,\mathbf {I} )}$

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

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

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

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

${\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}}}$

Reordenando arriba,

${\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})}{c_{J,J}}}&<\eta _{k}<&{\frac {b-(x_{J}\beta _{J}+\sum _{k=1}^{J-1}c_{J,k})}{c_{J,J}}}\\\end{array}}}$

Ahora, todo lo que se necesita hacer es extraer iterativamente de la distribución normal univariante truncada con los límites indicados anteriormente. Esto se puede hacer mediante el método de la función de distribución acumulativa inversa y teniendo en cuenta que la distribución normal truncada está 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 }})}}}$

Donde será un número entre 0 y 1 porque lo anterior es una función de distribución acumulativa. Esto sugiere generar extracciones aleatorias a partir de la distribución truncada que se debe resolver para obtener, ${\displaystyle u}$ ${\displaystyle x}$

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

donde y y es la CDF normal estándar. Con tales valores extraídos se puede reconstruir la por su ecuación simplificada utilizando la factorización de Cholesky. Estos valores extraídos estarán condicionados a los valores extraídos antes y utilizando las propiedades de las normales el producto de las PDF condicionales será la distribución conjunta de la , ${\displaystyle \alpha ={\frac {a-\mu }{\sigma }}}$ ${\displaystyle \beta ={\frac {b-\mu }{\sigma }}}$ ${\displaystyle F}$ ${\displaystyle \mathbf {y_{i}^{*}} }$ ${\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 está la distribución normal multivariada? ${\displaystyle q(\cdot )}$

Como la condicional está restringida al conjunto por la configuración que utiliza la factorización de Cholesky, entonces sabemos que es una normal multivariada truncada. La función de distribución de una normal truncada es, ${\displaystyle y_{j}^{*}}$ ${\displaystyle y_{k},\ k<j}$ ${\displaystyle A}$ ${\displaystyle q(\cdot )}$

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

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

${\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 está el pdf normal estándar para elección ? ${\displaystyle \phi _{j}}$ ${\displaystyle j}$

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

Sea el denominador y el numerador donde es la función de densidad de probabilidad normal multivariada. ${\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}}$ ${\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 )}$ ${\displaystyle f_{N}(\cdot )}$

Volviendo al objetivo original, evaluar la

${\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}}}$

Usando 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}}$

Referencias

1. Hajivassiliou, Vassilis (1994). "MÉTODOS CLÁSICOS DE ESTIMACIÓN PARA MODELOS LDV USANDO SIMULACIÓN" (PDF) . Manual de econometría .
2. Train, Kenneth (2003). Métodos de elección discreta con simulación . Cambridge University Press.
3. Greene, William (2003). Análisis econométrico . Prentice Hall.