Predicting variance of heteroscedastic data

Question

I am trying to do a regression on heteroscedastic data where I'm trying to predict the error variances as well as the mean values in terms of a linear model. Something like this:

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

In words, the data consists of repeated measurements of $y(x,t)$ at various values of $x$ and $t$. I assume these measurements consist of a "true" mean value $\bar{y}(x,t)$ which is a linear function of $x$ and $t$, with additive Gaussian noise $\xi(x,t)$ whose standard deviation (or variance, I haven't decided) also depends linearly on $x,t$. (I could allow more complicated dependencies on $x$ and $t$ – there isn't a strong theoretical motivation for a linear form – but I'd rather not overcomplicate things at this stage.)

I know the search term here is "heteroscedasticity," but all I've been able to find so far are discussions of how to reduce/remove it to better predict $\bar{y}$, but nothing in terms of trying to predict $\sigma$ in terms of the independent variables. I would like to estimate $y_0, a, b, \sigma_0, c$ and $d$ with confidence intervals (or Bayesian equivalents), and if there is an easy way to do it in SPSS so much the better! What should I do? Thanks.

Answer 1

I think your first problem is that $N\left(0,\sigma\left(x,t\right)\right)$ is not longer a normal distribution, and how the data needs to be transformed to be homoscedastic depends on exactly what $\sigma\left(x,t\right)$ is. For example, if $\sigma\left(x,t\right)= ax+bt$, then the error is proportional type and the logarithm of the y data should be taken before regression, or, the regression adjusted from ordinary least squares (OLS) to weighted least squares with a $1/y^2$ weight (that changes the regression to minimized proportional type error). Similarly, if $\sigma\left(x,t\right)= e^{a x+b t}$, one would have to take the logarithm of the logarithm and regress that.

I think the reason why prediction of error types is poorly covered is that one first does any old regression (groan, typically ordinary least squares, OLS). And from the residual plot, i.e., $model-y$, one observes the residual shape, and one plots the frequency histogram of the data, and looks at that. Then, if the residuals are a fan beam opening to the right, one tries proportional data modeling, if the histogram looks like an exponential decay one might try reciprocation, $1/y$, and so on and so forth for square roots, squaring, exponentiation, taking exponential-y.

Now, that is only the short story. The longer version includes an awful lot more types of regression including Theil median regression, Deming bivariate regression, and regression for minimization of ill-posed problems' error that have no particular goodness-of-curve-fit relationship to the propagated error being minimized. That last one is a whopper, but, see this as an example. So that it makes a big difference what answers one is trying to obtain. Typically, if one wants to establish a relationship between variables, routine OLS is not the method of choice, and Theil regression would be a quick and dirty improvement on that. OLS only minimizes in the y-direction, so the slope is too shallow, and the intercept too large to establish what the underlying relationship between the variables is. To say this another way, OLS gives a least error estimate of a y given an x, it does not give an estimate of how x changes with y. When the r-values are very high (0.99999+) is makes little difference what regression one uses and OLS in y is approximately the same as OLS in x, but, when the r-values are low, OLS in y is very different from OLS in x.

In summary, a lot depends on exactly what the reasoning is that motivated doing the regression analysis in the first place. That dictates the numerical methods needed. After that choice is made, the residuals then have a structure that is related to the purpose of the regression, and need to be analyzed in that larger context.

Answer 2

The STATS BREUSCH PAGAN extension command can both test residuals for heteroscedasticity and estimate it as a function of some or all of the regressors.

Answer 3

The general approach to problems of this kind is to maximize the (regularized) likelihood of your data.

In your case, the log-likelihood would look like $$ 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.