Hi, how to properly model circular dependence in Stochastic Volatility model? The problem is that \epsilon_t and h_t depend on each other and J_t is not normal.
\begin{aligned}
r_t &= \mu + \exp(\omega + \sigma h_t) \epsilon_t \\
\epsilon_t &\sim \mathcal N(0,1) \\
\\
h_t &= \phi h_{t-1} + \alpha(J_t + \lambda \epsilon_t) \\
J_t &= \exp(z_t - 1/2) - 1 \\
z_t &\sim \mathcal N(0,1) \\
\end{aligned}
What is a proper way to define such model?
model {
eps ~ std_normal();
z ~ std_normal();
{
real h = h0;
for (t in 1:T) {
// Heavy tailed asymmetrical innovation, centered.
real J = exp(z[t] - 0.5) - 1;
h = phi * h + alpha * (J + lambda * eps[t];);
real s = exp(omega + sigma * h);
// How to express r[t] = mu + s * eps[t]?
target += normal_lpdf(r[t] | mu, s); // <= Problem
}
}
}
It could be solved if model changed so that lagged \epsilon_{t-1} used or J_t made normal. But that changes the structure of the model and make it less realistic.
1 Like
Don’t make eps a parameter, just compute it by solving your first equation for eps.
1 Like
But it cant be solved analytically, its circular dependency, it can be solved for lagged \epsilon_{t-1} not for same time \epsilon_t.
I understand, the solution is to use probabilistic link instead of deterministic right?
\begin{aligned}
r_t &= \mu + \exp(\omega + \sigma h_t) \epsilon_t \\
h_t &= \phi h_{t-1} + \alpha J_t + \gamma \eta_t \\
J_t &= \exp(z_t - 1/2) - 1 \quad z_t \sim \mathcal N(0,1) \\
(\eta_t, \epsilon_t) &\sim \mathcal N_2(0,1, \rho) \\
\end{aligned}
or equivalently, rewritten same way as will be implemented in stan
\begin{aligned}
r_t &= \mu + \exp(\omega + \sigma h_t)\bigl(\rho \eta_t + \sqrt{1-\rho^2}\,\xi_t\bigr) \\
h_t &= c + \phi h_{t-1} + \alpha J_t + \gamma \eta_t \\
J_t &= \exp(z_t - 1/2) - 1 \\
z_t &\sim \mathcal N(0,1), \qquad
\eta_t \sim \mathcal N(0,1), \qquad
\xi_t \sim \mathcal N(0,1)
\end{aligned}
and the stan code would be
model {
z ~ std_normal();
eta ~ std_normal();
real h = h0;
for (t in 1:T) {
real J = exp(z[t] - 0.5) - 1;
h = phi*h + alpha*J + beta*eta[t];
real s = exp(omega + sigma*h);
target += normal_lpdf(r[t] |
mu + s*rho*eta[t],
s*sqrt(1 - rho^2)
);
}
}
The performance cost - we have the additional serie of latent variable \eta_t \sim N(0, 1) for MCMC to handle, there’s no way to avoid that?
And complexity cost - one extra parameter \rho although it could be fixed by strong prior close to -1.
I think what @jsocolar implies is that your model description entails that given r_t, \mu, \omega, \sigma, h_t we have the deterministic relationship:
\epsilon_t = \frac{r_t - \mu}{\exp(\omega + \sigma h_t)}\\
So you should be able to have something like:
model {
z ~ std_normal();
{
real h = h0;
for (t in 1:T) {
// Heavy tailed asymmetrical innovation, centered.
real J = exp(z[t] - 0.5) - 1;
real s = exp(omega + sigma * h);
real eps_t = (r[t] - mu)/s
h = phi * h + alpha * (J + lambda * eps_t);
target += std_normal_lpdf(eps_t);
}
}
}
If there is a reason, you cannot do this, please elaborate on the constraints you are facing.
Your stan code is correct, but it’s a different model with lagged dependency h_t \sim \epsilon_{t-1} this:
\epsilon_t=\frac{r_t-\mu}{exp(\omega + \sigma (h_t = \phi h_{t-1}+\alpha(J_t + \lambda \epsilon_{t-1})))}
The problem is to use same day dependency h_t \sim \epsilon_{t} I don’t know how to calculate it analytically, note the \epsilon_t on both sides of the equation:
\epsilon_t=\frac{r_t-\mu}{exp(\omega + \sigma (h_t = \phi h_{t-1}+\alpha(J_t + \lambda \epsilon_{t})))}
These are different models:
- Lagged - negative return on day t causes higher volatility on the next days.
- Same day - negative return on day t cause volatility increase on same day and next days.
If you plot left and right tails for same day model on LogLog plot the negative returns will be slightly larger. Lagged model can’t do that, left and right returns will be the same.
Oh I see, I missed that the equation in the first post is different. I think you should still be able to solve for \epsilon_t
Substituting for h_t in r_t = \mu + \exp(\omega + \sigma h_t) \epsilon_t gives
r_t = \mu + \exp(\omega + \sigma \phi h_{t-1} + \alpha(J_t + \lambda \epsilon_t)) \epsilon_t
Taking c = \omega + \sigma \phi h_{t-1} + \alpha J_t and d = \alpha\lambda we have
r_t = \mu + \exp(c + d \epsilon_t) \epsilon_t
And my friend Wolfram tells me this has solution in terms of the Lambert W_0 function, which is available in Stan (though it might cost you some performance)
\epsilon_t = \frac{W_0(d \exp(-c) (r-u))}{d}
(there are some additional conditions where you may get zero or two solutions - I assume those are impossible due to the way you set this up, but see the Wolfram output for more details).
Thank you, interesting idea. One thing that worries me - there’s the lookahead. The observation r_t used to calculate \epsilon_t which is used in latent state h_t which is then used to calculate the likelihood. Shouldn’t the likelihood be adjusted to account for lookahead?
I think you are highlighting a problem with the math formulation of your model. The solution I proposed is just what your math model entails (i.e. that once you know h_t and the other parameters, r_t is deterministic and vice versa). If the solution feels uncomfortable, something in the original math has to be wrong. I’ve never worked with any model related to Stochastic Volatility so can’t really help with the interpretation.
I also forgot to add a Jacobian adjustment, which would be needed due to the transform.
