Statistical inference for CIR process (cf .  
              Exercise 6.17 and Section 12.1)  dX (t) = 
      b*(a -  X (t)) dt + 
        X (t)^(1/2) dB (t)

        ----------------------------------------------------------------------------------------------------------------

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

The transition density function for the CIR process   pCIR[a, 
          b, x, y, t]  =  
      P {X (t) ∈ dy | X (0) = 
            x} dy   can be
found at various places in the literature :

Clear[c1]; c1[b_, t_] := 2*b/(1 - Exp[-b*t]); Clear[u1]; 
    u1[b_, x_, t_] := c1[b, t]*x*Exp[-b*t]; Clear[v1]; 
    v1[b_, y_, t_] := c1[b, t]*y; Clear[pCIR]; 
    pCIR[a_, b_, x_, y_, t_] := 
      c1[b, t]*Exp[-u1[b, x, t] - 
            v1[b, y, t]]*(v1[b, y, t]/u1[b, x, t])^(a*b - 1/2)*
        BesselI[2*a*b - 1, 2*Sqrt[v1[b, y, t]*u1[b, x, t]]];

The stationary distribution of the CIR process is given by (6.64) :

Exp[Integrate[2*b*(a - y)/y, {y, 1, x}, Assumptions -> x > 1]]/
      x

Integrate[%, {x, 0, Infinity}, 
      Assumptions -> b > 0 && a > 0]

Clear[PDFStat]; 
    PDFStat[a_, b_, 
        x_] := ((E^(2 (b - b x + a b Log[x])))/x)/((4^(-a b)) (b^(-2 a b)) (E^(2 b)) Gamma[2 a b]);

Check that    PDFStat   satisfies the Forward Equation (6.63)    (D_x)^2 (x*p)/2 + 
        D_x (b*(a - x)*p)  = 0 :

Simplify[
      D[D[x*PDFStat[a, b, x], x], x]/2 - 
        D[b*(a - x)*PDFStat[a, b, x], x]]

Simulation of an OU process is straightforward by the Euler method :

a0 = 5; b0 = 5; StepsPerTimeUnit = 100; T = 100; 
    CIRData = {x /. 
          FindRoot[
            PDFStat[a0, b0, x] == Random[UniformDistribution[0, 1]], {x, 
              1}]}; For[i = 2, i ≤ StepsPerTimeUnit*T, i++, 
      AppendTo[CIRData, 
        CIRData[[i - 1]] + 
          b0*(a0 - CIRData[[i - 1]])/StepsPerTimeUnit + 
          Sqrt[CIRData[[i - 1]]]*
            Random[NormalDistribution[0, 
                1/Sqrt[StepsPerTimeUnit]]]]];

ListPlot[CIRData, PlotJoined -> True];

Method 1  :  
      Estimation of   mu   by Likelihood Ratio according to Equation (10.42)
-------------------------------------------------------------------------

{a1 , 
        b1} = {a, 
          b} /. Solve[ {b*
                    Sum[CIRData[[i - 
                                1]]^(-1)*(CIRData[[i]] - 
                            CIRData[[i - 1]]), {i, 2, 
                        StepsPerTimeUnit*T}] - 
                  b^2*a*Sum[
                      CIRData[[i]]^(-1)/StepsPerTimeUnit, {i, 1, 
                        StepsPerTimeUnit*T}] + b^2*T == 0, 
              a*Sum[CIRData[[i - 1]]^(-1)*(CIRData[[i]] - 
                            CIRData[[i - 1]]), {i, 2, 
                        StepsPerTimeUnit*T}] - 
                  Sum[CIRData[[i]] - CIRData[[i - 1]], {i, 2, 
                      StepsPerTimeUnit*T}] - 
                  b*a^2*Sum[
                      CIRData[[i]]^(-1)/StepsPerTimeUnit, {i, 1, 
                        StepsPerTimeUnit*T}] + 2*b*a*T - 
                  b*Sum[CIRData[[i]]/StepsPerTimeUnit, {i, 1, 
                        StepsPerTimeUnit*T}] == 0}, {a, 
              b}][[2]]

Method 2  :  
      Estimation of   mu   by minimizing Kolmogorov - 
        Smirnov Distance between empirical CDF and stationary CDF
----------------------------------------------------------------------------------------------------------

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

Simplify[Integrate[PDFStat[a, b, y], {y, 0, x}], 
      Assumptions -> a > 0 && b > 0 && x > 0]

Clear[CDFStat]; 
    CDFStat[a_, b_, x_] := 
      1 - Gamma[2 a b, 2 b x]/Gamma[2 a b];

SampleDistance = 50; NData = StepsPerTimeUnit*T/SampleDistance; 
    NewCIRData = Table[CIRData[[i*SampleDistance]], {i, 1, NData}]; 
    SortNewCIRData = 
      Sort[NewCIRData]; {a2, 
        b2} = {a, 
          b} /. NMinimize[{Max[
                Table[Max[{Abs[
                        CDFStat[a, b, SortNewCIRData[[i]]] - (i - 1)/
                            NData], 
                      Abs[CDFStat[a, b, SortNewCIRData[[i]]] - 
                          i/NData]}], {i, 1, NData}]], a > 0, b > 0}, {a, 
              b}][[2]]

Show[{Plot[
            FEmp[NewCIRData, x], {x, SortNewCIRData[[1]], 
              SortNewCIRData[[NData]]}], 
          Plot[CDFStat[a2, b2, x], {x, SortNewCIRData[[1]], 
              SortNewCIRData[[NData]]}]}];

Method 3  :  
      Estimation of   mu   by Maximum Likelihood method based on    pCIR[a, b, x, y, t]
-------------------------------------------------------------------------------

{a3, 
        b3} = {a, 
          b} /. NMaximize[{Sum[
                Log[pCIR[a, b, NewCIRData[[i - 1]], NewCIRData[[i]], 
                    T/NData]], {i, 2, NData}], a > 0, b > 0}, {a, 
              b}][[2]]