18  Multiequationals linear models

library(tidyverse)
library(mvtnorm)
library(backports)
library(latex2exp)

Let’s consider a multivariate linear model, i.e. with $p>1$ in (Equation 14.3), then the model in matrix notation reads: $$\underset{n\times p}{\mathit{Y}}=\underset{n\times 1}{{\mathbf{J}}_{n,1}}\underset{1\times p}{{\mathit{a}}^{\mathrm{\top }}}+\underset{n\times k}{\mathit{X}}\underset{k\times p}{{\mathit{b}}^{\mathrm{\top }}}+\underset{n\times p}{\mathit{e}}$$

18.1 OLS estimate

As in the uni-variate case the optimal parameters are computed as: $$\begin{array}{rl}& {\mathbf{b}}^{\text{OLS}}=\mathbb{C}v(\mathit{Y},\mathit{X})\mathbb{C}v(\mathit{X}{)}^{-1}\\ & {\mathit{\alpha }}^{\text{OLS}}=\mathbb{E}\{\mathit{Y}\}-{\mathbf{b}}^{\text{OLS}}\mathbb{E}\{\mathit{X}\}\end{array}$$ And the variance covariance matrix of the residuals is computed as: $$\mathbf{\Sigma }=\mathbb{C}v(\mathit{e})=\mathbb{C}v(\mathit{Y})-{\mathbf{b}}^{\text{OLS}}\mathbb{C}v(\mathit{Y},\mathit{X})$$

18.1.1 Example

Let’s consider $n$-simulated observations of the explicative variables $\mathbf{X}$ drown from a multivariate normal, i.e. $\mathbf{X}\sim \mathcal{N}(\mathbb{E}\{X\},\mathbb{C}v\{X\})$, with parameters $$\mathbb{E}\{X\}=\left(\begin{array}{c}0.5\\ 0.5\\ 0.5\end{array}\right)\mathbb{C}v\{X\}=\left(\begin{array}{ccc}0.5& 0.2& 0.1\\ 0.2& 1.2& 0.1\\ 0.1& 0.1& 0.3\end{array}\right)$$ Let’s consider two dependent variables, hence $p=2$ and $k=3$. Let’s now simulate the $p\times k=6$ slopes parameter drown from a standard normal, i.e. for $j=1,\dots ,6$, $b_j\sim \mathcal{N}(0,1)$. The intercept parameters $\mathbf{a}$ are simulated drown from a uniform distribution in [0,1]. In the multivariate case $\mathbf{a}$ and $\mathbf{b}$ became matrices, i.e.  $$\underset{p\times k}{\mathbf{b}}=\left(\begin{array}{ccc}{b}_{1,1}& b_{1,2}& b_{1,3}\\ b_{2,1}& b_{2,2}& b_{2,3}\end{array}\right)\underset{p\times 1}{\mathit{\alpha }}=\left(\begin{array}{c}{\alpha }_1\\ {\alpha }_2\end{array}\right)$$

