PyMacLab Tutorial Series

Tutorial 1 - Getting started

Introduction

PyMacLab’s strength or orginal design goal has been that of providing users with a rich and flexible DSGE data structure (also called class in object-oriented programming speak) which allows them to do lots of interesting things with DSGE models. Don’t forget that another stated goal of PyMacLab is to permit researchers to work with DSGE models who know how to and enjoy to use the programming language Python. This tutorial assumes that you already know how to use Python, if you are completely new to it I suggest you consult one of the many tutorials available on the internet to get started with the language. Don’t forget that Python is free as in free as freedom, no proprietary software such as Matlab is tainting the freedom of the environment you will be working in when using PyMacLab.

The easiest and best way to get started and explore all of the available features of PyMacLab, or for that matter any other Python library, is to work with IPython, which is an abbreviation for Interactive Python. With this program, Python users who are scientifically inclined, can more or less emulate the behaviour of programs such as Matlab, Octave or Scilab. This means that sessions can be launched in which commands entered into a command prompt are executed immediately and resulting output is printed to screen. It also means that generated objects resulting from such computations carried out at the command prompt (such as matrices or DSGE models) begin to exist within the session’s variable scope and become ready for further manipulation, i.e. they persist. Before being able to work with PyMacLab we have to launch an IPython session and import the library into it’s scope:

# Import the pymaclab module into its namespace
In [1]: import pymaclab as pm

# Get the version and author's name
In [2]: pm.__version__
'0.95.1'

# Get the library's author's name
In [3]: pm.__author__
'Eric M. Scheffel'

Here we simply have imported the PyMacLab module and inspected some of its attributes, such as the current version numbering as well as them module’s author’s name. Let’s look deeper into the recesses of the module though to better understand who it is organized

# Import the pymaclab module into its namespace
In [1]: import pymaclab as pm

# Use the dir() command to view all available attributes and method calls,
# this command returns a list
In [2]: dir(pm)
['OPS',
'__builtins__',
'__doc__',
'__file__',
'__name__',
'__package__',
'__path__',
'dattrans',
'db_graph',
'dsge',
'explain',
'favar',
'filters',
'ldbs',
'linalg',
'lmods',
'lvars',
'macrolab',
'make_modfile',
'modedit',
'modfiles',
'modinfo',
'modsolve',
'newDB',
'newFAVAR',
'newMOD',
'newVAR',
'pyed',
'stats',
'sys',
'texedit',
'var']

As you can see the module contains a quite a few attributes, many of which are still experimental and perhaps best not called at this stage. The most mature and arguable most interesting method call is that called pm.newMOD, which allows users to instantiate a DSGE model instance, which would be done like so:

# Import the pymaclab module into its namespace, also import os module
In [1]: import pymaclab as pm
In [2]: import os

# Define the relative path to your modfiles
In [3]: modpath = "../pymaclab/modfiles/models/stable"

# Instantiate a new DSGE model instance like so
In [4]: rbc1 = pm.newMOD(os.path.join(modpath,"rbc1_res.txt"))

# As an example, check the models computed steady stated
In [5]: print rbc1.sstate
{'betta': 0.99009900990099009,
'c_bar': 2.7560505909330626,
'k_bar': 38.160700489842398,
'y_bar': 3.7100681031791227}

Alternatively, you can also test some of the DSGE model files which come supplied with PyMacLab’s standard installation. For this to work all you have to do is to import a provided handler module, pymaclab.modfiles.models, which contains all of the DSGE models’ names and their correspoding full file paths. Notice however that the models themselves are further classified into three categories, models.stable, models.testing and models.development which helps to distinguish between models which are in the process of being added and such which are known to work correctly:

 # Import the pymaclab module into its namespace, also import os module
 In [1]: import pymaclab as pm
 # Import the DSGE models' filepath handle
 In [2]: from pymaclab.modfiles import models

 #Check all of the available models in the stable branch
 In [3]: dir(models.stable)
 ['__builtins__',
 '__doc__',
 '__file__',
 '__name__',
 '__package__',
 '__path__',
 'jermann98',
 'jermann98_ext',
 'merz',
 'prog',
 'rbc1_cf',
 'rbc1_ext',
 'rbc1_extss',
 'rbc1_focs',
 'rbc1_num',
 'rbc1_res',
 'rbc1_sug',
 'rbc2',
 'sims']

 #Check all of the available models in the development branch
 In [4]: dir(models.development)
['RBC_Romer',
 '__builtins__',
 '__doc__',
 '__file__',
 '__name__',
 '__package__',
 '__path__',
 'max1',
 'max2',
 'mbc1',
 'model2',
 'model3',
 'nk_nocapital',
 'nkm',
 'nkm_nocapital']

 #Check all of the available models in the testing branch
 In [5]: dir(models.testing)
 ['__builtins__',
 '__doc__',
 '__file__',
 '__name__',
 '__package__',
 '__path__',
 'cee']


 # The DSGE models objects in pymaclab.modfiles.models
 # are just references to full file paths, i.e.

 In [6]: models.stable.rbc1_res
 '/usr/lib/python2.7/site-packages/pymaclab/modfiles/rbc1_res.txt'

 #Instantiate a new DSGE model instance like so
 In [7]: rbc1 = pm.newMOD(models.stable.rbc1_res)

 #As an example, check the models computed steady stated
 In [8]: print rbc1.sstate
 {'betta': 0.99009900990099009,
 'c_bar': 2.7560505909330626,
 'k_bar': 38.160700489842398,
 'y_bar': 3.7100681031791227}

Now we have already seen some of the power and simplicity we can leverage by using PyMacLab. Before learning some of its additional power, we do however still need to take a quick detour to study the model file rbc1.txt which we had to pass as an argument to the pm.newMOD method call, as its structure is determined by a set of conventions which are important to adhere to in order to enable PyMacLab to parse and employ the information contained therein correctly and efficiently.

Tutorial 2 - The PyMacLab model file

The PyMacLab DSGE model file

In order to be able to load or instantiate your first DSGE model and work with it, you have to make sure to first fill in a so-called PyMacLab DSGE model file. The idea behind this is the same as with the Dynare model file which typically ends with a .mod file suffix. PyMacLab already comes provided with a number of such files pre-compiled for you to experiment with. For instance the most basic real business cycle model is described in the model file rbc1_num.txt (and in many other similar files which demonstrate the different ways of computing the steady state). This file looks as follows:

%Model Description+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
This is just a standard RBC model, as you can see.


%Model Information+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Name = Standard RBC Model, NUM-SS;


%Parameters++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
rho       = 0.36;
delta     = 0.025;
betta     = 1.0/1.01;
eta       = 2.0;
psi       = 0.95;
z_bar     = 1.0;
sigma_eps = 0.052;


%Variable Vectors+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[1]  k(t):capital{endo}[log,bk]
[2]  c(t):consumption{con}[log,bk]
[4]  y(t):output{con}[log,bk]
[5]  z(t):eps(t):productivity{exo}[log,bk]
[6]  @inv(t):investment[log,bk]
[7]  @R(t):rrate

%Boundary Conditions++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
None


%Variable Substitution Non-Linear System++++++++++++++++++++++++++++++++++++++
[1]   @inv(t)   = k(t)-(1-delta)*k(t-1);
[2]   @inv_bar  = SS{@inv(t)};
[2]   @F(t)     = z(t)*k(t-1)**rho;
[2]   @Fk(t)    = DIFF{@F(t),k(t-1)};
[2]   @Fk_bar   = SS{@Fk(t)};
[2]   @F_bar    = SS{@F(t)};
[3]   @R(t)     = 1+DIFF{@F(t),k(t-1)}-delta;
[4]   @R_bar    = SS{@R(t)};
[3]   @R(t+1)   = FF_1{@R(t)};
[4]   @U(t)     = @I{eta!=1.0}{c(t)**(1-eta)/(1-eta)}+@I{eta==1.0}{LOG(c(t))};
[5]   @MU(t)    = DIFF{@U(t),c(t)};
[5]   @MU_bar   = SS{@U(t)};
[6]   @MU(t+1)  = FF_1{@MU(t)};



%Non-Linear First-Order Conditions++++++++++++++++++++++++++++++++++++++++++++
# Insert here the non-linear FOCs in format g(x)=0

[1]   @F(t)-@inv(t)-c(t) = 0;
[2]   betta*(@MU(t+1)/@MU(t))*@R(t+1)-1 = 0;
[3]   @F(t)-y(t) = 0;
[4]   LOG(E(t)|z(t+1))-psi*LOG(z(t)) = 0;


%Steady States [Closed Form]+++++++++++++++++++++++++++++++++++++++++++++++++++
None


%Steady State Non-Linear System [Manual]+++++++++++++++++++++++++++++++++++++++
[1]   @F_bar-@inv_bar-c_bar = 0;
[2]   y_bar-@F_bar = 0;
[3]   betta*@R_bar-1 = 0;
[4]   betta*R_bar-1 = 0;

[1]   c_bar = 1.0;
[2]   y_bar = 1.0;
[2]   k_bar = 1.0;
[3]   R_bar = 1.01;

%Log-Linearized Model Equations++++++++++++++++++++++++++++++++++++++++++++++++
None


%Variance-Covariance Matrix++++++++++++++++++++++++++++++++++++++++++++++++++++
Sigma = [sigma_eps**2];


%End Of Model File+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

So what does this file mean, and in particular, what is the meaning and purpose of the individual sections? These and related questions are addressed in the sections to follow below. They mostly discuss the syntax conventions model builders have to adhere to in order to use PyMacLab correctly.

A Description of the model file’s individual sections

Model Description Section

In the model description section of the model file you can use plain text in order to described more verbosely the type of the model summarized in the file, perhaps added with references to important academic journal articles in which the model appeared first.

Information Section

This section allows you to add more succinct model properties, including a shorter denominator given by Name= qualifier. These shorter attributes will then be attached to the model instance where they help to uniquely identify the model. In contrast to the information contained in the previous section these qualifiers should be short. You can also add yet another item using the Desc= qualifier. So an example of this would be:

%Model Information+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Name = Standard RBC Model, NUM-SS;
Desc = A fairly canonical RBC model with endogenous labour and physical capital;

Parameters Section

As the name suggests, this section provides space for writing down the model’s deep and presumably invariable parameters which are important as they appear in functionals such as the household’s utility or the firm’s production function. Don’t forget to close each declaration with a semi-colon, as this is one of the text parser’s conventions. It is also important, at least for the time being, to adhere to the convention of employing only floats in this section and to refrain from using integers. So this is discouraged:

%Parameters++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
rho        = 1/3;
ss_labour  = 1/3;

And instead one should use:

%Parameters++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
rho        = 1.0/3.0;
ss_labour  = 1.0/3.0;

Variable Vectors Section

This section is very important as it contains a summary of all of the (time-subscripted) variables of the model. The general format of this section for each variable is:

[1] x(t):var_name{endo|con|exo}[log,hp|bk|cf]

The first element is a descriptor of how the time-subscripted variable will appear in the system of nonlinear equations. The second descriptor is a more revealing but still short name, such as capital or consumption. It is preferable to write longer variable names with an underscore, such as for example physical_capital or human_capital. Thirdly, the descriptor in curly brackets allows you to specifically mark of each variable as either, control variable, endogenous state or exogenous state variable, using optimal control theory language. These are inserted in abbreviated style using either con, endo or exo.

Finally, the last option given enclosed in squared brackets allows for two additional options to be specified. Supplying the keyword log means that the approximation of the model should be formed about the log of the variable, while the last option allows to supply a filtering option which is applied to the computation of results based on simulations of the solved model. Currently available choices are either hp for the HP-Filter, bk for the Baxter-King-Filter or cf for the Christiano-Fitzgerald filter. Notice that for exogenous variables you also have to specify the name of the iid shock:

[7] x(t):eps(t):var_name{endo|con|exo}[log,hp|bk|cf]

Boundary Conditions Section

This section is currently not in use but has been included for future compatibility with solution methods which are not based on the perturbation paradigm.

Variable Substitution Non-Linear System

This is perhaps one of the most useful and convenient sections of the model file. In the section right after this one users are asked to insert the DSGE model’s firs-order conditions of optimality which can often be quite tedious and long algebraically. One way of giving users a more convenient and intuitive way of writing down the model’s FOCs is to work with a subsitution system which can be declared in this section. So for example if one wanted to write down the expression for output or the Euler equation for physical capital, one could resort to the following useful replacement definitions:

[1]   @inv(t)   = k(t)-(1-delta)*k(t-1);
[2]   @F(t)     = z(t)*k(t-1)**rho;
[3]   @F_bar    = SS{@F(t)};
[4]   @R(t)     = 1+DIFF{@F(t),k(t-1)}-delta;
[5]   @R(t+1)   = FF_1{@R(t)};
[6]   @U(t)     = @I{eta!=1.0}{c(t)**(1-eta)/(1-eta)}+@I{eta==1.0}{LOG(c(t))};
[7]   @MU(t)    = DIFF{@U(t),c(t)};
[8]   @MU(t+1)  = FF_1{@MU(t)};

These can then be used in the following section instead of having to work with the full expressions instead. Additionally, convience operators are accessible, given by:

DIFF{EXPRESSION,x(t)}     # replaced by first derivate if expression w.r.t. x(t)

SS{EXPRESSION}            # expression is converted to its steady state equivalent

