Perturbation AIM

Perturbation AIM is an algorithm that solves a linear or nonlinear dynamic, discrete-time system of rational expectations equations to arbitrarily high order, in the sense of a Taylor series approximation. The algorithm requires that your system of equations has at least one steady-state point around which to approximate a solution, that a solution exists in at least a neighborhood of that steady-state point, and that the solution and original system of equations be analytic (or at least sufficiently smooth) in at least a neighborhood of that steady-state point. The solution does not necessarily need to be first-order stable around the steady state: unit roots are handled automatically, and even explosive first-order roots are acceptable but require the modeler to specify which explosive roots are to be allowed (this requires the user to make minor modifications to the code).

An important feature of the algorithm is that, so long as the solution to your problem is analytic, Perturbation AIM provides a global solution in the sense that it is guaranteed to converge to the true solution over the whole domain of convergence of the Taylor series, and this domain is often very large. Moreover, the convergence is uniform on compact subsets, which implies that there exists a finite order of approximation n which achieves any desired level of accuracy over any given (possibly large) compact set. Note the contrast here with fixed-order approximation methods such as log-linearization or second-order approximation, which provide only local solutions.

For additional details about the class of admissible problems and the solution algorithm, see the following paper:

We have implemented the Perturbation AIM algorithm in Mathematica.  While second-order approximations can be done fairly easily in any computer language, generalization to higher orders is dramatically simpler in Mathematica, and Mathematica’s large library of built-in routines for symbolic differentiation and solving linear systems of equations has the following advantages:

• Allows for much simpler and intuitive code (less than 200 lines), making the implementation of the algorithm more transparent and more likely to be free of human coding error (bugs).
• Allows the algorithm to manipulate symbolic expressions and thus proceed in the same way and produce the same output as if the user were performing it by hand, making the implementation of the algorithm and the generated output more intuitive.
• Eliminates the need for user-written differentiation and linear solution routines, which are likely to be more amateur than built-in Mathematica routines and hence more likely to suffer from bugs and numerical inaccuracies.
• Allows the use of symbolic differentiation, improving numerical accuracy, particularly at higher orders.
• Allows the computation of all coefficients to arbitrarily high precision, eliminating inaccuracies from machine-precision computations that cumulate with each recursive step of the algorithm.

In addition to the internal advantages of PerturbationAIM described above, there are a number of advantages from the economic modeler's point of view:

• No need to categorize variables as “predetermined” or “non-predetermined.” Knowing that all variables dated t–1 or earlier are observed and all variables dated t are to be solved is sufficient for the computer to determine which linear combinations of variables are “predetermined” and which “non-predetermined.” Anderson and Moore (1985), Sims (2000), and Collard and Juillard (2001) also make this observation. This feature of the code becomes increasingly convenient as the number of variables in the model increases.
• No need for the modeler to substitute out identities. The computer is completely capable of handling these substitutions.
• Ability to solve models with arbitrarily long lags and leads (including shocks dated later than t or t+1), with no need for the modeler to transform the model by hand. Note also that the usual trick for linear models of defining an auxiliary variable z_t ≡ E_tx_t+1, and then considering z_t+1, fails for nonlinear models because F(E_tx_t+1) is unequal to E_tF(x_t+1) in general. Perturbation AIM handles multiple leads correctly in the nonlinear as well as the linear case.
• Ability to solve models with non-normally distributed shocks. Perturbation AIM allows the modeler to substitute in any set of moments for the stochastic shocks of the model, allowing consideration of models with essentially any shock distribution (for which the moments are finite).

Latest Version of the Code (Perturbation AIM 2.8.3)

The current version of the code is Perturbation AIM 2.8.3. Relative to version 2.8.2, the current version fixes some warnings that were appearing in Mathematica 10.

Relative to version 2.7, Perturbation AIM 2.8 includes:

• an improved algorithm that computes the steady state of the model to greater accuracy
• some minor numerical improvements that increase the performance and accuracy of arbitrary-precision solutions to a model
• a workaround for a bug in Mathematica that sometimes caused difficulty solving for the coefficients on the stochastic terms of a model, and
• some small improvements that should eliminate the occurrence of variable-name conflicts with user-defined variables (x, n, i, etc.)

Relative to version 2.6, Perturbation AIM 2.7 and 2.8 contain an internal optimization that speeds up the solution time.

Relative to version 2.5, Perturbation AIM 2.6 and later contain some numerical improvements (better precision and numerical stability) and some user interface enhancements (catch a few more user input errors, report estimated precision).

Relative to version 2.4, Perturbation AIM 2.5 and later feature:

