# Fitting a theoretical density function with log-return sample

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.

-
Welcome to Mathematica.SE, Geoffrey!. FindDist.Param does not look like Mathematica syntax. Are you sure you are not getting your software mixed up? Are you using RLink? Please add some more information about what you have tried for estimating this model, i.e. all your estimation code, not just the specific functions used. – Verbeia Mar 15 '13 at 5:39
@Verbeia I think the OP is a bit lazy and is arbitrarily shortening function names. He is referring to FindDistributionParameters, ExponentialDistribution, LaplaceDistribution. @Geoffrey: You have to be less sloppy with your questions if you like to get first class answers (or any answers at all). – Sjoerd C. de Vries Mar 15 '13 at 7:17
Please post the exact code you used and describe your problems you encounter more accurately than "without any success". If you don't improve your question I will have to vote to close. – Sjoerd C. de Vries Mar 15 '13 at 11:51

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.

-