FF_X{EXPRESSION}          # replaced with expression forwarded in time by X periods.
                          # Timing of the information set for expectations
                          # is unchanged!

BB_X{EXPRESSION}          # replaced with expression lagged in time by X periods.
                          # Timing if the information set for expectations
                          # is unchanged!

@ALL{EXPRESSION,[0-1],SS} # short-hand way of declaring many items in one go.
                          # See further below for detailed explanation.

@DISCOUNT                 # Special reserved keyword to define the discount factor.
                          # See further below for detailed explanation.

@I{CONDITION}{EXPRESSION} # An indicator function which inserts the expression whenever
                          # the condition evaluates to TRUE.

When declaring replacement items in this section make sure to adhere to the syntax of always naming them beginning with a @. Also, within this section substitutions within substitutions are permitted. Replacement items for steady-state calculations in the subsequent sections can also be supplied here, but have to be of the form such as:

[1]   @F_bar   = z_bar*k_bar**rho;

In PyMacLab steady state expressions of variables strictly have to adhere to the x_bar naming convention, i.e. be expressed by the stem variable name abbreviation followed by and underscore and the word bar. Finally, the DIFF{EXPRESSION,x(t)} is smart enough to differentiate across different time periods. So as an example with habit persistence in consumption our utility function depends on current and past consumption:

[1]   @DISCOUNT = betta;
[2]   @U(t)     = LOG(c(t)-B*c(t-1));
[3]   @Uc(t)    = DIFF{@U(t),c(t)};

Here the differentiation operator is smart enough to forward the expression by one period before taking the derivative w.r.t to c(t). In fact, internally the above will be replaced with:

[1]   @DISCOUNT = betta;
[2]   @U(t)     = LOG(c(t)-B*c(t-1));
[3]   @Uc(t)    = DIFF{LOG(c(t)-B*c(t-1))+betta*LOG(E(t)|c(t+1)-B*c(t)),c(t)};

This feature only works if the special reserved keyword @DISCOUNT is defined at the top of the list. This tells PyMacLab which discount rate to apply to future (or past) expressions. Finally, as of version 0.95.1 PyMacLab also supports another keyword which works as a short-cut to declare a large number of possible derivatives using only one command. This feature would work as follows:

%Variable Substitution Non-Linear System+++++++++++++
# The utility function and its derivatives
[1]   @MU(t)     = LOG(c(t))+em(t-1)**(1-1/ups)/(1-1/ups);
[2]   @ALL{@MU(t),[0-1],SS};

This command takes all of the partial derivatives (but no cross-partials!) of the supplied function @MU(t) both for the current and the future period, i.e period running from [0-1]. One could also specify this as a list like [0,1]. If the additional optional argument SS is also supplied then the steady state versions of both the original function and the derivatives would be declared. Essentially, the above is just a short-hand for the following manually declared version:

%Variable Substitution Non-Linear System+++++++++++++
# The utility function and its derivatives
[1]   @MU(t)     = LOG(c(t))+em(t-1)**(1-1/ups)/(1-1/ups);
[2]   @MU_bar    = SS{@MU(t)};
[3]   @MUc(t)    = DIFF{@MU(t),c(t)};
[4]   @MUc_bar   = SS{@MUc(t)};
[5]   @MUem(t)   = DIFF{@MU(t),em(t-1)};
[6]   @MUem_bar  = SS{@MUem(t)};
[7]   @MU(t+1)   = FF_1{@MU(t)};
[8]   @MUc(t+1)  = DIFF{@MU(t+1),E(t)|c(t+1)};
[9]   @MUem(t+1) = DIFF{@MU(t+1),em(t)};

Obviously, for reasons of brevity using the @ALL command is a much better option, in particular if the derivatives and steady state expressions one works with are kind of standard and flow naturally from the functional forms of utlity and production functions, for instance.

Note

The whole point of having the subsitutions section present in the library as a functionality to draw on is to reduce systems to a lower dimensionality without having to string together algebraic fragments into enormous mathematical expressions which are hard to read and understand by somebody who has not been involved in designing the model. This approach also reduces the likelihood of introducing mistakes. With the substitution systems everything looks clean and the intuition is immediately discernable from the simplified first-order conditions containing the substitution declarations.

Non-Linear First-Order Conditions Section

In this section users can supply the model’s first order conditions of optimality which are passed to PyMacLab for differentiation and evaluation. So to use the example from the RBC1 example file given above, filling in this section would look as follows:

[1]   @F(t)-@inv(t)-c(t) = 0;
[2]   betta*(@MU(t+1)/@MU(t))*@R(t+1)-1 = 0;
[3]   @F(t)-y(t) = 0;
[4]   LOG(E(t)|z(t+1))-psi*LOG(z(t)) = 0;

where we have made ample use of the convenient substitution definitions declared in the previous section. Expressions, such as the law of motion for the productivity shock, can be supplied in logs for the sake of readability, but otherwise could also alternatively be written as:

[4]   E(t)|z(t+1)/(z(t)**psi) = 0;

Deprecated since version 0.85: In previous versions of PyMacLab it was possible to write down the law of motion of exogenous states without expectations, i.e. z(t)/(z(t-1)**psi) = 0;. This behaviour is now deprecated and no longer supported.

Steady States [Closed Form] Section

For relatively simple models, closed form solutions for the steady state may exist and can be entered here as follows:

betta   = 1.0/R_bar;
k_bar   = ((rho*z_bar)/(R_bar - 1 + delta))**(1.0/(1 - rho));
y_bar   = (z_bar*k_bar)**rho;
c_bar   = y_bar - delta*k_bar;

Note that not only steady-state variables like x_bar can be supplied here, but indeed any variable who’s steady-state value has to be determined endogenously withing the model. Sometimes, depending on the model builder’s assumptions taken, this could also involve the’ determination of a parameter such as betta.

Sometimes the model’s full steady-state can be best determined using a combination of closed form expressions AND the additional numerical solution of a system on nonlinear equations, as is the case in the model file provided as rbc1_res.txt. Notice that here one set of steady state variables are calculated in closed from, given the knowledge of a set of other steady state variables, while these in turn are first solved for in the section using the nonlinear root-finding algorithm. This make sense as for many DSGE models a core set of steady state variables in physical capital and marginal utlity - as an example - can be computed using the non-linear root finder, while all of the other variables’ steady states follow immediately residually from this.

Steady State Non-Linear System [Manual] Section

In this section a partial list of or the entire model’s variables’ steady states can be determined numerically here using good starting values and a Newton-like root-finder algorithm. So this section would something like this:

%Steady State Non-Linear System [Manual]+++++++++++++
[1]   z_bar*k_bar**(rho)-delta*k_bar-c_bar = 0;
[2]   rho*z_bar*k_bar**(rho-1)+(1-delta)-R_bar = 0;
[3]   (betta*R_bar)-1 = 0;
[4]   z_bar*k_bar**(rho)-y_bar = 0;

[1]  c_bar = 1.0;
[2]  k_bar = 1.0;
[3]  y_bar = 1.0;
[4]  betta = 1.0;

Very often, this section is simply a restatement of the first order conditions of optimality but with time subscripts removed and instead replaced with the steady state x_bar notation. This section and the previous can often be the most difficult ones to specify well, as many more complex DSGE models’ steady states are not easy to determine and often require some good judegement, experience and good starting values for the root-finding algorithm.

As of version 0.95.1 of the PyMacLab library users can also use symbolic expressions in the starting values subsection following the non-linear system of equations, such as for instance:

%Steady State Non-Linear System [Manual]+++++++++++++
[1]   z_bar*k_bar**(rho)-delta*k_bar-c_bar = 0;
[2]   rho*z_bar*k_bar**(rho-1)+(1-delta)-R_bar = 0;
[3]   (betta*R_bar)-1 = 0;
[4]   z_bar*k_bar**(rho)-y_bar = 0;

[1]  k_bar = 30.0;
[2]  y_bar = k_bar**alpha;
[3]  c_bar = 2.0;
[4]  betta = 1.0;

Finally, again as of version 0.95.1, users can instead declare in this section the following:

%Steady State Non-Linear System [Manual]+++++++++++++
USE_FOCS=[0,1,2,3];

[1]  k_bar = 30.0;
[2]  y_bar = k_bar**alpha;
[3]  c_bar = 2.0;
[4]  betta = 1.0;

When using this USE_FOCS command, users are instructing the DSGE model instance to automatically form steady state versions of the non-linear system of equations, but doing this only for the equation numbers provided in the passed vector, i.e. [0,1,2,3], which instructs PyMacLab to pick equations 1,2,3,4 out of the system of FOCs declared before this section. Python uses 0-indexed vectors, that is why the list starts with 0 and not 1. If the FOCs are ordered differently, one can also employ different orderings, such as [0,2,3,4]. The point here is to have a way of disregarding certain equations we may not want to include in the non-linear root finding algorithm, such as certain exogenous laws of motion for which we may have calibrated steady state values and do not have to look for them.

Log-Linearized Model Equations Section

In this section you could theoretically also supply the first-order log-linearized equations manually, such as was necessary in Harald Uhlig’s toolbox. But this feature is perhaps best relegated to compatibility tests and proof-of-concept experiments to show that PyMacLab’s computed solutions based on automatic differentiation are identical with the ones computed from this section. An example would be:

# foc consumption
[1]   (1/C_bar)**Theta*X_bar**(Psi*(1-Theta))*x(t)...
     -(1/C_bar)**Theta*X_bar**(Psi*(1-Theta))*c(t)=...
       LAM_bar*lam(t)+A_bar*MU_bar*mu(t);
# foc leisure
[2]   (1-Theta)*c(t)+(Psi*(1-Theta)-1)*x(t)=lam(t)+...
       z(t)+(1-alpha)*k(t-1)-(1-alpha)*l(t);

In this case all variables already have to be interpreted as percentage deviations from steady state. Both in this and in the nonlinear FOCs section, model equations DO NOT necessarily have to be expressed as g(x)=0, but can also be written as f(x)=g(x). In this case the PyMacLab parser simply internally generates f(x)-g(x) = 0 and works with this expression instead.

Variance-Covariance Matrix Section

The standard way of supplying information on the variance-convariance structure of the iid shocks hitting the laws of motions of the exogenous state variables. So this section would look something like this:

Sigma = [sigma_eps**2];

or for more elaborate models like this:

Sigma = [sigma_eps**2   0;
         0    sigma_xi**2];

All sections

If in any of the lines of one of the sections the keyword None is inserted, even in a section which has otherwise been declared in the correct way as described above, then the entire section will be ignored and treated as empty, such as for instance:

%Log-Linearized Model Equations++++++++++++++++++++++++++++
None

If alebraic expression become to long, one can also employ a line-breaking syntax using the elipsis, such as:

[1]   (1-Theta)*c(t)+(Psi*(1-Theta)-1)*x(t)=lam(t)+...
       z(t)+(1-alpha)*k(t-1)-(1-alpha)*l(t);

Finally, as is customary from other programming languages, comments can also be inserted into DSGE model files. However, in contrast to other languages conventions, such as Python itself, at the moment the library will only parse model files correctly if the comments are on a line of their own, and not intermingled with model description items. As usual comments are identified by beginning a new line with the hash symbol #.

Finally, in all sections where it may be applicable, the operators LOG(x) and EXP(x) can be employed, where the former takes the natural logarithm of expression x while the latter raises e to the power x. An example of this would be:

[1]   @U(t)   = LOG(c(t));

More than one way to feed in model properties

As of PyMacLab version 0.95.1, there now exists more than one way to populate a DSGE model instance with information about the properties/features which comprise the model and dictate its ultimate behaviour. These changes have been implemented in order to make PyMacLab’s feature set more compatible with a programming paradigm often called “Meta-programming” or “Template programming” which encapsulates the idea of allowing programs to change their own “source code” or otherwise usually assumed fixed features during runtime.

Or at a more basic level, it simply offers a comfortable way for users of the library to change essential features of DSGE models or alternatively swap features between them while a program is running. This makes PyMacLab far more powerful in principle than for instance Dynare. So besides reading in a conformable DSGE model file from your computer’s file system, which other ways are on offer to populate a DSGE model instance?

Instead of passing the model file’s full path as a string to the DSGE model at instantiation time, we could have also alternatively passed the actual model file itself as a big triple-quoted string to the DSGE class generating instances. This could be defined inside a Python (batch) script and could for instance be done like this:

