How can you handle unstable $\beta$ estimates in linear regression with high multi-collinearity without throwing out variables? Beta stability in linear regression with high multi-collinearity?
Let's say in a linear regression, the variables $x_1$ and $x_2$ has high multi-collinearity (correlation is around 0.9).
We are concerned about the $\beta$ coefficient stability so we have to treat the multi-collinearity.
The textbook solution would be to just throw away one of variables. 
But we don't want to lose useful information by simply throwing away variables.
Any suggestions? 
 A: I like both of the answers given thus far.  Let me add a few things.  
Another option is that you can also combine the variables.  This is done by standardizing both (i.e., turning them into z-scores), averaging them, and then fitting your model with only the composite variable.  This would be a good approach when you believe they are two different measures of the same underlying construct.  In that case, you have two measurements that are contaminated with error.  The most likely true value for the variable you really care about is in between them, thus averaging them gives a more accurate estimate.  You standardize them first to put them on the same scale, so that nominal issues don't contaminate the result (e.g., you wouldn't want to average several temperature measurements if some are Fahrenheit and some are Celsius).  Of course, if they are already on the same scale (e.g., several highly-correlated public opinion polls), you can skip that step.  If you think one of your variables might be more accurate than the other, you could do a weighted average (perhaps using the reciprocals of the measurement errors).  
If your variables are just different measures of the same construct, and are sufficiently highly correlated, you really could just throw one out without losing much information.   As an example, I was actually in a situation once, where I wanted to use a covariate to absorb some of the error variance and boost power, but where I didn't care about that covariate--it wasn't germane substantively.  I had several options available and they were all correlated with each other $r>.98$.  I basically picked one at random and moved on, and it worked fine.  I suspect I would have lost power burning two extra degrees of freedom if I had included the others as well by using some other strategy.  Of course, I could have combined them, but why bother?  However, this depends critically on the fact that your variables are correlated because they are two different versions of the same thing; if there's a different reason they are correlated, this could be totally inappropriate. 
As that implies, I suggest you think about what lies behind your correlated variables.  That is, you need a theory of why they're so highly correlated to do the best job of picking which strategy to use.  In addition to different measures of the same latent variable, some other possibilities are a causal chain (i.e., $X_1\rightarrow X_2\rightarrow Y$) and more complicated situations in which your variables are the result of multiple causal forces, some of which are the same for both.  Perhaps the most extreme case is that of a suppressor variable, which @whuber describes in his comment below.  @Macro's suggestion, for instance, assumes that you are primarily interested in $X$ and wonder about the additional contribution of $Z$ after having accounted for $X$'s contribution.  Thus, thinking about why your variables are correlated and what you want to know will help you decide which (i.e., $x_1$ or $x_2$) should be treated as $X$ and which $Z$.  The key is to use theoretical insight to inform your choice.  
I agree that ridge regression is arguably better, because it allows you to use the variables you had originally intended and is likely to yield betas that are very close to their true values (although they will be biased--see here or here for more information).  Nonetheless, I think is also has two potential downsides:  It is more complicated (requiring more statistical sophistication), and the resulting model is more difficult to interpret, in my opinion.  
I gather that perhaps the ultimate approach would be to fit a structural equation model.  That's because it would allow you to formulate the exact set of relationships you believe to be operative, including latent variables.  However, I don't know SEM well enough to say anything about it here, other than to mention the possibility.  (I also suspect it would be overkill in the situation you describe with just two covariates.)  
A: Well, there is one ad hoc method that I've used before. I'm not sure if this procedure has a name but it makes sense intuitively. 
Suppose your goal is to fit the model
$$ Y_i = \beta_0 + \beta_1 X_i + \beta_2 Z_i + \varepsilon_i $$ 
where the two predictors - $X_i, Z_i$ - are highly correlated. As you've pointed out, using them both in the same model can do strange things to the coefficient estimates and $p$-values. An alternative is to fit the model 
$$ Z_i = \alpha_0 + \alpha_1 X_i + \eta_i $$ 
Then the residual $\eta_i$ will be uncorrelated with $X_i$ and can, in some sense, be thought of as the part of $Z_i$ that is not subsumed by its linear relationship with $X_i$. Then, you can proceed to fit the model 
$$ Y_i = \theta_0 + \theta_1 X_i + \theta_2 \eta_i + \nu_i $$
which will capture all of the effects of the first model (and will, indeed, have the exact same $R^2$ as the first model) but the predictors are no longer collinear. 
Edit: The OP has asked for an explanation of why the residuals do not, definitionally, have a sample correlation of zero with the predictor when you omit intercept like they do when the intercept is included. This is too long to post in the comments so I made an edit here. This derivation is not particularly enlightening (unfortunately I couldn't come up with a reasonable intuitive argument) but it does show what the OP requested: 
When the intercept is omitted in simple linear regression, $\hat \beta = \frac{ \sum x_i y_i}{\sum x_i^2}$, so $e_i = y_i - x_i \frac{ \sum x_i y_i}{\sum x_i^2}$. The sample correlation between $x_i$ and $e_i$ is proportional to $$\overline{xe} - \overline{x}\overline{e}$$ where $\overline{\cdot}$ denotes the sample mean of the quantity under the bar. I'll now show this is not necessarily equal to zero. 
First we have
$$\overline{xe} = \frac{1}{n} \left( \sum x_i y_i - x_{i}^2 \cdot \frac{ \sum x_i y_i}{\sum x_i^2} \right) = \overline{xy} \left( 1 - \frac{ \sum x_{i}^2}{ \sum x_{i}^2 } \right) = 0$$
but 
$$\overline{x} \overline{e} = \overline{x} \left( \overline{y} - \frac{\overline{x} \cdot \overline{xy}}{\overline{x^2}} \right) = \overline{x}\overline{y} - \frac{\overline{x}^2 \cdot \overline{xy}}{\overline{x^2}}$$
so in order for the $e_i$ and $x_i$ to have a sample correlation of exactly 0, we need $\overline{x}\overline{e}$ to be $0$. That is, we need $$ \overline{y} = \frac{ \overline{x} \cdot \overline{xy}}{\overline{x^2}} $$
which does not hold in general for two arbitrary sets of data $x, y$. 
A: You can try ridge regression approach in the case when correlation matrix is close to singular (i.e. variables have high correlations). It will provide you with a robust estimate of $\beta$. 
The only question is how to choose regularization parameter $\lambda$. It is not a simple problem, though I suggest to try different values.
Hope this helps!
