code: def lag_matrix(blk, max_lag=None):
    
    if max_lag is None:
      max_lag = len(blk) - 1
    elif max_lag >= len(blk):
      raise ValueError("Block length should be higher than order")
  
    return [[sum(blk[n - i] * blk[n - j] for n in xrange(max_lag, len(blk))
                ) for i in xrange(max_lag + 1)
            ] for j in xrange(max_lag + 1)]
docs: Finds the lag matrix for a given 1-D block sequence.
  
    Parameters
    ----------
    blk :
      An iterable with well-defined length. Don't use this function with Stream
      objects!
    max_lag :
      The size of the result, the lags you'd need. Defaults to ``len(blk) - 1``,
      the maximum lag that doesn't create fully zeroed matrices.
  
    Returns
    -------
    The covariance matrix as a list of lists. Each cell (i, j) contains the sum
    of ``blk[n - i] * blk[n - j]`` elements for all n that allows such without
    padding the given block.