# Performing 3 multivariate linear regressions at once

I have 3 variables X, Y and Z. I want to perform 3 OLS regressions: X dependent on Y and Z, Y dependent on X and Z, and Z dependent on X and Y.

Instead of doing the 3 of them sepparately, I want to know how can I do them in a single go via:

$$A = B \beta + \epsilon$$

Where $$A$$ is an n x 3 matrix (n observations of the 3 variables)

I asked ChatGPT and it suggested creating a $$B$$ matrix with shape (n, 6) (or 7 if we account for the intercept) where the columns would be the variables Y, Z, X, Z, X, Y and then regress it to get a parameters matrix beta with shape (6, 3) (or (7, 3) with the intercept). But after that I'm not sure how to interpret this resulting matrix or if it even makes sense (CGPT has definitely stopped making sense when talking about this resulting matrix and how to interpret it).

Another strange thing is that even if this makes sense and I'm doing it correctly (which I doubt), I'm not getting the results that I'd expect. If I run a regression of X on Y and Z only, I get the parameters

array([1.84477116, 0.74949417, 0.46818174])


But if I run all the regressions at the same time I get the matrix:

array([[ 1.32178002e-09,  9.02019792e-10, -2.18881269e-09],
[ 9.83568782e-11,  5.00000000e-01, -1.72848402e-10],
[ 5.98987526e-12,  2.35829134e-11,  5.00000000e-01],
[ 5.00000000e-01, -9.13527032e-11,  3.39086093e-10],
[ 1.90638616e-11,  3.31752403e-11,  5.00000000e-01],
[ 5.00000000e-01, -4.81339413e-11,  3.72506470e-10],
[ 1.19651844e-10,  5.00000000e-01, -1.60753189e-10]])


Which doesn't show the results of the single OLS anywhere.

Here is the code I'm using:

y = df.to_numpy() #df is a dataframe with 3 columns and n rows (observations)
x = add_constant(pd.DataFrame(data=[df.iloc[:,1], df.iloc[:,2], df.iloc[:,0], df.iloc[:,2], df.iloc[:,0], df.iloc[:,1]]).T).to_numpy() #it adds a new column of ones for the intercept
cov = np.dot(x.T, x)
inv = np.linalg.pinv(cov)
H = np.dot(inv, x.T)
betas = np.dot(H, y)

• What do you aim to achieve with a single equation that can't be achieved with three? Your question does not mention anything in particular that would benefit from having a single equation Mar 5, 2023 at 15:14
• Hi Firebug, it's 70% simply for learning, and 30% to streamline the process of calculating these parameters. In this case it's only 3 variables but I plan to use it for 10-20 variables, that would mean running 10-20 OLS regressions instead of just running it once and selecting the right elements from the output matrix Mar 5, 2023 at 15:20
• It would be invalid to apply standard regression methods here due to the inherent inconsistency: each regression views the explanatory variables as being known or measured without error and attributes all residual variation in the response variable to randomness. This inherent contradiction makes it difficult, or not impossible, to determine what you need to accomplish. Since you say it's for learning, I have two strong suggestions: (1) avoid ChatGPT because it will deceive you and (2) investigate the conditions needed to apply regression models.
– whuber
Mar 5, 2023 at 16:00
• @Hiperfly it doesn't simplify anything. If anything, you could concatenate the coefficients and call it a single equation, achieving the same Mar 5, 2023 at 17:11
• Also, you are not doing three multivariate regression, but three multiple regressions Mar 7, 2023 at 10:23

Assume $$X,Y,Z$$ are centered and organized as follows: $$A=\begin{pmatrix} \mid &\mid &\mid \\ X &Y &Z \\ \mid &\mid &\mid \end{pmatrix}$$

And

$$A = A B + E$$

Where

$$B = \begin{pmatrix} 0 &\beta_{X\to Y} &\beta_{X\to Z}\\ \beta_{Y\to X} &0 & \beta_{Y\to Z}\\ \beta_{Z\to X} &\beta_{Z\to Y} & 0 \end{pmatrix}$$

The general problem statement is

$$\cases{ \min_B\|E\|_2^{2}=\min_B\|A\cdot (\mathbb{I}-B)\|_2^{2}\\ \mathrm{diag}(B)=0 }$$

Using Lagrange multipliers:

$$f = \|A\cdot (\mathbb{I}-B)\|_2^{2}-\lambda^\top \cdot \mathrm{diag}(B)$$

$$\frac{\partial f}{\partial B} = -(2\cdot A^\top \cdot A\cdot (\mathbb{I}-B)+\mathrm{diag}(\lambda))=0$$

$$2\cdot A^\top \cdot A\cdot (B-\mathbb{I})=\mathrm{diag}(\lambda)$$