modstr='''
%Model Description+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
This is just a standard RBC model, as you can see.


%Model Information+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Name = Standard RBC Model, RES-SS;


%Parameters++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
rho       = 0.36;
delta     = 0.025;
R_bar     = 1.01;
betta     = 1.0/R_bar;
eta       = 2.0;
psi       = 0.95;
z_bar     = 1.0;
sigma_eps = 0.052;


%Variable Vectors+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[1]  k(t):capital{endo}[log,bk]
[2]  c(t):consumption{con}[log,bk]
[4]  y(t):output{con}[log,bk]
[5]  z(t):eps(t):productivity{exo}[log,bk]
[6]  @inv(t):investment[log,bk]
[7]  @R(t):rrate

%Boundary Conditions++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
None


%Variable Substitution Non-Linear System++++++++++++++++++++++++++++++++++++++
# Special discount variable
[1]   @DISCOUNT = betta;
[1]   @inv(t)   = k(t)-(1-delta)*k(t-1);
[2]   @inv_bar  = SS{@inv(t)};
[2]   @F(t)     = z(t)*k(t-1)**rho;
[2]   @Fk(t)    = DIFF{@F(t),k(t-1)};
[2]   @Fk_bar   = SS{@Fk(t)};
[2]   @F_bar    = SS{@F(t)};
[3]   @R(t)     = 1+DIFF{@F(t),k(t-1)}-delta;
[4]   @R_bar    = SS{@R(t)};
[3]   @R(t+1)   = FF_1{@R(t)};
[4]   @U(t)     = c(t)**(1-eta)/(1-eta);
[5]   @MU(t)    = DIFF{@U(t),c(t)};
[5]   @MU_bar   = SS{@U(t)};
[6]   @MU(t+1)  = FF_1{@MU(t)};



%Non-Linear First-Order Conditions+++++++++++++++++++++++++++++++++++++++++++++
# Insert here the non-linear FOCs in format g(x)=0

[1]   @F(t)-@inv(t)-c(t) = 0;
[2]   betta*(@MU(t+1)/@MU(t))*@R(t+1)-1 = 0;
[3]   @F(t)-y(t) = 0;
[4]   LOG(E(t)|z(t+1))-psi*LOG(z(t)) = 0;


%Steady States [Closed Form]+++++++++++++++++++++++++++++++++++++++++++++++++++
[1]   y_bar = @F_bar;


%Steady State Non-Linear System [Manual]+++++++++++++++++++++++++++++++++++++++
[1]   @F_bar-@inv_bar-c_bar = 0;
[2]   betta*@R_bar-1 = 0;
[3]   betta*R_bar-1 = 0;

[1]   c_bar = 1.0;
[2]   k_bar = 1.0;
[3]   betta = 0.9;

%Log-Linearized Model Equations++++++++++++++++++++++++++++++++++++++++++++++++
None


%Variance-Covariance Matrix++++++++++++++++++++++++++++++++++++++++++++++++++++
Sigma = [sigma_eps**2];


%End Of Model File+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
'''

import pymaclab as pm

rbc = pm.newMOD(modstr,mesg=True,ncpus='auto')

As you can see, the declared Python object modstr is just a string which holds a standard PyMacLab model file in its entirety (with line breaks!) This is then passed to the DSGE class to instantiate a new model and internally PyMacLab recognizes this not as a full path pointer to a physical model file existing in your computer’s file system but instead as the contents of the file itself ready for direct processing.

Finally, yet one more way open to users to instantiate and populate a new DSGE model with its characteristic features is closely related to the one described immediately above. This second way uses a Python templating library called wheezy.template which allows conformable PyMacLab model files to be generated on the fly from within a running Python script using a standard Python dictionary of DSGE model properties. Such a dictionary is always created by default and then attached to each DSGE model instance whenever they are created and is held inside the object model.template_paramdic. For a simple RBC model this would look like:

# Load the library and the models branch
In [1]: import pymaclab as pm
In [2]: import pymaclab.modfiles.models as models
In [3]: import pymaclab.modfiles.templates.wheezy_template as template

# Now instantiate the model
In [3]: rbc = pm.newMOD(models.stable.rbc1_num, mesg=True)

# Check contents of the template dictionary
In [4]: rbc.template_paramdic.keys()

# These are only the keys of the dictionary, but check the contents yourself
# to see that they are all standard Python data structures describing the model
Out[1]:
['use_focs',
 'vardic',
 'sigma',
 'mod_desc',
 'subs_list',
 'focs_list',
 'manss_sys',
 'mod_name',
 'llsys_list',
 'paramdic',
 'ssidic',
 'ssys_list']

 # Now use the template to automatically generate
 # a conformable PyMacLab model file string
 In [5]: modstr = template.render(rbc.template_paramdic)

 # Now print the modstr and check what it looks like
 In [6]: print modstr

 Out[2]:
 '''
 %Model Description++++++++++++++++++++++++++++++++++
 None


 %Model Information++++++++++++++++++++++++++++++++++
 # Short model name
 Name = Standard RBC Model, NUM-SS;
 # Short model description



 %Parameters+++++++++++++++++++++++++++++++++++++++++
 [1]   z_bar = 1.0;
 [2]   psi = 0.95;
 [3]   sigma_eps = 0.052;
 [4]   betta = 0.990099009901;
 [5]   eta = 2.0;
 [6]   rho = 0.36;
 [7]   delta = 0.025;



 %Variable Vectors++++++++++++++++++++++++++++++++++++
 [1]   k(t):capital{endo}[log,bk]
 [1]   c(t):consumption{con}[log,bk]
 [2]   y(t):output{con}[log,bk]
 [1]   z(t):eps(t):productivity{exo}[log,bk]
 [1]   @inv(t):investment [log,bk]
 [2]   @R(t):rrate


 %Boundary Conditions+++++++++++++++++++++++++++++++++
 None


 %Variable Substitution Non-Linear System+++++++++++++
 [1]   @inv(t) = k(t)-(1-delta)*k(t-1);
 [2]   @inv_bar = SS{@inv(t)};
 [3]   @F(t) = z(t)*k(t-1)**rho;
 [4]   @Fk(t) = DIFF{@F(t),k(t-1)};
 [5]   @Fk_bar = SS{@Fk(t)};
 [6]   @F_bar = SS{@F(t)};
 [7]   @R(t) = 1+DIFF{@F(t),k(t-1)}-delta;
 [8]   @R_bar = SS{@R(t)};
 [9]   @R(t+1) = FF_1{@R(t)};
 [10]   @U(t) = c(t)**(1-eta)/(1-eta);
 [11]   @MU(t) = DIFF{@U(t),c(t)};
 [12]   @MU_bar = SS{@U(t)};
 [13]   @MU(t+1) = FF_1{@MU(t)};


 %Non-Linear First-Order Conditions+++++++++++++++++++
 [1]   @F(t)-@inv(t)-c(t) = 0;
 [2]   betta*(@MU(t+1)/@MU(t))*@R(t+1)-1 = 0;
 [3]   @F(t)-y(t) = 0;
 [4]   LOG(E(t)|z(t+1))-psi*LOG(z(t)) = 0;


 %Steady States [Closed Form]++++++++++++++++++++++++++
 None


 %Steady State Non-Linear System [Manual]+++++++++++++
 [1]   @F_bar-@inv_bar-c_bar = 0;
 [2]   y_bar-@F_bar = 0;
 [3]   betta*@R_bar-1 = 0;
 [4]   betta*R_bar-1 = 0;

 [1]   c_bar = 1.0;
 [2]   k_bar = 1.0;
 [3]   y_bar = 1.0;
 [4]   R_bar = 1.01;



 %Log-Linearized Model Equations++++++++++++++++++++++
 None


 %Variance-Covariance Matrix++++++++++++++++++++++++++
 Sigma = [ 0.002704 ];


 %End Of Model File+++++++++++++++++++++++++++++++++++
 '''

 # You could now check if the model also loads with the generated modfile string
 In [7]: rbc_alt = pm.newMOD(modstr, mesg=True)

