code: def multidimensional_discretization(rho, sigma, N=3, method='rouwenhorst', m=2):
      
  
      # rho is assumed to be a scalar
      # sigma is a positive symmetric matrix
      # N number of points in each non-degenerate dimension
      # m : standard deviations to approximate
  
      import scipy.linalg
      from itertools import product
  
      d = sigma.shape[1]
  
      sigma = sigma.copy()
  
      zero_columns = np.where(sigma.sum(axis=0)==0)[0]
      for i in zero_columns:
          sigma[i,i] = 1
  
      L = scipy.linalg.cholesky(sigma)
  
      N = int(N)
  
      if method=='tauchen':
          [nodes_1d, probas_1d] = tauchen(N, 0, rho, 1, m=m)
      elif method=='rouwenhorst':
          [nodes_1d, probas_1d] = rouwenhorst(rho, 1, N)
  
      markov_nodes = np.array( list( product( *([nodes_1d]*d))) ).T
      markov_indices = np.array( list( product( *([range(N)]*d)  )  ) ).T
  
      markov_nodes = np.dot(L, markov_nodes)
  
  
      transition_matrix = 1
      for i in range(d):
          transition_matrix = np.kron(transition_matrix, probas_1d)
  
      markov_nodes = np.ascontiguousarray(markov_nodes.T)
  
      for i in zero_columns:
          markov_nodes[:,i] = 0
  
      return [markov_nodes, transition_matrix]
docs: Discretize an VAR(1) into a markov chain. The autoregression matrix is supposed to be a scalar.
  
      :param rho:
      :param sigma:
      :param N:
      :param method:
      :param m:
      :return: