\documentclass[article,nojss]{jss}

\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{float}
\usepackage{enumitem}

\hyphenation{bimets}

\author{Andrea Luciani\\Bank of Italy\thanks{Disclaimer: \emph{The views and opinions expressed in these pages are those of the author and do not necessarily reflect the official policy or position of the Bank of Italy. Examples of analysis performed within these pages are only examples. They should not be utilized in real-world analytic products as they are based only on very limited and dated open source information. Assumptions made within the analysis are not reflective of the position of the Bank of Italy.}}}
\Plainauthor{Andrea Luciani}

\title{The U.S. Federal Reserve quarterly model in \proglang{R} with \pkg{bimets}}
\Plaintitle{The U.S. Federal Reserve quarterly model in R with bimets}
\Shorttitle{The U.S. Federal Reserve quarterly model in R with bimets}

\Abstract{

The US Federal Reserve's econometric model for the US economy (i.e., FRB/US) is publicly available at \href{https://www.federalreserve.gov/econres/us-models-about.htm}{federalreserve.gov}. The website states, \emph{"FRB/US is a large-scale estimated general equilibrium model of the US economy that was developed at the Federal Reserve Board, where it has been in use since 1996 for forecasting, analysis of policy options, and research projects."}

\bigskip

FRB/US is a quarterly model with hundreds of equations and variables. The model definition and time series data are available for download on the Federal Reserve website, as is the source code, which allows users to perform several econometric exercises. 

\bigskip

However, the Federal Reserve publicly distributes source codes only for \proglang{EViews}\textsuperscript{\textregistered} and \proglang{python}. 

\bigskip

\pkg{bimets} is a software framework developed by using \proglang{R} language and designed for time series analysis and econometric modeling. In these pages, we will show how to use \pkg{bimets} capabilities to load the FRB/US model and perform in \proglang{R} the same econometric exercises provided by the Federal Reserve. For each exercise, we will compare the numerical results in \proglang{R} and \proglang{python}. 
}

\Keywords{\proglang{R}, bimets, system of simultaneous equations, 
  Federal Reserve quarterly model, FRB/US,
  model simulation, forecasting, endogenous targeting, 
  stochastic simulation, rational expectations}
\Plainkeywords{R, bimets, system of simultaneous equations, 
  Federal Reserve quarterly model, FRB/US,
  model simulation, forecasting, endogenous targeting, 
  stochastic simulation, rational expectations}
  
\Address{
   Andrea Luciani\\
   Bank of Italy\\
   Directorate General for Economics, Statistics and Research\\
   Via Nazionale, 91\\
   00184, Rome - Italy\\
   E-mail: \email{andrea.luciani@bancaditalia.it}\\
}

  
\begin{document}
\SweaveOpts{concordance=TRUE}
%\VignetteIndexEntry{US Federal Reserve quarterly model (FRB/US) in R with bimets}
%\VignetteKeywords{R, system of simultaneous equations,
% estimation, ols, instrumental variables, error autocorrelation,
% pdl, simulation, multipliers, renormalization, forecasting, rational expectations}
%\VignettePackage{bimets}
<<echo=FALSE>>=
options( prompt = "R> ", continue = "   " )
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The FRB/US model}
The Federal Reserve website states, \emph{"FRB/US is a large-scale estimated general equilibrium model of the US economy that was developed at the Federal Reserve Board, where it has been in use since 1996 for forecasting, analysis of policy options, and research projects. ... Compared with DSGE models, however, FRB/US applies optimization theory more flexibly, which permits its equations to better capture patterns in historical data and facilitates modeling the economy in greater detail. ... A distinctive feature of FRB/US is its ability to switch between alternative assumptions about how economic agents form expectations. Under the VAR-based option, expectations are derived from the average historical dynamics of the economy as manifested in the predictions of estimated VAR models. Under model-consistent (MC), agents are assumed to form accurate expectations of future outcomes as generated by simulations of FRB/US itself."}

\bigskip

FRB/US is a quarterly model, and counts 284 equations and 365 variables (Feb. 2024 version). The XML model definition is available for download on the Federal Reserve website, and contains, for each endogenous variable, the following information: the variable name, the variable definition with a short description, the economic sector the variable belongs to, the related equation in both \proglang{EViews}\textsuperscript{\textregistered} and \proglang{python} format, coefficients and exogenous variables involved in the equation.

\bigskip

64 endogenous variables are marked as stochastic and, during the stochastic simulation exercise (see section \ref{ssec:5}), will be transformed by applying sequences of shocks as drawn randomly from their historical residuals. 

\bigskip

14 endogenous variables belong to the MCE group (i.e., Model-Consistent Expectations) and have an alternative equation that contains forward-looking references (see section \ref{ssec:2}).

\bigskip

Finally, at the end of the XML model definition, users can find additional information on economic sectors and exogenous variables involved in the model definition.
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The pyfrbus python package}

\pkg{pyfrbus} is a \proglang{python}-based simulation platform for the Federal Reserve Board's FRB/US model, which depends on the \code{SuiteSparse} code, a widely used set of sparse-matrix-related optimized algorithms.

\bigskip

As stated in the reference manual, to load the model in \proglang{python}, create a new \code{Frbus} object using the constructor, pointing it at the provided model XML:

\begin{footnotesize}
\code{
from pyfrbus.frbus import Frbus \\
frbus = Frbus("model.xml")
}
\end{footnotesize}

\bigskip

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

\bigskip

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

\bigskip

The CSV FRB/US dataset \code{LONGBASE.TXT} can be loaded from the data-only-package with the \code{load_data} function:

\begin{footnotesize}
\code{
from pyfrbus.load_data import load_data \\
data = load_data("LONGBASE.TXT")
}
\end{footnotesize}

\bigskip

The model requires that series exist in the input \code{pandas} DataFrame for all endogenous and exogenous variables. From here, we can initialize the tracking residuals for the model using the \code{init_trac} function:

\begin{footnotesize}
\code{
start = "2019Q4" \\
end = "2030Q4" \\
baseline_with_adds = frbus.init_trac(start, end, data) 
}
\end{footnotesize}


\bigskip

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

\bigskip

Model simulations are run by calling the \code{solve} function (see section \ref{ssec:1}):

\begin{footnotesize}
\code{
sim = frbus.solve(start, end, baseline_with_adds)
}
\end{footnotesize}

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. 

\bigskip

The solver will automatically detect whether the model is forward-looking and will switch to an algorithm suited for solving such models (see section \ref{ssec:2}):

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

\bigskip

The \code{mcontrol} algorithm is a trajectory-matching control procedure (i.e., endogenous targeting) which adjusts the value of instrument variables such that target variables are forced to specified trajectories, as mediated by the model's dynamics.

\bigskip

\code{mcontrol} takes three lists of model variables as input: \code{targ} is the list of endogenous model variables to be forced; \code{traj} is the list of trajectories to force \code{targ} variables to; and \code{inst} is the list of exogenous model variables that will be moved freely to control the \code{targ} variables (see section \ref{ssec:4}).

\bigskip

The \code{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 \code{stochastic_type} tag in the model) have a shock applied from a particular quarter in the residual period (see section \ref{ssec:5}).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Moving to R}


