Nothing
# ############################################################################
# This file is a part of gEcon. #
# #
# (c) Chancellery of the Prime Minister of the Republic of Poland 2012-2015 #
# (c) Grzegorz Klima, Karol Podemski, Kaja Retkiewicz-Wijtiwiak 2015-2018 #
# License terms can be found in the file 'LICENCE' #
# #
# Authors: Grzegorz Klima, Karol Podemski, Kaja Retkiewicz-Wijtiwiak #
# ############################################################################
# Finding the steady state (equilibrium)
# ############################################################################
# ############################################################################
# Function steady_state computes the steady state
# of the dynamic model or equilibrium of the static model
# ############################################################################
# Input
# model - an object of the gecon_model class
# solver - the name of non-linear equations solver; if NULL the default solver
# is used; currently only the solver from the "nleqslv" package is
# supported
# use_jac - option to use Jacobian generated by symbolic library
# (it can resolve plenty of numerical problems), if FALSE numerical
# derivatives are computed
# calibration - if TRUE, calibrating equations are taken into account while
# solving for the steady state (equilibrium) and parameters
# are calibrated; otherwise calibrating equations are dropped
# options_list - nleqslv solver specific settings; the following options are
# available:
# -> method - a character string specifying the method used
# for finding a solution; it can be set to "Newton"
# or "Broyden"; default option is "Newton"
# -> global - a character string specifying global search strategy
# to use; it can be set to "dbldog", "pwldog", "qline",
# "gline", "none"; the default option is "qline"
# -> xscalm - a character string specifying the type of init_val
# scaling to be applied; it can be set to: "fixed", "auto";
# the default option is "fixed"
# -> max_iter - a numeric value, the maximal number of iterations;
# the default value is 150
# -> tol - a numeric value setting a numeric tolerance for a solution,
# the default value is 1e-6
# -> xtol - a numeric value setting a relative numeric tolerance for
# a solver step, the default value is 1e-6
# solver_status - if TRUE, solver status is printed
# Output
# This function sets the "steady" slot to a computed value.
# Solver status is also updated.
# ############################################################################
steady_state <- function(model,
solver = "nleqslv",
use_jac = TRUE,
calibration = TRUE,
options_list = NULL,
solver_status = FALSE)
{
if (!is.gecon_model(model)) {
stop("model argument should be of gecon_model class")
}
prepstat <- tryCatch(prep <- ss_prep(model, use_jac, calibration),
error = function(e) e)
if (inherits(prepstat, "error")) {
stop("evaluation of steady-state equations\' residuals has not been successful: ",
prepstat)
}
model@init_residual_vector <- prep$resid
slvstat <- tryCatch(res <- ss_find(ss_guess = prep$ss_guess,
ss_eval = prep$ss_eval,
ss_jac = prep$ss_jac,
solver = solver,
options_list = options_list),
error = function(w) w)
if (inherits(slvstat, "error")) {
stop(slvstat)
}
# clearing slots
model@loglin_var <- logical(length = 0)
model@re_solved <- FALSE
model@eig_vals <- matrix(nrow = 0, ncol = 0)
model@solution <- list(P = NULL, Q = NULL, R = NULL, S = NULL)
model@state_var_indices <- numeric(0)
model@solver_exit_info <- character(0)
model@solution_resid <- list(NULL)
model@corr_mat <- matrix(nrow = 0, ncol = 0)
model@corr_computed <- FALSE
model@autocorr_mat <- matrix(nrow = 0, ncol = 0)
model@ref_var_corr_mat <- matrix(nrow = 0, ncol = 0)
model@ref_var_idx <- 0L
model@var_dec <- matrix(nrow = 0, ncol = 0)
model@sdev <- matrix(nrow = 0, ncol = 0)
# results
tol <- 1e-6
tol_index <- which(names(options_list) %in% c("tol"))
if (length(tol_index) > 0) {
tol <- options_list[tol_index]
}
model@solver_status <- res$solver_status
solution <- round2zero(res$solution, tol)
resid <- prep$ss_eval(solution)
model@residual_vector <- resid
if (calibration && (length(model@parameters_calibr) != 0)) {
nv <- length(model@variables)
np <- length(model@parameters_calibr)
nvp <- nv + np
model@variables_ss_val <- solution[1:nv]
model@parameters_calibr_val <- solution[(nv + 1):nvp]
model@parameters_val[model@map_calibr_into_params] <- solution[(nv + 1):nvp]
} else {
model@variables_ss_val <- solution
}
if (any(is.na(model@residual_vector)) || any(!is.finite(model@residual_vector))) {
norm1 <- Inf
} else {
norm1 <- max(abs(model@residual_vector))
}
if (norm1 < tol) {
model@ss_solved <- TRUE
if (model@is_dynamic) { sseq <- "Steady state"; } else { sseq <- "Equilibrium"; }
cat(paste0(sseq, " has been FOUND\n"))
if (solver_status) {
cat(paste0("\nSolver exit info:\n", model@solver_status, "\n"))
}
} else {
if (calibration && any(is.na(model@parameters_calibr_val))) {
warning(paste0("initial values for some calibrated parameters ",
"have not been supplied (default value of 0.5 has been used."))
}
if (any(is.na(model@variables_ss_val))) {
warning(paste0("initial values for some variables ",
"have not been supplied (default value of 0.9 has been used."))
}
if (model@is_dynamic) { sseq <- "steady state"; } else { sseq <- "equilibrium"; }
mes <- paste0("The ", sseq, " has not been found, 1-norm of residuals is: ", norm1,
".\nIt is more than requested precision.\n")
if (solver_status) {
mes <- paste0(mes, "Solver exit info:\n", model@solver_status, "\n")
}
mes <- paste0(mes, "Change initial values and check if the ",
sseq, " can be found.")
warning(mes)
}
# returning
return(model)
}
# ############################################################################
# The get_residuals function prints the residuals of the
# steady state (equilibrium) equations before the first and
# after the last solver iteration
# ############################################################################
# Input
# model - an object of gecon_model class
# largest - the no. of equations with the largest errors to be printed
# calibration - logical. Should calibrating equations be taken into account?
# Output
# The function returns a list of two elements: initial and final.
# ############################################################################
get_residuals <- function(model, largest = 5, calibration = TRUE)
{
if (!is.gecon_model(model)) {
stop("model argument should be of gecon_model class")
}
output <- list()
ss_eq_no <- length(model@equations)
calibr_eq_no <- length(model@calibr_equations)
if (length(model@init_residual_vector)) {
irv <- model@init_residual_vector
if (calibration && (length(irv) < ss_eq_no + calibr_eq_no)) {
warning(paste0("solution was searched for without taking calibrating equations into account, ",
"ignoring option \'calibration = TRUE\'"))
}
if (!calibration && (length(irv) > ss_eq_no)) {
warning(paste0("solution was searched for with calibrating equations taken into account, ",
"ignoring option \'calibration = FALSE\'"))
}
} else {
prepstat <- tryCatch(prep <- ss_prep(model, use_jac = FALSE, calibration),
error = function(w) w)
if (inherits(prepstat, "error")) {
stop("evaluation of steady-state equations\' residuals has not been successful:", prepstat)
}
if (inherits(prepstat, "warning")) {
warning(prepstat)
}
irv <- prep$resid
}
eqnames <- paste0("Eq. ", 1:ss_eq_no)
if (length(irv) > ss_eq_no) {
eqnames <- c(eqnames, paste0("Calibr Eq. ", 1:calibr_eq_no))
}
names(irv) <- eqnames
cat("Initial residuals:\n")
print(round(irv, digits = 3))
tops <- min(length(irv), largest)
cat("\nEquations with the largest initial residuals:\n")
irvo <- eqnames[as.numeric(order(abs(irv), decreasing = TRUE)[1:tops])]
cat(list2str(irvo, ""))
cat("\n")
output$initial <- irv
if (length(model@residual_vector)) {
rv <- model@residual_vector
names(rv) <- eqnames
cat("\nFinal residuals:\n")
print(round(rv, digits = 3))
cat("\nEquations with the largest final residuals:\n")
rvo <- eqnames[as.numeric(order(abs(rv), decreasing = TRUE)[1:tops])]
cat(list2str(rvo, ""))
cat("\n")
output$final <- rv
}
return(invisible(output))
}
# ############################################################################
# Function ss_prep performs initial checks and constructs inputs for ss_find
# ############################################################################
# Input
# model - an object of gecon_model class
# use_jac - option to use Jacobian generated by symbolic library,
# if FALSE numerical derivatives are computed
# calibration - if TRUE, calibrating equations are taken
# into account while solving for the steady state
# or equilibrium and parameters are calibrated.
# Output
# 4-element list with inputs to ss_find and initial residuals
# ############################################################################
ss_prep <- function(model, use_jac = TRUE, calibration = TRUE)
{
if (!calibration &&
(length(model@parameters_free) != length(model@parameters_calibr)))
model@is_calibrated <- FALSE
if (any(is.nan(model@parameters_free_val) |
is.infinite(model@parameters_free_val))) {
nan_ind <- which(is.nan(model@parameters_free_val))
inf_ind <- which(is.infinite(model@parameters_free_val))
nans <- model@parameters_free[nan_ind]
infs <- model@parameters_free[inf_ind]
misspec <- c(nans, infs)
len <- (length(nans) + length(infs))
labs <- rep("Inf", len)
if (length(nans)) labs[1:length(nans)] <- "NaN"
stop("values of the following free parameters have not been specified correctly: ",
paste(paste0("\"", misspec, "\"", paste0(" (=", labs, ")")),
collapse = ", "))
}
if (any(is.na(model@parameters_free_val))) {
if (model@is_dynamic) { sseq <- "steady state"; } else { sseq <- "equilibrium"; }
mes <- paste0("values of the following free parameters have not been specified ",
"(all values have to be set before searching for ", sseq, "): ",
list2str(model@parameters_free[which(is.na(model@parameters_free_val))], "\""))
stop(mes)
}
if (!calibration && any(is.na(model@parameters_calibr_val))) {
mes <- paste0("calibration is turned off but values of the following calibrated parameters ",
"have not been specified (through a call to \'initval_calibr_par\'): ",
list2str(model@parameters_calibr[which(is.na(model@parameters_calibr_val))], "\""))
stop(mes)
}
if (!calibration && length(model@parameters_calibr)) {
mes <- paste0("calibration is turned off, initial values of calibrated parameters ",
"will be treated as their fixed values")
warning(mes)
}
# initial values
calibr_pars <- model@parameters_calibr_val
calibr_pars[which(is.na(calibr_pars))] <- 0.5
vars <- model@variables_ss_val
vars[which(is.na(vars))] <- 0.9
nv <- length(model@variables)
np <- length(model@parameters_calibr)
nvp <- nv + np
pars <- model@parameters_free_val
# construct initial vector and ss_eval
if (calibration) {
ss_guess <- c(vars, calibr_pars)
ss_eval <- function (x)
{
v <- x[1:nv]
cp <- x[(nv + 1):nvp]
return (c(model@ss_function(v, cp, pars),
model@calibr_function(v, cp, pars)))
}
} else {
ss_guess <- c(vars)
ss_eval <- function (x)
{
return (model@ss_function(x, calibr_pars, pars))
}
}
# Jacobian
if (use_jac) {
if (is.null(model@ss_calibr_jac_function)) {
warning(paste0("Jacobian will not be used to find steady state ",
"(Jacobian derivation was disabled in the .gcn file)"))
ss_jac <- NULL
} else {
if (calibration) {
ss_jac <- function(x)
{
v <- x[1:nv]
cp <- x[(nv + 1):nvp]
return (model@ss_calibr_jac_function(v, cp, pars))
}
} else {
ss_jac <- function(x)
{
return (model@ss_calibr_jac_function(x, calibr_pars, pars)[1:nv, 1:nv])
}
}
}
} else {
ss_jac <- NULL
}
# trial evaluation
residstat <- tryCatch(resid <- ss_eval(ss_guess),
error = function(w) w)
if (inherits(residstat, "error")) {
stop("indeterminate values")
}
if (any(!is.finite(resid))) {
eqnames <- paste0("Eq. ", 1:nv)
if (length(resid) > nv) {
eqnames <- c(eqnames, paste0("Calibr Eq. ", 1:np))
}
stop(paste0("non-finite residuals in equation(s): ",
list2str(eqnames[which(!is.finite(resid))], "")))
}
return (list(ss_guess = ss_guess, ss_eval = ss_eval, ss_jac = ss_jac, resid = resid))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.