pyFRB/US User Guide

This user guide describes basic usage of the FRB/US model with the pyFRB/US package. Users can observe the documented operations in action by running the bundled demo programs.

Installing

The pyFRB/US package can be installed by running pip3 install -e . or pip install -e . from the root directory of this package. Python dependencies should be automatically installed and are listed in setup.py. If scikit-umfpack is installed, it will be used to speed up solution times under certain conditions, but it is not required. pyFRB/US is only compatible with Python version 3.9 or greater. Compatibility has been tested up to Python 3.14 (inclusive), but cannot be guaranteed for later Python versions.

Basic Usage

Initialization

To load the model, create a new Frbus object using the constructor, pointing it at the provided model XML:

from pyfrbus.frbus import Frbus

frbus = Frbus("model.xml")

The equations of the model are loaded from the XML file’s python_equation tags, and each equation is modified by adding a tracking residual named with the suffix _trac.

By default, the backward-looking VAR expectations version of FRB/US is loaded. The constructor takes an optional argument mce which allows users to select which version of the forward-looking rational expectations model to load. Valid MCE types are all, mcap, wp, and mcap+wp.

The FRB/US dataset can be loaded from the data_only_package with the load_data function:

from pyfrbus.load_data import load_data

data = load_data("LONGBASE.TXT")

The model requires that series exist in the input DataFrame for all endogenous and exogenous variables.

From here, we can initialize the tracking residuals for the model using the init_trac function:

start = "2019Q4"
end = "2030Q4"
baseline_with_adds = frbus.init_trac(start, end, data)

init_trac returns a new DataFrame where the _trac variables have shocks filled in such that the model will solve to the specified baseline values.

Solving

Model simulations are run by calling the solve function:

sim = frbus.solve(start, end, baseline_with_adds)

This gives a new DataFrame where the series for endogenous variables take on values consistent with the exogenous series, lags, and model equations over the specified period. The returned DataFrame is identical to the original outside of the simulation period. The solver accepts some options, which are detailed in the Frbus Object API, and solver algorithms are described below in the Technical Appendix.

Model-Consistent Expectations

The solver will automatically detect whether the model is forward-looking and will switch to an algorithm suited for solving such models:

# MCE alert!
frbus = Frbus("model.xml", mce="mcap+wp")
...
# Returns a solution with forward-looking expectations
sim = frbus.solve(start, end, baseline_with_adds)

No further user configuration is required; additional options are detailed in the Frbus Object API. If you create a forward-looking model by appending equations with leads to a backward-looking model, this will be detected automatically.

The MCE type passed to the mce argument determines which expectations equations are replaced by their forward-looking versions. This is determined by the mce_group tag present in the model XML. mce="mcap" means that only equations with an mce_group of mcap will be replaced. mce="wp" means that only equations with an mce_group of mcwp will be replaced. mce="mcap+wp" will load both the mcap and mcwp variables, and mce="all" will replace all variables that have an MCE equation (mcap, mcwp, and mcother).

mcontrol

The mcontrol algorithm is a trajectory-matching control procedure which adjusts the value of instrument variables such that target variables are forced to specified trajectories, as mediated by the model’s dynamics.

mcontrol takes three lists of model variables as input: targ is the list of endogenous model variables to be forced; traj is the list of trajectories to force targ variables to; and inst is the list of exogenous model variables that will be moved freely to control the targ variables.

All three are specified as strings, which must take on the same names as those variables in the input DataFrame. targ and traj are specified in the same order - the first targ variable is forced to the first traj trajectory, and so on:

# Assemble trajectories
baseline_with_adds.loc[simstart:simend, "xgdp_t"] = numpy.nan
baseline_with_adds.loc[start_force:simend, "xgdp_t"] = (1.001 * data.loc[start_force:simend, "xgdp"])
baseline_with_adds.loc[simstart:simend, "ynicpn_t"] = numpy.nan
baseline_with_adds.loc[start_force:simend, "ynicpn_t"] = (1.001 * data.loc[start_force:simend, "ynicpn"])

