Hi, when I wanted to use lrvar() I discovered it's only for single series. So I ported some code to gretl to estimate the long-run covariance matrix for multivariate series (along with a helper function autocovar()). I haven't tested it yet, probably some bugs remaining. Actually I wanted to ask whether it would make sense to allow the multivariate case in built-in lrvar()? thanks, sven ----------- function matrix autocovar(matrix indatamat, int lag[0], bool demeaned[0]) /* computes the autocovariance at lag "lag" for a multivariate time series returns a NxN covariance matrix */ N = cols(indatamat) T = rows(indatamat) if demeaned=0 matrix datademeaned = cdemean(indatamat) else matrix datademeaned = indatamat endif matrix result = transp(datademeaned[lag+1:T,])*datademeaned[1:T-lag,] matrix result /= T return result end function function matrix longrunvar(matrix indatamat, bool demeaned[0], lagtrunc[4]) /* indatamat is a TxN time series matrix Bartlett window returns an NxN matrix (spectral density at frequency zero) */ N = cols(indatamat) matrix result = zeros(N,N) loop for tau=1..lagtrunc matrix Gamma = autocovar(indatamat,tau,demeaned) # positive and negative tau range together: matrix result += (1 -tau/(lagtrunc+1)) * (Gamma+Gamma') endloop # add the tau=0 part: matrix result += autocovar(indatamat,0,demeaned) return result end function