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

## Usage

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}$$.

## Value

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 $$U_n = X_n\beta + \epsilon_n,$$ 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 $$\Sigma \sim W(\kappa,E)$$ 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))
#>      [,1] [,2] [,3]
#> [1,]  1.0  0.5  0.2
#> [2,]  0.5  1.0  0.2
#> [3,]  0.2  0.2  2.0
### 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)
#>           [,1]        [,2]        [,3]
#> [1,] 1.2433933 0.599333881 0.303754609
#> [2,] 0.5993339 1.126533916 0.002703635
#> [3,] 0.3037546 0.002703635 2.054583290
apply(Sigma_draws, 1:2, stats::sd)
#>           [,1]      [,2]      [,3]
#> [1,] 0.1910950 0.1469450 0.1866532
#> [2,] 0.1469450 0.1481329 0.1588351
#> [3,] 0.1866532 0.1588351 0.3048582