Skip to contents

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