\pkg{bimets} is a software framework developed by using \proglang{R} language and designed for time series analysis and econometric modeling. More details about the package are available at the \href{https://cran.r-project.org/package=bimets/vignettes/bimets.pdf}{"Getting started with bimets"} vignette. FRB/US model definition has been translated into a \pkg{bimets} compliant syntax and is available to \proglang{R} users as the \code{FRB__MODEL} dataset, which contains the textual definition of the model, using \pkg{bimets} \code{MDL} compliant syntax, i.e., Model Description Language:

\begin{footnotesize}
<<>>=
#load bimets
library(bimets)
@
<<echo=FALSE>>=
sim_plot <- function(model,TSRANGE,plotidx)
{
  #define layout 
  par(mfrow = c(2, 2))
  
  if (plotidx==1)
    TSRANGE=c(normalizeYP(TSRANGE[1:2]-c(0,6),4),TSRANGE[3:4])
  
  if (plotidx==2)
    TSRANGE=c(normalizeYP(TSRANGE[1:2]-c(0,2),4),TSRANGE[3:4])
  
  if (plotidx==3)
    TSRANGE=c(normalizeYP(TSRANGE[1:2]-c(0,6),4),TSRANGE[3:4])
  
  if (plotidx==4)
    TSRANGE=c(normalizeYP(TSRANGE[1:2]-c(0,2),4),TSRANGE[3:4])
  
  xlim=c(TSRANGE[1]+(TSRANGE[2]-1)/4,TSRANGE[3]+(TSRANGE[4]-1)/4)
  
  #plot1
  series1=TSDELTAP(model$simulation$xgdp,4)
  series2=TSPROJECT(TSDELTAP(model$modelData$xgdp,4),
                    TSRANGE=TSRANGE)
  min=min(series1,series2)
  max=max(series1,series2)
  range=max-min
  plot(series2,
       font.main=1,
       col='blue',
       main='Real GDP Growth, Quarterly Annualized',
       ylab='Percent',
       xlab=NULL,
       ylim=c(min-0.05*range,max+0.05*range),
       xlim=xlim,
       yaxt='n',
       xaxt='n')
  axis(side=2,las=2)
  axis(side=1,at=time(series2),labels=as.character(as.yearqtr(time(series2))))
  lines(series1,
        lty='dashed',
        col='red')
   
  #plot2
  series1=TSPROJECT(model$simulation$lur,
                 TSRANGE=TSRANGE)
  series2=TSPROJECT(model$modelData$lur,
                  TSRANGE=TSRANGE)
  min=min(series1,series2)
  max=max(series1,series2)
  range=max-min
  plot(series2,
       font.main=1,
       col='blue',
       main='Unemployment Rate',
       ylab='Percent',
       xlab=NULL,
       ylim=c(min-0.05*range,max+0.05*range),
       xlim=xlim,
       yaxt='n',
       xaxt='n')
  axis(side=2,las=2)
  axis(side=1,at=time(series2),labels=as.character(as.yearqtr(time(series2))))
  lines(series1,
        lty='dashed',
        col='red')
  
  #plot3
  series1=TSDELTAP(model$simulation$pcxfe,4) 
  series2=TSPROJECT(TSDELTAP(model$modelData$pcxfe,4),
                  TSRANGE=TSRANGE) 
  min=min(series1,series2)
  max=max(series1,series2)
  range=max-min
  plot(series2,
       font.main=1,
       col='blue',
       main='Core PCE Inflation, Quarterly Annualized',
       ylab=ifelse(plotidx==2,'','Percent'),
       xlab=NULL,
       ylim=c(min-0.05*range,max+0.05*range),
       xlim=xlim,
       yaxt='n',
       xaxt='n')
  axis(side=2,las=2)
  axis(side=1,at=time(series2),labels=as.character(as.yearqtr(time(series2))))
  lines(series1,
        lty='dashed',
        col='red')
  
  #plot4
  series1=TSPROJECT(model$simulation$rff,
                 TSRANGE=TSRANGE)
  series2=TSPROJECT(model$modelData$rff,
                  TSRANGE=TSRANGE)
  min=min(series1,series2)
  max=max(series1,series2)
  range=max-min
  plot(series2,
       font.main=1,
       col='blue',
       main='Federal Funds Rate',
       ylab='Percent',
       xlab=NULL,
       ylim=c(min-0.05*range,max+0.05*range),
       xlim=xlim,
       yaxt='n',
       xaxt='n')
  axis(side=2,las=2)
  axis(side=1,at=time(series2),labels=as.character(as.yearqtr(time(series2))))
  lines(series1,
        lty='dashed',
        col='red')
}

stochsim_plot <- function(model,TSRANGE)
{
  TSRANGE[1:2]=normalizeYP(TSRANGE[1:2]-c(0,6),4)
    
  #define layout 
  par(mfrow = c(2, 2))
  xlim=c(TSRANGE[1]+(TSRANGE[2]-1)/4,TSRANGE[3]+(TSRANGE[4]-1)/4)
 
  #init stuff
  seriesName='xgdp'
  baseStochMatrix=model$simulation_MM[[seriesName]]
  repl=dim(baseStochMatrix)[2]-1
  simPeriods=dim(baseStochMatrix)[1]
  simStart=start(model$stochastic_simulation[[seriesName]]$mean)
  
  #plot1
  seriesName='xgdp'
  baseStochMatrix=model$simulation_MM[[seriesName]]
  baseStochMatrix=baseStochMatrix[,1+1:repl]
  historicalMatrix=matrix(TSPROJECT(model$modelData[[seriesName]],
                                    TSRANGE=c(normalizeYP(c(simStart[1],simStart[2]-4),4),
                                              normalizeYP(c(simStart[1],simStart[2]-1),4))
                                    ),nrow=4,ncol=repl)
  fullMatrix=rbind(historicalMatrix,baseStochMatrix)
  deltap4StochMatrix=100*((fullMatrix[5:(simPeriods+4),]-fullMatrix[1:(simPeriods),])/
                            fullMatrix[1:(simPeriods),])
  rowMeans=rowMeans(deltap4StochMatrix)
  rowSd=c()
  for (idx in 1:(simPeriods)) rowSd[idx]=sd(deltap4StochMatrix[idx,])
  series1=TSPROJECT(TSDELTAP(model$modelData[[seriesName]],4),
                    TSRANGE=TSRANGE)
  #90% confidence 1.644854
  series2=TSMERGE(TSERIES(rowMeans-1.644854*rowSd,START=simStart,FREQ=4),series1)
  series3=TSMERGE(TSERIES(rowMeans+1.644854*rowSd,START=simStart,FREQ=4),series1)
  
  #70% confidence 1.036433
  series4=TSMERGE(TSERIES(rowMeans-1.036433*rowSd,START=simStart,FREQ=4),series1)
  series5=TSMERGE(TSERIES(rowMeans+1.036433*rowSd,START=simStart,FREQ=4),series1)
  min=min(series1,series2,series3,series4,series5)
  max=max(series1,series2,series3,series4,series5)
  range=max-min
  plot(series1,
       col='blue',
       main='Real GDP Growth, Quarterly Annualized',
       ylab='',
       xlab=NULL,
       ylim=c(min-0.05*range,max+0.05*range),
       xlim=xlim,
       yaxt='n',
       xaxt='n')
  axis(side=2,las=2)
  axis(side=1,at=time(series1),labels=as.character(as.yearqtr(time(series1))))
  polygon(c(time(series1),rev(time(series1))),c(series3,rev(series2)),col='gray',border = F)
  polygon(c(time(series1),rev(time(series1))),c(series5,rev(series4)),col='lightgray',border = F)
  lines(series2,
        lwd=1,
        col='black')
  lines(series3,
        lwd=1,
        col='black')
  lines(series4,
        lwd=2,
        col='darkgray')
  lines(series5,
        lwd=2,
        col='darkgray')
  lines(series1,
        lwd=1,
        col='blue')
  
  #plot2
  seriesName='lur'
  series1=TSPROJECT(model$modelData[[seriesName]],
                    TSRANGE=TSRANGE)
  #90% confidence 1.644854
  series2=TSMERGE(TSERIES(model$stochastic_simulation[[seriesName]]$mean-1.644854*model$stochastic_simulation[[seriesName]]$sd,START=simStart,FREQ=4),series1)
  series3=TSMERGE(TSERIES(model$stochastic_simulation[[seriesName]]$mean+1.644854*model$stochastic_simulation[[seriesName]]$sd,START=simStart,FREQ=4),series1)
  #70% confidence 1.036433
  series4=TSMERGE(TSERIES(model$stochastic_simulation[[seriesName]]$mean-1.036433*model$stochastic_simulation[[seriesName]]$sd,START=simStart,FREQ=4),series1)
  series5=TSMERGE(TSERIES(model$stochastic_simulation[[seriesName]]$mean+1.036433*model$stochastic_simulation[[seriesName]]$sd,START=simStart,FREQ=4),series1)
  min=min(series1,series2,series3,series4,series5)
  max=max(series1,series2,series3,series4,series5)
  range=max-min
  plot(series1,
       col='blue',
       main='Unemployment Rate',
       ylab='',
       xlab=NULL,
       ylim=c(min-0.05*range,max+0.05*range),
       xlim=xlim,
       yaxt='n',
       xaxt='n')
  axis(side=2,las=2)
  axis(side=1,at=time(series1),labels=as.character(as.yearqtr(time(series1))))
  polygon(c(time(series1),rev(time(series1))),c(series3,rev(series2)),col='gray',border = F)
  polygon(c(time(series1),rev(time(series1))),c(series5,rev(series4)),col='lightgray',border = F)
  lines(series2,
        lwd=1,
        col='black')
  lines(series3,
        lwd=1,
        col='black')
  lines(series4,
        lwd=2,
        col='darkgray')
  lines(series5,
        lwd=2,
        col='darkgray')
  lines(series1,
        lwd=1,
        col='blue')
  
  #plot3
  seriesName='pcxfe'
  baseStochMatrix=model$simulation_MM[[seriesName]]
  baseStochMatrix=baseStochMatrix[,1+1:repl]
  historicalMatrix=matrix(TSPROJECT(model$modelData[[seriesName]],
                                    TSRANGE=c(normalizeYP(c(simStart[1],simStart[2]-4),4),
                                              normalizeYP(c(simStart[1],simStart[2]-1),4))
  ),nrow=4,ncol=repl)
  fullMatrix=rbind(historicalMatrix,baseStochMatrix)
  deltap4StochMatrix=100*((fullMatrix[5:(simPeriods+4),]-fullMatrix[1:(simPeriods),])/
                            fullMatrix[1:(simPeriods),])
  rowMeans=rowMeans(deltap4StochMatrix)
  rowSd=c()
  for (idx in 1:(simPeriods)) rowSd[idx]=sd(deltap4StochMatrix[idx,])
  series1=TSPROJECT(TSDELTAP(model$modelData[[seriesName]],4),
                    TSRANGE=TSRANGE)
  #90% confidence 1.644854
  series2=TSMERGE(TSERIES(rowMeans-1.644854*rowSd,START=simStart,FREQ=4),series1)
  series3=TSMERGE(TSERIES(rowMeans+1.644854*rowSd,START=simStart,FREQ=4),series1)
  #70% confidence 1.036433
  series4=TSMERGE(TSERIES(rowMeans-1.036433*rowSd,START=simStart,FREQ=4),series1)
  series5=TSMERGE(TSERIES(rowMeans+1.036433*rowSd,START=simStart,FREQ=4),series1)
  min=min(series1,series2,series3,series4,series5)
  max=max(series1,series2,series3,series4,series5)
  range=max-min
  plot(series1,
       col='blue',
       main='Core PCE Inflation, Quarterly Annualized',
       ylab='',
       xlab=NULL,
       ylim=c(min-0.05*range,max+0.05*range),
       xlim=xlim,
       yaxt='n',
       xaxt='n')
  axis(side=2,las=2)
  axis(side=1,at=time(series1),labels=as.character(as.yearqtr(time(series1))))
  polygon(c(time(series1),rev(time(series1))),c(series3,rev(series2)),col='gray',border = F)
  polygon(c(time(series1),rev(time(series1))),c(series5,rev(series4)),col='lightgray',border = F)
  lines(series2,
        lwd=1,
        col='black')
  lines(series3,
        lwd=1,
        col='black')
  lines(series4,
        lwd=2,
        col='darkgray')
  lines(series5,
        lwd=2,
        col='darkgray')
  lines(series1,
        lwd=1,
        col='blue')
  
  #plot4
  seriesName='rff'
  series1=TSPROJECT(model$modelData[[seriesName]],
                    TSRANGE=TSRANGE)
  #90% confidence 1.644854
  series2=TSMERGE(TSERIES(pmax(0,model$stochastic_simulation[[seriesName]]$mean-1.644854*model$stochastic_simulation[[seriesName]]$sd),START=simStart,FREQ=4),series1)
  series3=TSMERGE(TSERIES(pmax(0,model$stochastic_simulation[[seriesName]]$mean+1.644854*model$stochastic_simulation[[seriesName]]$sd),START=simStart,FREQ=4),series1)
  #70% confidence 1.036433
  series4=TSMERGE(TSERIES(pmax(0,model$stochastic_simulation[[seriesName]]$mean-1.036433*model$stochastic_simulation[[seriesName]]$sd),START=simStart,FREQ=4),series1)
  series5=TSMERGE(TSERIES(pmax(0,model$stochastic_simulation[[seriesName]]$mean+1.036433*model$stochastic_simulation[[seriesName]]$sd),START=simStart,FREQ=4),series1)
  min=min(series1,series2,series3,series4,series5)
  max=max(series1,series2,series3,series4,series5)
  range=max-min
  plot(series1,
       col='blue',
       main='Federal Funds Rate',
       ylab='',
       xlab=NULL,
       ylim=c(min-0.05*range,max+0.05*range),
       xlim=xlim,
       yaxt='n',
       xaxt='n')
  axis(side=2,las=2)
  axis(side=1,at=time(series1),labels=as.character(as.yearqtr(time(series1))))
  polygon(c(time(series1),rev(time(series1))),c(series3,rev(series2)),col='gray',border = F)
  polygon(c(time(series1),rev(time(series1))),c(series5,rev(series4)),col='lightgray',border = F)
  lines(series2,
        lwd=1,
        col='black')
  lines(series3,
        lwd=1,
        col='black')
  lines(series4,
        lwd=2,
        col='darkgray')
  lines(series5,
        lwd=2,
        col='darkgray')
  lines(series1,
        lwd=1,
        col='blue')
  
  
}
@
<<>>=
#load FRB/US MDL definition
data(FRB__MODEL)
@
<<>>=
#print first equations in model definition
cat(substring(FRB__MODEL,1,1615))
@
\end{footnotesize}

\bigskip

\proglang{bimets} users can save the model definition into a text file, then inspect and modify it in \proglang{RStudio} or in any other text editor, by using the following commands:

\begin{footnotesize}
<<eval=FALSE>>=
#define file path
modelDefinitionFile <- file('~/FRB__MODEL.txt')

#save FRB definition in the text file
writeLines(FRB__MODEL,modelDefinitionFile)    

#close connection
close(modelDefinitionFile)
@
\end{footnotesize}

\bigskip

The \pkg{bimets} dataset \code{FRB__MCAP__WP__MODEL} contains the MCE version of the FRB/US model, wherein, as said before, 14 equations have been modified accounting for rational expectations (see section \ref{ssec:2}).

\bigskip

Original FRB/US time series are provided in a CSV text file containing quarterly data from 1962 up to 2173. These data are available to \proglang{R} users in the \code{LONGBASE} dataset, which includes the list of all endogenous and exogenous time series.

\begin{footnotesize}
<<>>=
#load FRB/US model data
data(LONGBASE)
@
<<>>=
#print GDP in 2022-2024
TABIT(LONGBASE$xgdp,TSRANGE = c(2022,1,2024,1))
@
\end{footnotesize}
FRB/US \code{MDL} definition can be translated into a \pkg{bimets} model by using the \code{LOAD_MODEL} function, as usual:
\begin{footnotesize}
<<>>=
#create the bimets model
model <- LOAD_MODEL(modelText = FRB__MODEL)
@
<<>>=
#print a sample of endogenous variables
model$vendog[1:10]
@
<<>>=
#print a sample of exogenous variables
model$vexog[1:10]
@
<<>>=
#print GDP equation
model$identities$xgdp$eqFull
@
\end{footnotesize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The econometric excercises}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Dynamic simulation in a monetary policy shock} \label{ssec:1}
 
