# Does simulation time step have an impact on results (Vasicek model for instance)

I am trying to make some computations using Vasicek short rate model. Especially I am trying to compare exact expectation(obtained with the formula) and the expectation from Monte Carlo simulation.

### Exact computation

I use:

$E[r_t] = r_0 * \exp(-a*t) + (\theta/a)*(1-\exp(-a*t))$

public override double GetExpectation(double r0, double t)
{
double expectation = r0 * Math.Exp(-_a * t) + (_theta / _a) * (1 - Math.Exp(-_a * t));
return expectation;
}


### Monte Carlo simulations

I use the following method:

• I compute r from a time t to a time t+dt using:

public override double ComputeNextValue(double r0, double dt)
{
RandomVariableGenerator rvg = RandomVariableGenerator.GetInstance();

double randomGaussian = rvg.GetNextRandomGaussian();
double r_t_dt = (_theta - _a*r0)*dt + _sigma * Math.Sqrt(dt) * randomGaussian;

return r_t_dt;
}

• then a compute a path from 0 to t with dt as time step using:

public override double ComputeValue(double r0, double t, double dt)
{
double x = r0;

for(double slot = dt; slot <= t; slot += dt)
{
x = ComputeNextValue(x, dt);
}
return x;
}

• Then I compute the Monte Carlo Expectation using:

public override double ComputeMonteCarloExpectation(double r0, double t, double dt, int nreps)
{
double sum = 0.0;
double value;
for (int i = 0; i < nreps; i++)
{
value = ComputeValue(r0, t, dt);
sum += value;
}
return sum / nreps;
}


I use the following parameters:

double sigma = 0.03;
double r0 = 0.03;
double theta = 0.1;
double a = 0.3;

int nreps = 1000;
double t = 1;


For dt = 0.1:

Exact expectation: 0,108618473059879;
Monte Carlo Expectation: 0,0101464832161612


For dt = 1:

Exact expectation: 0,108618473059879;
Monte Carlo Expectation: 0,092058844704742


using dt = 1 leads to a result close to exact value while using dt = 0.1 seems to lead to a result having a 0.1 factor difference with exact one.

I think I am doing something wrong but I can't figure it out. Do you have an idea?

It would probably help if you sum up the Euler steps, either in

public override double ComputeNextValue(double r0, double dt)
{
...
return r_0+r_t_dt;
}


or in the integration loop

public override double ComputeValue(double r0, double t, double dt)
{
double x = r0;

for(double slot = dt; slot <= t; slot += dt)
{
x += ComputeNextValue(x, dt);
}
return x;
}


The SDE $dr_t=(θ-ar_t)dt+\sigma dW_t$ gets approximated in Euler fashion as $$\Delta r_t=(θ-ar_t)Δt+\sigma ΔW_t+O(Δt^{3/2})$$ where $\Delta r_t=r_{t+Δt}-r_t$ so that to compute the next value you need to add the increment to the current value.

One can also directly give an explicit formula for the paths, as $$r_t=e^{-at}(r_0+σ\widetilde W_{(e^{2at}-1)/(2a)})+\fracθa(1-e^{at})$$ where $\widetilde W$ is a different Wiener process.

• Thanks LutzL. Indeed the generator implements the singleton pattern so it always returns the same object. Concerning the sum, excuse me but I don't understand why it is necessary to sum. For me the program computes a path that will lead to r_t, going through many 0+dt. I mean the final value for me is r_t but you seem to say that r_t is the sum of rt_dt? Commented Dec 7, 2017 at 11:19
• Thanks I got it Commented Dec 7, 2017 at 13:19