Estimation of Multivariate Regression Models

Least Squares Estimation

• Ordinary Least Squares

• Covariance-Weighted Least Squares

• Error Covariance Estimation

• Feasible Generalized Least Squares

• Panel Corrected Standard Errors

Ordinary Least Squares

When you fit multivariate linear regression models using mvregress, you can use the optional name-value pair 'algorithm','cwls' to choose least squares estimation. In this case, by default, mvregress returns ordinary least squares (OLS) estimates using $$\Sigma =I_d$$. Alternatively, if you specify a covariance matrix for weighting, you can return covariance-weighted least squares (CWLS) estimates. If you combine OLS and CWLS, you can get feasible generalized least squares (FGLS) estimates.

The OLS estimate for the coefficient vector is the vector $$b$$ that minimizes

$${\displaystyle \sum _{i=1}^n{\left(y_i-X_{i}b\right)}^{\prime }}\left(y_i-X_{i}b\right).$$

Let $$y$$ denote the nd-by-1 vector of stacked d-dimensional responses, and $$X$$ denote the nd-by-K matrix of stacked design matrices. The K-by-1 vector of OLS regression coefficient estimates is

$$b_{OLS}={\left(X^{\prime }X\right)}^{-1}{X}^{\prime }y.$$

This is the first mvregress output.

Given $$\Sigma =I_d$$ (the mvregress OLS default), the variance-covariance matrix of the OLS estimates is

$$V(b_{OLS})={(X^{\prime }X)}^{-1}.$$

This is the fourth mvregress output. The standard errors of the OLS regression coefficients are the square root of the diagonal of this variance-covariance matrix.

If your data is not scaled such that $$\Sigma ={\sigma }^2{I}_d$$, then you can multiply the mvregress variance-covariance matrix by the mean squared error (MSE), an unbiased estimate of $${\sigma }^2$$. To compute the MSE, return the n-by-d matrix of residuals, $$E$$ (the third mvregress output). Then,

$$\text{MSE}=\frac{{\displaystyle \sum _{i=1}^n{e}_i{e}_i{}^{\prime }}}{n-K},$$

where $$e_i=(y_i-X_i\beta {)}^{\prime }$$ is the ith row of $$E$$.

Covariance-Weighted Least Squares

For most multivariate problems, an identity error covariance matrix is insufficient, and leads to inefficient or biased standard error estimates. You can specify a matrix for CWLS estimation using the optional name-value pair argument covar0, for example, an invertible d-by-d matrix named $$C_0$$. Usually, $$C_0$$ is a diagonal matrix such that the inverse matrix $$C_0^{-1}$$ contains weights for each dimension to model heteroscedasticity. However, $$C_0$$ can also be a nondiagonal matrix that models correlation.

Given $$C_0$$, the CWLS solution is the vector $$b$$ that minimizes

$${\displaystyle \sum _{i=1}^n{\left(y_i-X_{i}b\right)}^{\prime }}{C}_0\left(y_i-X_{i}b\right).$$

In this case, the K-by-1 vector of CWLS regression coefficient estimates is

$$b_{CWLS}={\left(X^{\prime }{\left(I_n\otimes C_0\right)}^{-1}X\right)}^{-1}{X}^{\prime }{\left(I_n\otimes C_0\right)}^{-1}y.$$

This is the first mvregress output.

If $$\Sigma =C_0$$, this is the generalized least squares (GLS) solution. The corresponding variance-covariance matrix of the CWLS estimates is