The first econometric exercise proposed by the Federal Reserve is a dynamic simulation of the FRB/US model under a monetary policy shock. The simulation is operated from 2040-Q1 to 2045-Q4, after the \code{rffintay} time series, defined as \emph{"Value of eff. federal funds rate given by the inertial Taylor rule"}, is shocked by 100 base points in 2040-Q1.

\bigskip

\proglang{python} code follows: 

\bigskip

\begin{footnotesize}
\code{
import pandas \\ \\
from pyfrbus.frbus import Frbus \\
from pyfrbus.sim_lib import sim_plot \\
from pyfrbus.load_data import load_data \\ \\
\# Load data \\
data = load_data("LONGBASE.TXT") \\ \\
\# Load model \\
frbus = Frbus("model.xml") \\ \\
\# Specify dates \\
start = pandas.Period("2040Q1") \\
end = start + 23 \\ \\
\# Standard configuration, use surplus ratio targeting \\
data.loc[start:end, "dfpdbt"] = 0 \\
data.loc[start:end, "dfpsrp"] = 1 \\ \\
\# Solve to baseline with adds \\
with_adds = frbus.init_trac(start, end, data) \\ \\
\# 100 bp monetary policy shock \\
with_adds.loc[start, "rffintay_aerr"] += 1 \\ \\
\# Solve \\
sim = frbus.solve(start, end, with_adds) \\ \\
\# View results \\
sim_plot(with_adds, sim, start, end) \\
}
\end{footnotesize}

