Maximum Likelihood Estimation

Questions
, ,
1

How come there is no discussion about MLE in the community?
Is it so straightforward to implement it in Mathematica?

I have been trying for a while to replicate the results from Barlow (2002) - “A diffusion model for electricity prices”, but I didn’t succeed.

Does anyone have any similar code, or any suggestions?

Letting DD be the vector of prices, here is what I have done:

code.jpg

And then, assuming some initial values, I use:

FindMaximum[{obslogL && \[Sigma] > 0 && \[Lambda] != 0 && \[Alpha] != 0}, {{\[Lambda], 10}, {a, 4}, {\[Sigma], 2}, {\[Alpha], 1.3}}]

However, I get the following error:

“FindMaximum::nrnum: The function value -False is not a real number at {\,a,\,\} = {9.0541,1.43219,-0.40076,0.265587}. >>”

and the hopeless result:
{-737.874, {\ → 9.99596, a → 4.00037, \ → 2.67958, \ → -4.17363}}

2

Sandu,

I’m not sure I completely understand what you are trying to do, but would EstimatedDistribution with ParameterEstimator → “MaximumLikelihood” be of any help?
see http://reference.wolfram.com/mathematica/ref/EstimatedDistribution.html section Details and Options

Peter

3

EstimateDistribution requires a specified distribution.

Mine, has the following PDF:

p(Si,Si-1) = |g'[Si| * 1/Sqrt[2 Pi Theta] * E^(-(g[Si] - rho*g[Si-1] - b)/(2 Theta)))^2

Where, theta, rho and b are known fucntions of my parameters.

It looks like a normal one, but it isn’t.
If could define a distribution from this PDF, then I guess I’m done.

4

Dear Sandu, we recommend avoiding using Subscripts in community post - they increase code volume and make it unreadable. We suggest to use, for example, instead of this

Subscript[g, ?]

simply

ga
  • plain English names for variables.
5

Sandu,

ProbabilityDistribution allows you to define your own Distribution from any given pdf or cdf, symbolic or numeric, while EmpiricalDistribution (and cousins) creates distributions from given data.

http://reference.wolfram.com/mathematica/ref/ProbabilityDistribution.html
http://reference.wolfram.com/mathematica/ref/EmpiricalDistribution.html

Hope this helps.

6

Hi,
a simple made up example to illustrate use of ProbabilityDistribution:

In[49]:= fun[
  x_] := (1 - c)*(1/(E^((-a + x - c*x)^2/(2*b^2))*(b*Sqrt[2*Pi])))

In[50]:= Integrate[fun[x], {x, -Infinity, Infinity}, 
 Assumptions -> b > 0 && 0 < c < 1]
Out[50]= 1

In[51]:= dist = 
  ProbabilityDistribution[fun[x], {x, -Infinity, Infinity}];

In[52]:= data = 
  RandomVariate[dist /. {a -> 0, b -> 1, c -> 1/3}, 10^3];

In[53]:= FindDistributionParameters[data, dist, {{a, .2}, {b, .8}, {c, .3}}, 
 ParameterEstimator -> "MaximumLikelihood"]
Out[53]= {a -> 0.0312187, b -> 0.900318, c -> 0.368235}

As for FindMaximum I would try to avoid != and use inequalities instead, maybe even split the region.

7

Thank you very much, Peter!

It does help.. somehow, but I still have problems with it.

7388D.jpg

Namely, I get different results if I start with a different “theta” for example.

8

There may not be “one theta to rule them all” :slight_smile:

Compare LogLikelihood values of the results:

In[1]:= fun[
  x_] := (1 - c)*(1/(E^((-a + x - c*x)^2/(2*b^2))*(b*Sqrt[2*Pi])))

In[2]:= dist = 
  ProbabilityDistribution[fun[x], {x, -Infinity, Infinity}];

In[3]:= data = 
  RandomVariate[dist /. {a -> 0, b -> 1, c -> 1/3}, 10^3, 
   WorkingPrecision -> 20];

In[4]:= res = 
 Table[Quiet@
   FindDistributionParameters[data, 
    dist, {{a, .2}, {b, .8}, {c, start}}, 
    ParameterEstimator -> "MaximumLikelihood", 
    WorkingPrecision -> 20], {start, .1, .9, .1}]

Out[4]= {{a -> -0.019932674343740414493, b -> 1.1096360920952880069, 
  c -> 0.26954206522212784697}, {a -> -0.018766451751580295519, 
  b -> 1.0447134039438880896, 
  c -> 0.31227976365858204554}, {a -> -0.016674487536403968415, 
  b -> 0.92825542430454599255, 
  c -> 0.38894242441693089468}, {a -> -0.015633567435080450467, 
  b -> 0.87030823228157333632, 
  c -> 0.42708825124255829112}, {a -> -0.014913446774718998486, 
  b -> 0.83021968892034345328, 
  c -> 0.45347797919258945110}, {a -> -0.014144317149532428768, 
  b -> 0.78740285539643202694, 
  c -> 0.48166370243199674812}, {a -> -0.013811227971511007229, 
  b -> 0.76886004700210827827, 
  c -> 0.49387017410287984506}, {a -> -0.015406220411774267890, 
  b -> 0.85765200449162658321, 
  c -> 0.43541966915346278263}, {a -> -0.015791913582145423886, 
  b -> 0.87912323840309263789, c -> 0.42128545588647653513}}

In[5]:= LogLikelihood[dist /. #, data] & /@ res

Out[5]= {-1837.0542834440804521, -1837.0542834440804521, \
-1837.0542834440804521, -1837.0542834440804521, \
-1837.0542834440804521, -1837.0542834440804521, \
-1837.0542834440804521, -1837.0542834440804521, -1837.0542834440804521}

(*compare to the "original" distribution*)

In[6]:= LogLikelihood[dist /. {a -> 0, b -> 1, c -> 1/3}, data]

Out[6]= -1837.3811539068860566

(*the original distribution may not be the best fit to the data, even \
though it was used to generate said data*)
9

I have a somewhat related question, which could also provide an answer to another unanswered question of mine: TimeSeriesModelFit // (G)ARCH estimating methods

Hopefully this thread will be revised by someone after so long time.

  • How do I maximize the MLE which contains many variables in a reasonable time?

Below I have uploaded a notebook where I try to find 103 parameters.

Excel’s Solver can do the maximization with less parameters relatively easy (I didn’t try with so many), but of course there is a precision concern. MaxBFGS from OxMetrics, on the other hand, does it reasonably fast.

FindMaximum in Mathematica, however, takes years…

10

Hi mam what about simulation of likelihood eatimates and their variances?