Previsão de variação de dados heterocedásticos

Estou tentando fazer uma regressão em dados heterocedásticos em que estou tentando prever as variações de erro, bem como os valores médios em termos de um modelo linear. Algo assim:

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

Em palavras, os dados consistem em medições repetidas de em vários valores de e . Suponho que essas medidas consistam em um valor médio "verdadeiro" que é uma função linear de e , com ruído gaussiano aditivo cujo desvio padrão (ou variância, Eu não decidi) também depende linearmente de . (Eu poderia permitir dependências mais complicadas de e - não há uma forte motivação teórica para uma forma linear - mas prefiro não complicar demais as coisas nesse estágio.)x t ˉ y ( x , t ) x t ξ ( x , t ) x , t x $y(x,t)$$x$$t$$\bar{y}(x,t)$$x$$t$$\xi(x,t)$$x,t$$x$$t$

Eu sei que o termo de pesquisa aqui é "heterocedasticidade", mas tudo o que consegui encontrar até agora são discussões sobre como reduzi-lo / removê-lo para prever melhor , mas nada em termos de tentativa de prever em termos de variáveis ​​independentes. Gostaria de estimar e com intervalos de confiança (ou equivalentes bayesianos), e se existe uma maneira fácil de fazer isso no SPSS, tanto melhor! O que devo fazer? Obrigado.$\bar{y}$ y 0 , a , b , σ 0 , c $\sigma$$y_0, a, b, \sigma_0, c$$d$

Acho que seu primeiro problema é que não é mais uma distribuição normal, e como os dados precisam ser transformados para serem homoscedásticos depende exatamente do que é. Por exemplo, se , o erro é do tipo proporcional e o logaritmo dos dados y deve ser tomado antes da regressão ou, a regressão ajustada a partir dos mínimos quadrados comuns (OLS ) para mínimos quadrados ponderados com um peso de (que altera a regressão para erro de tipo proporcional minimizado). Da mesma forma, se , seria necessário pegar o logaritmo do logaritmo e regredi-lo.$N\left(0,\sigma\left(x,t\right)\right)$$\sigma\left(x,t\right)$$\sigma\left(x,t\right)= ax+bt$$1/y^2$$\sigma\left(x,t\right)= e^{a x+b t}$

Penso que a razão pela qual a predição de tipos de erro é pouco abordada é que primeiro é feita qualquer regressão antiga (gemido, normalmente quadrados mínimos, OLS). E a partir da plotagem residual, isto é, , observa-se a forma residual, e plota-se o histograma de frequência dos dados e olha para isso. Então, se os resíduos são um feixe de ventilador que se abre à direita, tenta-se a modelagem proporcional de dados, se o histograma se parece com um decaimento exponencial, pode-se tentar reciprocidade, , e assim por diante, para raízes quadradas, quadratura, exponenciação , tomando exponencial-y.1 /$model-y$$1/y$

Agora, isso é apenas o conto. A versão mais longa inclui muito mais tipos de regressão, incluindo regressão mediana de Theil, regressão bivariada de Deming e regressão para minimizar o erro dos problemas incorretos 'que não têm nenhuma relação particular de ajuste de curva com o erro propagado sendo minimizado. Esse último é um enorme, mas, veja issocomo um exemplo. Para que faça uma grande diferença as respostas que se está tentando obter. Normalmente, se alguém deseja estabelecer um relacionamento entre variáveis, a rotina OLS não é o método de escolha e a regressão de Theil seria uma melhoria rápida e suja sobre isso. O OLS apenas minimiza na direção y, portanto a inclinação é muito rasa e a interceptação muito grande para estabelecer qual é o relacionamento subjacente entre as variáveis. Dito de outra maneira, o OLS fornece uma estimativa de erro mínima de ay dada um x, não fornece uma estimativa de como x muda com y. Quando os valores de r são muito altos (0,99999+), faz pouca diferença a regressão que se usa e OLS em y é aproximadamente o mesmo que OLS em x, mas, quando os valores de r são baixos, OLS em y é muito diferente de OLS em x.

Em resumo, muito depende exatamente de qual é o raciocínio que motivou a análise de regressão em primeiro lugar. Isso determina os métodos numéricos necessários. Após essa escolha, os resíduos têm uma estrutura relacionada ao objetivo da regressão e precisam ser analisados ​​nesse contexto maior.

O comando de extensão STATS BREUSCH PAGAN pode testar resíduos quanto à heterocedasticidade e estimar em função de alguns ou de todos os regressores.

A abordagem geral para problemas desse tipo é maximizar a probabilidade (regularizada) de seus dados.

No seu caso, a probabilidade do log seria semelhante a onde

$$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)$$ $$\phi(x, \mu, \sigma) = \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$$

Você pode codificar essa expressão em uma função no seu pacote estatístico favorito (eu preferiria Python, R ou Stata, pois nunca fiz programação no SPSS). Em seguida, você pode alimentá-lo com um otimizador numérico, que estimará o valor ideal de seus parâmetros .$\hat{\theta}$$\theta=(y_0, a, b, \sigma_0, c, d)$

Se você precisar de intervalos de confiança, esse otimizador também pode estimar a matriz Hessiana de (segundas derivadas) em torno do ideal. Teoria da estimativa de probabilidade máxima que diz para grande matriz de covariância de pode ser estimada como .$H$$\theta$$n$$\hat{\theta}$$H^{-1}$

Aqui está um exemplo de código em 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)

Observe que a formulação do seu problema pode produzir negativo , e eu tive que me defender dele pela substituição da força bruta de muito pequeno por .$\sigma$$\sigma$$10^{-10}$

O resultado (estimativas de parâmetros e seus erros padrão) produzidos pelo código é:

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

Você pode ver que as estimativas estão próximas de seus valores reais, o que confirma a exatidão dessa simulação.