\bigskip

\bigskip

\proglang{R} version of the same exercise follows:
\begin{footnotesize}
<<>>=
library(bimets)
@
<<>>=
# Load data
data(LONGBASE)
@
<<>>=
# Load model
data(FRB__MODEL)
model <- LOAD_MODEL(modelText = FRB__MODEL)
@
<<>>=
# Load data into model
model <- LOAD_MODEL_DATA(model, LONGBASE, quietly=TRUE)
@
<<>>=
# Specify dates
start <- c(2040,1)
end <- normalizeYP(start+c(0,23),4)
@
<<>>=
# Standard configuration, use surplus ratio targeting
model$modelData$dfpdbt[[start,end]] <- 0
model$modelData$dfpsrp[[start,end]] <- 1
@
<<>>=
# Solve to baseline with adds
model <- SIMULATE(model,
                  simType='RESCHECK',
                  TSRANGE=c(start,end),
                  ZeroErrorAC = TRUE,
                  quietly=TRUE)
@
<<>>=
# 100 bp monetary policy shock
trac <- model$ConstantAdjustmentRESCHECK
trac$rffintay[[start]] <- trac$rffintay[[start]]+1
@
<<>>=
# Solve
model <- SIMULATE(model,
                  simAlgo = 'NEWTON',
                  TSRANGE = c(start,end),
                  ConstantAdjustment = trac,
                  BackFill = 12,
                  quietly=TRUE)
@
<<eval=FALSE>>=
# View results
sim_plot(model,c(start,end),1)
@
\end{footnotesize}

\bigskip

\proglang{python} code produces the following charts:

\includegraphics[width=5in]{FRB_python_1.png}

\clearpage

On the other hand, \pkg{bimets} code produces very similar results:

\bigskip

<<fig=TRUE,echo=FALSE,width=7,height=7>>=
# View results
sim_plot(model,c(start,end),1)
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Rational expectations} \label{ssec:2}

The second econometric exercise proposed by the Federal Reserve is the simulation of a rational expectations variation of the FRB/US model (i.e, MCE, Model-Consistent Expectations) under a monetary policy shock. Again, the \code{rffintay} time series, is shocked by 100 base points.

\bigskip

14 equations of the original model (i.e. \code{zdivgr}, \code{zgap05}, \code{zgap10}, \code{zgap30}, \code{zpi10}, \code{zpi10f}, \code{zpib5}, \code{zpic30}, \code{zpic58}, \code{zpicxfe}, \code{zpieci}, \code{zrff10}, \code{zrff30}, \code{zrff5}), have been modified and present a forward-looking definition.

\bigskip

For example, in the MCE version of the model, the \code{zdivgr} equation presents a forward-looking reference to its future values: \\
\code{zdivgr=(0.009757264257434617)*TSLEAD(hgynid)+(0.9902427357425654)*TSLEAD(zdivgr)} 

\bigskip

Original simulation time range has been lowered from 60 years to 2 years in order to reduce calculation time in examples (see \emph{"Computational details"} on section \ref{ssec:6}). \\ \\

\proglang{python} code follows: \\ \\
\begin{footnotesize}
\code{
import pandas \\ \\
from pyfrbus.frbus import Frbus \\
from pyfrbus.sim_lib import sim_plot \\
from pyfrbus.load_data import load_data \\ \\
\# Load data \\
data = load_data("./LONGBASE.TXT") \\ \\
\# Load model \\
frbus = Frbus("./model.xml", mce="mcap+wp") \\ \\
\# Specify dates \\
start = pandas.Period("2040q1") \\
end = start + 8 \\ \\
\# Standard MCE configuration, use surplus ratio targeting, rstar endogenous in long run \\
data.loc[start:end, "dfpdbt"] = 0 \\
data.loc[start:end, "dfpsrp"] = 1 \\
data.loc[start:end, "drstar"] = 0 \\
data.loc[(start+4):end, "drstar"] = 1 \\ \\
\# Solve to baseline with adds \\
with_adds = frbus.init_trac(start, end, data) \\ \\
\# 100 bp monetary policy shock and solve \\
with_adds.loc[start, "rffintay_aerr"] += 1 \\ \\
\# Solve \\
sim = frbus.solve(start, end, with_adds) \\ \\
\# View results \\
sim_plot(with_adds, sim, start, end) 
}
\end{footnotesize}

\bigskip

\bigskip

