PerturbationDSGE.sample#

PerturbationDSGE.sample(theta0, *, n_particles=500, n_draws=2000, n_burn=500, thin=1, proposal_scale=None, seed=None)[source]#

Particle marginal Metropolis-Hastings under the specification’s prior.

Each iteration proposes in the unconstrained space, estimates the likelihood with a fresh particle filter, and accepts by the pseudo-marginal ratio, which targets the exact posterior whatever the particle count; the count only governs mixing.

Parameters:

  • theta0 (ArrayLike) – Starting structural parameters; a point estimate from fit() is the natural choice.

  • n_particles (int) – Particles per likelihood estimate. Pitt et al. (2012) suggest sizing it so the log-likelihood estimate’s standard deviation at the mode is near one.

  • n_draws (int) – Total iterations.

  • n_burn (int) – Burn-in, during which the proposal adapts.

  • thin (int) – Keep every thin-th post-burn draw.

  • proposal_scale (ArrayLike | None) – Initial proposal standard deviations in the unconstrained space, (d,), or a covariance (d, d); None for 0.1 per coordinate.

  • seed (int | None) – Seed.

Return type:

_PerturbationDSGEPosterior

Returns:

The _PerturbationDSGEPosterior with the kept draws, their log-likelihood and log-posterior values, the acceptance rate and convergence diagnostics.

Example

A short chain, far too short for inference, shows the shape of the surface:

>>> import numpy as np
>>> class Latent:
...     parameter_names, bounds = ("rho", "sigma"), ((0.0, 0.99), (0.01, 1.0))
...     n_states, n_controls = 1, 1
...     def equations(self, th, y_next, y, x_next, x):
...         return np.array([y[0] - x[0], x_next[0] - th[0] * x[0]])
...     def steady_state(self, th):
...         return np.zeros(1), np.zeros(1)
...     def shock_loading(self, th):
...         return np.array([[th[1]]])
...     def observation(self, th):
...         return np.array([[0.0, 1.0]]), np.zeros(1), np.array([[1e-4]])
...     def log_prior(self, th):
...         return 0.0
>>> rng = np.random.default_rng(0)
>>> x = np.zeros(120)
>>> for t in range(1, 120):
...     x[t] = 0.8 * x[t - 1] + 0.5 * rng.standard_normal()
>>> model = PerturbationDSGE(Latent(), x, order=2)
>>> res = model.fit([0.5, 0.3])
>>> post = model.sample(res.theta, n_particles=100, n_draws=40, n_burn=10, seed=0)
>>> post.n_kept, post.n_particles, bool(0.0 <= post.acceptance_rate <= 1.0)
(30, 100, True)