text: import statsmodels.api as sm
  import numpy as np
  
  ##################
  #Monte Carlo test#
  ##################
  modrand1 = np.random.RandomState(5676576)
  modrand2 = np.random.RandomState(1543543)
  modrand3 = np.random.RandomState(5738276)
  X = modrand1.uniform(0, 5, (1000, 4))
  X = sm.add_constant(X)
  beta = np.array([[10], [2], [3], [4], [5]])
  y = np.dot(X, beta)
  params = []
  for i in range(10000):
      yhat = y + modrand2.standard_normal((1000, 1))
      cens_times = 50 + (modrand3.standard_normal((1000, 1)) * 5)
      yhat_observed = np.minimum(yhat, cens_times)
      censors = np.int_(yhat < cens_times)
      model = sm.emplike.emplikeAFT(yhat_observed, X, censors)
      new_params = model.fit().params
      params.append(new_params)
  
  mc_est = np.mean(params, axis=0)  # Gives MC parameter estimate
  
  ##################
  #Koul replication#
  ##################
  
  koul_data = np.genfromtxt('/home/justin/rverify.csv', delimiter=';')
  # ^ Change path to where file is located.
  koul_y = np.log10(koul_data[:, 0])
  koul_x = sm.add_constant(koul_data[:, 2])
  koul_censors = koul_data[:, 1]
  koul_params = sm.emplike.emplikeAFT(koul_y, koul_x, koul_censors).fit().params
  