Fit of midday Nord Pool eletricity prices  year 2000 to CKLS equation     dX (t)  =  
      mu X (t) dt + sigma X (t)^gamma dB (t)

<< Statistics`DescriptiveStatistics`; << 
      Statistics`ContinuousDistributions`; $RecursionLimit = $IterationLimit = Infinity;

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

Clear[CDFStationary]; 
    CDFStationary[mu_, sigma_, gamma_, x_] := 
      Gamma[((1 - 2 gamma)/(2 - 2 gamma))]/
          Gamma[((1 - 2 gamma)/(2 gamma - 2))] -  (Gamma[((3 - 4 gamma)/(2 - 2 gamma))]/(2 gamma - 
                    1)  - (gamma Gamma[((1 - 2 gamma)/(2 - 2 gamma))] + (1 - 
                          2 gamma) Gamma[((1 - 2 gamma)/(2 - 2 gamma)), ((mu (x^(2 - 2 gamma)))/((gamma - 1) (sigma^2)))])/2)/
          Gamma[(1/(2 - 2 gamma))];

Clear[KSDistance]; 
    KSDistance[data_, mu_, sigma_, 
        gamma_] := {sdata = Sort[data], lt = Length[data], 
          Max[Table[
              Max[{Abs[
                    CDFStationary[mu, sigma, gamma, 
                        sdata[[i]]] - (i - 1)/lt], 
                  Abs[CDFStationary[mu, sigma, gamma, sdata[[i]]] - 
                      i/lt]}], {i, 1, lt}]]}[[3]]

\(NP = 
      Part[Transpose[
          Import["/users/math/palbin/user/courses/StokAnal/NP.csv", "CSV"]], 5]; lt = Length[NP]; 
    ListPlot[NP, PlotJoined -> True, 
      PlotLabel -> "NP prices 2000-2004"]; 
    NP2000 = Table[NP[\([i]\)], {i, 1, 350}]; 
    ListPlot[NP2000, PlotJoined -> True, 
      PlotLabel -> "NP prices 2000"];\)

NMinimize[{KSDistance[NP2000, mu, sigma/10000, gamma], mu > 0, 
        sigma > 0, gamma > 1}, {mu, sigma, gamma}]

{0.0904088, {gamma -> 2.90370, mu -> 0.37968, 
        sigma -> 1.23751}}

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

Simulation of CKLS model is a research level problem,  
    and cannot be discussed here in detail.

It turns out that it is suitable to simulate the transformed process  Y (t) = X (t)^(1 - gamma0),  
    rather that  X (t) itself,  
where,  
    by the Ito  formula,  
    dY (t)  =  
      phi (Y (t)) dt  +  (1 - gamma0) sigma0 dB (t),  
    with  phi (t)  as given below.

We use an implicit Euler scheme to simulate  Y (t),  
    and then get  X (t)  as  Y (t)^(1/(1 - 
                  gamma0)).

Clear[phi]; 
    phi[y_] := (1 - gamma0)*mu0*y + 
        gamma0*(gamma0 - 1)*sigma0^2/(2*y); step = 200; days = 350; 
    xi = Random[UniformDistribution[0, 1]]; 
    Y = {(y /. 
              FindRoot[
                CDFStationary[mu0, sigma0, gamma0, y] == xi, {y, 
                  80}])^(1 - gamma0)}; 
    For[i = 1, i ≤ step*days, i++, 
      AppendTo[Y, 
        y /. FindRoot[
            Y[[i]] + 
                phi[y]/step + (1 - gamma0)*sigma0*
                  Random[NormalDistribution[0, Sqrt[1/step]]] == y, {y, 
              SimNP[[i]]}]]]; 
    SimNP = (Table[
            Y[[1 + (i - 1)*step]], {i, 1, 
              days}])^(1/(1 - gamma0));

ListPlot[SimNP, PlotJoined -> True, 
        PlotLabel -> "Simulated CKLS Model Fitted to NP Prices"];