This page describes how to use Python libraries for macroeconomic modelling, allowing to extend Dynare functionalities.

Python ressources

* Formal calculus

• [http://code.google.com/p/sympy/ sympy ] python library to perform formal calculus.

• [http://www.sagemath.org/ sage ] (includes sympy) python wrapper around many opensource libraries which converts all of them into a single webbased interface. Is also able to delegate calculus to Maple or Mathematica.

* Numerical computing

• [http://numpy.scipy.org/ numpy] : contains an efficient implementation of matrix processing whose performance is close to matlab

• [http://www.scipy.org/ scipy] : set of scientific routines including optimization routines

• [http://www.cosc.canterbury.ac.nz/greg.ewing/python/Pyrex/ pyrex] and [http://psyco.sourceforge.net/ psyco ] (included in scipy) can speedup computations. The first allows to write C extensions in a python-like syntax, while the second performs just-in-time compilation of python functions.

* Other

• [http://epydoc.sourceforge.net/ epydoc ] generates online documentation from inline documentation

• [http://mlabwrap.sourceforge.net/ mlabwrap ] uses the Matlab Engine to interact with Matlab. It permits to send and retrieve data, as well as to call Matlab function from Python.

• [http://rpy.sourceforge.net/ rpy] : low-level access to statistical program R permitting to call every R functions as python functions.

• [http://www.py2exe.org/ py2exe] packs a python program with all its libraries in a single windows executable which includes the interpreter.

• [http://code.google.com/p/pymaclab/ pymaclab] another attempt to use python estimate DSGE models. It uses its own Modfile format which is then parsed by Python.

• [http://www.wiwi.uni-frankfurt.de/profs/nautz/schreibersoftware.html qz.py] : a wrapper around LAPACK's qz routine.

Writing/extending modfiles in Python

Here is the code of ramsey.py, which translates ramst.mod into python code. It is run by the standard Python interpreter. Line-by-line explication follows.

   1 ####################
   2 # Model Definition #
   3 ####################
   4 
   5 from daredare import *
   6 
   7 set_variables("k c")
   8 set_exovariables("z")
   9 set_shocks("x")
  10 set_parameters("alpha beta gamma delta aa")
  11 
  12 equations = [
  13     Equation(c + k , (1 - delta)*k(-1) + k(-1)**alpha*(1 + z)*aa),
  14     Equation(0 , -1/(1 + beta)*c(1)**(-gamma)*(1 - delta + k**(-1 + alpha)*(1 + z(1))*aa*alpha) + c**(-gamma)),
  15     Equation(z(1),x(1))
  16 ]
  17 
  18 parameters_values = {
  19     alpha : 0.5,
  20     gamma : 1*alpha,
  21     delta : 0.02,
  22     beta : 0.05,
  23     aa : 0.5
  24 }
  25 
  26 covariances = matrix([[0.01]])
  27 
  28 init_values = {
  29     k:(1/aa/alpha*(beta + delta))**(1/(-1 + alpha)),
  30     c:k**alpha*aa - delta*k
  31 }
  32 
  33 ####################
  34 # Model processing #
  35 ####################
  36 
  37 model = Model("ramsey",lookup=True)
  38 model.check_all('uhlig')
  39 
  40 solver = Solver(model,lookup=True)
  41 
  42 ss = solver.find_steady_state()
  43 print(ss)
  44 
  45 res = solver.solve_with_uhligs_toolkit()
  46 
  47 
  48 print('PP',res.PP)
  49 print('QQ',res.QQ)

Language extension (for model definition)

On the very first line of this code from daredare import * imports the library. The second block

   1 set_variables("k c")
   2 set_exovariables("z")
   3 set_shocks("x")
   4 set_parameters("alpha beta gamma delta aa")

creates formal objects "k,c,z" which all derive of the same class Variable. This class is a subclass of the standard Sympy class Symbol, which only adds the possibility to take lags or leads with a function call : for instance k(1) returns another Variable whose formal name is k_{+1}. k.lag returns the lag so that k(-k.lag))}} or equivalently {{{k.p() will always return a variable with lag 0. Parameters are of the class Parameter which doesn't have the same facility. It should be noted that, a priori, the formal object, and the Python reference name which is used by the compiler are to different things. For instance,

   1 somestrangename = Variable('x')
   2 anotherstrangename = somestrangename
   3 y = Variable('y')
   4 eq = somestrangename * anotherstrangename / y
   5 print(eq)

will print : x**2/y. In this respect the functions like set_variables are notable because instead of returning the created objects, they also add references to the interpreter context, with the same names as the variables. They also add the python names variables,shocks,exovariables,parameters. The next block

   1 equations = [
   2     Equation(c + k , (1 - delta)*k(-1) + k(-1)**alpha*(1 + z)*aa),
   3     Equation(0 , -1/(1 + beta)*c(1)**(-gamma)*(1 - delta + k**(-1 + alpha)*(1 + z(1))*aa*alpha) + c**(-gamma)),
   4     Equation(z(1),x)
   5 ]

defines equations. The object Equation is a child of Sympy's Equality object. It contains the expression as well as some attributes as (optionally) equation name, expectational nature, and other information that could be computed later. Three remarks :

• unfortunately we cannot write x+y == z because in python == must always return true or false (see Sympy mailing list).

• to convert from/to matlab, you only have to convert * from/to ^. For instance, if eq is a python formal expression, str(eq).replace('==','=').replace('**','^') will give the valid Dynare expression.

• every sympy operations on equalities are available

  • printing.latex(eq) (or printing.mathml(eq) will return the latex (or mathml) representation of the equality

  • preview(eq) will launch an viewer to preview the latexified output

  • eq.rhs.diff(c) will return the right side of the equation differentiated with respect to c

  • eq.rhs.subs(c(1),c) will replace c(1) by c in the right side

Processing model(s)

In the next lines we create the model and check its consistency :

   1 model = Model("ramsey",lookup=True)
   2 model.check_all('uhlig')

Note the option lookup=True. Normally, in the scope of a function, you can only access local variables and arguments. model function uses python's introspection routines to lookup in the stack for the lists variables,exovariables,shocks,parameters and equations and to attach them to the object model. We have already used this magic in the functions set_variables. Then, equations can be accessed with model.equations.

The second line checks for model consistency. So far the formal objects we have created are totally agnostic to modelling conventions which vary across software and users. The instance function check_all makes some basic checks to see if the model is consistent with one formulation. As for now, it accepts two arguments 'uhlig' and 'dynare'. Note that exovariables as it is defined in Uhlig's conventions don't exist in Dynare. shocks denote independent random variables and correspond to the varexo statement in Dynare.

   1 solver = Solver(model,lookup=True)
   2 res = solver.solve_with_uhligs_toolkit()

A solver object depends on a model and maps a set of parameters to a decision rule. Here, we have only one value for this set, so we can get them from the initial declaration using lookup=True. The second statement is useless when Dynare is the destination. It returns a dictionary containing the steady state. Many steps are involved :

• Initial values and initial parameters are computed from the formal declarations. This involve, solving a triangular system, which is easy.
• Formal equations are converted into a numeric python function which is added to interpreting context at runtime.
• This function is optimized to return the steady state

   1 res = solver.solve_with_uhligs_toolkit()
   2 print('PP',res.PP)
   3 print('QQ',res.QQ)

The last lines compute the decision rule. In the case of Uhlig, it returns two matrix representing the decision rules.

If the model were written following Dynare convention, we would have called :

   1 res = send_to_octave(solver,interactive=False,append="steady;\ncheck;\nstoch_simul(nograph);")
   2 M = res['M_']
   3 oo = res['oo_']
   4 options = res['options_']

to get the results from octave or

   1 res = send_to_matlab(solver,interactive=False,append="steady;\ncheck;\nstoch_simul(nograph);")
   2 print(res)
   3 M = res['M_']
   4 oo = res['oo_']
   5 options = res['options_']

to get the results from Matlab. For Octave, we write a modfile and send it to a new session while we use the Matlab engine to send it to Matlab with the added benefit that we can do many calculations without clearing the Matlab workspace each time.

Import and export