# Run mcontrol
targ = ["xgdp", "ynicpn"]
traj = ["xgdp_t", "ynicpn_t"]
inst = ["ecnia_aerr", "ynicpn_aerr"]
sim = frbus.mcontrol(simstart, simend, baseline_with_adds, targ, traj, inst)

User must set trajectories to numpy.nan when not in use. When nan is given, that trajectory is disabled for the specified quarter and the model’s endogenous behavior takes over. Trajectory matching can be done discontinuously by using nan during the simulation period.

Stochastic Simulations

The stochsim procedure performs a stochastic simulation by applying sequences of shocks to the model, drawn randomly from historical residuals.

The procedure begins by randomly drawing sequences of quarters from the residual period. In each quarter of a single replication, all stochastic variables (specified with a stochastic_type tag in the model) have a shock applied from a particular quarter in the residual period.

Before using this procedure, the DataFrame input_data must have residuals computed over both history and the simulation period with init_trac.

stochsim is run by specifing the simulation period, residual period, and the number of replications:

simstart = "2019Q4"
simend = "2030Q4"
residstart = "1975q1"
residend = "2018q4"
baseline_with_adds = frbus.init_trac(residstart, simend, data)
nrepl = 5000
frbus.stochsim(nrepl, baseline_with_adds, simstart, simend, residstart, residend)

Output is returned as a list of DataFrames, each corresponding to a single replication.

Other Features

Appending or Replacing Equations

Model equations can be added or modified at run-time with the append_replace command.

append_replace takes a dictionary, mapping endogenous variable names to their corresponding equations. Equations with the same name will be replaced in the model, while equations with new names will be appended:

new_eqs = {"anngr" : "anngr = (exp(xgdp)/exp(xgdp(-4)) - 1) * 100"}
frbus.append_replace(new_eqs)

If your equations contain new endogenous or exogenous variables, remember that those variables will each need a corresponding series in the DataFrame.

Exogenizing Variables

Endogenous variables can be exogenized with the use of the exogenize method:

frbus.exogenize(["rff", "rg5", "rg10", "rg30"])

During simulations, all specified variables will have their equations disabled and will instead take on the values given in the DataFrame.

Note that all all exogenized variables must be passed in at the same time; calling exogenize twice in a row will result in the loss of earlier exogenize settings:

# This set of modifications is dropped...
frbus.exogenize(["rff", "rg5"])
# when you make another change!
frbus.exogenize(["rg10"])

The correct usage requires all changes to be passed in as a single list.

Tips & Tricks

  • Many pyFRB/US functions takes a start and end date, which can be passed in either as strings or as pandas Period objects. Using pandas Period objects allows for easier date-shift operations in your programs - e.g. Period("2019Q4")+1 gives Period("2020Q1").

  • Simulations run under VAR expectations may fail to solve with the default algorithm - switching to the “newton” solver may improve convergence.

  • The sim_lib module has functions sim_plot and stochsim_plot, which can be used to visualize simulation and stochsim results. Standard usage can be seen in the demo programs.

  • An options dictionary can be passed to solve, stochsim, and mcontrol - documentation of usable options is listed under solve in the Frbus Object API.

  • If plots do not display, you may need to configure the default backend used by matplotlib.

Windows Installation

On Windows, it is suggested to install and run pyFRB/US via the Windows Subsystem for Linux (WSL). You must be on a recent version of Windows 10 or Windows 11 to use pyFRB/US via the WSL. First, install the WSL by following these instructions.

After opening a Linux terminal in the WSL, you can install the pyFRB/US package using pip (or pip3) and run example programs with python3.

Note that for most systems, WSL is not configured for graphical display by default so plots will not render if running from the command line. However, you can run pyFRB/US programs from Microsoft’s Visual Studio Code (VSCode) editor which will correctly display plots. Install VSCode on your Windows computer, then connect WSL with VSCode. Now you can open pyFRB/US demos in VSCode and run them with “Run Current File in Interactive Window” under the right-click menu.

Technical Appendix

This appendix includes technical information on the implementation of parts of the pyFRB/US platform. The solver accepts user-configurable options, some of which are detailed below. A full list is available in the Frbus Object API.

