R/class_model_ss.R

Defines functions steady_state get_residuals ss_prep

Documented in get_residuals steady_state

# ############################################################################
# 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))
}

Try the gEcon package in your browser

Any scripts or data that you put into this service are public.

gEcon documentation built on May 2, 2019, 6:52 p.m.