def var_acf(coefs, sig_u, nlags=None):
 (p, k, _) = coefs.shape
 if (nlags is None):
  nlags = p
 result = np.zeros(((nlags + 1), k, k))
 result[:p] = _var_acf(coefs, sig_u)
 for h in range(p, (nlags + 1)):
  for j in range(p):
   result[h] += np.dot(coefs[j], result[((h - j) - 1)])
 return result

def var_acf ( coefs , sig_u , nlags = None ) : ( p , k , _ ) = coefs . shape if ( nlags is None ) : nlags = p result = np . zeros ( ( ( nlags + 1 ) , k , k ) ) result [ : p ] = _var_acf ( coefs , sig_u ) for h in range ( p , ( nlags + 1 ) ) : for j in range ( p ) : result [ h ] += np . dot ( coefs [ j ] , result [ ( ( h - j ) - 1 ) ] ) return result
compute autocovariance function acf_y(h) up to nlags of stable var(p) process parameters coefs : ndarray coefficient matrices a_i sig_u : ndarray covariance of white noise process u_t nlags : int .