For $i=1,\dots ,n$, let’s consider a model of the form: $$\{\begin{array}{l}{Y}_{i,1}={\beta }_{0,1}+{\beta }_{1,1}{X}_{i,1}+{\beta }_{1,2}{X}_{i,2}+{\beta }_{1,k}{X}_{i,3}+u_{i,1}\\ Y_{i,2}={\beta }_{0,2}+{\beta }_{2,1}{X}_{i,1}+{\beta }_{2,2}{X}_{i,2}+{\beta }_{2,k}{X}_{i,3}+u_{i,2}\end{array}$$ where $u_{i,1}$ and $u_{i,2}$ are simulated from a multivariate normal random variables with true covariance matrix equal to: $$\mathbb{C}v\{\mathbf{u}\}=\left(\begin{array}{cc}0.55& 0.3\\ 0.3& 0.70\end{array}\right)$$ Hence, the procedure is structured as:

1. Simulate of the explanatory variables, the regression parameters and the residuals.
2. Simulate the perturbed $\tilde{\mathbf{Y}}$ (regression with errors).
3. Fit the regression parameters on the $\tilde{\mathbf{Y}}$.
   
4. Compute the fitted residuals from the prediction obtained with the parameters in step 3. and compute their variance covariance matrix.
5. Compare the results with the true parameters.

######################## Setup ########################
set.seed(1) # random seed 
n <- 500    # number of observations
p <- 2      # number of dependent variables 
k <- 3      # number of regressors 
# True regressor's mean 
true_e_x <- matrix(rep(0.5, k), ncol = 1)
# True regressor's covariance matrix 
true_cv_x <-  matrix(c(v_z1 = 0.5, cv_12 = 0.2, cv_13 = 0.1, 
                       cv_21 = 0.2, v_z2 = 1.2, cv_23 = 0.1, 
                       cv_31 = 0.1, cv_32 = 0.1, v_z3 = 0.3), 
                     nrow = k, byrow = FALSE)
# True covariance of the residuals 
true_cv_e <- matrix(c(0.55, 0.3, 0.3, 0.70), nrow = p)
##########################################################
# Generate a synthetic data set 
## Regressors  
X <- rmvnorm(n, true_e_x, true_cv_x) 
## Slope (Beta)
true_beta <- rnorm(p*k)
true_beta <- matrix(true_beta, ncol = k, byrow = TRUE) 
## Intercept (Alpha)
true_alpha <- runif(p, min = 0, max = 1)
true_alpha <- matrix(true_alpha, ncol = 1) 
## Matrix of 1 for matrix multiplication  
ones <- matrix(rep(1, n), ncol = 1)
## Fitted response variable 
Y <- ones %*% t(true_alpha) + X %*% t(true_beta)
## Simulated error 
eps <- rmvnorm(n, sigma = true_cv_e)
## Perturbed response variable 
Y_tilde <- Y + eps

# True Beta 
df_beta_true <- as_tibble(true_beta)
colnames(df_beta_true) <- paste0("$\\beta_", 1:ncol(df_beta_true), "$")
# Perturbed Beta (fitted)
fit_beta <- cov(Y_tilde, X) %*% solve(cov(X))
df_beta_pert <- as_tibble(fit_beta)
colnames(df_beta_pert) <- paste0("$\\beta_", 1:ncol(df_beta_pert), "$")
# True Alpha 
df_alpha_true <- as_tibble(t(true_alpha))
colnames(df_alpha_true) <- paste0("$\\alpha_", 1:ncol(df_alpha_true), "$")
# Perturbed Alpha (fitted)
## Perturbed mean  
e_y <- matrix(apply(Y_tilde, 2, mean), ncol = 1)
e_x <- matrix(apply(X, 2, mean), ncol = 1)
## Estimated Alpha (on perturbed data)
fit_alpha <- e_y - cov(Y_tilde, X) %*% solve(cov(X)) %*% e_x
df_alpha_pert <- as_tibble(t(fit_alpha))
colnames(df_alpha_pert) <- paste0("$\\alpha_", 1:ncol(df_alpha_pert), "$")

dplyr::bind_rows(
  dplyr::bind_cols(Type = "True", df_beta_true),
  dplyr::bind_cols(Type = "Fitted", df_beta_pert)
  ) %>%
  dplyr::mutate_if(is.numeric, format, digits = 4, scientific = FALSE) %>%
  knitr::kable(escape = FALSE)

dplyr::bind_rows(
  dplyr::bind_cols(Type = "True", df_alpha_true),
  dplyr::bind_cols(Type = "Fitted", df_alpha_pert)
  ) %>%
  dplyr::mutate_if(is.numeric, format, digits = 4, scientific = FALSE) %>%
  knitr::kable(escape = FALSE)

Table 18.1: Fitted parameters

Type	${\beta }_1$	${\beta }_2$	${\beta }_3$
True	0.8500	-0.9253	0.8936
True	-0.9410	0.5390	-0.1820
Fitted	0.8457	-0.8699	0.9396
Fitted	-0.9532	0.5518	-0.1804

Type	${\alpha }_1$	${\alpha }_2$
True	0.8137	0.8068
Fitted	0.7942	0.7423