PerturbationDSGE#

class cultivars.multivariate.structural.perturbation.PerturbationDSGE(specification, data, *, order=2)[source]#

A perturbation-solved DSGE model, estimated by likelihood.

The specification supplies the equilibrium conditions \(\mathbb{E}_t f(y_{t+1}, y_t, x_{t+1}, x_t; \theta) = 0\), the steady state, the shock loading, the measurement equation and a log prior, all as functions of the structural parameters; the engine solves the model at any \(\theta\) to first or second order by Schmitt-Grohé-Uribe perturbation with numerical derivatives, emits the solution as a state space on the substrate the order calls for, and estimates by whichever likelihood that substrate offers. The constructor is inherited from _PerturbationModel and validates the order, the data and the specification’s bounds; the parameters travel through the optimizer in an unconstrained transform of those bounds.

Parameters:

  • specification (_PerturbationModelSpecification) – An object with parameter_names, bounds, n_states, n_controls and the methods equations, steady_state, shock_loading, observation and log_prior; the protocol is _PerturbationModelSpecification in _internals.

  • data (ArrayLike) – The (nobs, p) observables, numpy.nan rows missing.

  • order (int) – Perturbation order, 1 or 2.

Raises:
Variables:

  • _spec – The specification.

  • _data – The validated (nobs, p) observables.

  • _order – The perturbation order.

  • _bounds – The (n_params, 2) parameter bounds.

See also

References

Fernández-Villaverde, J., & Rubio-Ramírez, J. F. (2007). Estimating macroeconomic models: A likelihood approach. Review of Economic Studies, 74(4), 1059-1087.

An, S., & Schorfheide, F. (2007). Bayesian analysis of DSGE models. Econometric Reviews, 26(2-4), 113-172.

Example

A stochastic growth model with log utility and full depreciation, whose exact policy is k' = alpha beta z k**alpha. Consumption is simulated from that policy, the second-order model is fitted through the unscented filter, and the solution’s first-order block carries the capital-accumulation and productivity dynamics:

>>> import numpy as np
>>> class Growth:
...     parameter_names = ("alpha", "rho", "sigma")
...     bounds = ((0.1, 0.6), (0.0, 0.999), (1e-4, 0.2))
...     n_states, n_controls, beta = 2, 1, 0.99
...     def equations(self, th, y_next, y, x_next, x):
...         a, rho, _ = th
...         (c,), (c_next,), (k, z), (k_next, z_next) = y, y_next, x, x_next
...         return np.array([
...             1 / c - self.beta / c_next * a * z_next * k_next ** (a - 1),
...             c + k_next - z * k ** a,
...             np.log(z_next) - rho * np.log(z),
...         ])
...     def steady_state(self, th):
...         k = (th[0] * self.beta) ** (1 / (1 - th[0]))
...         return np.array([k, 1.0]), np.array([k ** th[0] - k])
...     def shock_loading(self, th):
...         return np.array([[0.0], [th[2]]])
...     def observation(self, th):
...         return np.array([[0.0, 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)
>>> alpha, rho, sigma, beta = 0.33, 0.9, 0.03, 0.99
>>> k, z, consumption = (alpha * beta) ** (1 / (1 - alpha)), 1.0, []
>>> for t in range(300):
...     z = np.exp(rho * np.log(z) + sigma * rng.standard_normal())
...     output = z * k ** alpha
...     k = alpha * beta * output
...     consumption.append(output - k)
>>> model = PerturbationDSGE(Growth(), np.array(consumption[100:]), order=2)
>>> res = model.fit([0.3, 0.8, 0.05], filter="unscented")
>>> {name: round(value, 2) for name, value in res.params.items()}
{'alpha': 0.35, 'rho': 0.9, 'sigma': 0.02}
>>> own = np.diag(res.solution.h_x)
>>> bool(np.allclose(own, [res.params["alpha"], res.params["rho"]])), res.solution.is_stable
(True, True)

Attributes#

data

The observables.

n_params

Structural parameters.

order

The perturbation order.

parameter_names

Structural parameter labels.

Methods#

fit([theta0, filter, n_particles, seed])

Maximize one filter's likelihood from a starting point.

loglikelihood(theta, *[, filter, ...])

The log-likelihood (or its particle estimate) at one parameter point.

sample(theta0, *[, n_particles, n_draws, ...])

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

solve(theta)

Solve the model at a parameter point.

state_space(theta)

The solution at theta on the substrate its order calls for.