I need to perform manually two-stage Least Squares(to illustrate its advantages), where the first stage is repeated median estimate and the second stage should be weighted least squares, where weights are obtained(as far, as I understand) from polynomial regression of first-stage residuals on regressors.

Suppose I have generated the following heteroscedastic model:

Y_i = b_0+b_1 X_i + \epsilon_i

where error depends on regressor:

\epsilon_i \sim N(0,(X_i-1)^2)

b<-c(12,7.25) ## my coefficients
num<-50 ## number of observations

raw_x<-runif(num,min=0,max=2) ## regressors

     rnorm(num,mean=0,sd=(raw_x-1)^2)) ## observations

l<-lm(my_y~raw_x) ## let's create linear model

plot(fitted(l),residuals(l)) ## we see heteroskedasticity
## we got to higher values, our residuals explode

title("Residual vs Fit. value");

Residual vs fit value

So I perform repeated median regression (formula in the Introduction):

## Generating first model using repeated median

## slope 

fij = function(i,j)

bij<-outer(1:num,1:num,fij) ##NaN's were produced on the diagonal    

rowmeds <- apply(bij, 1, median,na.rm=TRUE)

# colmeds <- apply(bij, 2, median,na.rm=TRUE) ## column medians are the same
# b_med3<-median(colmeds)

## Intercept    

## med(y_i - b*x_i)


The fit is extremely accurate! In this example a_med is 11.97634 and b_med equals 7.27022.

Now I perform 2nd order polynomial regression of residuals:

\hat{\epsilon_i}=a_0 + a_1X_i+a_2X_i^2+\delta_i

\hat{\epsilon} = \begin{pmatrix} 1 & X_1 & X_1^2 \\ \vdots \\ 1 & X_m & X_m^2 \end{pmatrix}\begin{pmatrix}a_0 \\ a_1 \\ a_2 \end{pmatrix}+\delta

so that (X here is m x 3 matrix):

\hat{a} = [X^TX]^{-1}X^T\hat{\epsilon}

I was told that as long as residual variances can be roughly estimated from only one observation, actual fit from this model can be used; residual variances = coefficients for the weighted least squares:


## Obtaining 2nd order polynomial estimator for residuals

Xmatr = t(rbind(rep(1,num),raw_x,(raw_x)^2))

## coef_var = (X^T*X)^(-1)*X^T^e

## Obtaining sigma (residual deviation) esitmate


## performing regression with sigma-weights


## final plot

plot(raw_x,my_y, pch=20,col=alpha("salmon",0.6))

abline(b[1],b[2], col="black") ## real line
abline(a_med,b_med,col="blue") ## repeated median fit
abline(b_wls, col="magenta")
legend('bottomright', c("Real","Repeat median","Two-Level LS") , 
       lty=1, col=c('black', 'blue','magenta'), bty='n', cex=.75)

The resulting fit is always worse (and sometimes turns into complete garbage). Please, can you explain me what I'm doing wrong? I need to obtain the result where two-level LS is better than repeated median fit(provided the error depends on regressors as shown before).

enter image description here enter image description here

EDIT: Using R function lm seems to produce the same picture:

Xmatr = t(rbind(rep(1,num),raw_x,(raw_x)^2))
polycoef<-lm(residual ~ poly(raw_x, 2, raw=TRUE))$coefficients

EDIT2: Checking the quality of the residual fit (as far as I understand what's going on):

lines(sort(raw_x),(sort(raw_x)-1)^2) ## real residuals
lines(sort(raw_x), my_sigma[order(raw_x)],col = "magenta") ## fitted residuals
legend('topleft', c("Real res","Fitted res") , 
       lty=1, col=c('black','magenta'), bty='n', cex=.75)

enter image description here

EDIT3: As long as my standard deviation is $(X_i-1)^2$, so the variance has power 4 and maybe I should do for the residual fit


this makes residual fit by the quadratic function better, but the regression line moves even more far than before. enter image description here


Particular solution: looks like "weights" parameter in lm() is the diagonal of $W^{-1}$ in the formula of weighted LS: $b_{wls} = (X^T W^{-1} X)^{-1}W^{-1}y$, $W = diag(\sigma_1,...\sigma_m)$. So


greatly improved the accuracy, but it is still not as good as repeated median in general case. If max value is big in


The picture crashes again.


Not the answer you're looking for? Browse other questions tagged or ask your own question.