• A new AIMPrecision environment variable, which allows the user to easily specify the internal numerical precision that Mathematica should use to solve the model.  Although Perturbation AIM has always had the capability of using aribtrary-precision arthmetic, the new environment variable and other internal enhancements make this feature of the code much more user-friendly.
• Greater numerical precision and robustness. For example, when steady-state or first-order coefficients are found to be essentially zero (less than the AIMZeroTol environment variable), then Perturbation AIM now imposes exact zeros for those coefficients and re-solves with the hard zero constraint imposed. This increases the numerical precision of the solution.
• The ability to call Perturbation AIM with model equations and parameter values as separate arguments. This makes parameter estimation with Perturbation AIM faster, since now the derivatives of the model only need to be computed once and then the various parameter vectors can be substituted into those derivatives repeatedly.
• Improved compatibility with Mathematica 6.0. The previous version of Perturbation AIM generated occasional warnings regarding deprecated function calls. Those function calls have been updated in the current version of the code.

Some examples that utilize these new features of the code are provided in the “Using the Code” section below.

Obtaining the Code

Click on one of the following links to download Perturbation AIM (if your browser loads it in as a text file, then just save the text file to your hard drive):

For serious research work, you will want to use the most recent version of the code (version 2.8). The previous versions (especially 2.4 and the “simple” version 2.4) are provided for those who wish to understand how the algorithm works, and to provide backward-compatible versions of the algorithm that work on older versions of Mathematica. The “simple” version has about 200 lines of code (plus 100 lines of comments) and has been designed to be more transparent than the full version, hence easier to understand. The “full” version 2.4 is functionally identical to the “simple” version 2.4 but is a little more complicated and less transparent due to the implementation of some additional internal tricks that speed up computation and reduce the likelihood of numerical inaccuracies, particularly for larger models and higher orders of approximation. Versions 2.5 and later have even more internal tricks, optimizations, and enhancements that improve the numerical robustness and flexibility of the code.

No matter which version of the code you use, you are strongly encouraged to use its arbitrary-precision features when computing your final results, as experience has shown that numerical problems are surpisingly likely to occur, in any higher-order code written in any computer language. Iteratively solving huge systems of linear equations is inherently going to lose precision at every iterative step, not to mention potential condition number problems, near-singularities, and numerical inaccuracies that are increasingly likely to arise as one considers more and more curved models. Solutions computed using the arbitrary-precision features of the code are much more likely to be free of these kinds of numerical problems.

All versions of the code above require Mathematica version 5 or later—while it is not difficult to modify the routines to work in earlier versions of Mathematica, those earlier version are much slower at solving systems of linear equations. Perturbation AIM 2.5, 2.6, and 2.7 were designed to work with Mathematica 6 and 7 and may not work with earlier versions. Perturbation AIM 2.8 has only been tested on Mathematica 7. Note that Mathematica 7 introduces parallel computing procedures that I will try to incorporate into the Perturbation AIM code at some point in the near future.

Using the Code

To run the code, you must first start up Mathematica, and then type from a command line (including the “<<”)

<<perturbationAIM.m

If you are using Windows rather than Unix, you may first have to tell Mathematica which directory or folder you are working in. For example, on my laptop, the command:

SetDirectory["C:\Documents and Settings\Eric T. Swanson\My Documents\papers\pert"]

changes the current working directory to the one in quotes. Alternatively, you can use Mathematica’s pulldown menus to input the perturbationAIM.m file. This will define all of the functions you need to solve rational expectations models to arbitrarily high order. Of  course, you need to define the equations and parameters of your model. The following example files (particularly the first one) should serve as useful templates and include many comments, tips, and hints:

To run the example file etsSG.m, simply type:

<<etsSG.m

at a Mathematica command prompt (or use the pulldown menus). In order to play with the output, substitute in numerical values for the moments of your shocks, and so on, you will probably want to become familiar with basic Mathematica commands and syntax. The Mathematica Documentation Center provides complete and searchable documentation for all of Mathematica.

Some Caveats

Using numerical computation methods is always a bit of an art, and the Perturbation AIM code is no exception. With standard, machine-precision computation, numerical difficulties can easily arise, and in fact our experience has been that they occur with a surprising degree of frequency—for example, even in the simple stochastic growth model above, a few of the machine-precision coefficients in the third-order expansion are incorrect, as you can verify by checking them against their arbitrary-precision counterparts. Moreover, our experience has been that including representative-agent welfare as one of the equations of the model (a fairly common application of second-order methods) has a tendency to increase numerical problems, partly because steady-state welfare is typically 100 times larger than the steady-state values of the other variables in the model.

The simplest, most reliable way to avoid numerical inaccuracies in your answers is to use the aribtrary-precision feature of the code, and to experiment a little with the level of precision (using the AIMPrecision environment variable) to verify that the solution does not change meaningfully. The user should also experiment with varying the AIMZeroTol environment variable, which is the smallest numerical value that the code treats as being economically meaningful, as opposed to being simply an artifact of finite-precision arithmetic. The default value of AIMZeroTol is 10^-10. Note that, as with all machine-precision or finite-precision computations, a smaller value for this number does not necessarily lead to more accurate results, and can even lead to less accurate results.

Thus, users of Perturbation AIM are strongly encouraged to compute their final results using the arbitrary-precision feature of the code. While this comes at some cost in terms of computation time, the cost is small compared to the bottleneck in the algorithm, which is the evaluation of large symbolic derivative expressions at steady state.