# Vasicek Short rate simulation - analytical formula vs discretization

I've been using two approaches to simulate Vasicek short rate paths and I'm wondering if one of them is more correct than the other.

The first approach is based on the analytical formula (see code below, "OU_processes"), whereby every next simulated rate along a path is calculated from the previously simulated rate one timestep ago as follows:

def OU_processes(years, timestep, num_sims, startRate, kappa, theta, sigma):
"""
timestep has to be defined in years or a fraction of years
e.g. 0.1 => 1/10th of a year; 2 => 2 years
"""
times = np.arange(0,years+timestep,timestep)
epsilon = np.random.normal(0, 1, (num_sims, len(times)-1))
elt = 0.5 / kappa * (1.0 - np.exp (-2.0 * kappa * timestep))
V = elt * sigma ** 2
sqrt_V = np.sqrt(V)
ou = np.zeros((num_sims, len(times)))
ou[:,0] = startRate
ou[:, 1:] = np.kron(sqrt_V, np.ones((num_sims, 1))) * epsilon
for i in range(1, ou.shape[1]):
ou[:, i] += theta * (1 - np.exp(-kappa * timestep))
ou[:, i] += np.exp (-kappa * timestep) * ou[:, i-1]
ou = pd.DataFrame(np.transpose(ou))
ou.index = times
return ou


The second approach uses Euler discretization (see code below, "discretized_OU"):

def discretized_OU(years, timestep, num_sims, startRate, kappa, theta, sigma):
times = np.arange(0,years+timestep,timestep)
epsilon = np.random.normal(0, 1, (num_sims, len(times)-1))
ou=np.zeros((num_sims,len(times)))
ou[:,0] = startRate
for step in np.arange(1,int(years/timestep)+1,1):
ou[:,step]=ou[:,step-1]+kappa*(theta-ou[:,step-1])*timestep+sigma*np.sqrt(timestep)*epsilon[:,step-1]
ou = pd.DataFrame(np.transpose(ou))
ou.index = times
return ou


Could anyone tell me what the difference is between both methods? The paths simulated using the analytical method lead to a ZCB price that is very close to the price found using the analytical Vasicek ZCB formula, while the paths simulated using the discretized approach lead to a price that is somewhat further from the analytical price (although still very close to it) and with a seemingly larger standard deviation:

Moreover, is it correct that only the second approach is available for simulating CIR paths?