Jacobian and Supported Functions

In pyFRB/US, the model Jacobian is computed symbolically using SymPy and SymEngine; there is no option for a numerically-approximated Jacobian. As a result, the platform supports a limited number of functions to be used in model equations.

log, exp, and abs are supported. min and max are supported and take varargs.

Additionally, you can use the function ind_ltezero as an indicator function. ind_ltezero takes a single argument. If the input is less than or equal to 0, the function returns 1, and otherwise returns 0.

VAR Solver Algorithm

The default solver for the VAR expectations model uses the dense Powell hybrid method from scipy.optimize.root. The model itself is decomposed into purely simultaneous and non-simultaneous blocks so that each call to the solver is as small as possible. Each time period in the simulation range is solved iteratively in sequence. The initial starting “guess” for the algorithm is taken from that quarter in the input data.

The VAR expectations model can also be used with the sparse Newton and trust-region algorithms detailed under MCE Solver Algorithm by passing in the options dictionary. Users should try these methods if the model fails to solve with the default algorithm. In this case the model is still solved by iterating over each time period.

MCE Solver Algorithm

The MCE solver assembles the model into a stacked-time system. Model equations are duplicated for each period in the simulation range and terminal values are substituted for forward-looking variables past the end of that range. The model is then solved for all periods simultaneously.

The terminal conditions used for the stacked-time system are drawn from the input data at the given number of quarters past the end of the simulation range. For example, an equation containing the term xgdp(2) will, in the final period of the simulation, take its value from the series xgdp two quarters past the end of the simulation range. The initial starting “guess” for the algorithm is taken from each quarter in the input data.

The MCE solver provides two solution algorithms: sparse Newton and sparse trust-region. The default is Newton. Both algorithms use the sparse solver from either SuperLU or UMFPACK, and both algorithms use a diagonal preconditioner to standardize scaling across model variables. The preconditioner can be disabled by passing False as the precond option in the options dictionary.

Newton

The Newton algorithm includes a damping procedure to avoid domain violations in the model; if the Newton step would produce nan in a model equation, the stepsize is scaled down until the domain is satisfied. A scaling factor of less than 1e-5 will result in a ConvergenceError and a failed solution.

The solver terminates when the 2-norm of the next step is less than xtol, which can be specified by the user and defaults to 1e-4. If the solver iterates for more than maxiter iterations without reaching the tolerance level, the procedure will return a ConvergenceError. Difficult problems may require more iterations; maxiter defaults to 50 and can be specified by the user.

If xtol is satisfied but the remaining residual is larger than rtol, the procedure will return a ConvergenceError. A large residual with a small stepsize indicates that the algorithm has converged to a “point of no interest” that does not actually satisfy the model equations. rtol can be specified by the user and defaults to 5e-4.

Jacobian re-use is enabled by default – the algorithm will attempt to use the same Jacobian for multiple steps in a row to reduce computation time. If the Newton step based on a stale Jacobian does not sufficiently improve the residual, the Jacobian will be recomputed. This behavior can be disabled with the option force_recompute, which will force the Jacobian to be computed at each step. Additionally, the Jacobian equations are not checked for domain violations by default. If an invalid Jacobian is encountered, boundary checking can be enabled with the option check_jac.

Trust-Region

The trust-region algorithm is the standard Powell Newton dogleg method, which can be seen in Numerical Optimization by Nocedal and Wright. This procedure limits the iteration step to be within a certain radius where the model is known to be approximately linear. The size of this region is adjusted dynamically based on how close the actual improvement at each step is to the theoretical improvement in a linearized system. Additionally, the procedure takes “dogleg” steps where the step direction is taken to be the optimal linear combination of the Newton step and the gradient descent step which maximizes the residual improvement.

The trust-region solver takes an option trust_radius, which defaults to 1000000. trust_radius specifies the maximum possible size of the trust region; the starting size of the trust region is taken as trust_radius / 2. The trust-region algorithm also takes maxiter, xtol, and rtol options as specified in the Newton section.