Message Boards Message Boards

0
|
23405 Views
|
9 Replies
|
6 Total Likes
View groups...
Share
Share this post:

Maximum Likelihood Estimation

Posted 11 years ago
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:



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}}
POSTED BY: Sandu Ursu
9 Replies
Posted 10 years ago

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...

Attachments:
POSTED BY: Sandu Ursu
There may not be "one theta to rule them all"  emoticon

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*)
POSTED BY: Gosia Konwerska
Posted 10 years ago

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

POSTED BY: Maha haroon
Posted 11 years ago
Thank you very much, Peter!

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



Namely, I get different results if I start with a different "theta" for example.
POSTED BY: Sandu Ursu
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.
POSTED BY: Gosia Konwerska
Posted 11 years ago
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.
POSTED BY: Peter Fleck
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.
POSTED BY: Moderation Team
Posted 11 years ago
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.
POSTED BY: Sandu Ursu
Posted 11 years ago
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
POSTED BY: Peter Fleck
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract