Fitting a theoretical density function with log-return sample

Question

I am working on a exercise to fit a theoretical density function,presented below, with log-return series of a stock.

$\frac{1}{2}a e^{a(x-\mu)} \quad \text{if x } < \mu $

or

$ \frac{1}{2}be^{-b(x-\mu)} \quad \text{if x } \geq \mu \\$

I am now stuck on estimating the parameters for the function. I have tried different approaches, such as creating a probability distribution and use FindDist.Param. and also tried to fit it with a ExponentialDistr(or Laplace), multiplying neg. returns with -1 for case $(x < \mu)$, however without any success. Given I have the moments of the sample and at least think I have worked out the correct moments (mean, var, skew, kurt) I would also like to try a method of moments, but not sure how to do proceed with the above function? As I am relatively new to Mathematica and keep getting lost but enjoying the endless possibilities, any suggestions or recommendations on procedure is greatly appreciated.

For importing my log-returns I have used ImportString["...", "List"].

---

Many thanks for the feedback guys! Alas, my original post was written 4 am in the morning. Anyway, I'll be more thorough in coming posts. The following is the procedure I have tried,

1. logretun = ImportString["r1, r2, .., rT", "List"] Where $r_i$ are daily log-returns of a stock, copied from an .csv spreadsheet.

2. Define the above density as a function,

   f[x_] = 1/2*a*Exp[a*(x - m)]*Boole[x < m] + 
     1/2*b*Exp[-b*(x - m)]*Boole[x >= m]
   

   (is the use of Boole correct?)

3. Define a function from (2) as distribution

   fdist = ProbabilityDistribution[f, {x, -Infinity, Infinity}]
   

4. Run

   FindDistributionParameter[logreturn, fdist] 
   

   With no luck. (Just runs and runs with no answer - even if I take a smaller sample, so I do not think its me just being impatient..).

I have gotten some values for $a$ and $b$ (87.76922562 and -82.00007391, respectively) through the following method; however I am not sure it is justifiable.

Step (3) piecewise, for a ExponentialDistribution[λ]:

1. Sort the return series, select the ones that are $x < \mu$ and multiply by -1 and import as "minuslogreturns"

2. Run step 3 as

   FindDistributionParameter[minuslogreturns, ExponentialDistribution[a]]
   

3. Run again for the positive part, with $b$ as the parameter for ExponentialDistribution.

I am not confident that this is a feasible approach. I have tested my $a$ and $b$ by moments ($\mu$ = sample median), getting relative close mean, variance and skewness, however I cannot seem to get the kurtosis close. (sample kurtosis 9.494758128 vs estimated 15.33348985). On a side note, it might be very well that the expressions for expected return, variance, skewness and kurtosis also be off, as they are quite long expressions. I can post these if deemed necessary?

Hope this make my questions a bit more clearer.

Answer 1

If m is a known value, then the method of moments only requires the sample mean and sample variance to be equated to the mean and variance of the distribution (and one can have weights 1/2 by using w as the weight for the values of x < m:

(* Weight *)
w = 1/2

(* Mean *)
mean = FullSimplify[
  Integrate[w a x Exp[a (x - m)], {x, -Infinity, m}, 
    Assumptions -> a > 0] + 
   Integrate[(1 - w) b x Exp[-b (x - m)], {x, m, Infinity}, 
    Assumptions -> b > 0]]

(* Variance *)
variance = 
 FullSimplify[
  Integrate[w a x^2 Exp[a (x - m)], {x, -Infinity, m}, 
    Assumptions -> a > 0] + 
   Integrate[(1 - w) b x^2 Exp[-b (x - m)], {x, m, Infinity}, 
    Assumptions -> b > 0] - mean^2]

(* Solve for a and b *)
FullSimplify[Solve[{xbar == mean, sd^2 == variance}, {a, b}]]

with the following results

Out[228]= -(1/(2 a)) + 1/(2 b) + m

Out[229]= (3 a^2 + 2 a b + 3 b^2)/(4 a^2 b^2)

Out[230]= ((a - b) (7 a^2 + 10 a b + 7 b^2))/(4 a^3 b^3)

Out[231]= (3 (39 a^4 + 20 a^3 b + 10 a^2 b^2 + 20 a b^3 + 
   39 b^4))/(16 a^4 b^4)

In[233]:= FullSimplify[Solve[{xbar == mean, sd^2 == variance}, {a, b}]]

Out[233]= {{a -> 1/(
   m + Sqrt[-(m - sd - xbar) (m + sd - xbar)]/Sqrt[2] - xbar), 
  b -> 1/(-m + Sqrt[-(m - sd - xbar) (m + sd - xbar)]/Sqrt[2] + 
    xbar)}, {a -> 1/(
   m - Sqrt[-(m - sd - xbar) (m + sd - xbar)]/Sqrt[2] - xbar), 
  b -> 1/(-m - Sqrt[-(m - sd - xbar) (m + sd - xbar)]/Sqrt[2] + 
    xbar)}}

If m (or w) is not known, then you'll need to calculate the skewness and/or kurtosis and probably only get a numeric rather than symbolic solution.