```

This lecture. uses Bayesian methods offered by [pymc](https://www.pymc.io/projects/docs/en/stable/) and [numpyro](https://num.pyro.ai/en/stable/) to make statistical inferences about two parameters of a univariate first-order autoregression.

This lecture uses Bayesian methods offered by [pymc](https://www.pymc.io/projects/docs/en/stable/) and [numpyro](https://num.pyro.ai/en/stable/) to make statistical inferences about two parameters of a univariate first-order autoregression.

The model is a good laboratory for illustrating

consequences of alternative ways of modeling the distribution of the initial $y_0$:

- As a fixed number.

- As a fixed number

- As a random variable drawn from the stationary distribution of the $\{y_t\}$ stochastic process

The second component of the statistical model is

and

$$

y_0 \sim {\cal N}(\mu_0, \sigma_0^2)

$$ (eq:themodel_2)

First we'll use **pymc4**.

## `PyMC` Implementation

## PyMC Implementation

For a normal distribution in `pymc`,

$var = 1/\tau = \sigma^{2}$.

with AR1_model:

#start with priors

rho = pmc.Uniform('rho',lower=-1.,upper=1.) #assume stable rho

# Start with priors

rho = pmc.Uniform('rho', lower=-1., upper=1.) # Assume stable rho

sigma = pmc.HalfNormal('sigma', sigma = np.sqrt(10))

# Expected value of y at the next period (rho * y)

yhat = rho*y[:-1]

yhat = rho * y[:-1]

# Likelihood of the actual realization.

y_like = pmc.Normal('y_obs', mu = yhat, sigma=sigma, observed=y[1:])

# Likelihood of the actual realization

y_like = pmc.Normal('y_obs', mu=yhat, sigma=sigma, observed=y[1:])

```

```{code-cell} ipython3

The Hurwicz bias is worse the smaller is the sample (see {cite}`Orcutt_Winokur_69`.)

Be that as it may, here is more information about the posterior

Be that as it may, here is more information about the posterior.

```{code-cell} ipython3

with AR1_model:

summary = az.summary(trace, round_to=4)

# print

summary

```

with AR1_model_y0:

#start with priors

rho = pmc.Uniform('rho',lower=-1.,upper=1.) #assume stable rho

# Start with priors

rho = pmc.Uniform('rho', lower=-1., upper=1.) # Assume stable rho

sigma = pmc.HalfNormal('sigma', sigma=np.sqrt(10))

#standard deviation of ergodic y

y_sd = sigma/np.sqrt(1-rho**2)

# Standard deviation of ergodic y

y_sd = sigma / np.sqrt(1 - rho**2)

#yhat

yhat = rho*y[:-1]

y_data = pmc.Normal('y_obs', mu = yhat, sigma=sigma, observed=y[1:])

y0_data = pmc.Normal('y0_obs', mu = 0., sigma=y_sd, observed=y[0])

# yhat

yhat = rho * y[:-1]

y_data = pmc.Normal('y_obs', mu=yhat, sigma=sigma, observed=y[1:])

y0_data = pmc.Normal('y0_obs', mu=0., sigma=y_sd, observed=y[0])

```

```{code-cell} ipython3

with AR1_model_y0:

trace_y0 = pmc.sample(50000, tune=10000, return_inferencedata=True)

# grey vertical lines are the cases of divergence

# Grey vertical lines are the cases of divergence

```

```{code-cell} ipython3

with AR1_model:

summary_y0 = az.summary(trace_y0, round_to=4)

# print

summary_y0

```

We'll now repeat the calculations using `numpyro`.

## `Numpyro` Implementation

## Numpyro Implementation

