Multilinear regression The following code fits multi-linear regression model and produces coefficients that I don't quite know where they come from.
n <- 10000         # use large sample size to get precise estimates
x1 <- rnorm(n, 0, 1) #generated 10000 random numbers, mean =0, sd =1
x2 <- x1 + rnorm(n, 0, 0.1) #error random numbers, mean =0, sd =0.1
y <- 2 * x1 - x2 + rnorm(n, 0, 1)
x3 <- y + rnorm(n, 0, 1)

# Fit different regression models:
coef(lm(y ~ x1))           # beta1=1 
coef(lm(y ~ x2))           # beta2=1 
coef(lm(y ~ x3))           # beta3=2/3
coef(lm(y ~ x1 + x2))      # beta1=2, beta2=-1
coef(lm(y ~ x1 + x2 + x3)) # beta1=1, beta2=-1/2, beta3=1/2

The beta3 is 2/3, why is that? Can someone explicitly show me the calculation process?
The purpose of this exercise is to illustrate that $\beta_k$ is the "effect" of $x_k$ on Y when all other variables in the model are held constant. This "effect" should not be interpreted as a causal effect. When comparing two observations i and j, where $x_{ik} = x_{jk} + 1$ (i.e., values for $x_k$ differ by one) and $x_{ir} = x_{jr}$ for all other $r$ in $\{1,..,p\}$ (i.e., values for all other variables $x_r$ are identical), then $E(Y_i) = E(Y_j) + \beta_k$.
I do understand the purpose of this example perfectly, I just wonder what calculation is behind the 2/3 value
 A: Imagine the data are plotted with $y$ at y-axis and $x3$ at x-axis. The rnorm() term in "x3 <- y + rnorm(n, 0, 1)" smears the data point horizontally to the left and right.
Recalling that the slope, $\beta_1$, is computed by:
$$\beta_1 = \frac{\Sigma(x_i - \bar{x})(y_i - \bar{y})}{\Sigma(x_i - \bar{x})^2}$$
Holding $y_i$ unchanged, the more extreme the individual point, $x_i$, spread out, the denominator would increase faster than the numerator would, hence as you dial up the SD of in the "rnorm()", the beta moves closer to 0. For instance, try "rnorm(n, 0, 0.1)" and "rnorm(n, 0, 10)" and check it yourself.


Or more explicitly, why is it 2 in the numerator and 3 in the
  denominator?

Top part can be perceived as a function of the covariance of $y$ and $x3$ (sans the divide by (n-1) part but the ratio is not changing here), and the denominator can be perceived as the variance of $x3$ alone. Try:
n <- 10000      
x1 <- rnorm(n, 0, 1)
x2 <- x1 + rnorm(n, 0, 0.1) 
y <- 2 * x1 - x2 + rnorm(n, 0, 1)
x3 <- y + rnorm(n, 0, 1)

cov(x3,y)
var(x3)

Then you should see:
> cov(x3,y)
[1] 1.964835
> var(x3)
[1] 2.966569

Which is roughly 2/3.