As you can see, with the power of templating engines [#f1]_ such as wheezy.template we can generate PyMacLab-conformable DSGE model files on-the-fly by passing simple Python data structures to the template and calling its render() method. In the above script, the DSGE models called rbc and rbc_alt will be identical save for small numerical discrepancies introduced because of floating-point arithmatics imprecision. In the near future PyMacLab will include another template which will allow the automatic generation of Dynare-conformable model files, allowing users to compare and contrast results computed in both environments.

Footnotes

[1]There exist far more popular templating engines than wheezy.template. One of such, perhaps the most popular, is an engine library called jinja2 which is often used by programmers to design dynamic webpages. The other candidate is a library called cheetah. In spite of wheezy.template‘s lesser popularity, it was chosen for PyMacLab because it claims to be the fastest template engine of all of the above mentioned candidates.

Tutorial 3 - The PyMacLab DSGE instance

Introduction

As already stated in the introduction of the introductory basic tutorial, PyMacLab’s strength or orginal design goal has been that of providing users with a rich and flexible DSGE data structure (also called Class in object-oriented programming speak) which allows them to do lots of interesting things with DSGE models and to treat them as if they were some kind of primitive data type in their own right. While the previous tutorial described some basics as well as the all-important DSGE model file structure and syntax conventions, in this section I am going to stress some of the object-oriented programming features of PyMacLab, in particular the structure of a PyMacLab DSGE model instance or data structure.

Readers with a background in modern programming languages supporting the object-oriented programming (OOP) paradigm will easily relate to concepts in this sections, while for others it may appear more cryptic at first sight. But to convey these concepts to researchers is important, as it stresses many particular advantages of PyMacLab over other programs, and in particular its flexibility, transparency, consistency, persistence and enormous scope for extensibility. All example code fragments provided here assume that you are replicating them within an IPyton interactive session, but they could also be called from a Python program “batch” file.

Understanding the PyMacLab DSGE model class and its instances

PyMacLab has been written in the Python programming language which supports object-oriented programming. This means, that more than 80% of PyMacLab’s code is devoted to the definition of data fields and methods of the DSGE_model Class, which forms the basis for all DSGE models users can load, spawn or instantiate and interact with once they have imported they PyMacLab library into their programs. As already explained elsewhere, the basis of all DSGE model instances is the DSGE model’s model text file in which it is defined in terms of its specific characteristics, such as its parameters and first-order conditions of optimality. We recall that this process of loading or instantiating a DSGE model worked as follows:

# Import the pymaclab module into its namespace, also import os module
In [1]: import pymaclab as pm
In [2]: from pymaclab.modfiles import models

# Instantiate a new DSGE model instance like so
In [4]: rbc1 = pm.newMOD(models.stable.rbc1)

After executing these lines of code in an interactive environment such as that provided by IPython, which emulates well the feel and behaviour of the Matlab interactive environment, the DSGE instance or data object going by the name of rbc1 now exists in the namespace of the running program and can be put to further use. But before looking at these various ways possible to make effective use of this DSGE model instance, let’s first trace the various steps the programs goes through when a DSGE model get instantiated. So what happens internally when the above last line in the code fragment is called:

  1. The empty shell DSGE model instance gets instantiated
  2. The DSGE model instance reads the model file provided to it and any other arguments and saves them by attaching them to itself.
  3. Instantiation Step 1: The files get read in and a method defined on the instance simply splits the file into its individual sections and saves these raw sections by attaching them to itself.
  4. Instantiation Step 2: A parser method is called which disaggregates the raw information provided in each section of the model file and begins to extract meaningful information from it, each time saving this information by attaching it to itself as data fields. Also, the DSGE model instance is prepared and set up in order to attempt to solve for the steady state of it manually at the command line, instead of doing it automatically. If you want the model instance to do ONLY this next step and stop there for you to explore further interactively, you must call the command with and extra argument like this:
# Import the pymaclab module into its namespace, also import os module
In [1]: import pymaclab as pm
In [2]: from pymaclab.modfiles import models

# Instantiate a new DSGE model instance like so, but adding initlev=0 as extra argument
In [3]: rbc1 = pm.newMOD(models.stable.rbc1,initlev=0)
  1. Instantiation Step 3: The information is used in order to attempt to compute the numerical steady-state of the model. If you want the model instance to do ONLY this next step and stop there for you to explore further interactively, you must call the command with and extra argument like this:
# Import the pymaclab module into its namespace, also import os module
In [1]: import pymaclab as pm
In [2]: from pymaclab.modfiles import models

# Instantiate a new DSGE model instance like so, but adding initlev=1 as extra argument
In [3]: rbc1 = pm.newMOD(models.stable.rbc1,initlev=1)
  1. Instantiation Step 4: If the steady state was computed successfully then the model’s analytical and numerical Jacobian and Hessian are computed. Finally a preferred dynamic solution method is called which solves the model for its policy function and other laws of motion.

To give users a choice of “solution depths” at DSGE object instantiation time is important and useful, especially in the initial experimentation phase during which the DSGE model file gets populated. That way researchers can first carefully solve one part of the problem (i.e. looking for the steady state) and indeed choose to do so manually on the IPython interactive command shell, allowing them to immediately inspect any errors.

Instantiation options for DSGE model instances

There are a couple of instance invocation or instantiation arguments one should be aware of. At the time of writing these lines there are in total 5 other arguments (besides the DSGE model template file path) which can be passed to the pymaclab.newMOD function out of which 1 is currently not (yet) supported and not advisable to employ. The other 4 options determine the initiation level of the DSGE model (i.e. how far it should be solved if at all), whether diagnosis messages should be printed to screen during instantiation, how many CPU cores to employ when building the Jacobian and Hessian of the model, and finally whether the expensive-to-compute Hessian should be computed at all. Remember that the last option is useful as many researchers often - at least initially - want to explore the solution to their model to a first order of approximation before taking things further. So here are the options again in summary with their default values:

Option with default value Description
pm.newMOD(mpath,initlev=2) Initlev=0 only parses and prepares for manual steady state calculation
Initlev=1 does initlev=0 and attempts to solve for the model’s steady state automatically
Initlev=2 does initlev=0, initlev=1 and generates Jacobian and Hessian and solves model dynamically
pm.newMOD(mpath,mesg=False) Prints very useful runtime instantiation messages to the screen for users to follow progress
pm.newMOD(mpath,ncpus=1) CPU cores to be used in expensive computation of model’s derivatives, ‘auto’ for auto-detection
pm.newMOD(mpath,mk_hessian=True) Should Hessian be computed at all, as is expensive?
pm.newMOD(mpath,use_focs=False) Should only the model’s FOCs be used to computed the steady state? Accepts Python list or tuple
pm.newMOD(mpath,ssidic=None) Use in conjunction with previous argument to specify initial starting values as Python dictionary
pm.newMOD(mpath,sstate=None) Specify steady state values as Python dictionary and supply here. No SS computation in instance

Needless to say, all of the options can be and usually are called in combination, they are only shown separately here for sake of expositional clarity. Medium-sized to large-sized models can take considerable time to compute the Jacobian alone, let alone the Hessian. On the other hand passing more (real as opposed to virtual) CPU cores to the instantiation process can significantly cut down computation time. In this case, the FOCs nonlinear equations are distributed to individual cores for analytical differentiation as opposed to doing this serially on one CPU core.

Working with DSGE model instances

The most useful feature is to call the model with the option initlev=0, because this will allow you more control over the steady-state computation of the model by permitting a closer interactive inspection of the DSGE model instance as created thus far. Let’s demonstrate this here:

# Import the pymaclab module into its namespace, also import os module
In [1]: import pymaclab as pm
In [2]: from pymaclab.modfiles import models

# Instantiate a new DSGE model instance like so, but adding initlev=0 as extra argument
In [3]: rbc1 = pm.newMOD(models.stable.rbc1,initlev=0)

# This datafield contains the original nonlinear system expressed as g(x)=0
In [4]: rbc1.sssolvers.fsolve.ssm
['z_bar*k_bar**(rho)-delta*k_bar-c_bar',
 'rho*z_bar*k_bar**(rho-1)+(1-delta)-R_bar',
 '(betta*R_bar)-1',
 'z_bar*k_bar**(rho)-y_bar']

# This datafield contains the initial values supplied to the rootfinder algorithm
In [5]: rbc1.sssolvers.fsolve.ssi
{'betta': 1.0, 'c_bar': 1.0, 'k_bar': 1.0, 'y_bar': 1.0}

# Instead of letting the model during instantiation solve the model all the way through,
# we can solve for the steady state by hand, manually
In [6]: rbc1.sssolvers.fsolve.solve()

# And then inspect the solution and some message returned by the rootfinder
In [6]: rbc1.sssolvers.fsolve.fsout
{'betta': 0.9900990099009901,
 'c_bar': 2.7560505909330626,
 'k_bar': 38.1607004898424,
 'y_bar': 3.7100681031791227}

In [7]: rbc1.sssolvers.fsolve.mesg
'The solution has converged.'

Another useful lesson to take away from this example is that a DSGE model instance is like a many-branch tree structure, just like the Windows File Explorer so many people are familiar with, where individual “nodes” represent either data fields or methods (function calls) which equip the model instance with some functionality. This kind of approach of structuring and programming a solution to the problem of designing a program which handles the solution-finding of DSGE models offers enormous scope for experimentation and extensibility. After a DSGE model has been instantiated without passing the initlev argument, you can inspect this structure like so:

# Import the pymaclab module into its namespace, also import os module
In [1]: import pymaclab as pm
In [2]: from pymaclab.modfiles import models

# Instantiate a new DSGE model instance like so
In [3]: rbc1 = pm.newMOD(models.stable.rbc1)

# Inspect the data fields and methods of the DSGE model instance
In [4]: dir(rbc1)
['__class__',
 '__delattr__',
 '__dict__',
 '__doc__',
 '__format__',
 '__getattribute__',
 '__hash__',
 '__init__',
 '__module__',
 '__new__',
 '__reduce__',
 '__reduce_ex__',
 '__repr__',
 '__setattr__',
 '__sizeof__',
 '__str__',
 '__subclasshook__',
 '__weakref__',
 '_initlev',
 'audic',
 'author',
 'ccv',
 'dbase',
 'deltex',
 'getdata',
 'info',
 'init2',
 'manss_sys',
 'mkeigv',
 'mkjahe',
 'mkjahen',
 'mkjahenmat',
 'mkjahepp',
 'mkjaheppn',
 'mod_name',
 'modfile',
 'nall',
 'ncon',
 'nendo',
 'nexo',
 'nlsubs',
 'nlsubs_list',
 'nlsubs_raw1',
 'nlsubs_raw2',
 'nother',
 'nstat',
 'numssdic',
 'paramdic',
 'pdf',
 'setauthor',
 'ssidic',
 'sssolvers',
 'sstate',
 'ssys_list',
 'subs_vars',
 'switches',
 'texed',
 'txted',
 'txtpars',
 'updf',
 'updm',
 'vardic',
 'vreg']

As you can see, the attributes exposed at the root of the instance are plenty and can be acccessed in the usual way:

# Import the pymaclab module into its namespace, also import os module
In [1]: import pymaclab as pm
In [2]: from pymaclab.modfiles import models

# Instantiate a new DSGE model instance like so
In [3]: rbc1 = pm.newMOD(models.stable.rbc1)

# Access one of the model's fields
In [4]: rbc1.ssys_list
['z_bar*k_bar**(rho)-delta*k_bar-c_bar',
'rho*z_bar*k_bar**(rho-1)+(1-delta)-R_bar',
'(betta*R_bar)-1',
'z_bar*k_bar**(rho)-y_bar']

So one can observe that the data field rbc1.ssys_list simply summarizes the system of nonlinear equations which has been described in the relevant section of the DSGE model file. Now you know how to explore the DSGE model instance and understand its general structure, and we conclude this short tutorial by inviting you to do so. Don’t forget that some nodes at the root possess further sub-nodes, as was the case when cascading down the rbc1.sssolvers branch. To help your search, the only other node with many more sub-nodes is the rbc1.modsolvers branch, which we will explore more in the next section to this tutorial series.

DSGE modelling made intuitive

Before concluding this tutorial, we will demonstrate how PyMacLab’s DSGE data structure (or instance) approach allows researchers to implement ideas very intuitively, such as for instance “looping” over a DSGE model instance in order to explore how incremental changes to the parameter space alter the steady state of the model. Leaving our usual interactive IPyton shell, consider the following Python program file:

# Import the pymaclab module into its namespace
# Also import Numpy for array handling and Matplotlib for plotting
import pymaclab as pm
from pymaclab.modfiles import models
import numpy as np
from matplotlib import pyplot as plt

# Instantiate a new DSGE model instance like so
rbc1 = pm.newMOD(models.stable.rbc1)

# Create an array representing a finely-spaced range of possible impatience values
# Then convert to corresponding steady state gross real interest rate values
betarr = np.arange(0.8,0.99,0.001)
betarr = 1.0/betarr

# Loop over the RBC DSGE model, each time re-computing for new R_bar
ss_capital = []
for betar in betarr:
    rbc1.updaters.paramdic['R_bar'] = betar # assign new R_bar to model
    rbc1.sssolvers.fsolve.solve() # re-compute steady stae
    ss_capital.append(rbc1.sssolvers.fsolve.fsout['k_bar']) # fetch and store k_bar

# Create a nice figure
fig1 = plt.figure()
plt.grid()
plt.title('Plot of steady state physical capital against R\_bar')
plt.xlabel(r'Steady state gross real interest rate')
plt.ylabel(r'Steady State of physical capital')
plt.plot(betarr,ss_capital,'k-')
plt.show()

Anybody who has done some DSGE modelling in the past will easily be able to intuitively grasp the purpose of the above code snippet. All we are doing here is to loop over the same RBC model, each time feeding it with a slightly different steady state groos real interest rate value and re-computing the steady state of the model. This gives rise to the following nice plot exhibting the steady state relationship between the interest rate and the level of physical capital prevailing in steady state:

(Source code, png, hires.png, pdf)

../_images/test2.png

That was nice and simple, wasn’t it? So with the power and flexibility of PyMacLab DSGE model instances we can relatively painlessly explore simple questions such as how differing deep parameter specifications for the impatience factor \beta can affect the steady state level of physical capital. And indeed, as intuition would suggest, less patient consumers are less thrifty and more spend-thrifty thus causing a lower steady state level of physical capital in the economy. This last example also serves to make another important point. PyMacLab is not a program such as Dynare, but instead an add-in library for Python prividing an advanced DSGE model data structure in form of a DSGE model class which can be used in conjunction with any other library available in Python.

Tutorial 4 - The Python DSGE instance updater methods

Introduction

In the previous tutorial we described the general structure and some of the behaviour of PyMacLab’s DSGE model’s object-oriented design and what kind of advantages model builders can derive from this. In particular, at the end of the last tutorial, we saw how the decision to design DSGE model’s instances essentially in from of an object-oriented advanced data structure with well-defined data fields and methods allowed us to “loop over” a DSGE model instance many times over, each time feeding the instance with slightly a different value for the impatience factor \beta and obtaining a new value for the corresponding steady state value of physical capital.

To be able to do something like this simply and intuitively is one of the many advantages that PyMacLab has over traditional DSGE model solving programs such as Dynare. While Dynare is a complete program allowing users to solve for specific DSGE models with little added post-solution programmatic flexibility, PyMacLab is from ground up designed to be a library providing an object-oriented Class-defined advanced data structure with data fields and behaviour in form of methods. This makes it incredibly easy to build programs which treat DSGE models as if they were simple data structures, such as integers, floats, lists or any other data structure you are familiar with from using a number of different programming languages. This makes PyMacLab incredibly flexible for within- and post-solution program interventions, where in this section we will define more clearly what we meant by program interventions. The main point to take away from all of this is that Dynare is a canned program to get specific solutions, while PyMacLab is a Python plugin library, whose only purpose is to make available an advanced DSGE model class suitable for carrying out a variety of tasks.

While in programs like Dynare the most important aspect from the point of view of model builders to obtaining the solution of DSGE models lies in the specification of so-called model files, the proper way to understand the more flexible operation of PyMacLab is to view the model template only as a specific point of departure from which to start your analysis from. The only point of making use of model template files is to initialize or instantiate a DSGE model instance, the real power in using PyMacLab lies in the methods made available to researchers which become available after model files have been read in and the DSGE model instances become available inside the Python interactive environment. It is this post-instantiation scope for activity and alteration which makes using PyMacLab so much fun for researchers. As already illustrated in the previous tutorial, this design decision opens up the possibility of easily and intuitively obtaining a large permutation of possible solution outcomes quickly. There are two convenience methods or avenues open to the researcher to intelligently update DSGE model instances dynamically at runtime which we will describe in some detail next.

One-off alteration of one specific model property

Once a DSGE model is instantiated using a model template files as a point of departure and a specific source of information from which model properties get parsed and attached as data fields to the DSGE model instance, all we have done is to initialize or read into memory a specific state of our DSGE model. You may also recall that the nature of this state differs crucially depending upon what kind of initlev parameter value was passed at DSGE model instantiation. So we recall that at initlev=0 only information gets parsed in but nothing at all gets solved, not even the steady state. Other values for initlev gave rise to instantiated model instances which were solved all down to ever deeper levels including possibly including the steady state and dynamic elasticities (policy functions). In this section I will briefly describe methods allowing a functionality akin to comparative statics analysis in PyMacLab in which models can get easily be re-initialized with different model information injected into existing DSGE model instances dynamically at runtime.

One way of using PyMacLab to carry out a kind of comparative statics analysis has already been touched upon for one specific case in the previous tutorial, using the example of the rbc1.updaters.paramdic attribute. In general you will discover that navigating into the rbc1.updaters. branch of the model exposes a number of data fields which in exact name are also available at the root rbc1. of the model. So how then does rbc1.updaters.paramdic differ from its counterpart available at the root rbc1.paramdic?

The answer to this question is that while rbc1.paramdic is just the parameter dictionary which get populated at model instantiation time, rbc1.updaters.paramdic is a “wrapped” version of the same dictionary with the added behaviour that if a new values gets inserted into it using the relevant dictionary methods, such as rbc1.updaters.paramdic['R_bar']=1.02 or rbc1.updates.paramdic.update(dic_alt) where dic_alt is some other dictionary containing updated values of one or more keyed item in the original paramdic, then an internal function gets triggered at assignment time which re-initializes the model using the updated values for paramdic.

If the newly assigned value is exactly identical to the value which was already stored in the original paramdic, then the model will not get updated as its state has remained unaltered. The following smart one-off updaters are available all of which possess the above described behaviour. Notice also that the DSGE model will only get updated down to the level first specified in initlev at model instantiation time. Also, when any of these wrapped updater objects just gets called in the normal way they behave exactly as their non-wrapped counterparts by returning the values stored in them without any additional behaviour updating the model.

Wrapped updater object Description
model.updaters.paramdic Inserting updated values re-initializes the model with new set of parameter values
model.updaters.vardic Inserting new values updates/changes values and attributes defined for variables used in the model
model.updaters.nlsubsdic Inserting new values updates the dictionary of variable substitution items beginning with @
model.updaters.foceqs Inserting new string equations updates the set of non-linear first-order conditions of optimality
model.updaters.manss_sys In this list users can insert new equations as strings into the closed form steady state system
model.updaters.ssys_list In this list users can insert new equations as strings into the numerical steady state system
model.updaters.sigma A wrapped matrix which updates the model when changed values are inserted in the varcovar-matrix

Notice how in this branch rbc1.updaters. changing wrapped objects by assigning new values will trigger automatic updating of the DSGE model immediately upon assignment. This behaviour may not always be desirable whenever a series of changes need to be made before updating of the model can be considered. Whenever such situations occur an alternative route needs to be taken which we will explore next.

Altering many model properties and queued processing

At times researchers may want to load or instantiate a particular DSGE model instance using a corresponding template file but then perhaps plan to radically modify the model dynamically at runtime, by combining such actions as introducing new time-subscripted variables, altering the deep parameter space and adding new or augmenting existing equations in the system of non-linear FOCs. Whenever such radical alterations are considered, they will often have to happen in combindation before the model gets updated using the new information passed to it. In this case users will use the same wrapped objects already described above but instead use them in the rbc1.updaters_queued. branch.

Here, first a number of changes can be made to objects such as rbc1.updaters_queued.paramdic or rbc1.updaters_queued.foceqs, etc. which by themselves will not trigger an automatic model updating functionality. Instead all changes will be put into a queue which will then have to be processed manually by calling the method rbc1.updaters_queued.process_queue() after all desired changes have been made. This addscenormous flexibility to model builders’ options, as they can essentially build a completely new model at runtime dynamically starting from a simple model instantiated at the outset of their Python scripts/batch files. Therefore, this functionality allows users to dynamically update all information at runtime which was first parsed from the model template file, each time re-computing the DSGE model’s new state given the changes made after the call to the queue processing method has been made.

Wrapped updater object Description
model.updaters_queued.paramdic Inserting updated values re-initializes the model with new set of parameter values
model.updaters_queued.vardic Inserting new values updates/changes values and attributes defined for variables used in the model
model.updaters_queued.nlsubsdic Inserting new values updates the dictionary of variable substitution items beginning with @
model.updaters_queued.foceqs Inserting new string equations updates the set of non-linear first-order conditions of optimality
model.updaters_queued.manss_sys In this list users can insert new equations as strings into the closed form steady state system
model.updaters_queued.ssys_list In this list users can insert new equations as strings into the numerical steady state system
model.updaters_queued.sigma A wrapped matrix which updates the model when changed values are inserted in the varcovar-matrix
model.updaters_queued.queue The actual queue. Here objects which have been altered will be stored as strings
model.updaters_queued.process_queue The queue processing method which finally updates the queued objects in the right order

Tutorial 5 - Steady State Solution Methods

Introduction

In the previous tutorial we learnt how a PyMacLab DSGE model instance possesses the capability to intelligently upate its properties following the re-declaration at runtime of attached data fields such as the parameter space or the set of non-linear first-order conditions of optimality. In this section we will learn an important component to PyMacLab’s DSGE models which provides users with a large number of options available for solving models’ steady state solution. The great number of possible avenues to take here is quite deliberate; it would be reasonable to argue that for medium- to large-sized models the most difficult part to finding the general dynamic solution based on the approximation method of perturbations is to first obtain the steady state solution around which the approximations are computed. In total we are going to explore 7 different variants suitable for seeking to compute the steady state. So let’s get started.

Option 1: Using the model’s declared FOCs and passing arguments at model instantiation

Choosing option 1 allows users to leave the numerical as well as closed form steady state sections in the model template files entirely empty or unused indicated by the “None” keyword inserted into any line in these sections. In this case, the library has to rely on the time- subscripted non-linear first-order conditions of optimality, convert them to steady state equivalents and somehow discover the required set of initial guesses for the variables’ steady states to be searched for using the non-linear root-finding algorithm. This is accomplised in the following way. Consider first the following simple example of a DSGE model file:

%Model Description+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
This is just a standard RBC model, as you can see.


%Model Information+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Name = Standard RBC Model;


%Parameters++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
rho       = 0.36;
delta     = 0.025;
R_bar     = 1.01;
betta     = 1/R_bar;
eta       = 2.0;
psi       = 0.95;
z_bar     = 1.0;
sigma_eps = 0.052;


%Variable Vectors+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[1]  k(t):capital{endo}[log,bk]
[2]  c(t):consumption{con}[log,bk]
[4]  y(t):output{con}[log,bk]
[5]  z(t):eps(t):productivity{exo}[log,bk]
[6]  @inv(t):investment[log,bk]
[7]  @R(t):rrate

%Boundary Conditions++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
None


%Variable Substitution Non-Linear System++++++++++++++++++++++++++++++++++++++++++++++++
[1]   @inv(t)   = k(t)-(1-delta)*k(t-1);
[2]   @inv_bar  = SS{@inv(t)};
[3]   @F(t)     = z(t)*k(t-1)**rho;
[4]   @F_bar    = SS{@F(t)};
[5]   @R(t)     = 1+DIFF{@F(t),k(t-1)}-delta;
[6]   @R_bar    = SS{@R(t)};
[7]   @R(t+1)   = FF_1{@R(t)};
[8]   @U(t)     = c(t)**(1-eta)/(1-eta);
[9]   @MU(t)    = DIFF{@U(t),c(t)};
[10]  @MU(t+1)  = FF_1{@MU(t)};



%Non-Linear First-Order Conditions++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Insert here the non-linear FOCs in format g(x)=0

[1]   @F(t)-@inv(t)-c(t) = 0;
[2]   betta*(@MU(t+1)/@MU(t))*@R(t+1)-1 = 0;
[3]   @F(t)-y(t) = 0;
[4]   LOG(E(t)|z(t+1))-psi*LOG(z(t)) = 0;


%Steady States [Closed Form]+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
None


%Steady State Non-Linear System [Manual]+++++++++++++++++++++++++++++++++++++++++++++++++
None


%Log-Linearized Model Equations++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
None


%Variance-Covariance Matrix++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Sigma = [sigma_eps**2];


%End Of Model File+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Notice how we have left the usual sections employed to supply information useful for finding the steady state unused indicated by inserting the keyword “None”. As you can see by inspecting the system of non-linear first order conditions, a steady state could be obtained by passing a steady state version of the FOCs to the non-linear root-finding algorithm, with the additional qualifier that in this particular case we would ideally like to omit passing the last line which is just a declaration of the own-lagged law of motion of the exogenous state [#f1]_ productivity shock. This would lead to a 3 equation system in c_bar, k_bar and y_bar. Further more, we would have to let the model somehow know the set of intial guesses for these three variables, which we often tend to set to some generic values, such as 1.0 for all three of them. How is all of this accomplished? By passing the relevant variables directly to the DSGE model at instantiation time like so:

# Import the pymaclab module into its namespace, also import os module
In [1]: import pymaclab as pm
In [2]: from pymaclab.modfiles import models

# Define the ssidic of initial guesses or starting values
In [3]: ssidic = {}
In [4]: ssidic['c_bar'] = 1.0
In [5]: ssidic['k_bar'] = 1.0
In [6]: ssidic['y_bar'] = 1.0

# Instantiate a new DSGE model instance like so
In [7]: rbc1 = pm.newMOD(models.stable.rbc1_ext,use_focs=[0,1,2],ssidic=ssidic)

The default value passed to the DSGE model instance’s argument “use_focs” is False, the alternative value is a zero-indexed Python list (or tuple) indicating the equations of the declared system of FOCs to use in finding the steady state numerically. In the case of the model file given here, we don’t want to use the last line of 4 equations and thus set the list equal to [0,1,2]. We also define a dictionary of initial starting values or guesses for the three steady state values we wish to search for and pass this as a value to the argument ssidic. This method has the added advantage that steady state initial starting values can be determined intelligently at runtime external to the model file.

Option 2: Using the model’s declared FOCs and passing arguments inside the model file

Choosing option 2 is essentially mimicking the same method used in option 1, with the only difference that everything happens inside the model file itself and nothing has to be passed using arguments to the DSGE model instance at instantiation time externally. Instead, the list of FOC equations to be used in the calculation of the steady states is passed inside the numerical steady states section, as shown in model files rbc1_focs.txt, as follows:

%Model Description+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
This is just a standard RBC model, as you can see.


%Model Information+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Name = Standard RBC Model, USE_FOCS;


%Parameters++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
rho       = 0.36;
delta     = 0.025;
betta     = 1.0/1.01;
eta   = 2.0;
psi   = 0.95;
z_bar     = 1.0;
sigma_eps = 0.052;


%Variable Vectors+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[1]  k(t):capital{endo}[log,bk]
[2]  c(t):consumption{con}[log,bk]
[3]  y(t):output{con}[log,bk]
[4]  R(t):rrate{con}[log,bk]
[5]  z(t):eps(t):productivity{exo}[log,bk]
[6]  @inv(t):investment[log,bk]

%Boundary Conditions++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
None


%Variable Substitution Non-Linear System++++++++++++++++++++++++++++++++++++++++++++++++
[1]   @inv(t)   = k(t)-(1-delta)*k(t-1);
[2]   @inv_bar  = SS{@inv(t)};
[2]   @F(t)     = z(t)*k(t-1)**rho;
[2]   @Fk(t)    = DIFF{@F(t),k(t-1)};
[2]   @Fk_bar   = SS{@Fk(t)};
[2]   @F_bar    = SS{@F(t)};
[3]   @R(t)     = 1+DIFF{@F(t),k(t-1)}-delta;
[4]   @R_bar    = SS{@R(t)};
[3]   @R(t+1)   = FF_1{@R(t)};
[4]   @U(t)     = c(t)**(1-eta)/(1-eta);
[5]   @MU(t)    = DIFF{@U(t),c(t)};
[5]   @MU_bar   = SS{@U(t)};
[6]   @MU(t+1)  = FF_1{@MU(t)};



%Non-Linear First-Order Conditions++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Insert here the non-linear FOCs in format g(x)=0

[1]   y(t) - @inv(t) - c(t) = 0;
[2]   betta * (@MU(t+1)/@MU(t)) * E(t)|R(t+1) - 1 = 0;
[3]   R(t) - @R(t) = 0;
[4]   y(t) - @F(t) = 0;
[5]   LOG(E(t)|z(t+1)) - psi*LOG(z(t)) = 0;


%Steady States [Closed Form]+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
None


%Steady State Non-Linear System [Manual]+++++++++++++++++++++++++++++++++++++++++++++++++
USE_FOCS=[0,1,2,3];

[1]   c_bar = 2.0;
[2]   k_bar = 30.0;
[3]   k_bar = k_bar**alpha;
[4]   R_bar = 1.01;

%Log-Linearized Model Equations++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
None


%Variance-Covariance Matrix++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Sigma = [sigma_eps**2];


%End Of Model File+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Option 3: Supplying the non-linear steady state system in the model file

Yet another way available for finding the model’s steady state is similar to the one in option one in that it uses a system of non-linear equations specified in this case directly inside the model template file. The reason why one would want to prefer this option over option one has to do with the fact that the steady state version of the non-linear first-order conditions of optimality can often collapse to much easier to work with and succincter equations which the model builder would want to write down explicitly inside the model file. So this example would be exemplified by the following model template file:

%Model Description+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
This is just a standard RBC model, as you can see.


%Model Information+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Name = Standard RBC Model;


%Parameters++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
rho       = 0.36;
delta     = 0.025;
R_bar     = 1.01;
eta       = 2.0;
psi       = 0.95;
z_bar     = 1.0;
sigma_eps = 0.052;


%Variable Vectors+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[1]  k(t):capital{endo}[log,bk]
[2]  c(t):consumption{con}[log,bk]
[4]  y(t):output{con}[log,bk]
[5]  z(t):eps(t):productivity{exo}[log,bk]
[6]  @inv(t):investment[log,bk]
[7]  @R(t):rrate

%Boundary Conditions++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
None


%Variable Substitution Non-Linear System++++++++++++++++++++++++++++++++++++++++++++++++
[1]   @inv(t)   = k(t)-(1-delta)*k(t-1);
[2]   @inv_bar  = SS{@inv(t)};
[2]   @F(t)     = z(t)*k(t-1)**rho;
[2]   @Fk(t)    = DIFF{@F(t),k(t-1)};
[2]   @Fk_bar   = SS{@Fk(t)};
[2]   @F_bar    = SS{@F(t)};
[3]   @R(t)     = 1+DIFF{@F(t),k(t-1)}-delta;
[4]   @R_bar    = SS{@R(t)};
[3]   @R(t+1)   = FF_1{@R(t)};
[4]   @U(t)     = c(t)**(1-eta)/(1-eta);
[5]   @MU(t)    = DIFF{@U(t),c(t)};
[5]   @MU_bar   = SS{@U(t)};
[6]   @MU(t+1)  = FF_1{@MU(t)};



%Non-Linear First-Order Conditions++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Insert here the non-linear FOCs in format g(x)=0

[1]   @F(t)-@inv(t)-c(t) = 0;
[2]   betta*(@MU(t+1)/@MU(t))*@R(t+1)-1 = 0;
[3]   @F(t)-y(t) = 0;
[4]   LOG(E(t)|z(t+1))-psi*LOG(z(t)) = 0;


%Steady States [Closed Form]+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
None


%Steady State Non-Linear System [Manual]+++++++++++++++++++++++++++++++++++++++++++++++++
[1]   @F_bar-@inv_bar-c_bar = 0;
[2]   betta*@R_bar-1 = 0;
[3]   betta*R_bar-1 = 0;
[4]   y_bar-@F_bar = 0;

[1]   c_bar = 1.0;
[2]   k_bar = 1.0;
[3]   y_bar = 1.0;
[4]   betta = 0.9;

%Log-Linearized Model Equations++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
None


%Variance-Covariance Matrix++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Sigma = [sigma_eps**2];


%End Of Model File+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

As one can see easily in this case, we are instructing the model to solve the 4 equation system in the four variables c_bar, k_bar, y_bar and betta. This is also a very common option to choose in order to obtain the model’s steady state efficiently and conveniently.

Option 4: Use the numerical root finder to solve for some steady states and get remaining ones residually

Option 4 is perhaps one of the most useful ways one can employ in order to obtain a DSGE model’s steady state solution as it focuses the numerical non-linear root-finding algorithm on a very small set of equations and unknown steady state variables, leaving the computation of the remaining steady state variables to be done separately and residually after the small set of steady state variables have been solved for. So using again a slightly tweaked version of the model file given in option 3 we could write this as:

%Model Description+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
This is just a standard RBC model, as you can see.


%Model Information+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Name = Standard RBC Model;


%Parameters++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
rho       = 0.36;
delta     = 0.025;
R_bar     = 1.01;
eta       = 2.0;
psi       = 0.95;
z_bar     = 1.0;
sigma_eps = 0.052;


%Variable Vectors+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[1]  k(t):capital{endo}[log,bk]
[2]  c(t):consumption{con}[log,bk]
[4]  y(t):output{con}[log,bk]
[5]  z(t):eps(t):productivity{exo}[log,bk]
[6]  @inv(t):investment[log,bk]
[7]  @R(t):rrate

%Boundary Conditions++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
None


%Variable Substitution Non-Linear System++++++++++++++++++++++++++++++++++++++++++++++++
[1]   @inv(t)   = k(t)-(1-delta)*k(t-1);
[2]   @inv_bar  = SS{@inv(t)};
[2]   @F(t)     = z(t)*k(t-1)**rho;
[2]   @Fk(t)    = DIFF{@F(t),k(t-1)};
[2]   @Fk_bar   = SS{@Fk(t)};
[2]   @F_bar    = SS{@F(t)};
[3]   @R(t)     = 1+DIFF{@F(t),k(t-1)}-delta;
[4]   @R_bar    = SS{@R(t)};
[3]   @R(t+1)   = FF_1{@R(t)};
[4]   @U(t)     = c(t)**(1-eta)/(1-eta);
[5]   @MU(t)    = DIFF{@U(t),c(t)};
[5]   @MU_bar   = SS{@U(t)};
[6]   @MU(t+1)  = FF_1{@MU(t)};



%Non-Linear First-Order Conditions++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Insert here the non-linear FOCs in format g(x)=0

[1]   @F(t)-@inv(t)-c(t) = 0;
[2]   betta*(@MU(t+1)/@MU(t))*@R(t+1)-1 = 0;
[3]   @F(t)-y(t) = 0;
[4]   LOG(E(t)|z(t+1))-psi*LOG(z(t)) = 0;


%Steady States [Closed Form]+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[1]   y_bar = @F_bar;


%Steady State Non-Linear System [Manual]+++++++++++++++++++++++++++++++++++++++++++++++++
[1]   @F_bar-@inv_bar-c_bar = 0;
[2]   betta*@R_bar-1 = 0;
[3]   betta*R_bar-1 = 0;

[1]   c_bar = 1.0;
[2]   k_bar = 1.0;
[3]   betta = 0.9;

%Log-Linearized Model Equations++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
None


%Variance-Covariance Matrix++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Sigma = [sigma_eps**2];


%End Of Model File+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

In this case we have simply taken the equation for y_bar outside of the section passed on to the non-linear root-finder and instead included it into the section for closed form steady state expressions. Whenever a model is instantiate like this, it first attempts to solve the smaller steady state system in the Manual section, before turning to the Closed Form section in which remaining steady states are computed residually based on the subset of steady states already solved numerically in the first step.

This is an extremely useful way of splitting down the problem, as many complex DSGE models often possess a large number of such residually determinable steady state values, while the core system on non-linear equations in a subset of steady states can be kept small in dimension and thus easier to solve. This really keeps the iteration burden on the non-linear solver to a minimum and often also allows the researcher to be less judicious in his choice of starting values leaving them at the generic default values. As a general rule, passing ever more complex and larger-dimensioned non-linear systems to the root-finding algorithm will decrease the chances of finding a solution easily, especially when simple generic starting values are employed. The issue of starting values take us straight to the next available option available to PyMacLab users.

Option 5: Use the numerical root finder to solve for steady states with pre-computed starting values

It is often useful and sometimes even outright necessary to supply the root-finding algorithm with pre-computed “intelligently” chosen initial starting values which are better than the generic choice of just passing a bunch of 1.0s to the system. To this end, whenever the list of generic starting values given in the numerical Manual section is a subset of the list of variable declarations in the closed form section, then the generic starting values automatically get replaced by the computed suggestions found in the Closed Form section. So an example of this would be:

%Model Description+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
This is just a standard RBC model, as you can see.


%Model Information+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Name = Standard RBC Model;


%Parameters++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
rho       = 0.36;
delta     = 0.025;
R_bar     = 1.01;
eta       = 2.0;
psi       = 0.95;
z_bar     = 1.0;
sigma_eps = 0.052;


%Variable Vectors+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[1]  k(t):capital{endo}[log,bk]
[2]  c(t):consumption{con}[log,bk]
[4]  y(t):output{con}[log,bk]
[5]  z(t):eps(t):productivity{exo}[log,bk]
[6]  @inv(t):investment[log,bk]
[7]  @R(t):rrate

%Boundary Conditions++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
None


%Variable Substitution Non-Linear System++++++++++++++++++++++++++++++++++++++++++++++++
[1]   @inv(t)   = k(t)-(1-delta)*k(t-1);
[2]   @inv_bar  = SS{@inv(t)};
[2]   @F(t)     = z(t)*k(t-1)**rho;
[2]   @Fk(t)    = DIFF{@F(t),k(t-1)};
[2]   @Fk_bar   = SS{@Fk(t)};
[2]   @F_bar    = SS{@F(t)};
[3]   @R(t)     = 1+DIFF{@F(t),k(t-1)}-delta;
[4]   @R_bar    = SS{@R(t)};
[3]   @R(t+1)   = FF_1{@R(t)};
[4]   @U(t)     = c(t)**(1-eta)/(1-eta);
[5]   @MU(t)    = DIFF{@U(t),c(t)};
[5]   @MU_bar   = SS{@U(t)};
[6]   @MU(t+1)  = FF_1{@MU(t)};



%Non-Linear First-Order Conditions++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Insert here the non-linear FOCs in format g(x)=0

[1]   @F(t)-@inv(t)-c(t) = 0;
[2]   betta*(@MU(t+1)/@MU(t))*@R(t+1)-1 = 0;
[3]   @F(t)-y(t) = 0;
[4]   LOG(E(t)|z(t+1))-psi*LOG(z(t)) = 0;


%Steady States [Closed Form]+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[1]   k_bar = 10.0;
[2]   y_bar = @F_bar;
[3]   c_bar = y_bar - delta*k_bar;
[4]   betta = 1/(1+@Fk_bar-delta);


%Steady State Non-Linear System [Manual]+++++++++++++++++++++++++++++++++++++++++++++++++
[1]   @F_bar-@inv_bar-c_bar = 0;
[2]   betta*@R_bar-1 = 0;
[3]   betta*R_bar-1 = 0;
[4]   y_bar-@F_bar = 0;

[1]   c_bar = 1.0;
[2]   k_bar = 1.0;
[3]   y_bar = 1.0;
[3]   betta = 0.9;

%Log-Linearized Model Equations++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
None


%Variance-Covariance Matrix++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Sigma = [sigma_eps**2];


%End Of Model File+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

As is apparent, in this case the suggested values for the steady states given in the closed form section exactly mirror or overlap with the steady variables to be searched for using the non-linear root finder specified in the Manual section in the model file (which means that the latter is a subset of the former). Whenever this overlap is of the described subset-type, the values in the Closed Form section will always be interpreted as suggested starting values passed on to the non-linear root finder. Notice that in this case it is also possible to omit the additional specification of the generic starting values in the Manual section alltogether. However it is advisable to leave them there to give the program a better way of checking the overlap of the two sets of variables. Whenever they are omitted, this specific case of computing the steady state is triggered whenever the number of suggested starting values in the Closed Form section is exactly equal to the number of non-linear equations in the Manual section.

Option 6: Computing the steady state values dictionary entirely outside of PyMacLab and passing it

This is the most straightforward but at the same time possibly less “encapsulated” method of obtining the steady state of a model. In this case, we ignore everything inside the model file for the purpose of computing the steady state and instead do everything outside’ of the DSGE instance externally. When done we plug the values into a dictionary and pass it to the DSGE instance at instantiation time using the keyword pm.newMOD(modelfile,sstate=sstatedic). This method may only be necessary for extremely large models in which obtaining the steady state is a task so difficult that it may have to be dealt with in a separate programming block.

Option 7: Finding the steady state by only supplying information in the Closed Form section

This is the most straightforward but at the same time possibly also least-used method for finding a steady state and will not be explained in greater depth here. In this variant, the Manual section is marked as unused employing the “None” keyword and only information in the Closed Form section is provided. Since only the most simple DSGE models afford this option of finding the steady state, we will not discuss this option any further.

More about the steady state and its relevance to perturbation methods - A Detour

Perturbation methods discussed in the literature relevant to the discussion of solution methods for dynamic stochastic general equilibrium (DSGE) models are similar in spirit to Taylor series approximations of non-linear functions which are being taught to students in introductory/intermediate quantitative methods courses. Based on a functions derivatives, approximations are taken around a specific point or value of the parameter domain in order to approximate the fuction’s corresponding value when these values are perturbed away by small degrees from that point around which the approximation was taken.

This approach lends itself very well to the approximation of the exact solution of non-linear systems of expectatonal difference equations (i.e. DSGE models) chiefly because of two reasons. First of all, the assumption that shocks (perturbations) are not too large in some some in a real-world economic system is a reasonable one for most of the times. Secondly, while the general discussion of perturbation methods does not prescribe a preferred point around which the approximation should be taken, dynamic economic systems exhibting a steady state (equilibrium fixed point in the absence of any shocks) exhibit a perfect candidate value around which the approximation should be taken.

../_images/riskyss.jpg

In the most recent literature [2] this point has however been given additional attention in that researchers typically consider two possible candidate steady state equilibrium points [3] around which to form the approximation. One of them is the traditional deterministic steady state to which the economic system gravitates when (all moments of) future shocks are assumed to be zero, while the other steady state, suitably called the risky or stochastic steady state is the one the system comes to rest to when agents populating the system know that future shocks will continue to occur based on some known (non-Knightian, if you like) distribution of those shocks.

Footnotes

[2]In other more complicated cases the law of motion of some exogenous shock process may depend on other endogenous states of the system. In this case we would probably want to pass the line to the non-linear root finder as its specification would influence the steady state value of other steady state variables.
[3]Actually, there are more than 2 possible definitions of suitable steady state points. Here we want to draw attention to a stochastic point which has been defined as follows: “The risky steady-state is the point where agents choose to stay at a given date if they expect future risk and if the realization of shocks is 0 at this date.” (Coeurdacier, Rey and Winant, 2011)

Tutorial 6 - Dynamic Solution Methods

Introduction

In the previous tutorial we discovered the general structure of the PyMacLab DSGE model instance and saw how this programming approach lent itself well to the idea of inspecting and exploring the instantiated models’ current state, summarized by its data fields and supplied instance methods equipping them with functionality. This section finally discusses how the library equips DSGE model instances with methods which make use of the models’ computed Jacobian and Hessian, which are evaluated at the models’ numerical steady state. As a short reminder, we may recall here that it is often this step of obtaining a steady state can prove difficult in the case of a number of well-known models. Notwithstanding, for the remainder of this section we will assume that a steady state has successfully been attained and that the model’s Jacobian and Hessian have been computed. Let’s first start with our usual setup of lines of code entered into an IPython shell:

# Import the pymaclab module into its namespace
In [1]: import pymaclab as pm
In [2]: from pymaclab.modfiles import models

# Instantiate a new DSGE model instance like so
In [4]: rbc1 = pm.newMOD(models.stable.rbc1)

Since we have not supplied the instantiation call with the initlev argument, the model has been solved out completely, which includes the computation of a preferred dynamic solution, which in the current library’s version is Paul Klein’s 1st-order accurate method based on the Schur Decomposition. This means that one particular solution - if existing and found - has already been computed by default. Readers familar with this solution method will know that this method requires a partitioning of the model’s evaluated Jacobian into two separate matrices, which for simplicity we denote A and B. They are computed by default and are for instance inspectable at rbc1.jAA and rbc1.jBB.

The Jacobian and Hessian: A Detour

One of the great advantages of PyMacLab is that the Jacobian and Hessian are NOT computed or approximated numerically using the method of finite differences, but are calculated in exact fashion analytically using the special-purpose Python library Sympy, which is a CAS - or computer algebra system, similar in functionality to Mathematica and Maple. You can inspect the analytical counterparts to the exact numerical Jacobian and Hessian which have not yet been evaluated numerically at rbc1.jdicc and rbc1.hdicc, where the latter reference actually refers to the 3-dimensional analytical Hessian. If the model is composed of n equations describing equilibrium and optimality conditions, then the Jacobian is made up of n\times\left(n\times 2\right) elements and has dimension \left(n,n\times 2\right), because the derivatives are formed not only w.r.t. to current-period, but also future-period variables. Equally, the Hessian is a 3-dimensional matrix of dimension \left(n,n\times 2,n\times 2\right) [#f1]_.

Notice that for the numerical counterpart to the 3-dimensional Hessian, PyMacLab instead uses an alternatively dimensioned version of dimension \left(n\times n\times 2,n\times 2\right), which is the Magnus and Neudecker definition of a Hessian and is useful when one wishes to avoid using matrices of dimension larger than 2 and the corresponding Tensor notation. Again, you can exploit PyMacLab’s DSGE model instance’s design in order to inspect the derivatives contained in the Jacobian and Hessian inside an IPython interactive shell environment, such as follows:

# Import the pymaclab module into its namespace, also import os module
In [1]: import pymaclab as pm
In [2]: from pymaclab.modfiles import models

# Instantiate a new DSGE model instance like so
In [4]: rbc1 = pm.newMOD(models.stable.rbc1)


# Check some elements of the analytical Jacobian and Hessian
# Jacobian, equation 3, derivative w.r.t. k(t-1). Notice that Python arrays are 0-indexed.
In [5]: rbc1.jdicc[2]['k(t-1)']
'rho*exp(z(t))*exp(k(t-1)*rho)'

# You could also retrieve the evaluated equivalent of the above expression.
# Here you need to know the position of k(t-1) in the ordered list of derivatives, you can check
# that k(t-1) is on position 5 by inspecting rbc1.var_li
In [6]: rbc1.numj[2,5]
1.3356245171444843

# Hessian, equation 3, 2nd derivative w.r.t k(t-1).
In [6]: rbc1.hdicc[2]['k(t-1)']['k(t-1)']
'rho**2*exp(z(t))*exp(k(t-1)*rho)'

# The numerical evaluated equivalent can be retrieved as well.
# We are not retrieving the above value via rbc1.numh[2,5,5] as we are not working with
# the usual 3D notation of Hessians, but with the Magnus & Neudecker 2D definition of it.
# The result is correct, as the 2nd derivative is just rho(=0.36) times the first derivate.
In [7]: rbc1.numh[21,5]
0.48082482617201433

Now we have explored the ins and outs of PyMacLab’s way of handling the computation of a DSGE model’s Jacobian and Hessian. Equipped with these building blocks, it is now time to move on to a discussion of the actual solution methods which PyMacLab provides by default.

Dynamic Solution Methods - Nth-order Perturbation

Solving nonlinear rational expectations DSGE models via the method of perturbation represents an approximate solution around the computed steady state of the model. Since this approach is not too dissimilar from a Taylor Series expansion of a function around some point students learn about in some 101 maths course, it should come as no surprise that here too 1st and higher-order derivatives are needed in order to arrive at solutions which increasinly reflect the true exact solution of the system.

PyMacLab offers methods suitable for computing such approximated solutions based on linearization techniques which can either be 1st- or 2nd-order accurate. In order to obtain these solutions, we make use of Paul Klein’s solution algorithms, which are available on the internet and have been incorporated into PyMacLab. Needless to say, Klein’s 1st-order accurate method using the Schur Decomposition has been around for a while and only requires knowledge of the models Jacobian, while his latest paper (co-authored with Paul Gomme) spelling out the solution of the 2nd-order accurate approximation, also requires knowledge of the Hessian.

Choosing the degree of approximation

At the time of writing these words, PyMacLab includes full support for both of these two methods, where the first method has been made available by binding Klein’s original Fortran code into PyMacLab and making it accessible via the node rbc1.modsolvers.forkleind which provides the solution method callable via rbc1.modsolvers.forkleind.solve(). Once this method has been called and a solution has been found, it is essentially encapsulated in the matrices available at rbc1.modsolvers.forkleind.P and rbc1.modsolvers.forkleind.F, which represent matrices of dynamic elasticities summarizing the optimal laws of motion for the set of endogenous state and control variables, respectively. Since this method is actually internally calling a compiled Fortran dynamically linked library, its name is prefixed with for.

Klein & Gomme’s 2nd-order accurate method uses the solution from the 1st-order accurate method as a starting point but in addition also makes use of the model’s Hessian and the information provided by the model’s shocks variance-covariance matrix, in order to produce solutions which are risk-adjusted in some loosely defined sense. This solution method therefore no longer displays the well-known property of certainty equivalence for which first-order approximations are so well known for. At the moment, this solution method is completely implemented in the Python language itself and is callable at rbc1.modsolvers.pyklein2d.solve(). As already mentioned, the method makes use of rbc1.modsolvers.forkleind.P and rbc1.modsolvers.forkleind.F, the variance-covariance matrix rbc1.modsolvers.pyklein2d.ssigma, and the Magnus & Neudecker definition of the Hessian rbc1.modsolvers.pyklein2d.hes. It’s solution is encapsulated in the following objects:

Object Description
rbc1.modsolvers.forkleind.P Matrix of elasticities describing optimal law of motion for endog. state variables, 1st-order part
rbc1.modsolvers.forkleind.F Matrix of elasticities describing optimal law of motion for endog. state variables, 1st-order part
rbc1.modsolvers.forkleind.G Matrix of elasticities describing optimal law of motion for endog. state variables, 2nd-order part
rbc1.modsolvers.forkleind.E Matrix of elasticities describing optimal law of motion for endog. state variables, 2nd-order part
rbc1.modsolvers.forkleind.KX Array of risk-adjustment values for steady states of endogenous state variables
rbc1.modsolvers.forkleind.KY Array of risk-adjustment values for steady states of control variables

Once the above mentioned matrices are calculated, the solutions to either the 1st-order or 2nd-order accurate approximations are available and can be used by researchers to compute (filtered) simulations as well as impulse-response functions in order to either plot them or generate summary statistics from them. Luckily, neither of this has to be done by hand, as simulation- and IRF-generating methods are already supplied and convenience plotting functions are also readily available. But this will be the topic of our next tutorial in the tutorial series for PyMacLab.

Footnotes

[4]Obviously, we are abusing clearly defined mathematical definitions here to some extent, as a classical Jacobian would be of dimension \left(n,n\right) and a classical Hessian of dimension \left(n,n,n\right). The only reason why here computed dimensions tend to be twice as large has to do with the fact that for Klein’s first-order accurate solution method we require knowledge of derivatives w.r.t. both current and future (expected) versions of the set of all variables.

Tutorial 7 - Simulating DSGE models

Introduction

In the previous tutorial we discovered the general structure of the PyMacLab DSGE model instance and saw how this programming approach lent itself well to the idea of inspecting and exploring the instantiated models’ current state, summarized by its data fields and supplied instance methods equipping them with functionality. This section finally discusses how the library equips DSGE model instances with methods which make use of the models’ computed Jacobian and Hessian, which are evaluated at the models’ numerical steady state. As a short reminder, we may recall here that it is often this step of obtaining a steady state can prove difficult in the case of a number of well-known models. Notwithstanding, for the remainder of this section we will assume that a steady state has successfully been attained and that the model’s Jacobian and Hessian have been computed. Let’s first start with our usual setup of lines of code entered into an IPython shell:

# Import the pymaclab module into its namespace, also import os module
In [1]: import pymaclab as pm
In [2]: from pymaclab.modfiles import models

# Instantiate a new DSGE model instance like so
In [4]: rbc1 = pm.newMOD(models.stable.rbc1_res)

Simulating the model

Recall that instantiating a DSGE model without any additional parameters means that an automatic attempt of finding the steady state is being made and if successfully found the model is already solved dynamically using a preferred (1st-order approximate) method. This means that you now have all the information you need to simulate the model as well as to generate impulse-response functions and plot them. Let’s focus on the simulations first, they are being generate using the following command:

# Import the pymaclab module into its namespace, also import os module
In [1]: import pymaclab as pm
In [2]: from pymaclab.modfiles import models

# Instantiate a new DSGE model instance like so
In [3]: rbc1 = pm.newMOD(models.stable.rbc1_res)

# Now simulate the model
In [4]: rbc1.modsolvers.forkleind.sim(200,('productivity'))

This would simulate the model for 200 time periods using only iid shocks hitting the law of motion of the total factor productivity term in the model. Notice that here RBC1 is a very simple model only containing one structural shock, but more complicated models may possess more than one exogenous state variable. In that case, if you called instead:

# Import the pymaclab module into its namespace, also import os module
In [1]: import pymaclab as pm
In [2]: from pymaclab.modfiles import models

# Instantiate a new DSGE model instance like so
In [3]: rbc1 = pm.newMOD(models.stable.rbc1_res)

# Now simulate the model
In [4]: rbc1.modsolvers.forkleind.sim(200)

The model would be simulated using all shocks (exogenous state variables) specified in the model. However, since rbc1 only contains one shock, the two variants shown here of simulating the model would yield the same results as “productivity” is is the only exogenous state available here anyway. We can then also graph the simulation to get a better understanding of the model by running the command:

# Import the pymaclab module into its namespace, also import os module
In [1]: import pymaclab as pm
In [2]: from pymaclab.modfiles import models

# Also import matplotlib.pyplot for showing the graph
In [3]: from matplotlib import pyplot as plt

# Instantiate a new DSGE model instance like so
In [4]: rbc1 = pm.newMOD(models.stable.rbc1_res)

# Now solve and simulate the model
In [5]: rbc1.modsolvers.forkleind.solve()
In [6]: rbc1.modsolvers.forkleind.sim(200)

# Plot the simulation and show it on screen
In [7]: rbc1.modsolvers.forkleind.show_sim(('output','consumption'))
In [8]: plt.show()

This produces the following nice graph. Notice that you must specify the variables to be graphed and all simulated data is filtered according to the argument passed to each variable in the model file. So the key “hp” produces hp-filtered data, the key “bk” results in Baxter-King-filtered data while the key “cf” leads to cycles extraced using the Christiano-Fitzgerald filter.

(Source code, png, hires.png, pdf)

../_images/test4.png

Cross-correlation tables

Notice that filtered simulations are always stored in data fields which means that statistics such as correlations at leads and lags can easily be computed as well. Specifically, the simlulated data corresponding to the above graph can be retrieved from the object rbc1.modsolver.forkleind.insim [#f1]_. There already exist a number of simple convenience functions allowing users to generate cross-correlation tables for simulated data. The functions can be used as follows:

# Import the pymaclab module into its namespace, also import os module
In [1]: import pymaclab as pm
In [2]: from pymaclab.modfiles import models

# Also import matplotlib.pyplot for showing the graph
In [3]: from matplotlib import pyplot as plt
In [4]: from copy import deepcopy

# Instantiate a new DSGE model instance like so
In [5]: rbc1 = pm.newMOD(models.stable.rbc1_res)

# Now solve and simulate the model
In [6]: rbc1.modsolvers.forkleind.solve()
In [7]: rbc1.modsolvers.forkleind.sim(200)

# Generate the cross-correlation table and show it
# Produce table with 4 lags and 4 leads using output as baseline
In [8]: rbc1.modsolvers.forkleind.mkact('output',(4,4))
In [9]: rbc1.modsolvers.forkleind.show_act()

Autocorrelation table, current output
=================================================================
productivity  |-0.016  0.109  0.335  0.663  0.997  0.619  0.264  0.034 -0.084
capital       |-0.433 -0.429 -0.381 -0.258 -0.024  0.318  0.522  0.599  0.596
consumption   |-0.134 -0.009  0.228  0.587  0.98   0.699  0.404  0.198  0.08
output        |-0.049  0.077  0.308  0.647  1.     0.646  0.305  0.08  -0.039

If users wish to obtain the data of the above table directly in order to import them into a different environment more suitable for producing publication-quality tables, the cross-correlation data can be accesssed at rbc1.modsolvers.forkleind.actm which is a matrix object of cross-correlations at the leads and lags specified in the previous calling function generating that table data.

Simulating while keeping random shocks fixed

Yet another useful feature to know about is that after each call to rbc1.modsolvers.forkleind.sim() the vector of randomly drawn iid shocks gets saved into object rbc1.modsolver.forkleind.shockvec. This is useful because when calling the simulation function, we can also pass an existing pre-computed vector of shocks as an argument instead of allowing the call to generate a new draw of random shocks. That way we can keep the random shocks fixed from model run to model run. So this would be accomplished as follows:

# Import the pymaclab module into its namespace, also import os module
In [1]: import pymaclab as pm
In [2]: from pymaclab.modfiles import models

# Also import matplotlib.pyplot for showing the graph
In [3]: from matplotlib import pyplot as plt
In [4]: from copy import deepcopy

# Instantiate a new DSGE model instance like so
In [5]: rbc1 = pm.newMOD(models.stable.rbc1_res)

# Now solve and simulate the model
In [6]: rbc1.modsolvers.forkleind.solve()
In [7]: rbc1.modsolvers.forkleind.sim(200)

# Plot the simulation and show it on screen
In [8]: rbc1.modsolvers.forkleind.show_sim(('output','consumption'))
In [9]: plt.show()

# Now save the shocks, by saving a clone or copy, instead of a reference
In [10]: shockv = deepcopy(rbc1.modsolvers.forkleind.shockvec)

# Now we could run the simulation again, this time passing the randomly drawn shocks
In [11]: rbc1.modsolvers.forkleind.sim(200,shockvec=shockv)

# Plot the simulation and show it on screen
In [12]: rbc1.modsolvers.forkleind.show_sim(('output','consumption'))
In [13]: plt.show()

Notice that in this script the graphs plotted to screen using the plt.show() command will produce identical graphs as the random draw of shocks only occurs in the first call to sim() while in the second it gets passed as an argument with a value retrieved and retained from the first simulation run. The reason why this feature is so useful has to do with the fact that sometimes we wish to produce summary statistics from simulation runs of one version of a model, then tweak the model’s properties dynamically at runtime and re-compute the very same summary statistics, under the assumption of holding the iid errors fixed, so that we can observe the pure net effect from changing the model’s properties elimiting any unwanted variation from “sampling variation”. As an example of this we demonstrate a script in which simulations are run and plotted under different filtering assumption.

# Import the pymaclab module into its namespace, also import os module
In [1]: import pymaclab as pm
In [2]: from pymaclab.modfiles import models

# Also import matplotlib.pyplot for showing the graph
In [3]: from matplotlib import pyplot as plt
In [4]: from copy import deepcopy

# Instantiate a new DSGE model instance like so
In [5]: rbc1 = pm.newMOD(models.stable.rbc1_res,mk_hessian=False)

# Now solve and simulate the model
In [6]: rbc1.modsolvers.forkleind.solve()
In [7]: rbc1.modsolvers.forkleind.sim(200)

# Plot the simulation and show it on screen
In [8]: rbc1.modsolvers.forkleind.show_sim(('output','consumption'))
In [9]: plt.show()

# Now save the shocks, by saving a clone or copy, instead of a reference
In [10]: shockv = deepcopy(rbc1.modsolvers.forkleind.shockvec)

# Change the filtering assumption of output and consumption using the queued updater branch
In [11]: rbc1.updaters_queued.vardic['con']['mod'][0][1] = 'hp'
In [12]: rbc1.updaters_queued.vardic['con']['mod'][1][1] = 'hp'
In [13]: rbc1.updaters_queued.process_queue()

# Now we could run the simulation again, this time passing the randomly drawn shocks
In [14]: rbc1.modsolvers.forkleind.solve()
In [15]: rbc1.modsolvers.forkleind.sim(200,shockvec=shockv)

# Plot the simulation and show it on screen
In [16]: rbc1.modsolvers.forkleind.show_sim(('output','consumption'))
In [17]: plt.show()

# Change the filtering assumption of output and consumption using the queued updater branch
In [18]: rbc1.updaters_queued.vardic['con']['mod'][0][1] = 'cf'
In [19]: rbc1.updaters_queued.vardic['con']['mod'][1][1] = 'cf'
In [20]: rbc1.updaters_queued.process_queue()

# Now we could run the simulation again, this time passing the randomly drawn shocks
In [21]: rbc1.modsolvers.forkleind.solve()
In [22]: rbc1.modsolvers.forkleind.sim(200,shockvec=shockv)

# Plot the simulation and show it on screen
In [23]: rbc1.modsolvers.forkleind.show_sim(('output','consumption'))
In [24]: plt.show()

(Source code, png, hires.png, pdf)

../_images/test6_00_00.png

As is apparent from the three plots produced above, the simulated data is first filtered using the Baxter-King filter, then the more commonly used Hodrick-Prescott filter and finally the Christian-Fitzgerald asymmetric filter. Notice that the BK filter by default (or rather by specification) cuts off 6 time periods at the beginning and at the end of the simulated sample. The purpose for using any of the three filters is of course to make the simulated data stationary and to extract the cycle only.

Generating impulse-response functions

Dynamic solutions obtained to first-order approximated DSGE models using the method of perturbations have a great deal in common with standard Vector Autoregression (VAR) models commonly used in applied Macroeconometrics. This in turn implies that solved DSGE models can be described using so-called impulse-response functions (also abbreviated as IRFs) or impulse-response graphs which show how the solved model responds to a one-off shock to a particular exogenous state variable. In PyMacLab this can easily be achieved as follows:

# Import the pymaclab module into its namespace, also import os module
In [1]: import pymaclab as pm
In [2]: from pymaclab.modfiles import models

# Also import matplotlib.pyplot for showing the graph
In [3]: from matplotlib import pyplot as plt

# Instantiate a new DSGE model instance like so
In [4]: rbc1 = pm.newMOD(models.stable.rbc1_res)

# Now solve and simulate the model
In [5]: rbc1.modsolvers.forkleind.solve()
In [6]: rbc1.modsolvers.forkleind.irf(100,('productivity',))

# Plot the simulation and show it on screen
In [7]: rbc1.modsolvers.forkleind.show_irf(('output','consumption'))
In [8]: plt.show()

This produces the following nice graph. Notice that here the shock to total productivity has been normalized to 100%.

(Source code, png, hires.png, pdf)

../_images/test5.png

Notice that the IRF-generating function also accepts (a tuple or list of) values defining the sign and magnitude of the shock. So the above could also have been written as follows:

# Import the pymaclab module into its namespace, also import os module
In [1]: import pymaclab as pm
In [2]: from pymaclab.modfiles import models

# Also import matplotlib.pyplot for showing the graph
In [3]: from matplotlib import pyplot as plt

# Instantiate a new DSGE model instance like so
In [4]: rbc1 = pm.newMOD(models.stable.rbc1_res)

# Now solve and simulate the model
In [5]: rbc1.modsolvers.forkleind.solve()
In [6]: rbc1.modsolvers.forkleind.irf(100,('productivity',),(-2.0,))

# Plot the simulation and show it on screen
In [7]: rbc1.modsolvers.forkleind.show_irf(('output','consumption'))
In [8]: plt.show()

Footnotes

[5]The simulated time series data contained within rbc1.modsolvers.forkleind.insim is NOT filtered yet. In designing the library I have decided to delay filtering to the stage where the user calls rbc1.modsolvers.forkleind.show_sim() or when similar functions are called to for instance generate the cross-correlation table.

Tutorial 8 - PyMacLab-to-Dynare++ Converter

Introduction

PyMacLab comes shipped with a very useful dynare++ wrapper/translator which essentially possesses two types of functionality. First of all, it can translate a PyMacLab-formatted DSGE model files into a valid dynare++ model file. Secondly, it internally can pass this translated model file to dynare++, instruct it to solve the model, and then pass all of the solution objects back to PyMacLab where they get stored in the solution branch model.modsolvers.dynare..

Note

Notice that in order for PyMacLab’s dynare++ wrapper/translator to function correctly, in a standard *nix environment you need to have a functioning (compiled or otherwise obtained) dynare++ binary (executable) program file in your local PATH. This means that dynare++ on your system can be called from any location inside you shell. Also important is that the Python library Mako is installed which is a templating library used in generating the dynare++ conformable model file. This library should however be automatically installed when you installed PyMacLab on your system.

The most basic functionality of the wrapper/translator works as follows calling the method model.mk_dynare() without any additional arguments:

# Import the pymaclab module into its namespace, also import os module
In [1]: import pymaclab as pm
In [2]: from pymaclab.modfiles import models

# Instantiate a new DSGE model instance like so
In [4]: rbc1 = pm.newMOD(models.stable.rbc1_focs)

# Use the dynare++ translator/wrapper as follows, by calling
In [5]: rbc1.mk_dynare()

# If the model solved, then check the contents of the dynare branch
In [6]: dir(rbc1.modsolvers.dynare)
['FF',
 'PP',
 '__class__',
 '__delattr__',
 '__dict__',
 '__doc__',
 '__format__',
 '__getattribute__',
 '__hash__',
 '__init__',
 '__module__',
 '__new__',
 '__reduce__',
 '__reduce_ex__',
 '__repr__',
 '__setattr__',
 '__sizeof__',
 '__str__',
 '__subclasshook__',
 '__weakref__',
 'attdic',
 'dyn_g_0',
 'dyn_g_1',
 'dyn_i_R',
 'dyn_i_c',
 'dyn_i_eps',
 'dyn_i_k',
 'dyn_i_y',
 'dyn_i_z',
 'dyn_irfm_eps_mean',
 'dyn_irfm_eps_var',
 'dyn_irfp_eps_mean',
 'dyn_irfp_eps_var',
 'dyn_mean',
 'dyn_nboth',
 'dyn_nforw',
 'dyn_npred',
 'dyn_nstat',
 'dyn_shocks',
 'dyn_ss',
 'dyn_state_vars',
 'dyn_steady_states',
 'dyn_vars',
 'dyn_vcov',
 'dyn_vcov_exo',
 'dynorder',
 'dynsorder',
 'modfile',
 'sstate']

Notice that all of dynare++’s computed solution objects (which are mostly matrices containing various types of information characterizing the solution to the model) have been attached to this solution branch as standard numpy matrices/scalars, which can then be inspected inside a normal scientific Python stack environment. The only two objects which are being computed in a derived fashion using dynare++’s results are the matrices model.modsolvers.dynare.PP and model.modsolvers.dynare.FF which are the implied solutions from dynare++’s output in the solution format which get internally computed in PyMacLab using the model.modsolvers.forkleind.solve() method instead.

This is useful as they can be directly compared and so checked for consistency between the different solution methods. Also important is the attribute model.modsolvers.dynare.modfile which is a string containing the PyMacLab-to-dynare++ translated model file which was passed to dynare++ in order to find the solution. Notice that in test runs it has not been uncommon to encounter situations in which PyMacLab’s internal algorithms successfully computed a model’s steady state, while dynare++ failed in doing so.

Finally, bear in mind that dynare++’s standard solution method uses an approximation around the level of the model’s steady state, and not the logarithm! This means that when you wish to compare solutions from dynare++ with those computed internally using PyMacLab’s model.modsolvers.forkleind.solve() method, to take one specific example, then first make sure to omit the [log] qualifier in the variable section for each variable defined in the PyMacLab DSGE model file. That way you are instructing PyMacLab to also employ approximations around the level and not the logarithm of the model’s SS.

The mk_dynare() solution method’s options

Sometimes it may be useful to instruct PyMacLab not to solve the model using dynare++ but instead just to generate the dynare++ conformable file for inspection by the model-builder. This can be done as follows:

# Import the pymaclab module into its namespace, also import os module
In [1]: import pymaclab as pm
In [2]: from pymaclab.modfiles import models

# Instantiate a new DSGE model instance like so
In [4]: rbc1 = pm.newMOD(models.stable.rbc1_focs)

# Use the dynare++ translator/wrapper as follows, by calling
In [5]: rbc1.mk_dynare(fpath='./test.mod')

In this case, no finding of a model solution would be attempted, but instead the PyMacLab-to-dynare++ translated model file would be copied into the location passed to the fpath= option. In the above example we are saving the translated dynare++ file into the current directory using the filename test.mod. For the above model file the generated dynare++ model file would look as follows:

var k, c, y, R, z;
varexo eps;

parameters z_bar, psi, sigma_eps, betta, eta, rho, delta, c_bar, R_bar, k_bar, y_bar;

z_bar = 1.0;
psi = 0.95;
sigma_eps = 0.052;
betta = 0.990099009901;
eta = 2.0;
rho = 0.36;
delta = 0.025;
c_bar = 2.0;
R_bar = 1.01;
k_bar = 30.0;
y_bar = 3.40222985854;


model;
0=y - (k-(1-delta)*k(-1)) - c;
0=betta * ((((c(+1)^(-eta))))/((c^(-eta)))) * R(+1) - 1;
0=R - (1+(z(-1)*k(-1)^(-1 + rho)*rho)-delta);
0=y - (z(-1)*k(-1)^rho);
0=log(z) - psi*log(z(-1)) - eps;
end;


initval;
R = 1.01;
c = 2.0;
k = 30.0;
y = 3.40222985854;
z = 1.0;
z = 1.0;
end;


order = 1;

vcov = [ 0.002704 ];

Notice here that the starting values passed to dynare++’s non-linear steady state solver are already the solutions found by the non-linear solution algorithm used internally within the PyMacLab library itself. This particular function is useful in the sense that it may help model builders identify bugs or translation mistakes which may have occurred when converting the PyMacLab DSGE model to a dynare++ conformable format. Finally the full set of arguments which can be passed to the mk_dynare() method are given by model.mk_dynare(order=1,centralize=False,fpath=None,focli=None).

The first argument order= is clearly used to specify the order of approximation used, centralize= can be used to specify whether the approximation should be computed around the deterministic or the stochastic (or risky) steady state [#f1]_, while focli= accepts a list or tuple of index numbers in order to pick out a subset of specific equations from PyMacLab’s DSGE model file’s declaration of the nonlinear system of equations of the model’s first-order conditions of optimality. fpath= has already been discussed in the above and is an option which can be used in order to output the dynare++ conformable model file to a specific location in your computer’s file system.

Footnotes

[6]Dynare++ usually as default behaviour computes approximations around the stochastic or risk steady state. In order to be able to make comparisons between dynare++’s solutions and those obtained using most of the solution algorithms which are available internally to PyMacLab, it is important to pass the centralize=False option, which is the standard choice for mk_dynare‘s method call. That way dynare++ computes the solutions around the deterministic steady state of the model.

Table Of Contents

Previous topic

Macroeconomic Analysis in Python

Next topic

Release History

This Page