R/sfcr_baseline.R
sfcr_baseline.RdThe sfcr_baseline() function is used to simulate a SFC model.
sfcr_baseline( equations, external, periods, initial = NULL, hidden = NULL, max_iter = 350, .hidden_tol = 0.1, tol = 1e-08, method = "Broyden", rhtol = FALSE, ... )
| equations | A |
|---|---|
| external, initial | A |
| periods | A number specifying the total number of periods of the model to be simulated. |
| hidden | Named object that identify the two variables that make the hidden equality
in the SFC model, e.g., |
| max_iter | Maximum iterations allowed per period. |
| .hidden_tol | Error tolerance to accept the equality of the hidden equation. Defaults to 1.
In growth models, computational errors might buildup in the hidden equation, which renders any absolute
comparison inadequate. For such models, please turn |
| tol | Tolerance accepted to determine convergence. |
| method | The method to use to find a solution. Defaults to "Broyden". |
| rhtol | A logical argument that defines whether the a relative measure is used to evaluate
the hidden equation or not. Defaults to |
| ... | Extra arguments to pass to |
A sfcr_tbl.
The output of a sfcr_baseline() is a sfcr_tbl. The only difference between
a sfcr_tbl and a standard tbl_df is that the former has two extra attributes:
matrix and call. The matrix attribute, for example, can be accessed by
calling attributes(sfcr_sim_object)$matrix.
It is possible to see, in the matrix, the number of iterations required to calculate each
block of equations in the model.
The call attribute shows the blocks of equations and preserve the call that are used
internally.
The equations, exogenous, and parameters arguments must be written
with the R formula syntax, i.e., the left-hand side of each item is separated to the
right-hand side by a twiddle. Variables that represent lags of endogenous or exogenous
variables must be followed by [-1]. See examples for details on the syntax.
Before solving the system of equations, two consecutive depth-first searches identify
and order the blocks of independent equations in the system. The system is then solved
sequentially, i.e., the variables that depend only on lagged or exogenous values are evaluated
first, and then the variables that depends on these variables, etc. The solving algorithms
are only applied to the blocks of mutually dependent equations. The great igraph
package is used to implement the two consecutive depth-first searches.
Methods:
The sfcr package provides three algorithms to solve the blocks of cyclical equations:
the Gauss-Seidel algorithm, the Broyden algorithm, and the Newton-Raphson algorithm. The
default method is "Broyden" as it tends to be fastest one.
See (Kinsella and OShea 2010) for details on the Gauss-Seidel algorithm and (Peressini et al. 1988) for details on the Broyden and Newton-Raphson algorithms.
The "Broyden" algorithm uses the rootSolve::jacobian.full() function to get the
initial Jacobian matrix, and compiled code from RcppArmadillo to invert the
jacobians. See also https://www.math.usm.edu/lambers/mat419/lecture11.pdf.
The Gauss Seidel algorithm is implemented as described by (Kinsella and OShea 2010)
.
Finally, the "Newton" method uses the rootSolve::multiroot() function to solve the system.
Hidden equation:
One of the defining aspects of a SFC model is its water tight accounting. One way
to check whether the model was correctly defined is to see if the hidden (redundant)
equation is satisfied after the model is simulated. In stationary models, an absolute
comparison should suffice as the model converges to a stationary state. However,
growth models converge to a stable growth rate where stocks are perpetually increasing.
It is inadequate to use a absolute comparison in such models. In these cases, the
rhtol argument ("relative hidden tolerance") must be set to TRUE in order
to perform a relative comparison. The relative comparison evaluates the numerical
discrepancy in the hidden equation as a ratio of one of its elements. For example,
if hidden = c("Bbs" = "Bbd"), the hidden equation will be evaluated according to
the following steps:
d = (Bbs - Bbd)
isTRUE(d/Bbs < .hidden_tol)
In general, the .hidden_tol argument should be set to a small number (e.g. 1e-6).
The function will check that this proportion remains the same for all simulated periods.
Kinsella S, OShea T (2010). “Solution and Simulation of Large Stock Flow Consistent Monetary Production Models via the Gauss Seidel Algorithm.” SSRN Electronic Journal. doi: 10.2139/ssrn.1729205 . Peressini AL, Sullivan FE, Uhl JJ (1988). The Mathematics of Nonlinear Programming. Springer-Verlag, Berlin, Heidelberg. ISBN 0387966145.
João Macalós, joaomacalos@gmail.com
eqs <- sfcr_set( TXs ~ TXd, YD ~ W * Ns - TXs, Cd ~ alpha1 * YD + alpha2 * Hh[-1], Hh ~ YD - Cd + Hh[-1], Ns ~ Nd, Nd ~ Y / W, Cs ~ Cd, Gs ~ Gd, Y ~ Cs + Gs, TXd ~ theta * W * Ns, Hs ~ Gd - TXd + Hs[-1] ) external <- sfcr_set(Gd ~ 20, W ~ 1, alpha1 ~ 0.6, alpha2 ~ 0.4, theta ~ 0.2) # Periods is set to 10 to run faster. A usual model should run at # least 50 periods to find a steady state sfcr_baseline(equations = eqs, external = external, periods = 10)#> # A tibble: 10 × 17 #> period TXs YD Cd Hh Ns Nd Cs Gs #> * <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> #> 1 1 1 e-15 1 e-15 1 e-15 1 e-15 1 e-15 1 e-15 1 e-15 1e-15 #> 2 2 7.69e+ 0 3.08e+ 1 1.85e+ 1 1.23e+ 1 3.85e+ 1 3.85e+ 1 1.85e+ 1 2e+ 1 #> 3 3 9.59e+ 0 3.83e+ 1 2.79e+ 1 2.27e+ 1 4.79e+ 1 4.79e+ 1 2.79e+ 1 2e+ 1 #> 4 4 1.12e+ 1 4.48e+ 1 3.59e+ 1 3.15e+ 1 5.59e+ 1 5.59e+ 1 3.59e+ 1 2e+ 1 #> 5 5 1.25e+ 1 5.02e+ 1 4.27e+ 1 3.90e+ 1 6.27e+ 1 6.27e+ 1 4.27e+ 1 2e+ 1 #> 6 6 1.37e+ 1 5.48e+ 1 4.85e+ 1 4.53e+ 1 6.85e+ 1 6.85e+ 1 4.85e+ 1 2e+ 1 #> 7 7 1.47e+ 1 5.86e+ 1 5.33e+ 1 5.06e+ 1 7.33e+ 1 7.33e+ 1 5.33e+ 1 2e+ 1 #> 8 8 1.55e+ 1 6.19e+ 1 5.74e+ 1 5.52e+ 1 7.74e+ 1 7.74e+ 1 5.74e+ 1 2e+ 1 #> 9 9 1.62e+ 1 6.47e+ 1 6.09e+ 1 5.90e+ 1 8.09e+ 1 8.09e+ 1 6.09e+ 1 2e+ 1 #> 10 10 1.68e+ 1 6.71e+ 1 6.38e+ 1 6.22e+ 1 8.38e+ 1 8.38e+ 1 6.38e+ 1 2e+ 1 #> # … with 8 more variables: Y <dbl>, TXd <dbl>, Hs <dbl>, Gd <dbl>, W <dbl>, #> # alpha1 <dbl>, alpha2 <dbl>, theta <dbl>