Let's assume I have a set of samples of a random variable
$$ X = Y + Z \>, $$
where $Y$ is Gaussian (with a mean of zero and variance $\sigma^2$) and $Z$ has a symmetric $\alpha$-stable distribution with index of stability $\alpha$ and scale parameter $c$. I know $\alpha$, but I don't know $\sigma$ and $c$. How do I go about estimating them?
I am going to assume that you have observations from $Y$ because you are specifying distributional assumptions on $X$ and $Z$. Then $Y=Z-X$.
Given that $Y=Z-X$, we have that the density of $Y$ can be written as the convolution:
$$f_Y(y)=\int_{-\infty}^{\infty}\dfrac{1}{c}f_Z\left(\dfrac{y+x}{c}\right)\dfrac{1}{\sigma}f_X\left(\dfrac{x}{\sigma}\right)dx,$$
where $f_Z$ is the standard $\alpha -$stable density and $f_X$ represents the standard normal density. Note that, using a change of variable, we can rewrite this density as follows
$$f_Y(y)=\int_{-\infty}^{\infty}\dfrac{1}{c}f_Z\left(\dfrac{y+u\sigma}{c}\right)f_X\left(u\right)du = \int_{-\infty}^{\infty}\dfrac{1}{\sigma}f_X\left(\dfrac{cu-y}{\sigma}\right)f_Z\left(u\right)du,$$
It might be difficult (if feasible) to obtain this density in closed form but it can be approximated by simulating $(x_1,...,x_N)$ from a standard normal distribution and calculating the average
$$f_Y(y;c,\sigma)\approx \dfrac{1}{N}\sum_{j=1}^N \dfrac{1}{c}f_Z\left(\dfrac{y+x_j \sigma}{c}\right),$$
or by simulating from a standard $\alpha-$stable distribution $(z_1,...,z_N)$ and calculating the average
$$f_Y(y;c,\sigma)\approx \dfrac{1}{N}\sum_{j=1}^N \dfrac{1}{\sigma}f_X\left(\dfrac{c z_j -y}{\sigma}\right),$$
Note that the first approximation implies an easy simulation but a difficult evaluation of the density $f_Z$. The second approximation implies a difficult simulation but an easy evaluation of $f_X$. I also guess that the number of needed simulations for a good approximation is less in the first approximation. Either way, this might imply the use of intensive computation.
Now, if you have a sample $(y_1,...,y_n)$, then you can write the likelihood of $(c,\sigma)$ as
$${\mathcal L}(c,\sigma)\propto \prod_{j=1}^n f_Y(y_j;c,\sigma).$$
Using this, you can approximate the Maximum Likelihood Estimators (MLE) of $(c,\sigma)$ by maximising the corresponding likelihood using the approximations I described above.
Toy example in R
In this example, consider $(c,\sigma)=(0.5,1)$, $n=100$ and $N=1000$.
rm(list=ls())
library(stabledist)
# Values of the theoretical parameters
alpha0 = 1.75
sigma0 = 1
c0 = 0.5
set.seed(2)
# Simulated sample
y0 = rnorm(100) - rstable(n=100, alpha=alpha0 , beta=0, gamma = c0, delta = 0, pm = 0)
# A histogram of the sample
hist(y0)
# Second approximation of the density of Y
fy = function(y,c,sigma,ns){
z = rstable(n=ns, alpha=alpha0 , beta=0, gamma = 1, delta = 0, pm = 0)
return( mean(dnorm(c*z-y),mean=0,std=sigma ))
}
# -log likelihood
ll = function(par){
temp = rep(0,length(y0))
if(par[1]>0&par[2]>0){
for(j in 1:length(y0)) temp[j]=fy(y0[j],par[1],par[2],1000)
return(-sum(log(temp)))
}
else return(Inf)
}
# optimisation
optim(c(0.5,1),ll,control = list(maxit=500))
The estimators I got are $(\hat c,\hat\sigma)=(0.657,0.942)$, sort of close to the theoretical values.
You can easily play with this code to obtain the corresponding estimators for your sample.
I hope this helps.
A much slower (at least in my implementation) alternative to Procrastinator's answer is the brute-force method: form the likelihood as the convolution of a Gaussian and $\alpha$-stable distribution, calculated by numerical integration, then maximize it.
library(stabledist)
library(MASS)
library(stats4)
# True parameter values
alpha0 <- 1.75
sigma0 <- 1
c0 <- 0.5
set.seed(2)
# Simulated sample
Z <- rnorm(100) + rstable(n=100, alpha=alpha0 , beta=0, gamma = c0)
# -log likelihood
# (uses log of scale parameters as input, so range is (-Inf,Inf))
ll = function(lsigma, lsc) {
fconv <- function(x, zi, sigma, sc) {
dnorm(x, 0, sigma)*dstable(zi-x, alpha=alpha0, beta=0, gamma=sc)
}
sigma <- exp(lsigma)
sc <- exp(lsc)
f <- 0
for (zi in Z) {
f <- f + log(integrate(fconv, lower=-5*sigma, upper=5*sigma, zi=zi, sigma=sigma, sc=sc)$value)
}
-f
}
# optimisation
# Note: reltol should probably be set larger than the accuracy of integrate
# or you may have convergence problems
foo <- mle(ll, start=list(lsigma=0, lsc=log(0.5)),
control=list(reltol=4*(.Machine$double.eps^0.25)))
summary(foo)
... blah blah blah ...
Coefficients:
Estimate Std. Error
lsigma 0.1703237 0.9844097
lsc -0.6680191 2.0703820
-2 log L: 361.0518
> exp(foo@coef)
lsigma lsc
1.1856886 0.5127232
Runtime, however, is an issue; on my reasonably fast computer, this took about 25 minutes to run. Larger samples, or starting well away from the MLE, would no doubt take longer. You clearly wouldn't want to form confidence intervals by bootstrapping the procedure...