Tutorial: Using the Perturbation Helper

This document uses the object-oriented programming syntax of MATLAB. Necessary things are introduced, but refer to the MATLAB documentation more details.

Introduction

One should never differentiate large economic models by hand. Most economic toolboxes like DYNARE and RISE come with automation tools for differentiation. For example, DYNARE uses differentiation rules written in C++ to generate MATLAB codes to differentiate model expressions. However, the m-file generation methods are impossible to generalize to bigger problems where we need to differentiate (sparse) matrix-vector multiplications. Therefore, I wrote a MATLAB code to automate the differentiation MATLABAutoDiff . The ADHelper is written to help with data organization so that one would not need to deal directly with “raw” indices, and to help with organization of ordering as required by krylov_reduction and schur_solver as introduced in [Ahn et al., 2017] . Though this class is written mostly for heterogeneous agent models in continuous time, we will solve the 3-equation New Keynesian model in this tutorial as people are more familiar with that setup. [1]

Subclassing the ADHelper

Before we can get started with the main section, we need a slight detour with object-oriented programming syntax. ADHelper is a class that is intended to be subclassed, and I routinely subclass the object that represent “economy” with this class. In the MATLAB object-oriented programming syntax, this is done by

classdef economy < ADHelper
    properties
    end

    methods
    end
end

in a economy.m file. A class constructor is run when an object is instantiated. Hence, we can define all the relevant parameters with the constructor. For example, we can define set parameters taylor_on_inflation to 1.5 and taylor_on_gdp 0 by

classdef economy < ADHelper
    properties
    end
    methods
        function obj = economy()
            obj.taylor_on_inflation = 1.8;
            obj.taylor_on_gdp = 0.63;
        end
    end
end

However, if you run this, this will lead to an error because the parameters are not declared as properties. Hence, one needs to define them as properties

classdef economy < ADHelper
    properties
        taylor_on_gdp
        taylor_on_inflation
    end
    methods
        function obj = economy()
            obj.taylor_on_inflation = 1.5;
            obj.taylor_on_gdp = 1.5;
        end
    end
end

As the economy is a subclass of the ADHelper, it will have many properties and functions already available to it. Refer to the technical documentation (Perturbation Method ) for all available properties and methods. Lastly, one instantiates an instance of the class by

econ = economy;

Preparing Variables

Unlike the other differentiation methods as in DYNARE, the derivatives are computed at an given point. Hence, the point to differentiate needs to be decided before differentiation. Under the perturbation method in economics, that point is the steady-state. Hence, we need to start by defining steady-state values. This is done in ADHelper class using var2ss . var2ss is a struct that one need to fill with variable names and steady state values

econ.var2ss.variable_name = steady_state_value;

For example, if the steady-state value of inflation is zero, we set

econ.var2ss.inflation = 0;

For the DYNARE crowd, this can be considered as implementing the var and initval portion of the model file at the same time. One needs to do these for all endogenous and exogenous variables. Hence, for example, for the monetary policy shock needs to be coded as

econ.var2ss.mp_shock = 0;

In ADHelper, the exogenous variables are denoted with var2isshock . The default is false , so one only need to set the flag for the exogenous variables to true ,i.e.,

econ.var2isshock.mp_shock = true;

Again, for the DYNARE crowd, this is equivalent to declaring mp_shock as an varexo . Now, in DYNARE, there is no special treatment for predetermined vs decision variables. This is because there are algorithms, e.g., QZ-decomposition, without requiring the predeterminedness. However, for large scale problems arising with heterogeneous agent models, QZ-decomposition scales worse than Schur-decomposition, and for problems that require model-reductions (see [Ahn et al., 2017] ), it is not possible to use QZ-decomposition at all. Hence, we require labels for decision variables for large scale problems. [2] The decision variables are labeled with var2ischoice . For example,

econ.var2ischoice.consumption = true;