\proglang{R} version of the same exercise follows:
\begin{footnotesize}
<<>>=
library(bimets)
@
<<>>=
# Load data
data(LONGBASE)
@
<<>>=
# Load model
data(FRB__MCAP__WP__MODEL)
model <- LOAD_MODEL(modelText = FRB__MCAP__WP__MODEL)
@
<<>>=
# Load data into model
model <- LOAD_MODEL_DATA(model, LONGBASE, quietly=TRUE)
@
<<>>=
# Specify dates
start <- c(2040,1)
end <- normalizeYP(start+c(0,8),4)
@
<<>>=
# Standard MCE configuration, use surplus ratio targeting, rstar endogenous in long run
model$modelData$dfpdbt[[start,end]] <- 0
model$modelData$dfpsrp[[start,end]] <- 1
model$modelData$drstar[[start,end]] <- 0
model$modelData$drstar[[normalizeYP(start+c(0,4),4),end]] <- 1 
@
<<>>=
# Solve to baseline with adds
model <- SIMULATE(model,
                  simType = 'RESCHECK',
                  TSRANGE = c(start,end),
                  ZeroErrorAC = TRUE,
                  quietly=TRUE)
@
<<>>=
# 100 bp monetary policy shock
shock <- model$ConstantAdjustmentRESCHECK
shock$rffintay[[start]] <- shock$rffintay[[start]]+1
@
<<>>=
# Solve
model <- SIMULATE(model,
                  simAlgo = 'NEWTON',
                  TSRANGE = c(start,end),
                  ConstantAdjustment = shock,
                  BackFill = 12,
                  quietly=TRUE)
@
<<eval=FALSE>>=
# View results
sim_plot(model,c(start,end),2)
@
\end{footnotesize}

\clearpage

\proglang{python} code produces the following charts: \\

\includegraphics[width=5in]{FRB_python_2.png}

\clearpage

On the other hand, \pkg{bimets} code produces very similar results:

\bigskip

<<fig=TRUE,echo=FALSE,width=7,height=7>>=
# View results
sim_plot(model,c(start,end),2)
@


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Auto-correlation on tracking residuals} \label{ssec:3}
The 3rd econometric exercise proposed by the Federal Reserve is a dynamic simulation of the FRB/US model (backward-looking version) wherein historical tracking residuals have been forward extended by using an auto-correlation (i.e., \emph{persistence}) factor of \code{0.5}.

\bigskip

\proglang{python} code follows: \\ \\
\begin{footnotesize}
\code{
import pandas \\ \\
from pyfrbus.frbus import Frbus \\
from pyfrbus.sim_lib import sim_plot \\
from pyfrbus.load_data import load_data \\
from pyfrbus.time_series_data import TimeSeriesData \\ \\
\# Load data \\
data = load_data("./LONGBASE.TXT") \\ \\
\# Load model \\
frbus = Frbus("./model.xml") \\ \\
\# Specify dates \\
start = pandas.Period("2040Q1") \\
end = start + 24 \\ \\
\# Standard configuration, use surplus ratio targeting \\
data.loc[start:end, "dfpdbt"] = 0 \\
data.loc[start:end, "dfpsrp"] = 1 \\ \\
\# Use non-inertial Taylor rule \\
data.loc[start:end, "dmptay"] = 1 \\
data.loc[start:end, "dmpintay"] = 0 \\ \\
\# Enable thresholds \\
data.loc[start:end, "dmptrsh"] = 1 \\ \\
\# Arbitrary threshold values \\
data.loc[start:end, "lurtrsh"] = 6.0 \\
data.loc[start:end, "pitrsh"] = 3.0 \\ \\
\# Solve to baseline with adds \\
with_adds = frbus.init_trac(start, end, data) \\ \\
\# Zero tracking residuals for funds rate and thresholds \\
with_adds.loc[start:end, "rfftay_trac"] = 0 \\
with_adds.loc[start:end, "rffrule_trac"] = 0 \\
with_adds.loc[start:end, "rff_trac"] = 0 \\
with_adds.loc[start:end, "dmptpi_trac"] = 0 \\
with_adds.loc[start:end, "dmptlur_trac"] = 0 \\
with_adds.loc[start:end, "dmptmax_trac"] = 0 \\
with_adds.loc[start:end, "dmptr_trac"] = 0 \\ \\
\# Shocks vaguely derived from historical residuals \\
with_adds.loc[start:start + 3, "eco_aerr"] = [-0.002, -0.0016, -0.0070, -0.0045] \\
with_adds.loc[start:start + 3, "ecd_aerr"] = [-0.0319, -0.0154, -0.0412, -0.0838] \\
with_adds.loc[start:start + 3, "eh_aerr"] = [-0.0512, -0.0501, -0.0124, -0.0723] \\
with_adds.loc[start:start + 3, "rbbbp_aerr"] = [0.3999, 2.7032, 0.3391, -0.7759] \\
with_adds.loc[start:start + 8, "lhp_aerr"] = [-0.0029,-0.0048,-0.0119,-0.0085, ... \\
 -0.0074,-0.0061,-0.0077,-0.0033,-0.0042,] \\ \\
\# Set up time-series object \\
d = TimeSeriesData(with_adds) \\
\# Roll off residuals with 0.5 persistence \\
rho = 0.5 \\
\# Set range \\
d.range = pandas.period_range(start + 4, end) \\
d.eco_aerr = rho * d.eco_aerr(-1) \\
d.ecd_aerr = rho * d.ecd_aerr(-1) \\
d.eh_aerr = rho * d.eh_aerr(-1) \\
d.rbbbp_aerr = rho * d.rbbbp_aerr(-1) \\
d.range = pandas.period_range(start + 9, end) \\
d.lhp_aerr = rho * d.lhp_aerr(-1) \\ \\
\# Adds so that thresholds do not trigger before shocks are felt \\
with_adds.loc[start, "dmptr_aerr"] = -1 \\
with_adds.loc[start : start + 2, "dmptlur_aerr"] = -1 \\ \\
\# Solve \\
sim = frbus.solve(start, end, with_adds) \\ \\
\# View results, unemployment threshold binds \\
sim_plot(with_adds, sim, start, end) \\
}
\end{footnotesize}

\bigskip

\proglang{R} version of the same exercise follows:

\begin{footnotesize}
<<>>=
library(bimets)
@
<<>>=
# Load data
data(LONGBASE)
@
<<>>=
# Load model
data(FRB__MODEL)
model <- LOAD_MODEL(modelText = FRB__MODEL)
@
<<>>=
# Load data into model
model <- LOAD_MODEL_DATA(model, LONGBASE, quietly=TRUE)
@
<<>>=
# Specify dates
start <- c(2040,1)
end <- normalizeYP(start+c(0,24),4)
@
<<>>=
# Standard configuration, use surplus ratio targeting
model$modelData$dfpdbt[[start,end]] <- 0
model$modelData$dfpsrp[[start,end]] <- 1
@
<<>>=
# Use non-inertial Taylor rule
model$modelData$dmptay[[start,end]] <- 1
model$modelData$dmpintay[[start,end]] <- 0
@
<<>>=
# Enable thresholds
model$modelData$dmptrsh[[start,end]] <- 1
@
<<>>=
# Arbitrary threshold values
model$modelData$lurtrsh[[start,end]] <- 6
model$modelData$pitrsh[[start,end]] <- 3
@
<<>>=
# Solve to baseline with adds
model <- SIMULATE(model,
                  simType = 'RESCHECK',
                  TSRANGE = c(start,end),
                  ZeroErrorAC = TRUE,
                  quietly=TRUE)