```{code-cell} ipython3

"""

Plot trace and histogram

"""

# to np array

# To np array

rhos = sample['rho']

sigmas = sample['sigma']

rhos, sigmas, = np.array(rhos), np.array(sigmas)

fig, axs = plt.subplots(2,2, figsize=(17,6))

# plot trace

axs[0,0].plot(rhos) # rho

axs[1,0].plot(sigmas) # sigma

fig, axs = plt.subplots(2, 2, figsize=(17, 6))

# Plot trace

axs[0, 0].plot(rhos) # rho

axs[1, 0].plot(sigmas) # sigma

# plot posterior

axs[0,1].hist(rhos, bins=50, density=True, alpha=0.7)

axs[0,1].set_xlim([0,1])

axs[1,1].hist(sigmas, bins=50, density=True, alpha=0.7)

# Plot posterior

axs[0, 1].hist(rhos, bins=50, density=True, alpha=0.7)

axs[0, 1].set_xlim([0, 1])

axs[1, 1].hist(sigmas, bins=50, density=True, alpha=0.7)

axs[0,0].set_title("rho")

axs[0,1].set_title("rho")

axs[1,0].set_title("sigma")

axs[1,1].set_title("sigma")

axs[0, 0].set_title("rho")

axs[0, 1].set_title("rho")

axs[1, 0].set_title("sigma")

axs[1, 1].set_title("sigma")

plt.show()

```

sigma = numpyro.sample('sigma', dist.HalfNormal(scale=np.sqrt(10)))

# Expected value of y at the next period (rho * y)

yhat = rho*data[:-1]

yhat = rho * data[:-1]

# Likelihood of the actual realization.

y_data = numpyro.sample('y_obs', dist.Normal(loc=yhat, scale=sigma), obs=data[1:])

```{code-cell} ipython3

:tag: [hide-output]

# make jnp array

# Make jnp array

y = jnp.array(y)

# set NUTS kernal

# Set NUTS kernal

NUTS_kernel = numpyro.infer.NUTS(AR1_model)

# run MCMC

# Run MCMC

mcmc = numpyro.infer.MCMC(NUTS_kernel, num_samples=50000, num_warmup=10000, progress_bar=False)

mcmc.run(rng_key=random.PRNGKey(1), data=y)

```

```{code-cell} ipython3

def AR1_model_y0(data):

# set prior

# Set prior

rho = numpyro.sample('rho', dist.Uniform(low=-1., high=1.))

sigma = numpyro.sample('sigma', dist.HalfNormal(scale=np.sqrt(10)))

#standard deviation of ergodic y

y_sd = sigma/jnp.sqrt(1-rho**2)

# Standard deviation of ergodic y

y_sd = sigma / jnp.sqrt(1 - rho**2)

# Expected value of y at the next period (rho * y)

yhat = rho*data[:-1]

yhat = rho * data[:-1]

# Likelihood of the actual realization.

y_data = numpyro.sample('y_obs', dist.Normal(loc=yhat, scale=sigma), obs=data[1:])

```{code-cell} ipython3

:tag: [hide-output]

# make jnp array

# Make jnp array

y = jnp.array(y)

# set NUTS kernal

# Set NUTS kernal

NUTS_kernel = numpyro.infer.NUTS(AR1_model_y0)

# run MCMC

# Run MCMC

mcmc2 = numpyro.infer.MCMC(NUTS_kernel, num_samples=50000, num_warmup=10000, progress_bar=False)

mcmc2.run(rng_key=random.PRNGKey(1), data=y)

```

Look what happened to the posterior!

It has moved far from the true values of the parameters used to generate the data because of how Bayes Law (i.e., conditional probability)

is telling `numpyro` to explain what it interprets as "explosive" observations early in the sample

It has moved far from the true values of the parameters used to generate the data because of how Bayes' Law (i.e., conditional probability)

is telling `numpyro` to explain what it interprets as "explosive" observations early in the sample.

Bayes Law is able to generates a plausible likelihood for the first observation is by driving $\rho \rightarrow 1$ and $\sigma \uparrow$ in order to raise the variance of the stationary distribution.

Bayes' Law is able to generate a plausible likelihood for the first observation by driving $\rho \rightarrow 1$ and $\sigma \uparrow$ in order to raise the variance of the stationary distribution.

Our example illustrates the importance of what you assume about the distribution of initial conditions.

Our example illustrates the importance of what you assume about the distribution of initial conditions.