$$B=\frac12(A^\top \cdot A)^{-1}\cdot \mathrm{diag}(\lambda)+\mathbb{I}$$

Using our equality constraint:

$$\mathrm{diag}(B)=0=\mathrm{diag}\left(\frac12(A^\top \cdot A)^{-1}\cdot \mathrm{diag}(\lambda)+\mathbb{I}\right)=\\ \frac12\mathrm{diag}\left((A^\top \cdot A)^{-1}\right)\cdot\mathrm{diag}(\lambda)+\mathrm{diag}(\mathbb{I})\\ \therefore \lambda_i=\frac{-2}{(A^\top \cdot A)^{-1}_{ii}} \therefore \mathrm{diag}(\lambda)=-2 \left( \mathbb I \odot (A^\top \cdot A)^{-1}\right)^{-1}$$

Plugging it back:

$$B=\frac12(A^\top \cdot A)^{-1}\cdot \mathrm{diag}(\lambda)+\mathbb{I}\\ =\frac12(A^\top \cdot A)^{-1}\cdot -2 \left( \mathbb I \odot (A^\top \cdot A)^{-1}\right)^{-1}+\mathbb{I}\\ B=\mathbb{I}-(A^\top \cdot A)^{-1}\cdot \left( \mathbb I \odot (A^\top \cdot A)^{-1}\right)^{-1}\\$$

Call $$n\Sigma = A^\top \cdot A,\Omega = \Sigma^{-1}, D_\Omega = \mathrm{diag}(\Omega)$$

$$B = \mathbb{I} - \Omega\cdot D_\Omega^{-1}$$

Compare $$B$$ with the partial correlation matrix, $$R=2\mathbb{I}-D_\Omega^{-1/2}\Omega\cdot D_\Omega^{-1/2}$$ for some intuition

You can check it is true with the following R code:

# generate X, Y, Z as a 100x3 matrix
A <- matrix(rnorm(300), ncol=3)
A <- scale(A)
# generate a 3x3 mixing matrix M
M <- matrix(rnorm(9), ncol=3)

# generate a 100x3 matrix of observed data
B <- A %*% M

# perform three linear regressions, one for each column of B from the other two
# columns of B
lm1 = lm(B[,1] ~ 0 + B[,2] + B[,3])
lm2 = lm(B[,2] ~ 0 + B[,1] + B[,3])
lm3 = lm(B[,3] ~ 0 + B[,1] + B[,2])

S = solve(t(B) %*% B)
S = diag(1, 3, 3) - S %*% diag(1/diag(S))

coef(lm1)
#    B[, 2]     B[, 3]
#-0.7249205 -0.9375626
coef(lm2)
> coef(lm2)
#    B[, 1]    B[, 3]
# -1.351588 -1.280748
coef(lm3)
#     B[, 1]     B[, 2]
# -1.0485396 -0.7682355
S
#            [,1]      [,2]       [,3]
# [1,]  0.0000000 -1.351588 -1.0485396
# [2,] -0.7249205  0.000000 -0.7682355
# [3,] -0.9375626 -1.280748  0.0000000


I released a package that can help with this at https://github.com/bhvieira/avaols. You can install it simply doing (requires devtools)

# install.packages("devtools")
devtools::install_github("bhvieira/avaols")


Then you can simply do

library(avaols)
# generate X, Y, Z as a 100x3 matrix
A = matrix(rnorm(300), ncol=3)
# generate a 3x3 mixing matrix M
M <- matrix(rnorm(9), ncol=3)
# generate a 100x3 matrix of observed data
B <- data.frame(A %*% M)

# fit avaols object
obj = avaols(B)
coef(obj)
#                     X1         X2           X3
# Intercept -0.009905788 -0.1840545 2.097737e-02
# X1         0.000000000 -3.1097798 1.219195e+00
# X2        -0.192439394  0.0000000 2.372571e-01
# X3         0.782301262  2.4601173 1.110223e-16

# compare with three linear regressions
lm1 = lm(X1 ~ ., data = B)
lm2 = lm(X2 ~ ., data = B)
lm3 = lm(X3 ~ ., data = B)

coef(lm1)
#  (Intercept)           X2           X3
# -0.009905788 -0.192439394  0.782301262
coef(lm2)
# (Intercept)          X1          X3
#  -0.1840545  -3.1097798   2.4601173
coef(lm3)
# (Intercept)          X1          X2
#  0.02097737  1.21919453  0.23725706

• can this includethe intercept? Mar 14, 2023 at 20:37
• @onyambu it can be included afterwards, as in linear regression Mar 15, 2023 at 8:14
• @onyambu I've released a library in R that can help with that, here: github.com/bhvieira/avaols Mar 15, 2023 at 13:16

