Wnioskowanie w modelu liniowym z warunkową heteroskedastycznością

Załóżmy, że obserwuję wektory zmiennych niezależnych i i zmienną zależną . Chciałbym dopasować model o postaci: gdzie jest jakąś podwójnie różniczkowalną funkcją o dodatniej wartości, jest nieznanym parametrem skalowania, a jest średnią losową zmienną Gaussa o wariancji jednostkowej (zakładaną jako niezależną od i ). Jest to zasadniczo konfiguracja testu heteroskedastyczności Koenkera (przynajmniej tak dalece, jak to rozumiem).$\vec{x}$$\vec{z}$$y$

$$y = \vec{x}^{\top}\vec{\beta_1} + \sigma g\left(\vec{z}^{\top} \vec{\beta_2}\right) \epsilon,$$ $g$$\sigma$$\epsilon$$\vec{x}$$\vec{z}$

Mam obserwacji z i , i chciałbym oszacować i . Mam jednak kilka problemów:$n$$\vec{x}, \vec{z}$$y$$\vec{\beta_1}$$\vec{\beta_2}$

1. Nie jestem pewien, jak przedstawić problem oszacowania jako coś w rodzaju najmniejszych kwadratów (zakładam, że istnieje dobrze znana sztuczka). Moje pierwsze przypuszczenie byłoby takie jak $$min_{\vec{\beta_1}, \vec{\beta_2}} \left(\sum_{i=1}^n \frac{\left(y_i - \vec{x_i}^{\top}\vec{\beta_1}\right)^2}{g\left(\vec{z_i}^{\top}\vec{\beta_2}\right)^2}\right)\left(\sum_{i=1}^n \frac{1}{g\left(\vec{z_i}^{\top}\vec{\beta_2}\right)^2}\right)^{-1},$$ ale nie jestem pewien, jak rozwiązać to numerycznie (być może mogłaby to zrobić iteracyjna metoda quasi-Newtona).
2. Zakładając, że potrafię postawić problem w rozsądny sposób i znaleźć jakieś oszacowania $\hat{\beta}_1, \hat{\beta}_2$, Chciałbym poznać rozkład szacunków, aby np. Wykonać testy hipotez. Byłbym w porządku z osobnym testowaniem dwóch wektorów współczynników, ale wolałbym jakiś sposób przetestowania, np $H_0: \vec{w_1}^{\top} \vec{\beta_1} + \vec{w_2}^{\top} \vec{\beta_2} \le c$ za dane $\vec{w_1}, \vec{w_2}, c$.

W nieco bardziej ogólnym kontekście z $Y$ na $n$-wymiarowy wektor $y$- obserwacje (odpowiedzi lub zmienne zależne), $X$ na $n \times p$ macierz $x$- obserwacje (zmienne towarzyszące lub zmienne zależne) i $\theta = (\beta_1, \beta_2, \sigma)$ parametry takie, że $Y \sim N(X\beta_1, \Sigma(\beta_2, \sigma))$ wtedy prawdopodobieństwo minus-log jest

$$l(\beta_1, \beta_2, \sigma) = \frac{1}{2}(Y-X\beta_1)^T \Sigma(\beta_2, \sigma)^{-1} (Y-X\beta_1) + \frac{1}{2}\log |\Sigma(\beta_2, \sigma)|$$ W pytaniu PO $\Sigma(\beta_2, \sigma)$ jest przekątna z $$\Sigma(\beta_2, \sigma)_{ii} = \sigma^2 g(z_i^T \beta_2)^2$$ więc determinant staje się $\sigma^{2n} \prod_{i=1}^n g(z_i^T \beta_2)^2$ i wynikowe prawdopodobieństwo logarytmu ujemnego staje się $$\frac{1}{2\sigma^2} \sum_{i=1}^n \frac{(y_i-x_i^T\beta_1)^2}{ g(z_i^T \beta_2)^2} + n \log \sigma + \sum_{i=1}^n \log g(z_i^T \beta_2)$$ Istnieje kilka sposobów podejścia do minimalizacji tej funkcji (przy założeniu, że trzy parametry są niezależne od zmian).

