Linear regression

The ommitted variable problem

1 Trivariate Gaussian distribution

Let \((y_i,x_i,z_i)\) for \(i=1,\ldots,n\) be a random sample from a trivariate Gaussian distribution, as follows:

\[ (y_i,x_i,z_i)' \sim N\left( \begin{bmatrix}0 \\ 0 \\ 0\end{bmatrix}, \begin{bmatrix} 1 & \sigma_{yx} & \sigma_{yz} \\ \sigma_{yx} & 1 & \rho \\ \sigma_{yz} & \rho & 1\\ \end{bmatrix} \right). \]

1.1 Conditional distribution of \(y_i|x_i,z_i\)

Basic multivariate normal properties lead to the following distribution for \(y_i\) given \(x_i\) and \(y_i\): \[ y_i \mid x_i,z_i \sim N\left(\frac{(\sigma_{xy} - \sigma_{yz}\rho)}{1-\rho^2}x + \frac{(\sigma_{yz} - \sigma_{xy}\rho)}{1-\rho^2}z, \; \frac{1 - \rho - \sigma_{xy}^2 - \sigma_{yz}^2 + 2\sigma_{xy}\sigma_{yz}\rho}{1 - \rho^2} \right), \] or \[ y_i|x_i,z_i = \beta x_i + \gamma z_i + \varepsilon_i \qquad \varepsilon_i \sim N(0,\tau^2), \] for \(\beta\), \(\gamma\) and \(\tau^2\) defined in the previous expression.

1.2 Conditional distribution of \(y_i|x_i,z_i\) when \(\rho=0\)

Notice that, when \(cor(x_i,z_i)=\rho=0\), it follows that \[ y_i|x_i,z_i = \sigma_{xy} x_i + \sigma_{zy} z_i + \varepsilon_i \qquad \varepsilon_i \sim N(0,1-(\sigma_{xy}^2+\sigma_{zy}^2)), \]

#install.packages("MASS")
library(MASS)
set.seed(123)
  
n     = 30
beta  = 0.75
gamma = 1.20
N     = 1000
rhos = seq(0.0,0.95,by=0.05)
coefs = matrix(0,length(rhos),N)
coefs1 = matrix(0,length(rhos),N)

j = 0
for (rho in rhos){
  Sigma = matrix(c(1,rho,rho, 1),2,2)
  j     = j + 1
  for (i in 1:N){
    draw = mvrnorm(n,c(0,0),Sigma)
    x = draw[,1]
    z = draw[,2]
    y = beta*x+gamma*z + rnorm(n)
    coefs[j,i] = lm(y~x-1)$coef[1]
    coefs1[j,i] = lm(y~x+z-1)$coef[1]
  }
}

2 Fitted model: \(y_i = \beta x_i + \gamma z_i + \epsilon_i\)

par(mfrow=c(1,1))
boxplot(t(coefs1),outline=FALSE,xlab="corr(x,z)",ylab="OLS estimate of beta",names=rhos)
title("Fitted model: y=beta*x+gamma*z+error")
abline(h=beta,col=2,lwd=2)
legend("topleft",legend=paste("Data generating process: y=",beta,"*x+",gamma,"*z+epsilon",sep=""))

3 Fitted model: \(y_i = \beta x_i + \epsilon_i\)

par(mfrow=c(1,1))
boxplot(t(coefs),outline=FALSE,xlab="corr(x,z)",ylab="OLS estimate of beta",names=rhos)
title("Fitted model: y=beta*x+error")
abline(h=beta,col=2,lwd=2)
legend("topleft",legend=paste("Data generating process: y=",beta,"*x+",gamma,"*z+epsilon",sep=""))