<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def random_data(N=5000, K=3, unobservables=False, **kwargs):
 	""" 	Function that generates data according to one of two simple models that 	satisfies the unconfoundedness assumption. 	The covariates and error terms are generated according to 		X ~ N(mu, Sigma), epsilon ~ N(0, Gamma). 	The counterfactual outcomes are generated by 		Y0 = X*beta + epsilon_0, 		Y1 = delta + X*(beta+theta) + epsilon_1. 	Selection is done according to the following propensity score function: 		P(D=1|X) = Lambda(X*beta). 	Here Lambda is the standard logistic CDF. 	Parameters 	N: int 		Number of units to draw. Defaults to 5000. 	K: int 		Number of covariates. Defaults to 3. 	unobservables: bool 		Returns potential outcomes and true propensity score 		in addition to observed outcome and covariates if True. 		Defaults to False. 	mu, Sigma, Gamma, beta, delta, theta: NumPy ndarrays, optional 		Parameter values appearing in data generating process. 	Returns 	------- 	tuple 		A tuple in the form of (Y, D, X) or (Y, D, X, Y0, Y1) of 		observed outcomes, treatment indicators, covariate matrix, 		and potential outomces. 	"""


	mu = kwargs.get('mu', np.zeros(K))
	beta = kwargs.get('beta', np.ones(K))
	theta = kwargs.get('theta', np.ones(K))
	delta = kwargs.get('delta', 3)
	Sigma = kwargs.get('Sigma', np.identity(K))
	Gamma = kwargs.get('Gamma', np.identity(2))

	X = np.random.multivariate_normal(mean=mu, cov=Sigma, size=N)
	Xbeta = X.dot(beta)
	pscore = logistic.cdf(Xbeta)
	D = np.array([np.random.binomial(1, p, size=1) for p in pscore]).flatten()

	epsilon = np.random.multivariate_normal(mean=np.zeros(2), cov=Gamma, size=N)
	Y0 = Xbeta + epsilon[:,0]
	Y1 = delta + X.dot(beta+theta) + epsilon[:,1]
	Y = (1-D)*Y0 + D*Y1

	if unobservables:
		return Y, D, X, Y0, Y1, pscore
	else:
		return Y, D, X