If you generalize directly the covariance to three variables, i.e.

Cov(X1,X2,X3) = E[ (X-E[X1])·(X-E[X2])·(X-E[X3]) ]

you cannot interpret the result univocally. For instance, if Cov(X1,X2,X3)>0 were large, you wouldn't know what is happening, as different combinations of signs (like +·+·+, +·-·-) contribute to the same result. Thus, you cannot distinguish whether X2 and X3 are both large when X1 is large, or X2 and X3 are both small when X1 is large, or X2 and X3 are either both large or both small when X1 is large.

Another argument against using that generalization to more than two variables is the following. For two variables, you have Cov(X,X)=Var(X), so it is plausible to interpret covariance as being related to variability. But for more variables, Cov(X,X,X) and so on are related to higher moments of X which are not interpreted as variability, but as symmetry (third moment), kurtosis (fourth moment), and so on.

(Where `and so on' means `I have no idea what the interpretation of the fifth moment is.)

Hello,

Denote X as a vector of n Random Variables i.e.

X = transpose([ X1,X2,..,Xn])

The covariance matrix of dimensions nxn is defined as 

Cxx = Expectation { (X-mX)*transpose(X-mX) }

where mX is the mean of X i.e. mX = Expectation { X }.

That's the only definition :-)

Yes, there is a direct generalization of the covariance concept to any dimension formulated with an ( at least) semidefinite positive matrix- the covariance matrix- An explanation  easy to understand is given by the above response by A. Bazzi.

Yes, there is. I agree with Ahmad Bazzi

Thank you very much for the answer, but I still wonder though.

The covariance matrix is very much the presentation of pair-wise covariances, not of all variables. For time series data, we have the co-integration analysis, but I wonder for panel data, would it be possible? 

are you  talking about partial covariance!

You might be thinking of multiple correlation. See

https://en.wikipedia.org/wiki/Multiple_correlation

If you generalize directly the covariance to three variables, i.e.

Cov(X1,X2,X3) = E[ (X-E[X1])·(X-E[X2])·(X-E[X3]) ]

you cannot interpret the result univocally. For instance, if Cov(X1,X2,X3)>0 were large, you wouldn't know what is happening, as different combinations of signs (like +·+·+, +·-·-) contribute to the same result. Thus, you cannot distinguish whether X2 and X3 are both large when X1 is large, or X2 and X3 are both small when X1 is large, or X2 and X3 are either both large or both small when X1 is large.

Another argument against using that generalization to more than two variables is the following. For two variables, you have Cov(X,X)=Var(X), so it is plausible to interpret covariance as being related to variability. But for more variables, Cov(X,X,X) and so on are related to higher moments of X which are not interpreted as variability, but as symmetry (third moment), kurtosis (fourth moment), and so on.

(Where `and so on' means `I have no idea what the interpretation of the fifth moment is.)

I do not think it is reasonable to calculate covariance of more than two variables, because the relationship among many variables cannot be defined clearly. However, if we have just only two variables, it is easy to think about the relationship between them. Hope can help you.

I worked on this three variable covariance problem and so far, I have not been successful. The formulas that I tried led to contradictions. 

use python numpy library

assume you have data u represented by 3 variables

u = [ [ x1, x2, x3, ..., xn], [y1, y2, y3, ..., yn], [z1,z2,z3,...,zn]]

numpy.cov(u)

it will produce 3*3 matrix, the diagonal is the auto-covariance of the three variables, and off-diagonal are cross-covariance of the three variable

you will find the matrix is positive semi-definitive

The formula for the covariance of three or more variables (if it exists) should reflect the measurement unit of each variable.

Also, how do you define the correlation coefficient of three or more variables?. Will it be between -1 and 1 ?.

yes, you can you use principal component analysis techniques on the numpy.cov generated matrix to find how the variable related to each other. for example singular value decomposition can be used to decompose the matrix into its eigenvalues and eigen-vectors. analysis can be made on these eigen value and eign-vectors to examine dependency of the analysed variables. I use this technique in seismic data analysis recorded by a geophone having three orthogonal sensors.

This might give a simple demonstration; http://www.itl.nist.gov/div898/handbook/pmc/section5/pmc541.htm