$$V(b_{CWLS})={\left(X\text{'}{\left(I_n\otimes C_0\right)}^{-1}X\right)}^{-1}.$$

This is the fourth mvregress output. The standard errors of the CWLS regression coefficients are the square root of the diagonal of this variance-covariance matrix.

If you only know the error covariance matrix up to a proportion, that is, $$\Sigma ={\sigma }^2{C}_0$$, you can multiply the mvregress variance-covariance matrix by the MSE, as described in Ordinary Least Squares.

Error Covariance Estimation

Regardless of which least squares method you use, the estimate for the error variance-covariance matrix is

$$\widehat{\Sigma }=\left(\begin{array}{cccc}{\widehat{\sigma }}_1^2& {\widehat{\sigma }}_{12}& \cdots & {\widehat{\sigma }}_{1d}\\ {\widehat{\sigma }}_{12}& {\widehat{\sigma }}_2^2& \cdots & {\widehat{\sigma }}_{2d}\\ ⋮& ⋮& \ddots & ⋮\\ {\widehat{\sigma }}_{1d}& {\widehat{\sigma }}_{2d}& \cdots & {\widehat{\sigma }}_d^2\end{array}\right)=\frac{E^{\prime }E}{n},$$

where $$E$$ is the n-by-d matrix of residuals. The ith row of $$E$$ is $$e_i={\left(y_i-X_{i}b\right)}^{\prime }.$$

The error covariance estimate, $$\widehat{\Sigma }$$, is the second mvregress output, and the matrix of residuals, $$E$$, is the third output. If you specify the optional name-value pair 'covtype','diagonal', then mvregress returns $$\widehat{\Sigma }$$ with zeros in the off-diagonal entries,

$$\widehat{\Sigma }=\left(\begin{array}{ccc}{\widehat{\sigma }}_1^2& & 0\\ & \ddots & \\ 0& & {\widehat{\sigma }}_d^2\end{array}\right).$$

Feasible Generalized Least Squares

The generalized least squares estimate is the CWLS estimate with a known covariance matrix. That is, given $$\Sigma$$ is known, the GLS solution is

$$b_{GLS}={\left(X^{\prime }{(I_n\otimes \Sigma )}^{-1}X\right)}^{-1}{X}^{\prime }{(I_n\otimes \Sigma )}^{-1}y,$$

with variance-covariance matrix

$$V(b_{GLS})={\left(X^{\prime }{(I_n\otimes \Sigma )}^{-1}X\right)}^{-1}.$$

In most cases, the error covariance is unknown. The feasible generalized least squares (FGLS) estimate uses $$\widehat{\Sigma }$$ in place of $$\Sigma$$. You can obtain two-step FGLS estimates as follows:

1. Perform OLS regression, and return an estimate $$\widehat{\Sigma }$$.

2. Perform CWLS regression, using $$C_0=\widehat{\Sigma }$$.

You can also iterate between these two steps until convergence is reached.

For some data, the OLS estimate $$\widehat{\Sigma }$$ is positive semidefinite, and has no unique inverse. In this case, you cannot get the FGLS estimate using mvregress. As an alternative, you can use lscov, which uses a generalized inverse to return weighted least squares solutions for positive semidefinite covariance matrices.

Panel Corrected Standard Errors

An alternative to FGLS is to use OLS coefficient estimates (which are consistent) and make a standard error correction to improve efficiency. One such standard error adjustment—which does not require inversion of the covariance matrix—is panel corrected standard errors (PCSE) [1]. The panel corrected variance-covariance matrix for OLS estimates is

$$V_{pcse}(b_{OLS})={(X^{\prime }X)}^{-1}{X}^{\prime }(I_n\otimes \Sigma )X{(X^{\prime }X)}^{-1}.$$

The PCSE are the square root of the diagonal of this variance-covariance matrix. Fixed Effects Panel Model with Concurrent Correlation illustrates PCSE computation.

Maximum Likelihood Estimation

• Maximum Likelihood Estimates

• Standard Errors

Maximum Likelihood Estimates

The default estimation algorithm used by mvregress is maximum likelihood estimation (MLE). The loglikelihood function for the multivariate linear regression model is

$$\begin{array}{l}\log L(\beta ,\Sigma |y,X)=\frac{1}{2}nd\log (2\pi )+\frac{1}{2}n\log (\det (\Sigma ))\\ \text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}+\frac{1}{2}{\displaystyle \sum _{i=1}^n{\left(y_i-X_i\beta \right)}^{\prime }}{\Sigma }^{-1}\left(y_i-X_i\beta \right).\\ \text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\end{array}$$

The MLEs for $$\beta$$ and $$\Sigma$$ are the values that maximize the loglikelihood objective function.

mvregress finds the MLEs using an iterative two-stage algorithm. At iteration m + 1, the estimates are

$$b_{MLE}^{(m+1)}={\left(X^{\prime }{\left(I_n\otimes {\Sigma }^{(m)}\right)}^{-1}X\right)}^{-1}{X}^{\prime }{\left(I_n\otimes {\Sigma }^{(m)}\right)}^{-1}y$$

and

$${\widehat{\Sigma }}^{(m+1)}=\frac{1}{n}{\displaystyle \sum _{i=1}^n\left(y_i-X_i{b}_{MLE}^{(m+1)}\right)}{\left(y_i-X_i{b}_{MLE}^{(m+1)}\right)}^{\prime }.$$

The algorithm terminates when the changes in the coefficient estimates and loglikelihood objective function are less than a specified tolerance, or when the specified maximum number of iterations is reached. The optional name-value pair arguments for changing these convergence criteria are tolbeta, tolobj, and maxiter, respectively.

Standard Errors

The variance-covariance matrix of the MLEs is an optional mvregress output. By default, mvregress returns the variance-covariance matrix for only the regression coefficients, but you can also get the variance-covariance matrix of $$\widehat{\Sigma }$$ using the optional name-value pair 'vartype','full'. In this case, mvregress returns the variance-covariance matrix for all K regression coefficients, and d or d(d + 1)/2 covariance terms (depending on whether the error covariance is diagonal or full).

By default, the variance-covariance matrix is the inverse of the observed Fisher information matrix (the 'hessian' option). You can request the expected Fisher information matrix using the optional name-value pair 'vartype','fisher'. Provided there is no missing response data, the observed and expected Fisher information matrices are the same. If response data is missing, the observed Fisher information accounts for the added uncertainty due to the missing values, whereas the expected Fisher information matrix does not.

The variance-covariance matrix for the regression coefficient MLEs is

$$V(b_{MLE})={\left(X^{\prime }{(I_n\otimes \widehat{\Sigma })}^{-1}X\right)}^{-1},$$

evaluated at the MLE of the error covariance matrix. This is the fourth mvregress output. The standard errors of the MLEs are the square root of the diagonal of this variance-covariance matrix.

For $$\widehat{\Sigma }$$, let $$\theta$$ denote the vector of parameters in the estimated error variance-covariance matrix. For example, if d = 2, then:

• If the estimated covariance matrix is diagonal, then $$\theta =({\widehat{\sigma }}_1^2,{\widehat{\sigma }}_2^2)$$.

• If the estimated covariance matrix is full, then $$\theta =({\widehat{\sigma }}_1^2,{\widehat{\sigma }}_{12},{\widehat{\sigma }}_2^2)$$.

The Fisher information matrix for $$\theta$$, $$I(\theta )$$, has elements

$$I{(\theta )}_{u,v}=\frac{1}{2}tr\left({\widehat{\Sigma }}^{-1}\text{}\frac{\partial \widehat{\Sigma }}{\partial {\theta }_u}{\widehat{\Sigma }}^{-1}\text{}\frac{\partial \widehat{\Sigma }}{\partial {\theta }_v}\right),u,v=1,\dots ,n_{\theta },$$

where $$n_{\theta }$$ is the length of $$\theta$$ (either d or d(d + 1)/2). The resulting variance-covariance matrix is

$$V(\theta )=I{(\theta )}^{-1}.$$

When you request the full variance-covariance matrix, mvregress returns (as the fourth output) the block diagonal matrix

$$\left(\begin{array}{cc}V(b_{MLE})& 0\\ 0& V(\theta )\end{array}\right).$$

Missing Response Data

• Expectation/Conditional Maximization

• Observed Information Matrix

Expectation/Conditional Maximization

If any response values are missing, indicated by NaN, mvregress uses an expectation/conditional maximization (ECM) algorithm for estimation (if enough data is available). In this case, the algorithm is iterative for both least squares and maximum likelihood estimation. During each iteration, mvregress imputes missing response values using their conditional expectation.

Consider organizing the data so that the joint distribution of the missing and observed responses, denoted $$\tilde{y}$$ and $$y$$ respectively, can be written as

$$\left(\begin{array}{l}\tilde{y}\\ y\end{array}\right)\sim MVN\left\{\left(\begin{array}{l}\tilde{X}\beta \\ X\beta \end{array}\right),\left(\begin{array}{cc}{\Sigma }_{\tilde{y}}& {\Sigma }_{\tilde{y}}{}_y\\ {\Sigma }_{y\tilde{y}}& {\Sigma }_y\end{array}\right)\right\}.$$

Using properties of the multivariate normal distribution, the conditional expectation of the missing responses given the observed responses is

$$\text{E}\left(\tilde{y}|y\right)=\tilde{X}\beta +{\Sigma }_{\tilde{y}y}{\Sigma }_y^{-1}(y-X\beta ).$$

Also, the variance-covariance matrix of the conditional distribution is

$$\text{COV}(\tilde{y}|y)={\Sigma }_{\tilde{y}}-{\Sigma }_{\tilde{y}y}{\Sigma }_y^{-1}{\Sigma }_{y\tilde{y}}.$$

At each iteration of the ECM algorithm, mvregress uses the parameter values from the previous iteration to:

• Update the regression coefficients using the combined vector of observed responses and conditional expectations of missing responses.

• Update the variance-covariance matrix, adjusting for missing responses using the variance-covariance matrix of the conditional distribution.

Finally, the residuals that mvregress returns for missing responses are the difference between the conditional expectation and the fitted value, both evaluated at the final parameter estimates.

If you prefer to ignore any observations that have missing response values, use the name-value pair 'algorithm','mvn'. Note that mvregress always ignores observations that have missing predictor values.

Observed Information Matrix

By default, mvregress uses the observed Fisher information matrix (the 'hessian' option) to compute the variance-covariance matrix of the regression parameters. This accounts for the additional uncertainty due to missing response values.

The observed information matrix includes contributions from only the observed responses. That is, the observed Fisher information matrix for the parameters in the error variance-covariance matrix has elements

$$I{(\theta )}_{u,v}=\frac{1}{2}{\displaystyle \sum _{i=1}^{n}tr\left({\widehat{\Sigma }}_{{}^i\text{}}^{-1}\frac{\partial {\widehat{\Sigma }}_i}{\partial {\theta }_u}{\widehat{\Sigma }}_{{}^i\text{}}^{-1}\text{}\frac{\partial {\widehat{\Sigma }}_i}{\partial {\theta }_v}\right)},u,v=1,\dots ,n_{\theta },$$

where $${\widehat{\Sigma }}_i$$ is the subset of $$\widehat{\Sigma }$$ corresponding to the observed responses in $$y_i.$$

For example, if d = 3, but $$y_{i2}$$ is missing, then

$${\widehat{\Sigma }}_i=\left(\begin{array}{cc}{\widehat{\sigma }}_1^2& {\widehat{\sigma }}_{13}\\ {\widehat{\sigma }}_{13}& {\widehat{\sigma }}_3^2\end{array}\right).$$

The observed Fisher information for the regression coefficients has similar contributions from the design and covariance matrices.

References

[1] Beck, N. and J. N. Katz. "What to Do (and Not to Do) with Time-Series-Cross-Section Data in Comparative Politics." American Political Science Review, Vol. 89, No. 3, pp. 634–647, 1995.