• Możesz spróbować zminimalizować tę funkcję za pomocą standardowego algorytmu optymalizacji, pamiętając o tym ograniczeniu $\sigma > 0$.
• Możesz obliczyć profil minus-log-prawdopodobieństwo $(\beta_1, \beta_2)$ poprzez minimalizację $\sigma$ dla ustalonych $(\beta_1, \beta_2)$, a następnie podłącz wynikową funkcję do standardowego nieograniczonego algorytmu optymalizacji.
• Możesz na przemian optymalizować każdy z trzech parametrów osobno. Optymalizacja ponad$\sigma$ można to zrobić analitycznie, optymalizując ponad $\beta_1$ jest ważonym problemem regresji metodą najmniejszych kwadratów i optymalizacją $\beta_2$ odpowiada dopasowaniu do uogólnionego modelu liniowego gamma $g^2$ odwrotny link.

Ostatnia propozycja przemawia do mnie, ponieważ opiera się na rozwiązaniach, które już dobrze znam. Ponadto pierwsza iteracja jest czymś, co chciałbym rozważyć. To znaczy, najpierw obliczyć wstępne oszacowanie$\beta_1$ przez zwykłe najmniejsze kwadraty ignorując potencjalną heteroskedastyczność, a następnie dopasuj gamma glm do kwadratowych reszt, aby uzyskać wstępne oszacowanie $\beta_2$ $-$aby sprawdzić, czy bardziej skomplikowany model wydaje się opłacalny. Iteracje uwzględniające heteroskedastyczność w roztworze najmniejszych kwadratów, ponieważ wagi mogą następnie poprawić się po oszacowaniu.

Jeśli chodzi o drugą część pytania, prawdopodobnie rozważyłbym obliczenie przedziału ufności dla kombinacji liniowej $w_1^T\beta_1 + w_2^T\beta_2$ albo przez użycie standardowej asymptotyki MLE (sprawdzanie za pomocą symulacji, że asymptotyka działa) lub przez ładowanie.

Edycja: Przez standardowe asymptotyki MLE mam na myśli stosowanie wielowymiarowej normalnej aproksymacji do rozkładu MLE z macierzą kowariancji odwrotnej informacji Fishera. Informacja Fishera jest z definicji macierzą kowariancji gradientu$l$. To zależy ogólnie od parametrów. Jeśli możesz znaleźć wyrażenie analityczne dla tej ilości, możesz spróbować podłączyć MLE. Alternatywnie, możesz oszacować informacje Fishera na podstawie zaobserwowanej informacji Fishera, którą jest Hesjan$l$w MLE. Twój parametr będący przedmiotem zainteresowania to liniowa kombinacja parametrów w dwóch$\beta$-wektory, stąd w przybliżeniu wielowymiarowej normalnej MLE można znaleźć normalne przybliżenie rozkładu estymatorów, jak opisano tutaj . Daje to przybliżony błąd standardowy i można obliczyć przedziały ufności. Jest dobrze opisany w wielu (matematycznych) statystykach, ale dość przystępną prezentacją, którą mogę polecić, jest In All Likelihood Yudi Pawitan. W każdym razie formalne wyprowadzenie teorii asymptotycznej jest dość skomplikowane i opiera się na szeregu warunków prawidłowości i daje tylko prawidłowy asymptotycznydystrybucje. Dlatego w razie wątpliwości zawsze przeprowadzałbym niektóre symulacje z nowym modelem, aby sprawdzić, czy mogę zaufać wynikom w zakresie realistycznych parametrów i wielkości próbek. Proste, nieparametryczne ładowanie początkowe, w którym próbkujesz trzykrotnie$(y_i,x_i,z_i)$ z obserwowanego zestawu danych z wymianą może być użyteczną alternatywą, jeśli procedura dopasowania nie jest zbyt czasochłonna.