I want to estimate the volatility $\sigma$ of a process $(X_t)$ following an arithmetic Brownian motion, that is, for a constant time step $\Delta$, $X_{t+\Delta} = X_t + \sigma B_{\Delta}$ , where $B_t$ is a Brownian motion.
Problem: Results from a straightforward implementation with a realized $\sigma_{real} = 0.2$, an initial $\sigma^{(1)}=0.1$, 10000 iterations and a data vector length of $N=$1000, however, show an always increasing $\sigma$ diverging to $+\infty$.
In contrast, if we choose the transformation $g_2(X;\sigma) = \sigma X$ instead, the approach yields acceptable results. Also, including a Jacobian term $J$ along
\begin{equation}
p(Y|\sigma, X^{(k+1)} ) = p(X^{(k+1)}|\sigma, X^{(k+1)}) \; |J| = p(X^{(k+1)}|\sigma) \; |J| ,
\end{equation}
did not improve results.
Question: What is not right?
It could be that in the expression for $p(\sigma | X^{(k+1)} , Y)$ in the Dirac term, it should be $\sigma$ (which we want to sample) instead of $\sigma_k$.
If so, I do not see how to sample from this distribution. Any ideas?
Approach in Detail: Defining $dx_i = x_i - x_{i-1}$, the likelihood of observations $X= (x_1, ..., x_N)$ at times $(\Delta, 2\Delta, ..., N \Delta)$ conditional on $\sigma$ is:
$L(\sigma|X) = p(X|\sigma) = \prod_{i=2}^N n(\Delta x_i; 0,\sigma^2 \Delta) = (\sqrt{2 \pi \Delta} \, \sigma^{N-1})^{-1} \exp\left(-\frac{1}{2 \sigma^2 \Delta} \sum_{i=2}^N dx_i^2 \right) .$
We do not observe $X$ itself but only discrete transformed data $Y = (y_1, ..., y_N)$, which are related to $X$ by $Y = g_1(X; \sigma) = X / \sigma$. We simulate a path of $X$, compute the transformed data $Y$ and try to estimate the parameter $\sigma$.
A first candidate approach to do this would be a Gibbs sampler, which generates a chain of pairs $(X^{(k)}=(x^{(k)}_1, ..., x^{(k)}_N), \sigma^{(k)})$, where in each iteration $k$, new values of $X$ and $\sigma$ are drawn from the conditional distributions:
$$X^{(k+1)} \sim p(X | \sigma^{(k)}, Y) \quad \text{and} \\
\sigma^{(k+1)} \sim p(\sigma | X^{(k+1)}, Y) \,.$$
In the above set-up, $p(X | \sigma^{(k)}, Y) = \delta(X - g_1(Y;\sigma^{(k)})) = \delta(X-\sigma^{(k)}Y)$ with $\delta$ denoting the Dirac delta, which indicates that drawing $X^{(k+1)}$ from $p(X | \sigma^{(k)}, Y)$ can be achieved by computing $X^{(k+1)} = \sigma^{(k)} Y$.
For the second conditional distribution it looks like we could write
\begin{equation}
p(\sigma | X^{(k+1)} , Y) \propto p(Y|\sigma, X^{(k+1)} ) \; p(\sigma|X^{(k+1)} ) = \delta(X^{(k+1)} - Y \sigma^{(k)}) \; L(X^{(k+1)}|\sigma) = L(X^{(k+1)}|\sigma)
\end{equation}
since in the step before $X^{(k+1)}$ was computed as $ X^{(k+1)} = Y \sigma^{(k)}$; in this case, it can be seen that we can sample $s=\sigma^2$ from an inverse gamma distribution derived from the likelihood of $X$.
