Re: [stan-users] General Stochastic Volatility models (with correlated shocks and jumps)

[Bob Carpenter]

This is another old message that nobody answered.
So I'll take a stab at it below.


> On Nov 4, 2016, at 2:16 PM, nuffe <erik...@gmail.com> wrote:
>
> I was thinking about attempting to try out Stan on a set of Stochastic Volatility models (SVm), originally fitted manually with MCMC from a Journal of Finance (JF) paper, The Impact of Jumps in Volatility and Returns. I know the Stan reference, and The No-U-Turn Sampler has an example of a simple SVm, but these are not in the same class as the ones considered in the JF paper.
>
> That paper consider the general class of models adhering to (equation (2) page 11)
>
> Y_t = Y_{t-1} + mu + sqrt(V_t) * epsy_t + xiy_t * Jy_t
> V_t = V_{t-1} + k * (theta - V_{t-1}) + sigmav * sqrt(V_t) * epsv_t + xiv_t * Jv_t
>
>
> i.e. Y is the measurement equation (returns), and V is the process equation (variance). This can characterise four SVm's commonly abbreviated as: SV, SVJ, SVCJ, and CVIJ,
>
> where
> • epsy_t and epsv_t are standard normal with correlation rho
> • Jy and Jv are the jump arrivals (=1 indicates an arrival), which are Bernoulli with constant intensities lambday and lambdav. For SVCJ the arrivals are contemporaneous
> • xiy_t and xiv_t are the jumps sizes. For SVCJ: xiv ~ exp(muv) and xiy | xiv ~N(muy + rhoj*xiv, sigmay)
> Is this something Stan should be able to do? If that is the case, and it is not too much to ask, what would a Stan model for SVCJ, and the simplest, SV:
>
> Y_t = Y_{t-1} + mu + sqrt(V_t) * epsy_t
> V_t = V_{t-1} + k * (theta - V_{t-1}) + sigmav * sqrt(V_t) * epsv_t
>
> look like?


You just need to translate to sampling notation, and you can use slicing:

Y_t[2:T] ~ normal(Y[1:(T - 1)] + mu, sqrt(V_t));
V_t[2:T] ~ normal(V[1:(T - 1)] + k * (theta - V[1:(T - 1)]), sigmav * sqrt(V_t));

Then you can just add all the priors you want.

>
> If epsy and epsx had been independent this would look more similar to the SV model from the Stan reference, but here every time point is characterized by a bivariate normal with time changing mean and covariance matrix, which I think in pseudo notation would look like:
>
> [Y_t, V_t] ~ multi_normal(MUt, VCMt)
> #where
> MUt = [Y_{t-1} + mu, V_{t-1} + k(theta-V_{t-1}]
> VCMt[1,1] = V_{t-1}
> VCM[1, 2] = rho * sigmav * V_{t-1}
> VCM[2, 2] = sigmav^2 * V_{t-1}
>
> And since all these models are discrete time versions of continuous ones, I assume one would have to specify bounds on V to avoid the variance turning negative?

Yes, if you need V to be positive it needs to be declared with <lower = 0>.
And you need to have lower=-1, upper=1 for rho and lower=0 for sigmav.

Then the most efficient thing to do is this:

vector[2] YV[T];
for (t in 1:T)
YV[t] = to_vector({Y_t[t], V_t[t]});

YV ~ multi_normal(MU_t, VCM);

You need to make sure to set VCM[2, 1] = VCM[1, 2].
I didn't follow your covariance matrix construction otherwise, I'm
afraid, so you need to make sure that's what you really want.

- Bob