@
<<>>=
# Get tracking residuals        
trac <- model$ConstantAdjustmentRESCHECK
@
<<>>=
# Zero tracking residuals for funds rate and thresholds 
trac$rfftay[[start,end]] <- 0
trac$rffrule[[start,end]] <- 0
trac$rff[[start,end]] <- 0
trac$dmptpi[[start,end]] <- 0
trac$dmptlur[[start,end]] <- 0
trac$dmptmax[[start,end]] <- 0
trac$dmptr[[start,end]] <- 0
@
<<>>=
# Shocks vaguely derived from historical residuals
aerr <- list()
aerr$eco <- TSERIES(c(-0.002, -0.0016, -0.0070, -0.0045),START=start,FREQ=4)
aerr$ecd <- TSERIES(c(-0.0319, -0.0154, -0.0412, -0.0838),START=start,FREQ=4)
aerr$eh <- TSERIES(c(-0.0512, -0.0501, -0.0124, -0.0723),START=start,FREQ=4)
aerr$rbbbp <- TSERIES(c(0.3999, 2.7032, 0.3391, -0.7759),START=start,FREQ=4)
aerr$lhp <- TSERIES(c(-0.0029,-0.0048,-0.0119,-0.0085,-0.0074,-0.0061,-0.0077,-0.0033,-0.0042),
            START=start,FREQ=4)
@
<<>>=
# Roll off residuals with 0.5 persistence
rho <- 0.5
aerr$eco <- TSEXTEND(aerr$eco,UPTO=end,EXTMODE='MYRATE',FACTOR=rho)
aerr$ecd <- TSEXTEND(aerr$ecd,UPTO=end,EXTMODE='MYRATE',FACTOR=rho)
aerr$eh <- TSEXTEND(aerr$eh,UPTO=end,EXTMODE='MYRATE',FACTOR=rho)
aerr$rbbbp <- TSEXTEND(aerr$rbbbp,UPTO=end,EXTMODE='MYRATE',FACTOR=rho)
aerr$lhp <- TSEXTEND(aerr$lhp,UPTO=end,EXTMODE='MYRATE',FACTOR=rho)
@
<<>>=
# Adds so that thresholds do not trigger before shocks are felt
aerr$dmptr <- TSERIES(c(-1),START=start,FREQ=4)
aerr$dmptlur <- TSERIES(c(-1,-1,-1),START=start,FREQ=4)
@
<<>>=
# Create Constant Adjustments for SIMULATE op.
for (idx in names(aerr)) trac[[idx]] <- trac[[idx]]+aerr[[idx]]
@
<<>>=
# Solve
model <- SIMULATE(model,
                  simAlgo = 'NEWTON',
                  TSRANGE = c(start,end),
                  ConstantAdjustment = trac,
                  BackFill = 12,
                  quietly=TRUE)  
@
<<eval=FALSE>>=
# View results, unemployment threshold binds
sim_plot(model,c(start,end),3)
@
\end{footnotesize}

\clearpage

\proglang{python} code produces the following charts: \\

\includegraphics[width=5in]{FRB_python_3.png}

\clearpage

On the other hand, \pkg{bimets} code produces very similar results:

\bigskip

<<fig=TRUE,echo=FALSE,width=7,height=7>>=
# View results
sim_plot(model,c(start,end),3)
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Endogenous targeting} \label{ssec:4}

The 4th econometric exercise proposed by the Federal Reserve is an endogenous targeting exercise. The \code{mcontrol} procedure in the \pkg{pyfrbus} package calculates the constant adjustments values to be applied to 5 instruments variables (i.e., \code{eco}, \code{lhp}, \code{picxfe}, \code{rff}, \code{rg10p}, all in the \code{inst} list) such that the simulated values for the endogenous variables in the target list \code{targ} (i.e., \code{xgdp}, \code{lur}, \code{picxfe}, \code{rff}, \code{rg10}) are equal to the values of the arbitrary trajectories in the \code{traj} list.

\bigskip

In a similar way in \proglang{R}, the \pkg{bimets} \code{RENORM} procedure determines the values for the \code{INSTRUMENT} exogenous variables (i.e., the constant adjustments in the \code{inst} list in this exercise) that allow the objective \code{TARGET} endogenous 
values to be achieved, with respect to the constraints given by the model equations.

\bigskip

\proglang{python} code follows: \\ \\
\begin{footnotesize}
\code{
import pandas \\
from numpy import array, cumprod \\ \\
from pyfrbus.frbus import Frbus \\
from pyfrbus.sim_lib import sim_plot \\
from pyfrbus.load_data import load_data \\ \\
\# Load data \\
data = load_data("./LONGBASE.TXT") \\ \\
\# Load model \\
frbus = Frbus("./model.xml") \\ \\
\# Specify dates \\
start = pandas.Period("2021Q3") \\
end = "2022Q3" \\ \\
\# Standard configuration, use surplus ratio targeting \\
data.loc[start:end, "dfpdbt"] = 0 \\
data.loc[start:end, "dfpsrp"] = 1 \\ \\
\# Solve to baseline with adds \\
with_adds = frbus.init_trac(start, end, data) \\ \\
\# Scenario based on 2021Q3 Survey of Professional Forecasters \\
with_adds.loc[start:end, "lurnat"] = 3.78 \\ \\
\# Set up trajectories for mcontrol \\
with_adds.loc[start:end, "lur_t"] = [5.3, 4.9, 4.6, 4.4, 4.2] \\
with_adds.loc[start:end, "picxfe_t"] = [3.7, 2.2, 2.1, 2.1, 2.2] \\
with_adds.loc[start:end, "rff_t"] = [0.1, 0.1, 0.1, 0.1, 0.1] \\
with_adds.loc[start:end, "rg10_t"] = [1.4, 1.6, 1.6, 1.7, 1.9] \\ \\
\# Get GDP level as accumulated growth from initial period \\
gdp_growth = cumprod((array([6.8, 5.2, 4.5, 3.4, 2.7]) / 100 + 1) ** 0.25) \\
with_adds.loc[start:end, "xgdp_t"] = with_adds.loc[start - 1, "xgdp"] * gdp_growth \\ \\
targ = ["xgdp", "lur", "picxfe", "rff", "rg10"] \\
traj = ["xgdp_t", "lur_t", "picxfe_t", "rff_t", "rg10_t"] \\
inst = ["eco_aerr", "lhp_aerr", "picxfe_aerr", "rff_aerr", "rg10p_aerr"] \\ \\
\# Run mcontrol \\
sim = frbus.mcontrol(start, end, with_adds, targ, traj, inst) \\ \\
\# View results \\
sim_plot(with_adds, sim, start, end) 
}
\end{footnotesize}

\bigskip

\bigskip

\proglang{R} version of the same exercise follows:
\begin{footnotesize}
<<>>=
library(bimets)
@
<<>>=
# Load data
data(LONGBASE)
@
<<>>=
# Load model
data(FRB__MODEL)
model <- LOAD_MODEL(modelText = FRB__MODEL)
@
<<>>=
# Load data into model
model <- LOAD_MODEL_DATA(model, LONGBASE, quietly=TRUE)
@
<<>>=
# Specify dates
start <- c(2021,3)
end <- c(2022,3)
@
<<>>=
# Standard configuration, use surplus ratio targeting
model$modelData$dfpdbt[[start,end]] <- 0
model$modelData$dfpsrp[[start,end]] <- 1
@
<<>>=
# Solve to baseline with adds
model <- SIMULATE(model,
                  simType = 'RESCHECK',
                  TSRANGE = c(start,end),
                  ZeroErrorAC = TRUE
                  ,quietly=TRUE)
