Model API¶
For numerical purposes, models are essentially represented as a set of symbols, calibration and functions representing the various equation
types of the model. This data is held in a Model object whose API is described in this chapter. Models are usually created by writing a Yaml files as described in the the previous chapter, but as we will
see below, they can also be written directly as long as they satisfy the requirements detailed below.
Model Object¶
As previously, let's consider, the Real Business Cycle example, from the introduction. The model object can be created using the yaml file:
The object contains few meta-data:
display( model.name ) # -> "Real Business Cycles"
display( model.infos ) # -> {"name": "Real Business Cycle", "filename": "examples/models/rbc.yaml", "type": "dtcc"}
Calibration¶
Each models stores calibration information as model.calibration. It is a special dictionary-like object, which contains calibration information, that is values for parameters and initial values (or steady-state) for all other variables of the model.
It is possible to retrieve one or several variables calibrations:
display( model.calibration['k'] ) # -> 2.9
display( model.calibration['k', 'delta'] # -> [2.9, 0.08]
When a key coresponds to one of the symbols group, one gets one or several vectors of variables instead:
model.calibration['states'] # - > np.array([2.9, 0]) (values of states [z, k])
model.calibration['states', 'controls'] # -> [np.array([2.9, 0]), np.array([0.29, 1.0])]
To get regular dictionary mapping states groups and vectors, one can use the attributed .grouped
The values are vectors (1d numpy arrays) of values for each symbol group. For instance the following code will print the calibrated values of the parameters:
In order to get a (key,values) of all the values of the model, one can call model.calibration.flat.
Note
The calibration object can contain values that are not symbols of the model. These values can be used to calibrate model parameters and are also evaluated in the other yaml sections, using the supplied value.
One uses the model.recalibrate() routine to derive a new model with updated calibration values (models in Dolo are immutable):
# Recalibrate returns a new, independent Model instance
new_model = model.recalibrate({'delta': 0.01})
new_model = model.recalibrate(delta=0.08, k=2.8)
This method also understands symbolic expressions (as string) which makes it possible to define symbols as a function of other symbols:
new_model = model.recalibrate(beta='1/(1+delta)')
print(new_model.calibration['beta'])
new_model = model.recalibrate(delta=0.04)
print(new_model.calibration['beta', 'delta'])
Under the hood, the method resolves symbolic relations between symbols. In particular, the calibration order is irrelevant as long as all parameters can be deduced one from another.
Functions¶
A model of a specific type can feature various kinds of functions. For instance, a continuous-states-continuous-controls models, solved by iterating on the Euler equations may feature a transition equation \(g\) and an arbitrage equation \(f\). Their signature is respectively \(s_t=g(m_{t-1},s_{t-1},x_{t-1},m_t)\) and \(E_t[f(m_t,s_t,x_t,s_{t+1},x_{t+1},m_{t+1})]\), where \(s_t\), \(x_t\) and \(m_t\) respectively represent a vector of states, controls and exogenous shock. Implicitly, all functions are also assumed to depend on the vector of parameters \(p\).
These functions can be accessed by their type in the model.functions dictionary:
Let's call the arbitrage function on the steady-state value, to see the residuals at the deterministic steady-state:
m = model.calibration["exogenous"]
s = model.calibration["states"]
x = model.calibration["controls"]
p = model.calibration["parameters"]
res = f(m,s,x,m,s,x,p)
display(res)
The output (res) is two element vector, representing the residuals of the two arbitrage equations at the steady-state. It should be full of zero. Is it ? Great !
By inspecting the arbitrage function ( f? ), one can see that its call api is:
Since m, s and x are the short names for exogenous shocks, states and controls, their values at date \(t+1\) is denoted with S and X. This simple convention prevails in most of dolo source code: when possible, vectors at date t are denoted with lowercase, while future vectors are with upper case. We have already commented the presence of the parameter vector p.
Now, the generated functions also gives the option to perform in place computations, when an output vector is given:
It is also possible to compute derivatives of the function by setting diff=True. In that case, the residual and jacobians with respect to the various arguments are returned as a list:
Since there are two states and two controls, the variables r_s, r_x, r_S, r_X are all 2 by 2 matrices.
The generated functions also allow for efficient vectorized evaluation. In order to evaluate the residuals :math:N times, one needs to supply matrix arguments, instead of vectors, so that each line corresponds to one value to evaluate as in the following example:
N = 10000
vec_m = m[None,:].repeat(N, axis=0) # we repeat each line N times
vec_s = s[None,:].repeat(N, axis=0) # we repeat each line N times
vec_x = x[None,:].repeat(N, axis=0)
vec_X = X[None,:].repeat(N, axis=0)
vec_p = p[None,:].repeat(N, axis=0)
# actually, except for vec_s, the function repeat is not need since broadcast rules apply
vec_s[:,0] = linspace(2,4,N) # we provide various guesses for the steady-state capital
vec_S = vec_s
out = f(vec_m, vec_s,vec_x,vec_M, vec_S,vec_X,vec_p) # now a 10000 x 2 array
out, out_m, out_s, out_x, out_M, out_S, out_X = f(vec_m, vec_s,vec_x, vec_m, vec_S,vec_X,vec_p)
The vectorized evaluation is optimized so that it is much faster to make a vectorized call rather than iterate on each point.
Note
In the preceding example, the parameters are constant for all evaluations, yet they are repeated. This is not mandatory, and the call f(vec_m, vec_s, vec_x, vec_M, vec_S, vec_X, p) should work exactly as if p had been repeated along the first axis. We follow there numba's guvectorize conventions, even though they slightly differ from numpy's ones.
Exogenous shock¶
The exogenous field contains information about the driving process. To get its default, discretized version, one can call model.exogenous.discretize().
The model.options structure holds an information required by a particular solution method. For instance, for global methods, model.options['grid'] is supposed to hold the boundaries and the number nodes at which to interpolate.
Discretized Model and Transitions¶
For global nonlinear solution methods, continuous models are converted into a discrete numerical representation called a DiscretizedModel.
DiscretizedModel¶
A DiscretizedModel encapsulates the triplet \(\langle \text{Model}, \; \mathcal{G}, \; \tau \rangle\):
from dolo import yaml_import
from dolo.numeric import DiscretizedModel
model = yaml_import("examples/models/rbc.yaml")
dmodel = DiscretizedModel(model)
Attributes and methods:
- dmodel.grid: The computational domain grid \(\mathcal{G}\) (e.g. ProductGrid of discrete exogenous nodes and Cartesian endogenous nodes).
- dmodel.tau: The transition operator \(\tau\) generating next-state branching distributions.
- dmodel.initial_guess(): Returns a GArray of calibrated control values across all grid nodes.
- dmodel.lower_bounds(), dmodel.upper_bounds(): Returns GArray lower and upper bounds for controls.
- dmodel.transition_matrix(dr): Computes the full row-stochastic Markov transition matrix \(P\) on grid nodes given a solved decision rule dr.
- dmodel.ergodic_distribution(dr): Solves for the stationary distribution \(\mu P = \mu\).
Qualified Points and Batches (QualifiedPoint / QP)¶
To separate the computational location (address) from evaluated economic values, Dolo uses qualified points:
QualifiedPoint(aliasQP): Bundlesloc(e.g. regime index or discrete locator) andval(continuous economic coordinates \([m, s]\)).QualifiedBatch(aliasSliceQP): Bundlesloc,exo, andendofor vectorized operations over an entire continuous slice.
GArray (Grid-Bound Fields)¶
A GArray binds an array of data to a grid \(\mathcal{G}\):
- Supports geometric access: ga[i_m] gives the slice for discrete regime \(i_m\).
- Supports numerical solver duality: ga.ravel() provides a flat 1D view, and ga.unravel(flat_array) reconstructs the structured field.
- Provides vector arithmetic and maximum-norm calculation: ga.norm().