code: def block(
      rg: Generator,
      trend: str,
      sample_sizes: list[int],
      cross_section_size: int,
      simulations: int,
      idx: int,
  ) -> None:
      filename = f"engle-granger-results-trend-{trend}-{idx:04d}.npz"
      fullfile = os.path.join(ROOT, filename)
      if os.path.exists(fullfile):
          return
      max_k = cross_section_size
      max_t = max(sample_sizes) + 1
      max_blocks_size = max(MAX_MEMORY // (max_t * max_k * 8), 1)
      print(f"Max block size: {max_blocks_size}")
      remaining = simulations
      completed = 0
      tstats = np.empty((len(sample_sizes), simulations, max_k))
      start = dt.datetime.now()
      last_report = remaining + 1001
      while remaining:
          this_block = min(remaining, max_blocks_size)
          remaining -= this_block
          e = rg.standard_normal((this_block, max_t + 50, max_k))
          data = np.cumsum(e, 1)
          data = data[:, 50:, :]
          del e
          for j, sub_size in enumerate(sample_sizes):
              # Copy for needed for 'nc' trend
              sub = data[:, : sub_size + 1].copy()
              if trend not in ("c", "ct", "ctt", "nc"):
                  raise NotImplementedError
              if trend in ("c", "ct", "ctt"):
                  # Mean 0
                  mu = sub.mean(1)
                  mu.shape = (this_block, 1, cross_section_size)
                  sub -= mu
              if trend in ("ct", "ctt"):
                  # Orthogonalize to trend
                  tau = np.arange(float(sub_size + 1))
                  tau -= tau.mean()
                  tau.shape = (sub_size + 1, 1)
                  coefs = (tau * sub).sum(1) / (tau**2).sum()
                  coefs.shape = (this_block, 1, cross_section_size)
                  sub -= coefs * tau
              if trend == "ctt":
                  tau = np.arange(float(sub_size + 1))
                  tau -= tau.mean()
                  tau2 = np.arange(float(sub_size + 1)) ** 2
                  tau2 -= tau2.mean()
                  tau2 -= (tau * tau2).sum() / (tau**2).sum() * tau
                  tau2.shape = (sub_size + 1, 1)
                  coefs = (tau2 * sub).sum(1) / (tau2**2).sum()
                  coefs.shape = (this_block, 1, cross_section_size)
                  sub -= coefs * tau2
              errors = np.empty((this_block, sub_size + 1, max_k))
              errors[:, :, :1] = sub[:, :, :1]
              rhs = sub[:, :, 1:]
              ip = (
                  np.matmul(np.transpose(rhs, (0, 2, 1)), rhs)
                  / sub_size
              )
              ip_chol_inv = np.linalg.inv(np.linalg.cholesky(ip))
              orth = np.matmul(
                  rhs, np.transpose(ip_chol_inv, (0, 2, 1))
              )
              del ip_chol_inv, ip, rhs, sub
              for i in range(1, max_k):
                  y = errors[:, :, (i - 1) : i]
                  x = orth[:, :, (i - 1) : i]
                  xp = np.transpose(x, (0, 2, 1))
                  xpy = np.matmul(xp, y)
                  xpx = np.matmul(xp, x)
                  coefs = xpy / xpx
                  errors[:, :, i : (i + 1)] = y - coefs * x
              del orth, x, y
              # Compute ADF statistics for each of the error series
              err_lag = errors[:, :-1]
              derrors = errors[:, 1:] - err_lag
  
              ypy = (derrors * derrors).sum(1)
              xpy = (derrors * err_lag).sum(1)
              xpx = (err_lag * err_lag).sum(1)
              del errors, err_lag, derrors
  
              coefs = xpy / xpx
              sse = ypy - coefs**2 * xpx
              sigma2 = sse / sub_size
              se = np.sqrt(sigma2 / xpx)
              tstats[j, completed : (completed + this_block)] = (
                  coefs / se
              )
          completed += this_block
          elapsed = dt.datetime.now() - start
          time_per_iter = elapsed.total_seconds() / completed
          remaining_time = int(time_per_iter * remaining)
          if last_report - remaining > 1000:
              last_report = remaining
              print(
                  f"Index: {idx}, Trend: {trend}, Remaining: {remaining}"
              )
              print(
                  f"Est. time remaining: {str(dt.timedelta(seconds=remaining_time))}"
              )
  
      out = np.percentile(tstats, PERCENTILES, 1)
      np.savez(
          fullfile,
          quantiles=out,
          trend=np.array([trend]),
          sample_sizes=np.array(sample_sizes),
      )
  
file_path: arch/arch/unitroot/critical_values/simulation/engle_granger_simulation.py
imports: import argparse
  import datetime
  from functools import partial
  from itertools import product
  import os
  from random import shuffle
  import sys
  import colorama
  from joblib import Parallel, cpu_count, delayed
  import numpy
  from numpy.random import PCG64, Generator, SeedSequence
label_desc: Use image node_id 1 for calling a global function with example usage: block(rg, trend, sample_sizes, cross_section_size, simulations, idx) without return types