@
<<>>=         
# Scenario based on 2021Q3 Survey of Professional Forecasters
model$modelData$lurnat[[start,end]] <- 3.78
@
<<>>=
# Set up trajectories for mcontrol
targ <- list()
targ$lur <- TSERIES(c(5.3, 4.9, 4.6, 4.4, 4.2),START=start,FREQ=4)
targ$picxfe <- TSERIES(c(3.7, 2.2, 2.1, 2.1, 2.2),START=start,FREQ=4)
targ$rff <- TSERIES(c(0.1, 0.1, 0.1, 0.1, 0.1),START=start,FREQ=4)
targ$rg10 <- TSERIES(c(1.4, 1.6, 1.6, 1.7, 1.9),START=start,FREQ=4)
@
<<>>=
# Get GDP level as accumulated growth from initial period
gdp_growth <- model$modelData$xgdp[[2021,2]]*CUMPROD((c(6.8,5.2,4.5,3.4,2.7) / 100 + 1) ** 0.25)
targ$xgdp <- TSERIES(gdp_growth,START=start,FREQ=4)
@
<<>>=
# define INSTRUMENT
inst <- c("eco", "lhp", "picxfe", "rff", "rg10p")
@
<<>>=
# Run RENORM
model <- RENORM(model,
                simAlgo = 'NEWTON',
                TSRANGE=c(start,end),
                ConstantAdjustment = model$ConstantAdjustmentRESCHECK,
                TARGET=targ,
                INSTRUMENT=inst,
                BackFill = 8,
                quietly=TRUE)
@
<<eval=FALSE>>=
# View results
sim_plot(model,c(start,end),4)
@  
\end{footnotesize}
 
\proglang{python} code produces the following charts: \\

\includegraphics[width=5in]{FRB_python_4.png}

\clearpage

On the other hand, \pkg{bimets} code produces very similar results:

\bigskip

<<fig=TRUE,echo=FALSE,width=7,height=7>>=
# View results
sim_plot(model,c(start,end),4)
@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Stochastic simulation} \label{ssec:5}
The 5th econometric exercise proposed by the Federal Reserve is a stochastic simulation of the FRB/US model (backward-looking version). The \code{stochsim} procedure in the \pkg{pyfrbus} package performs a stochastic simulation by applying sequences of shocks to the model, as drawn randomly from historical residuals.

\bigskip

The \code{stochsim} procedure begins by drawing \code{nrepl} sequences of quarters over the periods \code{residstart} to \code{residend}, where the length of that sequence goes from \code{simstart} to \code{simend}. That is, for a particular replication, each quarter in the simulation period is randomly assigned a quarter from residual period. In that quarter of the simulation, all the 64 stochastic variables (specified with a \code{stochastic_type} tag in the XML model file) have a shock applied from a particular quarter in the residual period.

\bigskip

In a similar way in \proglang{R}, the \pkg{bimets} \code{STOCHSIMULATE} procedure allow users to shock, with arbitrary values, all the \code{StochReplica} replicas of variables that appear in the \code{StochStructure} list; the shocks, in this exercise, are the randomly sampled historical residuals in the time range from \code{residstart} to \code{residend}, and are added up to the 64 constant adjustment of the related endogenous variable listed as stochastic, thus replicating the \proglang{python} code. The initial mapping per period, between historical residuals and stochastic realizations, is stored in the \code{sampleHistoricalResidual} variable.

\bigskip

\proglang{python} code follows: \\ \\
\begin{footnotesize}
\code{
from pyfrbus.frbus import Frbus \\
from pyfrbus.sim_lib import stochsim_plot \\
from pyfrbus.load_data import load_data \\ \\
\# Load data \\
data = load_data("./LONGBASE.TXT") \\ \\
\# Load model \\
frbus = Frbus("./model.xml") \\ \\
\# Specify dates and other params \\
residstart = "1975q1" \\
residend = "2018q4" \\
simstart = "2040q1" \\
simend = "2045q4" \\
\# Number of replications \\
nrepl = 1000 \\
\# Run up to 5 extra replications, in case of failures \\
nextra = 5 \\ \\
\# Policy settings \\
data.loc[simstart:simend, "dfpdbt"] = 0 \\
data.loc[simstart:simend, "dfpsrp"] = 1 \\ \\
\# Compute add factors \\
\# Both for baseline tracking and over history, to be used as shocks \\
with_adds = frbus.init_trac(residstart, simend, data) \\ \\
\# Call FRBUS stochsim procedure \\
solutions = frbus.stochsim(nrepl, with_adds, simstart, simend, residstart, residend, nextra=nextra) \\ \\
\# View results \\
stochsim_plot(with_adds, solutions, simstart, simend)
}
\end{footnotesize}

\clearpage

\proglang{R} version of the same exercise follows:
\begin{footnotesize}
<<>>=
library(bimets)
@
<<>>=
# Load data
data(LONGBASE)
@
<<>>=
# Load model
data(FRB__MODEL)
model <- LOAD_MODEL(modelText = FRB__MODEL)
@
<<>>=
# Load data into model
model <- LOAD_MODEL_DATA(model, LONGBASE, quietly=TRUE)
@
<<>>=
# Specify dates and other params
residstart <- c(1975,1)
residend <- c(2018,4)
simstart <- c(2040,1)
simend <- c(2045,4)
@
<<>>=
# Number of replications
nrepl <- 1000
@
<<>>=
# Policy settings
model$modelData$dfpdbt[[simstart,simend]] <- 0
model$modelData$dfpsrp[[simstart,simend]] <- 1
@
<<>>=
# Compute add factors
# Both for baseline tracking and over history, to be used as shocks
model <- SIMULATE(model,
                  simType = 'RESCHECK',
                  TSRANGE = c(residstart,simend),
                  ZeroErrorAC = TRUE,
                  quietly=TRUE)
@
<<>>=
# Get tracking residuals
trac <- model$ConstantAdjustmentRESCHECK
@
<<>>=
# Set seed
set.seed(9)
@
<<>>=
# 64 Stochastic vars as listed in XML FRB/US model file
stochasticVars <- c("ebfi","ecd","ech","eco","egfe","egfen","egfet","egfl","egse",
                  "egsen","egset","egsl","eh","emo","emp","ex","fpxrr","fxgap",
                  "ugfsrp","gtn","gtr","gtrd","hmfpt","hqlfpr","hqlww","ki","leg",
                  "leo","lfpr","lhp","lurnat","lww","mfpt","pbfir","pcer",
                  "pcfr","pegfr","pegsr","phouse","phr","picxfe","pieci","pmo",
                  "poilr","pxr","rbbbp","rcar","rcgain","reqp","rfynic",
                  "rfynil","rg10p","rg30p","rg5p","rgfint","rme","tcin","tpn",
                  "trci","trp","trpt","uynicpnr","ynidn","ynirn")
@
<<>>=
# Pseudo random array that maps residuals range to simulation range for each replica
residualsLength <- NUMPERIOD(residstart,residend,4)+1
stochSimLength <- NUMPERIOD(simstart,simend,4)+1
sampleHistoricalResidual <- sample(1:residualsLength,stochSimLength*nrepl,replace=T)
@
<<>>=
# Create BIMETS stochastic structure  
modelStochStructure <- list()
for (tmpStochVar in stochasticVars)
{
    #see BIMETS reference manual for details on STOCHSIMULATE and StochStructure
    modelStochStructure[[tmpStochVar]] <- list()
    modelStochStructure[[tmpStochVar]]$TSRANGE <- TRUE
    modelStochStructure[[tmpStochVar]]$TYPE <- 'MATRIX'
    shockMatrix <- matrix(trac[[tmpStochVar]][sampleHistoricalResidual],
                          nrow=stochSimLength,ncol=nrepl)
    shockMatrix <- shockMatrix - ave(shockMatrix)
    modelStochStructure[[tmpStochVar]]$PARS <- shockMatrix
}
@
<<>>=
# Call BIMETS stoch sim procedure
model <- STOCHSIMULATE(model,
                       simAlgo = 'NEWTON',
                       TSRANGE = c(simstart,simend),
                       StochStructure = modelStochStructure,
                       StochReplica = nrepl,
                       ConstantAdjustment = trac,
                       quietly=TRUE)
