update_Sigma function

Update error term covariance matrix of multiple linear regression

This function updates the error term covariance matrix of a multiple linear regression.

update_Sigma(kappa, E, N, S)

Arguments

• kappa: The degrees of freedom (a natural number greater than J-1) of the Inverse Wishart prior for Sigma. Per default, kappa = J + 1.
• E: The scale matrix of dimension J-1 x J-1 of the Inverse Wishart prior for Sigma. Per default, E = diag(J - 1).
• N: The draw size.
• S: A matrix, the sum over the outer products of the residuals $(\epsilon_n)_{n=1,\dots,N}$.

Returns

A matrix, a draw from the Inverse Wishart posterior distribution of the error term covariance matrix in a multiple linear regression.

Details

This function draws from the posterior distribution of the covariance matrix $\Sigma$ in the linear utility equation

where $U_n$ is the (latent, but here assumed to be known) utility vector of decider $n = 1,\dots,N$, $X_n$

is the design matrix build from the choice characteristics faced by $n$, $\beta$ is the coefficient vector, and $\epsilon_n$ is the error term assumed to be normally distributed with mean $0$ and unknown covariance matrix $\Sigma$. A priori we assume the (conjugate) Inverse Wishart distribution

with $\kappa$ degrees of freedom and scale matrix $E$. The posterior for $\Sigma$ is the Inverted Wishart distribution with $\kappa + N$ degrees of freedom and scale matrix $E^{-1}+S$, where $S = \sum_{n=1}^{N} \epsilon_n \epsilon_n'$ is the sum over the outer products of the residuals $(\epsilon_n = U_n - X_n\beta)_n$.

Examples

### true error term covariance matrix
(Sigma_true <- matrix(c(1,0.5,0.2,0.5,1,0.2,0.2,0.2,2), ncol=3))
### coefficient vector
beta <- matrix(c(-1,1), ncol=1)
### draw data
N <- 100
X <- replicate(N, matrix(rnorm(6), ncol=2), simplify = FALSE)
eps <- replicate(N, rmvnorm(mu = c(0,0,0), Sigma = Sigma_true), simplify = FALSE)
U <- mapply(function(X, eps) X %*% beta + eps, X, eps, SIMPLIFY = FALSE)
### prior parameters for covariance matrix
kappa <- 4
E <- diag(3)
### draw from posterior of coefficient vector
outer_prod <- function(X, U) (U - X %*% beta) %*% t(U - X %*% beta)
S <- Reduce(`+`, mapply(outer_prod, X, U, SIMPLIFY = FALSE))
Sigma_draws <- replicate(100, update_Sigma(kappa, E, N, S))
apply(Sigma_draws, 1:2, mean)
apply(Sigma_draws, 1:2, stats::sd)