Predicción de la varianza de los datos heteroscedasticos

Estoy tratando de hacer una regresión en los datos heteroscedasticos donde trato de predecir las varianzas de error, así como los valores medios en términos de un modelo lineal. Algo como esto:

$$\begin{align}\\ y\left(x,t\right) &= \bar{y}\left(x,t\right)+\xi\left(x,t\right),\\ \xi\left(x,t\right) &\sim N\left(0,\sigma\left(x,t\right)\right),\\ \bar{y}\left(x,t\right) &= y_{0}+ax+bt,\\ \sigma\left(x,t\right) &= \sigma_{0}+cx+dt. \end{align}$$

En palabras, los datos consisten en mediciones repetidas de a diversos valores de x y t . Asumo estas mediciones consisten en un "verdadero" valor medio ˉ y ( x , t ) que es una función lineal de x y t , con ruido aditivo gaussiano ξ ( x , t ) cuya desviación estándar (o la varianza, no tengo decidido) también depende linealmente de x , t . (Podría permitir dependencias más complicadas en x y$y(x,t)$$x$$t$$\bar{y}(x,t)$$x$$t$$\xi(x,t)$$x,t$$x$ : no hay una fuerte motivación teórica para una forma lineal, pero prefiero no complicar demasiado las cosas en esta etapa).$t$

Sé que el término de búsqueda aquí es "heteroscedasticidad", pero todo lo que he podido encontrar hasta ahora son discusiones sobre cómo reducirlo / eliminarlo para predecir mejor , pero nada en términos de tratar de predecir σ en términos de variables independientes. Me gustaría estimar y 0 , un , b , σ 0 , c y d con intervalos de confianza (o equivalentes) bayesianos, y si hay una manera fácil de hacerlo en SPSS tanto mejor! ¿Qué tengo que hacer? Gracias.$\bar{y}$ $\sigma$$y_0, a, b, \sigma_0, c$$d$

Creo que su primer problema es que ya no es una distribución normal, y cómo los datos deben transformarse para ser homoscedastic depende exactamente de qué es σ ( x , t ) . Por ejemplo, si σ ( x , t ) = a x + b t , entonces el error es de tipo proporcional y el logaritmo de los datos y debe tomarse antes de la regresión, o, la regresión ajustada de mínimos cuadrados ordinarios (MCO) a ponderados mínimos cuadrados con un $N\left(0,\sigma\left(x,t\right)\right)$$\sigma\left(x,t\right)$$\sigma\left(x,t\right)= ax+bt$ ponderación (que cambia la regresión a un error de tipo proporcional minimizado). De manera similar, si σ ( x , t ) = e a x + b t , uno tendría que tomar el logaritmo del logaritmo y retrocederlo.$1/y^2$$\sigma\left(x,t\right)= e^{a x+b t}$

Creo que la razón por la cual la predicción de los tipos de error está mal cubierta es que primero se hace una regresión antigua (gemido, normalmente mínimos cuadrados ordinarios, MCO). Y a partir de la gráfica residual, es decir, , uno observa la forma residual y la gráfica del histograma de frecuencia de los datos, y lo observa. Luego, si los residuos son un haz de abanico que se abre a la derecha, se intenta el modelado de datos proporcionales, si el histograma se ve como una disminución exponencial, se puede intentar la reciprocidad, 1 / y , y así sucesivamente para raíces cuadradas, cuadratura, exponenciación , tomando exponencial-y.$model-y$$1/y$

Ahora, esa es solo la historia corta. La versión más larga incluye muchos más tipos de regresión, incluida la regresión mediana de Theil, la regresión bivariada de Deming y la regresión para minimizar el error de problemas mal planteados que no tienen una relación particular de bondad de ajuste de curva con el error propagado que se minimiza. Ese último es un whopper, pero, mira estocomo ejemplo. Para que haga una gran diferencia las respuestas que uno está tratando de obtener. Por lo general, si se desea establecer una relación entre variables, la OLS de rutina no es el método de elección, y la regresión de Theil sería una mejora rápida y sucia en eso. OLS solo se minimiza en la dirección y, por lo que la pendiente es demasiado superficial y la intersección demasiado grande para establecer cuál es la relación subyacente entre las variables. Para decir esto de otra manera, OLS da una estimación de error mínimo de ay dada una x, no da una estimación de cómo x cambia con y. Cuando los valores de r son muy altos (0.99999+), la regresión que se usa es mínima y la OLS en y es aproximadamente la misma que OLS en x, pero cuando los valores de r son bajos, la OLS en y es muy diferente de MCO en x.

En resumen, mucho depende exactamente de cuál sea el razonamiento que motivó a hacer el análisis de regresión en primer lugar. Eso dicta los métodos numéricos necesarios. Después de hacer esa elección, los residuos tienen una estructura relacionada con el propósito de la regresión y deben analizarse en ese contexto más amplio.

El comando de extensión STATS BREUSCH PAGAN puede probar los residuales para detectar heterocedasticidad y estimarlo en función de algunos o todos los regresores.

El enfoque general para los problemas de este tipo es maximizar la probabilidad (regularizada) de sus datos.

$$LL(y_0, a, b, \sigma_0, c, d) = \sum_{i=1}^n \log \phi(y_i, y_0 + a x_i + b t_i, \sigma_0 + c x_i + d t_i)$$ where $$\phi(x, \mu, \sigma) = \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$$

You can code this expression into a function in your favorite statistical package (I would prefer Python, R or Stata, for I never did programming in SPSS). Then you can feed it to a numerical optimizer, which will estimate optimal value $\hat{\theta}$ of your parameters $\theta=(y_0, a, b, \sigma_0, c, d)$.

If you need confidence intervals, this optimizer can also estimate Hessian matrix $H$ of $\theta$ (second derivatives) around the optimum. Theory of maximum likelihood estimation says that for large $n$ covariance matrix of $\hat{\theta}$ may be estimated as $H^{-1}$.

Here is an example code in Python:

import scipy
import numpy as np

# generate toy data for the problem
np.random.seed(1) # fix random seed
n = 1000 # fix problem size
x = np.random.normal(size=n)
t = np.random.normal(size=n)
mean = 1 + x * 2 + t * 3
std = 4 + x * 0.5 + t * 0.6
y = np.random.normal(size=n, loc=mean, scale=std)

# create negative log likelihood
def neg_log_lik(theta):
    est_mean = theta[0] + x * theta[1] + t * theta[2]
    est_std = np.maximum(theta[3] + x * theta[4] + t * theta[5], 1e-10)
    return -sum(scipy.stats.norm.logpdf(y, loc=est_mean, scale=est_std))

# maximize
initial = np.array([0,0,0,1,0,0])
result = scipy.optimize.minimize(neg_log_lik, initial)
# extract point estimation
param = result.x
print(param)
# extract standard error for confidence intervals
std_error = np.sqrt(np.diag(result.hess_inv))
print(std_error)

Notice that your problem formulation can produce negative $\sigma$, and I had to defend myself from it by brute force replacement of too small $\sigma$ with $10^{-10}$.

The result (parameter estimates and their standard errors) produced by the code is:

[ 0.8724218   1.75510897  2.87661843  3.88917283  0.63696726  0.5788625 ]
[ 0.15073344  0.07351353  0.09515104  0.08086239  0.08422978  0.0853192 ]

You can see that estimates are close to their true values, which confirms correctness of this simulation.