@
<<eval=FALSE>>=
# View results
stochsim_plot(model,c(simstart,simend))
@
\end{footnotesize}

\clearpage

\proglang{python} code produces the following charts: \\

\includegraphics[width=5in]{FRB_python_5.png}

\clearpage

On the other hand, \pkg{bimets} code produces similar results, despite the random numbers generator is different between \proglang{R} and \proglang{python}:

\bigskip

<<fig=TRUE,echo=FALSE,width=7,height=7>>=
# View results
stochsim_plot(model,c(simstart,simend))
@


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Computational details and capabilities: pyfrbus vs bimets} \label{ssec:6}

Both \proglang{R} (\proglang{R}-4.2.0 optimized with Intel\textsuperscript{\textregistered} MKL libraries) and \proglang{python} (Intel\textsuperscript{\textregistered} \proglang{python} 3.7.9) code have been tested on Red Hat Enterprise Linux 8.8 with Intel\textsuperscript{\textregistered} Xeon\textsuperscript{\textregistered} Gold 6248 CPU @ 2.50GHz 40 cores and 1.2TB RAM; all econometric exercises on the backward-looking version of the FRB/US model are fast in both \pkg{bimets} and \pkg{pyfrbus}, having just a few seconds of computation time for the simulation in both scenarios. Endogenous targeting (see \ref{ssec:4}) is faster in \pkg{pyfrbus}, while stochastic simulation (see \ref{ssec:5}) is faster in \pkg{bimets}, but generally speaking the \pkg{pyfrbus} package relies on compiled code and, in scenarios that require heavy computational resources, e.g., rational expectation models,  is faster than \pkg{bimets}, that is pure \proglang{R} code.
\bigskip

Indeed, \pkg{pyfrbus} is way faster than \pkg{bimets} in forward-looking large models with an extensive simulation time range. \pkg{pyfrbus} can simulate 60 years in FRB/US rational expectations exercise (see \ref{ssec:2}), with just a few minutes of execution time; in \pkg{bimets}, exclusively for that exercise, a FRB/US simulation with a time range greater than 10 years is not feasible. The \proglang{python} package relies on the \code{SuiteSparse} libraries, a suite of compiled \proglang{C} code that provides optimized algorithms for sparse matrices. Moreover, in \pkg{pyfrbus}, the model Jacobian is computed symbolically using \pkg{SymPy} and \pkg{SymEngine}, thus the computation time is very fast also in Netwon algorithms; on the other hand, there is no option for a numerically-approximated Jacobian. As a result, the platform supports a limited number of functions; thus, in \pkg{pyfrbus} only the following functions are supported in model equations: \code{log}, \code{exp}, \code{abs}, \code{min}, and \code{max}.

\bigskip

On the other hand, the \pkg{pyfrbus} package has been developed with a focus on the FRB/US model and presents the following limitations with respect to \pkg{bimets}, which is a general framework for econometric modeling:

\begin{itemize}
\item[--] \pkg{pyfrbus} only supports quarterly models, while \pkg{bimets} supports models of any frequency listed in the time series constructor \code{TIMESERIES} (e.g., annual, semiannual, quarterly, monthly, daily, weekly, etc.);

\item[--] in \pkg{pyfrbus}, no estimation nor structural stability procedures are available, while in \pkg{bimets} users can insert coefficients in model equations, then perform estimation and structural stability analysis accounting for linear restrictions on coefficients, polynomial distributed lags, and residuals auto-correlation;

\item[--] \pkg{pyfrbus} does not support conditional evaluation in equations, while in \pkg{bimets} users can use the \code{IF>} statements in model equations;

\item[--] as said, in \pkg{pyfrbus}, due to the symbolic Jacobian implementation, only \code{log}, \code{exp}, \code{abs}, \code{min}, and \code{max} functions are allowed in model equations, while \pkg{bimets} allows for more complex time series transformations to be inserted  directly into the model definition, e.g., \code{TSDELTA}, \code{TSDELTALOG}, \code{MOVAVG}, etc.; on the other hand, at the moment \pkg{bimets} does not support \code{MIN} and \code{MAX} functions in model equations, but users can use \code{IF>} statements in \code{\code{MDL}} model definitions in order to create an equivalent syntax, as in the FRB/US \code{MDL} version in these exercises; \code{MIN}, \code{MAX}, and other \code{MDL} functions will be available in the upcoming \pkg{bimets} versions.

\item[--] in \pkg{bimets}, detailed messages in error and in verbose operations can help users in identifying and solving grammatical and numerical issues;

\item[--] in time series operations, \pkg{pyfrbus} relies on \pkg{numpy} and \pkg{pandas} capabilities, while \pkg{bimets} allows for more complex time series transformations, e.g., (dis)aggregations, extension, merge, projection, etc.;

\item[--] optimal control for econometric models is a missing capability in \pkg{pyfrbus};

\item[--] stochastic simulation in \pkg{pyfrbus} only allows for shocks randomly selected from past residuals tracking, while in \pkg{bimets} users have more flexibility, i.e., uniform and normal stochastic shocks, and arbitrary matrices of shocks, in arbitrary periods; moreover in \pkg{pyfrbus} the list of stochastic variables is hard-coded into the model definition, while in \pkg{bimets} users can change stochastic variables interactively in code by changing the \code{StochStructure} variable in the procedure's call;

\item[--] in \pkg{bimets}, constant adjustments are function's arguments, therefore users can disable or change, in arbitrary periods, the add-factor impact in model equilibrium, interactively in code with no need to modify the model data set, as required in \pkg{pyfrbus};

\item[--] in \pkg{bimets}, users can exogenize an endogenous variable in arbitrary periods, while in \pkg{pyfrbus} exogenization is forced on the whole simulation time range;

\item[--] in \pkg{pyfrbus}, multiplier analysis is not available;

\item[--] \pkg{pyfrbus} is distributed only for Linux operating system, with instructions on how to run Linux code on a Microsoft\textsuperscript{\textregistered} Windows machine. Moreover, the package has a complex list of dependencies (e.g., \code{UMFPACK}, \code{SuiteSparse}, \code{scipy}, \code{sympy}, \code{symengine}, \code{scikit}, \code{numpy}, \code{pandas}, etc.) and its installation could be difficult; \pkg{bimets}, which depends only on base \proglang{R} and \code{xts}, does not require any compilation, and is easy to install in Linux, Microsoft\textsuperscript{\textregistered} Windows, and Apple\textsuperscript{\textregistered} OS-X.
\end{itemize}

More details about \pkg{bimets} are available at the \href{https://cran.r-project.org/package=bimets/vignettes/bimets.pdf}{"Getting started with bimets"} vignette.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{References}
\begin{enumerate}[label={[\arabic*]}]
\item Andrea Luciani, Bank of Italy - \textit{bimets: Time Series and Econometric Modeling}. R package version 4.0.1, \href{https://CRAN.R-project.org/package=bimets/}{https://CRAN.R-project.org/package=bimets/}, \\ GIT: \href{https://github.com/andrea-luciani/bimets}{https://github.com/andrea-luciani/bimets}
\item Board of governors of The Federal Reserve System - \textit{pyfrbus: a Python-based platform to run simulations with the FRB/US model.}. python package version 1.0.0, \href{https://www.federalreserve.gov/econres/us-models-about.htm}{https://www.federalreserve.gov/econres/us-models-about.htm}
\end{enumerate}
\end{document}