Just the way one works through var , varexo , and initval in DYNARE. One needs to do this for all variables (and do not have spurious variables in var2ss ). Just in case there is a DYNARE person wondering, since model equations will not go through a parser, one do not need to declare parameters . In fact, I usually just set it as a class property and use it everywhere as I have set with taylor_on_inflation above. One just need to do this before solve_perturbation is called. I usually create find_steady_state method within the subclass. In this example, I will put this in the constructor. Putting everything up to this point so far together, we have

classdef simple_nk_economy < ADHelper
    properties
        tau
        kappa
        taylor_on_inflation
        taylor_on_gdp
        beta_util
        inflation_ss
        interest_ss
        persistence_fiscal
        persistence_interest
        persistence_tfp
        std_monetary_policy_shock
        std_government_shock
        std_tfp_shock
    end

    methods
        function obj = simple_nk_economy()
            % Parameter Values
            obj.tau = 2.83;
            obj.kappa = 0.78;
            obj.taylor_on_inflation = 1.80;
            obj.taylor_on_gdp = 0.63;
            obj.beta_util = 1/(1+0.42/400);
            obj.inflation_ss = 1+3.3/400;
            obj.interest_ss = 1+0.52/400;
            obj.persistence_interest = 0.77;
            obj.persistence_fiscal = 0.98;
            obj.persistence_tfp = 0.88;
            obj.std_monetary_policy_shock = 0.0022;
            obj.std_government_shock = 0.0071;
            obj.std_tfp_shock = 0.0031;

            % Set Pre-determined Variables
            obj.var2ischoice.pi = true;
            obj.var2ischoice.y = true;

            % Set Exogenous Variables
            obj.var2isshock.g_shock = true;
            obj.var2isshock.z_shock = true;
            obj.var2isshock.mp_shock = true;

            % This problem is for discrete time
            obj.is_cont_time = false;
        end

        function find_steady_state(obj)
            % Set steady-state values
            obj.var2ss.y = 0;
            obj.var2ss.pi = 0;
            obj.var2ss.r = 0;
            obj.var2ss.g = 0;
            obj.var2ss.z = 0;

            obj.var2ss.g_shock = 0;
            obj.var2ss.z_shock = 0;
            obj.var2ss.mp_shock = 0;
        end
    end
end

Writing the setup_eqns Function

Though one can use other approaches that uses lower-level function calls, one can just rely on solve_perturbation for most applications. This function requires the model equations that define dynamics. This is set by write setup_eqns in the subclass, i.e.,

classdef economy < ADHelper
    ...
    methods
        ...
        function setup_eqns(obj)
            % Fill in here

            obj.eqns = getderivs(obj.eqns);
        end
    end
end

As ADHelper calls this functions after all the pre-processing is done. At this point, var2advars is populated and available to use which is a wrapper over MATLABAutoDiff so that one do not have to access things by raw indices. One can access endogenous variable, via obj.var2advars.variable_name . Also, at this point, the expression for \(\displaystyle \left(\frac{d}{dt} var\right)\) is populated by prepending d in front of the variable name for continuous time. In discrete time, one can consider d notation as one time period forward, i.e., (+1) in DYNARE. Given this notation, we can for example, write the IS curve

\[\widehat{y}_t = \mathbb{E}_t[\widehat{y}_{t+1}] - \frac{1}{\tau}\left(\widehat{R}_t - \mathbb{E}_t[\widehat{\pi}_{t+1}] - \mathbb{E}_t[\widehat{z}_{t+1}]\right) + \widehat{g}_t - \mathbb{E}_t[\widehat{g}_{t+1}]\]

as

obj.var2advars.y = obj.var2advars.dy - 1/obj.tau*(obj.var2advars.r ...
                   - obj.var2advars.dpi - obj.var2advars.dz) ...
                   + obj.var2advars.g - obj.var2advars.dg;

