How can I calculate standardized central cross-moments for 2 time-series?

The 4th-order standardized central moment, kurtosis, is;

 kurt = mean((x-mean(x)).^4)./mean(x.^2).^2; 

It should be something like the kurt, for x and y.

I also want to know the 3rd order cross-moment (something like skewness).

I am trying to write a code in MATLAB.

  • $\begingroup$ Is this an attempt at LaTeX? Or are you trying to enter code? If the latter, what language are you using? $\endgroup$
    – wolfies
    Oct 18, 2014 at 16:16
  • $\begingroup$ @wolfies I'm trying to write a code in Matlab. $\endgroup$
    – Rash
    Oct 18, 2014 at 16:22
  • $\begingroup$ You included the tag for [discrete-data]. Are your data discrete? (You don't mention anything about that in the question.) Is that important for your situation in some way? $\endgroup$ Oct 18, 2014 at 17:21
  • $\begingroup$ @gung, My mistake, My data is discrete, in matrix x, and I want to write the codes for calculation. $\endgroup$
    – Rash
    Oct 18, 2014 at 17:34

1 Answer 1


You are in essence looking for a multivariate measure of skew and kurtosis. There are many. I would start with the most establish ones, which are the multivariate skew and kurtosis measures of Mardia 1977 [0].

It seems to me you are more asking for an implementation than about the measures themselves. I don t know of any matlab implementation, but the R code below (from the R library psych) should be fairly easy to translate in matlab:

mardia <- function(x, na.rm=TRUE, plot=TRUE) {
  cl <- match.call()
  x <- as.matrix(x)     # in case it was a dataframe
  if(na.rm) x <- na.omit(x)
  n <- dim(x)[1]
  p <- dim(x)[2]
  x <- scale(x,scale=FALSE)  # zero center
  S <- cov(x)
  S.inv <- solve(S)
  D <- x %*% S.inv %*% t(x)
  b1p <- sum(D^3)/n^2
  b2p <- tr(D^2)/n 
  chi.df <- p*(p+1)*(p+2)/6
  k <- (p+1)*(n+1)*(n+3)/(n*((n+1)*(p+1) -6))

  small.skew <- n*k*b1p/6
  M.skew <- n*b1p/6
  M.kurt <- (b2p - p * (p+2))*sqrt(n/(8*p*(p+2)))
  p.skew <- 1-pchisq(M.skew,chi.df)
  p.small <- 1 - pchisq(small.skew,chi.df)
  p.kurt <- 2*(1- pnorm(abs(M.kurt)))
  d = sqrt(diag(D))
  if(plot) {qqnorm(d)
  results <- list(n.obs=n, n.var=p, b1p=b1p,b2p= b2p, skew=M.skew, 
                  small.skew=small.skew, p.skew=p.skew, p.small=p.small, 
                  kurtosis=M.kurt, p.kurt=p.kurt, d=d, Call=cl)
  class(results) <- c("psych", "mardia")
  • [0] K.V. Mardia (1970). Measures of multivariate skewness and kurtosis with applications. Biometrika, 57(3):pp. 519-30, 1970.
  • 2
    $\begingroup$ I couldn't understand much about R code, but the paper was very useful, Thanks. $\endgroup$
    – Rash
    Oct 18, 2014 at 16:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.