Fit of OU process  dX (t)  =  -mu X (t) dt + 
        sigma dB (t)  to SP500 volatility process

<< Statistics`DescriptiveStatistics`; << 
      Statistics`ContinuousDistributions`;

SP500 = 
      Part[Transpose[
          Reverse[Import["/users/math/palbin/user/courses/StokAnal/SP500.\
csv", "CSV"]]], 7]; lt = Length[SP500]; 
    LogIncrSP500 = 
      Table[Log[SP500[[i + 1]]] - Log[SP500[[i]]], {i, 1, lt - 1}]; 
    lt = Length[LogIncrSP500]; 
    Volatility = 
      Table[StandardDeviation[
          Table[LogIncrSP500[[i]], {i, j - 19, j}]], {j, 20, lt}]; 
    lt = Length[Volatility]; 
    ListPlot[Volatility, PlotJoined -> True, 
      PlotLabel -> "Volatility of SP500 log-increments"];

Clear[FEmp]; 
    FEmp[data_, 
        x_] := {lt = Length[data], 
          Sum[If[data[[i]] ≤ x, 1, 0], {i, 1, lt}]/
            lt}[[2]];

Clear[pOU]; 
    pOU[mu_, sigma_, x_, y_, t_] := 
      Exp[-(y - x*Exp[-mu*t])^2/(2*
                sigma^2*
(1 - Exp[-2*mu*t])/(2*
                      mu))]/(Sqrt[2*Pi]*sigma*
            Sqrt[1 - Exp[-2*mu*t]]/
              Sqrt[2*mu]);
Clear[CDFOU]; 
    CDFOU[mu_, sigma_, x_, z_, t_] := 
      (1/2) (1 + 
            Erf[(((Sqrt[(mu/(-1 + (E^(2 mu t))))]) (-x + (E^(mu t)) z))/sigma)]);
Clear[CDFStationary]; 
    CDFStationary[mu_, sigma_, x_] := Erf[(((Sqrt[mu]) x)/sigma)];

NMaximize[{Sum[
          Log[pOU[mu, sigma, Volatility[[i]], Volatility[[i + 1]], 
              1]], {i, 1, lt - 1}], mu > 0, sigma > 0}, {mu, 
        sigma}]

{mu0, sigma0} = {mu, sigma} /. %[[2]]; 
    Show[{Plot[
          CDFStationary[mu0, sigma0, x], {x, Min[Volatility], 
            Max[Volatility]}, PlotRange -> {0, 1}], 
        Plot[FEmp[Volatility, x], {x, Min[Volatility], Max[Volatility]}, 
          PlotRange -> {0, 1}]}, 
      PlotLabel -> "Empirical CDF vs fitted CDF"];

Simulation of OU model is straightforward by Euler method

SimVol = {Mean[Volatility]}; 
    For[i = 1, i <= lt - 1, i++, 
      AppendTo[SimVol, 
y /. {xi = 
                  Random[UniformDistribution[0, 1]], 
                FindRoot[
                  CDFOU[mu0, sigma0, 
SimVol[[i]], y, 
                      1] == xi, {y, 
                    SimVol[[1]]}]}[[2]][[1]]]];

ListPlot[SimVol, PlotJoined -> True, 
        PlotLabel -> "Simulated fitted OU-model"];