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.