There is another way to obtain these coefficients in a single regression. It consists in having six variables, one for each coefficient we want to estimate. If all you want are the OLS coefficients, then the idea is that we vectorize the response, and then place the correct values in each of the six columns so they correspond to the right problem. The rest of the elements should be zero, so they cancel the respective coefficient.

See the example below:

# generate X, Y, Z as a 100x3 matrix
A = matrix(rnorm(300), ncol=3)
# generate a 3x3 mixing matrix M
M <- matrix(rnorm(9), ncol=3)
# generate a 100x3 matrix of observed data
B <- data.frame(A %*% M)

mu = colMeans(B)
A = B
B = scale(B, center=TRUE, scale=FALSE)
vecA = Reduce(c, B)
vecA12 = vecA13 = vecA23 = vecA21 = vecA31 = vecA32 = vecA * 0
vecA12[101:200] = vecA[1:100]
vecA13[201:300] = vecA[1:100]
vecA23[201:300] = vecA[101:200]
vecA21[1:100] = vecA[101:200]
vecA31[1:100] = vecA[201:300]
vecA32[101:200] = vecA[201:300]
coef(lm(vecA ~ 0 + vecA12 + vecA13 + vecA23 + vecA21 + vecA31 + vecA32))
vecA12      vecA13      vecA23      vecA21      vecA31      vecA32
#  0.09034396  0.67427978 -1.03436202  0.17768726  0.45582021 -0.35552393


Compare with:

lm1 = lm(X1 ~ 0 + ., data = data.frame(B))
lm2 = lm(X2 ~ 0 + ., data = data.frame(B))
lm3 = lm(X3 ~ 0 + ., data = data.frame(B))
coef(lm1)
#       X2        X3
#0.1776873 0.4558202
coef(lm2)
#         X1          X3
# 0.09034396 -0.35552393
coef(lm3)
#        X1         X2
# 0.6742798 -1.0343620


Note that this only works for getting the coefficients: if you want the test-statistics, covariance matrix of coefficients, or confidence intervals, these will obviously differ from the individual regressions: the assumption of homoscedasticity is obviously violated.

In your solution, you have 7 regressors, which are Y, Z, X, Z, X, Y, which are all regressed onto X,Y,Z. In the final matrix, each regressor has a beta (= rows in the final array) onto each target (= column in the final array). This is clearly not what you want.

I would propose looking into regression theory and figuring this out on your own instead of relying on ChatGPT. The answer you got is very wrong and the code quality is also pretty bad. From the way you posed your question it also appears as if you have some basic misconceptions about how linear regression works, so diving into the theory might be a good idea. If you are a visual learner, I can highly recommend the ritvikmath YouTube channel, which is the place I usually go to if I need to get an overview over a topic. However, there is a plethora of high quality material on OLS, so this is just one suggestion.

Lastly, I would like to add one point: I have personally never heard of any scenario like the one you're describing and feel like it's a very uncommon use case to the extent where I would say it's not really valid. Could you maybe explain a bit about what you're trying to achieve?

• Also, it's probably easiest to streamline this process using a for loop over all your target variables and creating temporary DataFrames within the loop. I doubt the solution with one big X matrix works for all equations, at least I fail to see what this would look like. Mar 5, 2023 at 15:32
• Hi David, thanks for the suggestions, that's a great channel and I remember using it to learn about time series forecasting. In this case I have the prices of 3 assets X, Y and Z. I want to describe each one of them through the other 2 (for more clarification: log prices and assuming group cointegration, that's why I'm using prices and not returns). In the real case I will actually have around 10-20 assets, so I'd rather just compute the parameters once and select the right elements from the resulting matrix, instead of running 10-20 multivariate regressions. Mar 5, 2023 at 15:33
• Okay, I'm not familiar with this use case at all, sorry. However, you could probably start by doing something like for column in df:, which will iterate over the name of each column of your initial dataframe. Within the loop create a dataframe X = df.drop([column]) to create your X matrix and Y = df[column] for your Y. From there, you can just treat it as a simple univariate OLS problem and print the resulting matrix in each loop. :) Mar 5, 2023 at 15:42
• Also, would you mind explaining why do you think the code quality is pretty bad? I'm genuily asking, I thought the code is simple enough to not look bad at all. I just divided the steps of the $\beta$ function to make it more clear. Thanks! Mar 5, 2023 at 15:42
• There's two things that stood out to me from the code: The add_constant() line is very long, which reduces readability of the code and is somewhat frowned upon by Python aficionados. Looking at the code again, I don't see how you could improve it in your case though, so maybe it was more of a reflex from my side. Secondly, you explicitly compute the inverse, which is very costly and should be avoided if possible. np.linalg.lstsq should be an alternative that would probably be way more efficient (and require less manual code, which is also nice). Mar 5, 2023 at 15:51