Note that setup_eqns is executed completely within MATLAB without any parsing, so one can define intermediate variables instead of typing thing over and over again (without special syntax like # ). To reduce typing, one can also write

y = obj.var2advars.y;
y_1 = obj.var2advars.dy;
r = obj.var2advars.r;
pi_1 = obj.var2advars.dpi;
z_1 = obj.var2advars.dz;
g = obj.var2advars.g;
g_1 = obj.var2advars.dg;

y = y_1 - 1/obj.tau*(r - pi_1 - z_1) + g - g_1;

The latter is probably easier to understand as there is less clutter in equations. Again, setup_eqns is executed within MATLAB, so one is free to do anything without the restrictions of the model languages through a parser. Not heeding my own advise of reducing clutter, the function can be finished for the three equation model as

classdef simple_nk_economy < ADHelper
    ...
    methods
        ...
        function setup_eqns(obj)
            % IS Curve
            is_curve = -obj.var2advars.y + obj.var2advars.dy - 1/obj.tau*(obj.var2advars.r ...
                       - obj.var2advars.dpi - obj.var2advars.dz) + obj.var2advars.g ...
                       - obj.var2advars.dg;

            % Philips Curve
            phil = -obj.var2advars.pi + obj.beta_util*obj.var2advars.dpi ...
                + obj.kappa*(obj.var2advars.y - obj.var2advars.g);

            % Taylor Rule
            taylor = -obj.var2advars.dr + obj.persistence_interest*obj.var2advars.r ...
                + (1-obj.persistence_interest)*(obj.taylor_on_inflation*obj.var2advars.dpi ...
                + obj.taylor_on_gdp*(obj.var2advars.dy - obj.var2advars.dg)) ...
                + obj.std_monetary_policy_shock*obj.var2advars.mp_shock;

            % Fiscal Policy
            gov_process = obj.var2advars.dg - obj.persistence_fiscal*obj.var2advars.g ...
                - obj.std_government_shock*obj.var2advars.g_shock;

            % TFP Shock
            tfp_process = obj.var2advars.dz - obj.persistence_tfp*obj.var2advars.z ...
                - obj.std_tfp_shock*obj.var2advars.z_shock;

            % Combine Equations
            obj.eqns = getderivs([is_curve; phil; taylor; gov_process; tfp_process]);
        end
    end
end

Every other parts of the lines will be similar to the syntax of DYNARE model language by replacing (+1) with d except for the last line. var2advars is an myAD object, and eqns need to be set as the derivatives. We can get the derivatives by using getderivs function of MATLABAutoDiff , and the resulting is the last line. For more explanation of the MATLABAutoDiff , refer to the github page and the contained PDF documentation.

Solving the Model

Writing the setup_eqns is the only involved part, and we are ready to solve the model with the linear perturbation method. First, we instantiate the simple_nk_economy class and set steady-state by

economy = simple_nk_economy;
economy.find_steady_state();

Then, we can just call solve_perturbation function

economy.solve_perturbation();

Underneath the hood, the problem is solved with get_sspace_form and then using schur_solver . The dynamics matrix and the transfer function is saved as G1 and impact following the same variable names as that from Gensys. We can also just compute the impulse-response functions by using simulate_irf . At this stage, one can using shock2loc and var2loc to avoid raw-indices. For example, to plot the irf of GDP on the monetary policy shock, we can do

simulate_time = 20;
irf_mp_shock = economy.simulate_irf(simulate_time, simulate_time, economy.shock2loc.mp_shock);
plot(irf_mp_shock(economy.var2loc.y, :));

Warning

Lastly, though one can use the underlying codes with discrete time problems and I have checked the output of this exercise with those from the RISE toolbox , I actually recommend one to continue to use the DYNARE or RISE for small scale discrete time problems as those toolboxes were written for those applications and have been more. [3] This was just a tutorial to show the syntax of the ADHelper class. Of course, for large-scale problems like heterogeneous agent models, one can’t use DYNARE or RISE.

[1]	For small scale problems in discrete time, one can stick with pre-existing toolboxes like DYNARE and RISE as they have a lot more features.

[2]	In fact, for small scale problem, one can just set set var2ischoice randomly, and just use other solvers using get_gensys_form , but we will not elaborate on that here.

[3]	As much as people can try to avoid bugs, it is impossible to be 100% bug free without a lot of market test. In this sense, DYNARE has been through a extensive market test from its bigger user base.