How to code the stochastic volatility model in numpyro

numpyro
1

Hello,
I am trying to model a stochastic volatility model that follows
y_t = \exp{(h_t/2)} \epsilon_t \\ h_{t+1}= \mu + \phi (h_{t} - \mu) \eta_t
where \epsilon_t, \eta_t is independent for all t = 1,...,T. y_t is stock return and h_t is a latent variable.

I have created simulation data that follows this model and estimated in numpyro, but the parameters are converging to wrong values.
I would appreciate it if you could tell me how to fix it.

The code is as follows

import numpyro
from numpyro.infer import MCMC, NUTS, Predictive
from numpyro import sample
import numpyro.distributions as dist
import jax
import jax.numpy as jnp
from jax import random
from numpyro.contrib.control_flow import scan

import  matplotlib.pyplot as plt

import numpy as np
def make_data(T, mu, phi, sigma2):
    ht = make_ht(T, mu, phi, sigma2)
    yt = make_yt_from_ht(ht)
    return yt, ht

def make_yt_from_ht(ht):

    eps = np.random.randn(ht.shape[0])
    yt = (np.exp(ht/2))*eps
    return yt

def make_ht(T, mu, phi, sigma2, option=0):

    h0 = np.random.normal(mu, np.sqrt(sigma2/(1-phi**2)))
    eta = np.random.randn(T,)

    ht = np.zeros((T,))
    for t in range(ht.shape[0]):
        if t == 0:
            ht[t] = mu + phi*(h0 - mu) + np.sqrt(sigma2)*eta[t]
        else:
            ht[t] = mu + phi*(ht[t-1] - mu) + np.sqrt(sigma2)*eta[t]

    return ht

def sv_model(ys, T):

    phi = sample("phi", dist.Uniform(-1, 1))
    sigma = sample("sigma", dist.Exponential(50))
    mu = sample("mu", dist.Normal(-10, 10))

    def trainsition(carry, t):
        h0 = carry
        h1 = sample("h", dist.Normal(mu + phi*(h0 - mu), sigma))
        y1 = sample("y", dist.Normal(0, jnp.exp(h1)), obs=ys[t])
        return (h1), (h1, y1)

    # h0 = sample("h_0", dist.Normal(mu, sigma/jnp.sqrt(1-phi**2)))
    # init_state = h0

    init_state = ys[0]
    with numpyro.handlers.condition(data={"y": ys[1:]}):
        _, (hs, ys) = scan(trainsition, (init_state), jnp.arange(1, T))

mu, phi, sigma = -10, 0.9, 0.5
T = 500
yt, true_ht = make_data(T, mu, phi, sigma**2)
true_ht = jnp.asarray(true_ht)
yt = jnp.asarray(yt)

nuts_kernel = NUTS(sv_model)
mcmc = MCMC(nuts_kernel, num_warmup=5000, num_samples=5000)
rng_key = random.PRNGKey(0)
mcmc.run(rng_key, yt, T)

plt.plot(true_ht)
plt.plot(mcmc.get_samples()["h"].mean(axis=0))
plt.plot()

The true parameters are \mu, \phi, \sigma = -10, 0.9, 0.5. The estimated values are -5.12, 0.86, 0.27. So something is wrong.

Thanks.

2

Hi @kurara, I’m not sure about your specific model, but you might compare your model with Pyro’s stochastic volatility example. Good luck!

1 Like
3

@fritzo Thank you for your reply. I will try to work on it some more!

Actually, I am using the sv model to practice the state space model. So Time Series Forecasting example is closer to sv model. But this example is a bit complicated for me as a beginner. So a simpler example of a state space model would be great for me. But anyway thanks for your reply!

4

You might want to reparameterize this site. Kind of

h1 = mu + phi * (h0 - mu) + sigma * sample("h_base", dist.Normal(0, 1))
1 Like
5

@fehiepsi Thank you for your reply. But the reparametrization did not work.
I also tried the numpyro.handlers.reparam but it converges to the wrong value.

6

Hmm, I think the result is expected. Please consider using improper prior or weak prior rather than strong prior (like Exp(50))

1 Like