diff --git a/.Rbuildignore b/.Rbuildignore index 0db7381..442bcbb 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -9,3 +9,4 @@ ^to-do.md$ ^cran-comments.md$ ^CRAN-RELEASE$ +^development_test.R$ \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION index ee1d74c..3363496 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: sovereign Title: State-Dependent Empirical Analysis -Version: 1.1.0 +Version: 1.2.0 Authors@R: person(given = "Tyler J.", family = "Pike", @@ -19,7 +19,9 @@ Roxygen: list(markdown = TRUE) RoxygenNote: 7.1.1 Imports: broom, - dplyr, + dplyr, + future, + furrr, ggplot2, gridExtra, lmtest, diff --git a/NAMESPACE b/NAMESPACE index f3087a4..dcf1a6e 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,6 +1,9 @@ # Generated by roxygen2: do not edit by hand export("%>%") +export(FEVD) +export(HD) +export(IRF) export(LP) export(RLP) export(RVAR) diff --git a/NEWS.md b/NEWS.md index 43239dc..569de0b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -12,12 +12,18 @@ - no structural restrictions, - short-term restrictions (via Cholesky decomposition), - instrumental variable estimation - - combination of instrument variable and short-term restricions - - The choice run bootstrapping routines in parallel for greater computational efficiency + - combination of instrument variable and short-term restrictions + - Parallel processing for greater computational efficiency in computing bootstrapped confidence intervals + +- Minor updates + - IRF plots use darker blue for CI + - Analysis function wrappers implemented for IRFs, FEVDs, and HDs + - Updated model documentation with academic references and corrected minor errors - Bug fixes - plot_irf() plots all targets by default - - test-threvhold_var.R changed to test-threshold_var.R + - test-threvhold_var.R changed to test-regime_var.R + - Fixed IRF results format listed in documentation ## Version 1.1.0 (2021-06-01) @@ -33,7 +39,7 @@ - Clarify model specification of RVAR - Bug fixes - - Convert errorConditiont() to stop() + - Convert errorCondition() to stop() ## Version 1.0.0 (2021-04-02) diff --git a/R/analysis_wrappers.R b/R/analysis_wrappers.R new file mode 100644 index 0000000..e0658fc --- /dev/null +++ b/R/analysis_wrappers.R @@ -0,0 +1,264 @@ +#' Estimate impulse response functions +#' +#' See VAR, RVAR, LP, and RLP documentation for details +#' regarding models and structural errors. +#' +#' @param model VAR, RVAR, LP, or RLP class object +#' @param horizon int: number of periods +#' @param CI numeric vector: c(lower ci bound, upper ci bound) +#' @param bootstrap.type string: bootstrapping technique to use ('auto', 'standard', or 'wild'); if auto then wild is used for IV or IV-short, else standard is used +#' @param bootstrap.num int: number of bootstraps +#' @param bootstrap.parallel boolean: create IRF draws in parallel +#' @param bootstrap.cores int: number of cores to use in parallel processing; -1 detects and uses half the available cores +#' +#' @return data frame with columns `target`, `shock`, `horizon`, `response.lower`, `response`, `response.upper`; regime-based models return a list with a data frame per regime. +#' +#' @seealso [var_irf()] +#' @seealso [rvar_irf()] +#' @seealso [lp_irf()] +#' @seealso [rlp_irf()] +#' +#' @examples +#' \donttest{ +#' +#' # simple time series +#' AA = c(1:100) + rnorm(100) +#' BB = c(1:100) + rnorm(100) +#' CC = AA + BB + rnorm(100) +#' date = seq.Date(from = as.Date('2000-01-01'), by = 'month', length.out = 100) +#' Data = data.frame(date = date, AA, BB, CC) +#' +#' # estimate VAR +#' var = +#' sovereign::VAR( +#' data = Data, +#' horizon = 10, +#' freq = 'month', +#' lag.ic = 'BIC', +#' lag.max = 4 +#' ) +#' +#' # impulse response function +#' var.irf = sovereign::IRF(var) +#' +#' # local projection forecasts +#' lp = +#' sovereign::LP( +#' data = Data, +#' horizon = c(1:10), +#' lag.ic = 'AIC', +#' lag.max = 4, +#' type = 'both', +#' freq = 'month') +#' +#' # LP impulse response function +#' lp.irf = sovereign::IRF(lp) +#' +#' } +#' +#' @export + +IRF = function( + model, # VAR, RVAR, LP, or RLP class object + horizon = 10, # int: number of periods + CI = c(0.1, 0.9), # numeric vector: c(lower ci bound, upper ci bound) + bootstrap.type = 'auto', # string: bootstrapping technique to use ('auto', 'standard', or 'wild'); if auto then wild is used for IV or IV-short, else standard is used + bootstrap.num = 100, # int: number of bootstraps + bootstrap.parallel = FALSE, # boolean: create IRF draws in parallel + bootstrap.cores = -1 # int: number of cores to use in parallel processing; -1 detects and uses half the available cores +){ + + if(class(model) == 'VAR'){ + return( + var_irf( + var = model, + horizon = horizon, + CI = CI, + bootstrap.type = bootstrap.type, + bootstrap.num = bootstrap.num, + bootstrap.parallel = bootstrap.parallel, + bootstrap.cores = bootstrap.cores + ) + ) + }else if(class(model) == 'RVAR'){ + return( + rvar_irf( + rvar = model, + horizon = horizon, + CI = CI, + bootstrap.type = bootstrap.type, + bootstrap.num = bootstrap.num, + bootstrap.parallel = bootstrap.parallel, + bootstrap.cores = bootstrap.cores + ) + ) + }else if(class(model) == 'LP'){ + return( + lp_irf( + lp = model, + CI = CI + ) + ) + }else if(class(model) == 'RLP'){ + return( + rlp_irf( + rlp = model, + CI = CI + ) + ) + }else{ + stop('model must be a sovereign VAR, RVAR, LP, or RLP class object') + } + +} + +#' Estimate forecast error variance decomposition +#' +#' Estimate the forecast error variance decomposition for VARs with +#' either short or 'IV-short' structural errors. See VAR +#' and RVAR documentation for details regarding structural errors. +#' +#' @param model VAR or RVAR class object +#' @param horizon int: number of periods +#' @param scale boolean: scale variable contribution as percent of total error +#' +#' @return long-form data.frame +#' +#' @seealso [VAR()] +#' @seealso [var_fevd()] +#' @seealso [RVAR()] +#' @seealso [rvar_fevd()] +#' +#' @examples +#' \donttest{ +#' +#' # simple time series +#' AA = c(1:100) + rnorm(100) +#' BB = c(1:100) + rnorm(100) +#' CC = AA + BB + rnorm(100) +#' date = seq.Date(from = as.Date('2000-01-01'), by = 'month', length.out = 100) +#' Data = data.frame(date = date, AA, BB, CC) +#' +#' # estimate VAR +#' var = +#' sovereign::VAR( +#' data = Data, +#' horizon = 10, +#' freq = 'month', +#' lag.ic = 'BIC', +#' lag.max = 4) +#' +#' # impulse response functions +#' var.irf = sovereign::IRF(var) +#' +#' # forecast error variance decomposition +#' var.fevd = sovereign::FEVD(var) +#' +#' # historical shock decomposition +#' var.hd = sovereign::HD(var) +#' +#' } +#' +#' @export + +FEVD = function( + model, # VAR or RVAR class object + horizon = 10, # int: number of periods + scale = TRUE # boolean: scale variable contribution as percent of total error +){ + + if(class(model) == 'VAR'){ + return( + var_fevd( + var = model, + horizon = horizon, + scale = scale + ) + ) + }else if(class(model) == 'RVAR'){ + return( + rvar_fevd( + rvar = model, + horizon = horizon, + scale = scale + ) + ) + }else{ + stop('model must be a sovereign VAR or RVAR class object') + } + +} + + +#------------------------------------------------------- +# Function calculate historical decomposition of shocks +#------------------------------------------------------- + +#' Estimate historical decomposition +#' +#' Estimate the historical decomposition for VARs with +#' either 'short' or 'IV-short' structural errors. See VAR +#' and RVAR documentation for details regarding structural errors. +#' +#' @param model VAR or RVAR class object +#' +#' @return long-from data.frame +#' +#' @seealso [VAR()] +#' @seealso [var_hd()] +#' @seealso [RVAR()] +#' @seealso [rvar_hd()] +#' +#' @examples +#' \donttest{ +#' +#' # simple time series +#' AA = c(1:100) + rnorm(100) +#' BB = c(1:100) + rnorm(100) +#' CC = AA + BB + rnorm(100) +#' date = seq.Date(from = as.Date('2000-01-01'), by = 'month', length.out = 100) +#' Data = data.frame(date = date, AA, BB, CC) +#' +#' # estimate VAR +#' var = +#' sovereign::VAR( +#' data = Data, +#' horizon = 10, +#' freq = 'month', +#' lag.ic = 'BIC', +#' lag.max = 4) +#' +#' # impulse response functions +#' var.irf = sovereign::IRF(var) +#' +#' # forecast error variance decomposition +#' var.fevd = sovereign::FEVD(var) +#' +#' # historical shock decomposition +#' var.hd = sovereign::HD(var) +#' +#' } +#' +#' @export + +HD = function( + model # VAR or RVAR class object +){ + + if(class(model) == 'VAR'){ + return( + var_hd( + var = model + ) + ) + }else if(class(model) == 'RVAR'){ + return( + rvar_hd( + rvar = model + ) + ) + }else{ + stop('model must be a sovereign VAR or RVAR class object') + } + +} diff --git a/R/development_test.R b/R/development_test.R deleted file mode 100644 index aef04af..0000000 --- a/R/development_test.R +++ /dev/null @@ -1,803 +0,0 @@ - -TVAR = function( - data, # data.frame, matrix, ts, xts, zoo: Endogenous regressors - threshold, # string: variable to estimate threshold values - horizon = 10, # int: forecast horizons - freq = 'month', # string: frequency of data (day, week, month, quarter, year) - type = 'const', # string: type of deterministic terms to add ('none', 'const', 'trend', 'both') - p = 1, # int: lags - lag.ic = NULL, # string: information criterion to choose the optimal number of lags ('AIC' or 'BIC') - lag.max = NULL, # int: maximum number of lags to test in lag selection - structure = 'short', # string: type of structural identification strategy to use in model analysis (NULL, 'short', or 'IV') - instrument = NULL, # string: name of instrumental variable contained in the data matrix - instrumented = NULL # string: name of variable to be instrumented in IV and IV-short procedure; default is the first non-date variable in data -){ - - # function warnings - if(!is.numeric(p) | p %% 1 != 0){ - stop('p must be an integer') - } - if(!is.null(lag.ic)){ - if(!lag.ic %in% c('BIC','AIC')){ - stop("lag.ic must be either 'BIC', 'AIC', or NULL") - } - } - if(!is.null(lag.max)){ - if(lag.max %% 1 != 0){ - stop('lag.max must be an integer if IC-based lag selection is used') - } - } - if(!is.matrix(data) & !is.data.frame(data)){ - stop('data must be a matrix or data.frame') - } - if(!is.numeric(p) | p %% 1 != 0){ - stop('p must be an integer') - } - if(!is.numeric(horizon) | horizon %% 1 != 0 | horizon <= 0){ - stop('horizon must be a positive integer') - } - if(!freq %in% c('day','week','month','quarter','year')){ - stop("freq must be one of the following strings: 'day','week','month','quarter','year'") - } - if(!structure %in% c('short', 'IV', 'IV-short') & !is.null(structure)){ - stop("strucutre must be one of 'strucutre', 'IV', 'IV-short', or NULL.") - } - if(!is.null(instrument)){ - if(!instrument %in% colnames(data)){ - stop("instrument must be the name of a variable found in data.") - } - } - if(!is.null(instrumented)){ - if(!instrumented %in% colnames(data)){ - stop("instrumented must be the name of a variable found in data.") - } - } - - # cast as data frame if ts, xts, or zoo object - if(stats::is.ts(data) | xts::is.xts(data) | zoo::is.zoo(data)){ - data = data.frame(date = zoo::index(date), data) - } - - # set aside instruments - if(!is.null(instrument)){ - data.instrument = dplyr::select(data, date, dplyr::all_of(instrument)) - data = dplyr::select(data, -dplyr::all_of(instrument)) - }else{ - data.instrument = NULL - } - - # detect variable to be instrumented - if(is.null(instrumented)){ - var_to_instrument = colnames(dplyr::select(data, -date))[1] - }else{ - var_to_instrument = instrumented - } - - - for(t in unique(data[,c(threshold)])){ - - data.threshold = data %>% - mutate(t_regime = if_else(threshold >= t, 1, 0)) - - rvar = - RVAR( - data = data.threshold, - horizon = horizon, - freq = freq, - type = type, - p = p, - lag.ic = lag.ic, - lag.max = lag.max, - regime = t_regime, - ) - - rvar.ll = rvar$ll - - } - - -} - - -VAR = function( - data, # data.frame, matrix, ts, xts, zoo: Endogenous regressors - horizon = 10, # int: forecast horizons - freq = 'month', # string: frequency of data (day, week, month, quarter, year) - type = 'const', # string: type of deterministic terms to add ('none', 'const', 'trend', 'both') - p = 1, # int: lags - lag.ic = NULL, # string: information criterion to choose the optimal number of lags ('AIC' or 'BIC') - lag.max = NULL, # int: maximum number of lags to test in lag selection - structure = 'short', # string: type of structural identification strategy to use in model analysis (NULL, 'short', or 'IV') - instrument = NULL, # string: name of instrumental variable contained in the data matrix - instrumented = NULL # string: name of variable to be instrumented in IV and IV-short procedure; default is the first non-date variable in data -){ - - # function warnings - if(!is.numeric(p) | p %% 1 != 0){ - stop('p must be an integer') - } - if(!is.null(lag.ic)){ - if(!lag.ic %in% c('BIC','AIC')){ - stop("lag.ic must be either 'BIC', 'AIC', or NULL") - } - } - if(!is.null(lag.max)){ - if(lag.max %% 1 != 0){ - stop('lag.max must be an integer if IC-based lag selection is used') - } - } - if(!is.matrix(data) & !is.data.frame(data)){ - stop('data must be a matrix or data.frame') - } - if(!is.numeric(p) | p %% 1 != 0){ - stop('p must be an integer') - } - if(!is.numeric(horizon) | horizon %% 1 != 0 | horizon <= 0){ - stop('horizon must be a positive integer') - } - if(!freq %in% c('day','week','month','quarter','year')){ - stop("freq must be one of the following strings: 'day','week','month','quarter','year'") - } - if(!structure %in% c('short', 'IV', 'IV-short') & !is.null(structure)){ - stop("strucutre must be one of 'strucutre', 'IV', 'IV-short', or NULL.") - } - if(!is.null(instrument)){ - if(!instrument %in% colnames(data)){ - stop("instrument must be the name of a variable found in data.") - } - } - if(!is.null(instrumented)){ - if(!instrumented %in% colnames(data)){ - stop("instrumented must be the name of a variable found in data.") - } - } - - # cast as data frame if ts, xts, or zoo object - if(stats::is.ts(data) | xts::is.xts(data) | zoo::is.zoo(data)){ - data = data.frame(date = zoo::index(date), data) - } - - # set aside instruments - if(!is.null(instrument)){ - data.instrument = dplyr::select(data, date, dplyr::all_of(instrument)) - data = dplyr::select(data, -dplyr::all_of(instrument)) - }else{ - data.instrument = NULL - } - - # detect variable to be instrumented - if(is.null(instrumented)){ - var_to_instrument = colnames(dplyr::select(data, -date))[1] - }else{ - var_to_instrument = instrumented - } - - # VAR estimation - if(!is.null(lag.ic)){ - - ic.scores = vector(length = lag.max+1) - - models = c(1:lag.max) %>% - purrr::map(.f = function(p){ - - # estimate candidate model - model = - VAR_estimation( - data = data, - p = p, - horizon = horizon, - freq = freq, - type = type - ) - - # calculate IC - ic.score = - IC( - ic = lag.ic, - errors = model$residuals[[1]], - data = data, - p = p - ) - - ic.scores[p] = ic.score - - # return candidate model - return(model) - - }) - - # return IC minimizing VAR - min.ic = which.min(ic.scores) - model = models[[min.ic]] - - }else{ - - model = - VAR_estimation( - data = data, - p = p, - horizon = horizon, - freq = freq, - type = type - ) - - } - - # add structure - model$structure = structure - model$instrument = data.instrument - model$instrumented = var_to_instrument - - # assign class and return - class(model) = 'VAR' - return(model) - -} - -##################################################################### -# IRF FUNCTIONS - -solve_B = function(var, report_iv = TRUE){ - - if(is.null(var$structure) == TRUE){ - - # retrieve reduced-form residuals - data.residuals = var$residuals[[1]] - - # reduced form variance-covariance matrix - cov.matrix = stats::var(stats::na.omit(dplyr::select(data.residuals, -date, -forecast.date))) - - B = cov.matrix - - return(B) - - }else if(var$structure == 'short'){ - - # retrieve reduced-form residuals - data.residuals = var$residuals[[1]] - - # reduced form variance-covariance matrix - cov.matrix = stats::var(stats::na.omit(dplyr::select(data.residuals, -date, -forecast.date))) - - # take cholesky decomposition - B = t(chol(cov.matrix)) - - return(B) - - }else if(var$structure == 'IV'){ - - # retrieve reduced-form residuals - data.residuals = var$residuals[[1]] - covariates = colnames(dplyr::select(data.residuals, -date, -forecast.date)) - col_to_instrument = var$instrumented - - # retrieve instrument - data.instrument = var$instrument - instrument = colnames(dplyr::select(data.instrument, -date)) - - # combine data - data = - dplyr::inner_join(data.residuals, data.instrument, by = 'date') %>% - na.omit() - - # first stage least squares - model.first_stage = - lm(col_to_instrument ~., - data = data %>% - select(col_to_instrument = dplyr::all_of(col_to_instrument), all_of(instrument)) %>% - na.omit()) - p_hat = model.first_stage$fitted.values - - if(report_iv == TRUE){ - print(summary(model.first_stage)) - } - - # second stage least squares - # to estimate first column of B - # (automatically scales first entry to 1) - instrumented_shocks = covariates %>% - purrr::map_df(.f = function(X){ - - model.second_stage = lm(data[,X] ~ p_hat) - second_stage_beta = model.second_stage$coefficients[2] - return(second_stage_beta) - - }) %>% - as.vector() - - # scale size of the shock - # see Gertler and Karadi (2015) for background - # see Cesa-Bianchi's VAR-Toolbox for MATLAB implementation - X.demean = - select(data, -date, -forecast.date, -dplyr::all_of(instrument)) %>% - mutate_all(function(X){return(X - mean(X, na.rm = T))}) %>% - as.matrix() - - sigma_b = 1/(nrow(data) - ncol(var$model$coef) + 1) * t(X.demean) %*% X.demean - - s21s11 = instrumented_shocks[2:length(covariates),] %>% as.matrix() - S11 = sigma_b[1,1] %>% as.numeric() - S21 = sigma_b[2:nrow(sigma_b),1] %>%as.vector() - S22 = sigma_b[2:nrow(sigma_b),2:ncol(sigma_b)] - - Q = (s21s11 * S11) %*% t(s21s11) - (S21 %*% t(s21s11) + as.matrix(s21s11) %*% t(S21)) + S22 - sp = sqrt( S11 - t(S21 - as.matrix(s21s11) %*% S11) %*% as.matrix(solve(Q) %*% (S21 - s21s11 * S11)) ) - - scaled_instrumented_shocks = instrumented_shocks * as.numeric(sp) - - # prepare B matrix - B = matrix(0, ncol = (length(covariates) - 1), nrow = length(covariates)) - B = cbind(scaled_instrumented_shocks, B) - - return(B) - - }else if(var$structure == 'IV-short'){ - - # align instrument and residuals - valid_dates = - dplyr::inner_join( - dplyr::select(na.omit(var$residuals[[1]]), date), - dplyr::select(na.omit(var$instrument), date), - by = 'date' - ) - - # set residuals matrix - residuals = var$residuals[[1]] %>% - dplyr::inner_join(valid_dates, by = 'date') %>% - dplyr::select(-date, -forecast.date) %>% - as.matrix() - - residuals = -1*residuals - - # set instrument - instrument = var$instrument %>% - data.frame() %>% - dplyr::inner_join(valid_dates, by = 'date') %>% - dplyr::select(-date) - - instrument.mean = instrument %>% - dplyr::mutate_all(mean, na.rm = T) - - # number of observations - n.obs = nrow(var$data) - - # solve for IV implied impact - pi = solve(n.obs^(-1) * t(residuals) %*% residuals) %*% - (n.obs^(-1) * t(residuals) %*% as.matrix(instrument - instrument.mean)) - - phi_sq = - (n.obs^(-1) * t(instrument - instrument.mean) %*% residuals) %*% - solve( n.obs^(-1) * t(residuals) %*% residuals ) %*% - ( n.obs^(-1) * t(residuals) %*% as.matrix(instrument - instrument.mean) ) - - B1 = - n.obs^(-1) * - ( t(residuals) %*% as.matrix(instrument - instrument.mean) ) %*% - (phi_sq)^(-0.5) - - B = matrix(ncol = ncol(residuals), nrow = ncol(residuals)) - B[,1:ncol(instrument)] = B1 - rownames(B) = colnames(B) = colnames(residuals) - - # solve for orthogonalized structural shock - model.first_stage = lm(instrument[,1] ~ residuals) - orthogonal_instrument = instrument - model.first_stage$residuals - orthogonal_instrument = orthogonal_instrument[,] / sd(orthogonal_instrument[,], na.rm = T) - - shock_matrix = matrix(nrow = nrow(residuals), ncol = ncol(residuals)) - shock_matrix[,1] = orthogonal_instrument - - # reduced form variance-covariance matrix - sigma = stats::var(residuals) - - # solve additional entries - # with a lower triangular restriction - order_sequence = c(1:ncol(residuals)) - - for(i in order_sequence){ - - Y = residuals[,i] - X = shock_matrix[,c(1:i)] - model.second_stage = lm(Y~X) - - if(i != tail(order_sequence,1)){ - - B[i,] = - c( model.second_stage$coef[-1], - sd(model.second_stage$residuals), - rep(0, length(order_sequence) - ncol(instrument) - i) ) - - shock_matrix[,i+1] = model.second_stage$residuals / sd(model.second_stage$residuals) - - }else{ - - B[i,] = model.second_stage$coef[-1] - - } - - } - - # make diagonal entries positive - shock_ordering = data.frame(residuals) %>% - select(var$instrumented, dplyr::everything()) %>% - colnames() - B.sign = diag(sign(diag(B[shock_ordering,]))) - B = B %*% B.sign - - return(B) - - }else{ - - stop('structure must be set to either NULL, "short", or "IV".') - - } - -} - -IRF = function (Phi, B, lag, structure = NULL){ - - # cast data as matrices - if (!is.matrix(Phi)) - Phi = as.matrix(Phi) - if (!is.matrix(B)) - B = as.matrix(B) - - - # set dimensions - k = nrow(Phi) - m = ncol(Phi) - p = floor(m/k) - Si = diag(rep(1, k)) - wk = c(Si) - awk = c(wk) - acuwk = c(awk) - - # set lag - if (p < 1) - p = 1 - if (lag < 1) - lag = 1 - - # estimate IRFs - for (i in 1:lag) { - if (i <= p) { - idx = (i - 1) * k - tmp = Phi[, (idx + 1):(idx + k)] - } - else { - tmp = matrix(0, k, k) - } - jj = i - 1 - jp = min(jj, p) - if (jp > 0) { - for (j in 1:jp) { - jdx = (j - 1) * k - idx = (i - j) * k - w1 = Phi[, (jdx + 1):(jdx + k)] - w2 = Si[, (idx + 1):(idx + k)] - tmp = tmp + w1 %*% w2 - } - } - Si = cbind(Si, tmp) - wk = cbind(wk, c(tmp)) - awk = awk + c(tmp) - acuwk = cbind(acuwk, awk) - } - orSi = NULL - wk1 = NULL - awk1 = NULL - acuwk1 = NULL - - irf = Si - - if (!is.null(structure)) { - - P = B - wk1 = cbind(wk1, c(P)) - awk1 = wk1 - acuwk1 = wk1 - orSi = cbind(orSi, P) - for (i in 1:lag) { - idx = i * k - w1 = Si[, (idx + 1):(idx + k)] - w2 = w1 %*% P - orSi = cbind(orSi, w2) - wk1 = cbind(wk1, c(w2)) - awk1 = awk1 + c(w2) - acuwk1 = cbind(acuwk1, awk1) - } - - irf = orSi - - } - - # return results - return(irf) - -} - -var_irf = function( - var, # VAR output - horizon = 10, # int: number of periods - CI = c(0.1, 0.9), # numeric vector: c(lower ci bound, upper ci bound) - bootstrap.type = 'auto', # string: bootstrapping technique to use ('auto', 'standard', or 'wild'); if auto then wild is used for IV or IV-short, else standard is used - bootstrap.num = 100, # int: number of bootstraps - bootstrap.parallel = FALSE, # boolean: create IRF draws in parallel - bootstrap.cores = -1 # int: number of cores to use in parallel processing; -1 detects and uses half the available cores -){ - - # NOTES: scaled wild is used for IV, consistent with the proxy SVAR literature, - # while standard resample is used for others - - # function warnings - if(!is.numeric(bootstrap.num) | bootstrap.num %% 1 != 0){ - stop('bootstrap.num must be an integer') - } - if(!is.numeric(CI) | length(CI) != 2 | min(CI) < 0 | max(CI) > 1 | is.na(sum(CI))){ - stop('CI must be a two element numeric vector bound [0,1]') - } - if(!is.numeric(horizon) | horizon %% 1 != 0 | horizon <= 0){ - stop('horizon must be a positive integer') - } - - # function variables - y = forecast.date = shock = target = response = response.lower = response.upper = response.med = response.adjust = NULL - - # set data - coef = var$model$coef - residuals = var$residuals[[1]] - data = var$data - - # set model variables - p = var$model$p - freq = var$model$freq - structure = var$structure - - if('const' %in% colnames(coef) & 'trend' %in% colnames(coef)){ - type = 'both' - }else if('const' %in% colnames(coef)){ - type = 'const' - }else if('trend' %in% colnames(coef)){ - type = 'trend' - } - - regressors = colnames(dplyr::select(data, -date)) - if(!is.null(var$regime)){regressors = regressors[regressors != var$regime$regime]} - regressors.cols = paste0(regressors, '.l1') - - # set IRF variables - p.lower = CI[1] - p.upper = CI[2] - - ### calculate impulse responses -------------- - - # estimate error-covariance matrix or structural impact matrix - B = solve_B(var) - B = as.matrix(B) - - # estimate IRFs - phi = dplyr::select(coef, -y, -dplyr::contains('cosnt'), -dplyr::contains('trend')) - - irf = IRF(Phi = phi, - B = B, - lag = horizon, - structure = structure) - - # reorganize results - irf = data.frame(t(irf)) - - if(!is.null(structure)){ - if(structure == 'IV-short'){ - shocks = c(var$instrumented, regressors[!regressors %in% var$instrumented]) - irf$shock = rep(shocks, horizon + 1) - } - }else{ - irf$shock = rep(regressors, horizon + 1) - } - - irf$horizon = sort(rep(c(0:horizon), length(regressors))) - irf = dplyr::arrange(irf, shock, horizon) - rownames(irf) = NULL - colnames(irf) = c(regressors, 'shock', 'horizon') - - ### bootstrap IRF confidence intervals --------- - - # 0. set up parallel back-end - if(bootstrap.parallel == TRUE){ - if(bootstrap.cores == -1){ - n.cores = floor(future::availableCores() / 2) - }else{ - n.cores = bootstrap.parallel - } - future::plan(future::multisession, workers = n.cores) - }else{ - future::plan(future::sequential) - } - - # 1. create bootstrap time series - bagged.irf = as.list(1:bootstrap.num) %>% - furrr::future_map(.f = function(count){ - - # bootstrap residuals - if(bootstrap.type == 'wild' | - bootstrap.type == 'auto' & structure == 'IV' | - bootstrap.type == 'auto' & structure == 'IV-short'){ - - # 'wild' bootstrap technique for simple distribution - # using observed scaled residuals with random sign flip. - # See the Rademacher distribution. - U = residuals[,-c(1,2)] - r = sample(c(-1,1), size = nrow(U), replace = T) - U = sweep(U, MARGIN = 1, r, `*`) - - }else if(bootstrap.type == 'standard' | bootstrap.type == 'auto' & structure != 'IV'){ - - # standard bootstrap technique a al Lutkepohl (2005) - # draw residuals with replacement - U = na.omit(residuals) - U = U[sample(c(1:nrow(U)), - size = nrow(residuals), - replace = TRUE),] - U = U %>% - dplyr::select(-date, -forecast.date) %>% - dplyr::mutate_all(function(X){return(X-mean(X, na.rm = T))}) - } - - # recursively build synthetic data - Y = matrix(nrow = nrow(var$data), ncol = ncol(var$data)-1) - Y = data.frame(Y); colnames(Y) = regressors - Y[1:p, ] = var$data[1:p, -1] - - for(i in (p+1):nrow(var$data)){ - - X = Y[(i-p):(i-1),] - X = embed(as.matrix(X), dimension = p) - X = data.frame(X) - - if(type %in% c('const', 'both')){X$const = 1} - if(type %in% c('trend', 'both')){X$trend = c(1:nrow(X))} - - X.hat = as.matrix(coef[,-1]) %*% t(as.matrix(X)) - X.hat = t(X.hat) - - Y[i, ] = X.hat - U[i,] - - } - - Y$date = data$date - - # estimate VAR with synthetic data - var.bag = - VAR( - data = Y, - p = p, - horizon = 1, - freq = freq, - type = type, - structure = structure) - - # bootstrap instrument - if(!is.null(structure)){ - if(structure == 'IV' | structure == 'IV-short'){ - var.bag$instrument = var$instrument - var.bag$instrument[,-1] = sweep(var.bag$instrument[,-1], MARGIN = 1, r, `*`) - - var.bag$instrumented = var$instrumented - } - } - - # estimate error-covariance matrix or structural impact matrix - B.bag = solve_B(var.bag, report_iv = FALSE) - - # set bagged coef matrix - coef.bag = dplyr::select(var.bag$model$coef, -y, -dplyr::contains('const'), -dplyr::contains('trend')) - - # estimate IRFs - irf.bag = IRF(Phi = coef.bag, - B = B.bag, - lag = horizon, - structure = structure) - - # reorganize results - irf.bag = data.frame(t(irf.bag)) - - if(!is.null(structure)){ - if(structure == 'IV-short'){ - shocks = c(var$instrumented, regressors[!regressors %in% var$instrumented]) - irf.bag$shock = rep(shocks, horizon + 1) - } - }else{ - irf.bag$shock = rep(regressors, horizon + 1) - } - - irf.bag$horizon = sort(rep(c(0:horizon), length(regressors))) - irf.bag = dplyr::arrange(irf.bag, shock, horizon) - rownames(irf.bag) = NULL - colnames(irf.bag) = c(regressors, 'shock', 'horizon') - - return(irf.bag) - - }) - - # 3. calculate confidence intervals - - # collapse to data.frame - ci = bagged.irf %>% - purrr::reduce(dplyr::bind_rows) - - # remove unused shocks in the case of IV - if(!is.null(var$structure)){ - if(var$structure == 'IV'){ - ci = ci %>% - dplyr::filter(shock == regressors[1]) - } - } - - # correct shock names in the case of IV-short - if(!is.null(var$structure)){ - if(var$structure == 'IV-short'){ - - } - } - - # estimate bands - ci = ci %>% - tidyr::pivot_longer(cols = regressors, names_to = 'target', values_to = 'response') %>% - group_by(shock, horizon, target) %>% - summarize(response.lower = quantile(response, probs = p.lower, na.rm = T), - response.upper = quantile(response, probs = p.upper, na.rm = T), - response.mean = mean(response, na.rm = T) ) - - # combine point estimates and CI - irf = irf %>% - tidyr::pivot_longer(cols = regressors, names_to = 'target', values_to = 'response') %>% - dplyr::arrange(shock, horizon) - - irf = - full_join( - irf, ci, - by = c('shock', 'target', 'horizon') - ) - - # adjust for bias in CI - # (bias can be introduced in the bootstrapping if residuals are not actually mean zero) - irf = irf %>% - dplyr::mutate( - response.adjust = response - response.mean, - response.lower = response.lower + response.adjust, - response.upper = response.upper + response.adjust - ) %>% - dplyr::select(target, shock, horizon, response.lower, response, response.upper) %>% - dplyr::arrange(target, shock, horizon) - - ### return output -------------- - return(irf) -} - - -# ##################################################################### -# TEST - -library(tidyverse) -library(lubridate) - -# Data -data.macro = read_csv('/scratch/m1tjp01/Andrea/FinancialStability/Data/Intermediate/model_monthly_macro_covariates.csv') -data.macro = data.macro %>% filter(year(date) >= 1991) -data.macro = data.macro %>% select(date, PCE, UNEMP, PR, EBP, mp_shock) - -# MODEL -var = - VAR( - data = data.macro, - horizon = 1, - freq = 'month', - p = 1, - structure = 'IV-short', - instrument = 'mp_shock', - instrumented = 'PR' - ) - -irf = var_irf(var, horizon = 20, bootstrap.type = 'wild', bootstrap.num = 10, bootstrap.parallel = TRUE) -plot_irf(irf, shocks = 'PR') - diff --git a/R/helper.R b/R/helper.R new file mode 100644 index 0000000..2685d11 --- /dev/null +++ b/R/helper.R @@ -0,0 +1,98 @@ +#------------------------------------------ +# Data helper functions +#------------------------------------------ +# create n lags +n.lag = function( + Data, # data.frame: data frame of variables to lag and a 'date' column + lags, # int: number of lags to create + variables = NULL # string: vector of variable names to lag, default is all non-date variables +){ + + # function variables + date = NULL + + if(is.null(variables)){ + variables = names(dplyr::select(Data, -dplyr::contains('date'))) + } + + if(!'date' %in% colnames(Data)){ + no.date = TRUE + Data$date = c(1:nrow(Data)) + }else{ + no.date = FALSE + } + + Data = c(0:lags) %>% + purrr::map( + .f = function(n){ + + if(n == 0){return(Data)} + + X = Data %>% + dplyr::mutate_at(variables, dplyr::lag, n) + + names(X)[names(X) != 'date'] = paste0(names(X)[names(X) != 'date'], '.l', n) + + return(X) + } + ) %>% + purrr::reduce(dplyr::full_join, by = 'date') + + if(no.date == TRUE){ + Data = dplyr::select(Data, -date) + } + + return(Data) +} + +# adjust forecast dates +forecast_date = function( + forecast.date, + horizon, + freq +){ + + date = forecast.date + + if(freq == 'day'){ + date = forecast.date + horizon + }else if(freq == 'week'){ + lubridate::week(date) = lubridate::week(date) + horizon + }else if(freq == 'month'){ + lubridate::month(date) = lubridate::month(date) + horizon + }else if(freq == 'quarter'){ + lubridate::month(date) = lubridate::month(date) + horizon*3 + }else if(freq == 'year'){ + lubridate::year(date) = lubridate::year(date) + horizon + } + + return(date) +} + +#------------------------------------------ +# Information criterion +#------------------------------------------ +IC = function( + ic = 'AIC', + errors, + data, + p +){ + + cov.matrix = stats::var(stats::na.omit(dplyr::select(errors, -date))) + + regressors = colnames(dplyr::select(data, -date)) + regressors.cols = paste0(regressors, '.l1') + k = length(regressors.cols) + + sample.size = nrow(stats::na.omit(data)) + + if(ic == 'AIC'){ + score = log(sum(abs(cov.matrix))) + (2*k^2*p)/sample.size + }else if(ic == 'BIC'){ + score = log(sum(abs(cov.matrix))) + (log(sample.size)*k^2*p)/sample.size + } + + return(score) + +} diff --git a/R/lp.R b/R/lp.R index 153cdbd..10537c4 100644 --- a/R/lp.R +++ b/R/lp.R @@ -170,6 +170,9 @@ LP_estimate = function( #' @seealso [RLP()] #' @seealso [rlp_irf()] #' +#' @references +#' 1. Jorda, Oscar "[Estimation and Inference of Impulse Responses by Local Projections](https://www.aeaweb.org/articles?id=10.1257/0002828053828518)" 2005. +#' #' @examples #' \donttest{ #' @@ -530,6 +533,11 @@ RLP_estimate = function( #' Estimate regime-dependent local projections #' +#' Estimate a regime-dependent local projection (i.e. a state-dependent LP), with an exogenous state indicator, of the specification: +#' \deqn{Y_{t+h} = X_t \beta_{s_t} + \epsilon_t} +#' where *t* is the time index, and *s* is a mutually exclusive state of the world observed at time *t*. When the regime vector is +#' not supplied by the user, then a two-state regime series is estimated via random forest. +#' #' @param data data.frame, matrix, ts, xts, zoo: Endogenous regressors #' @param horizons int: forecast horizons #' @param freq string: frequency of data ('day', 'week', 'month', 'quarter', or 'year') @@ -544,13 +552,18 @@ RLP_estimate = function( #' @param regime.method string: regime assignment technique ('rf', 'kmeans', 'EM', 'BP') #' @param regime.n int: number of regimes to estimate (applies to kmeans and EM) #' -#' @return list object with elements `data`, `model`, `forecasts`, `residuals`; if there is more than one forecast horizon estimated, then `model`, `forecasts`, `residuals` will each be a list where each element corresponds to a single horizon +#' @return list of lists, one list per regime, each regime with objects with elements `data`, `model`, `forecasts`, `residuals`; +#' if there is more than one forecast horizon estimated, then `model`, `forecasts`, `residuals` +#' will each be a list where each element corresponds to a single horizon #' #' @seealso [LP()] #' @seealso [lp_irf()] #' @seealso [RLP()] #' @seealso [rlp_irf()] #' +#' @references +#' 1. Jorda, Oscar "[Estimation and Inference of Impulse Responses by Local Projections](https://www.aeaweb.org/articles?id=10.1257/0002828053828518)" 2005. +#' #' @examples #' \donttest{ #' diff --git a/R/plot_irf.R b/R/plot_irf.R index 1cf4c42..8acc1b8 100644 --- a/R/plot_irf.R +++ b/R/plot_irf.R @@ -34,8 +34,8 @@ plot_individual_irf = function( irf.plot <- plotdata %>% ggplot2::ggplot(ggplot2::aes(x=horizon, y=response, ymin=response.lower, ymax=response.upper)) + ggplot2::geom_hline(yintercept = 0, color="black") + - ggplot2::geom_ribbon(fill="lightblue", alpha=0.2) + - ggplot2::geom_line(color = 'darkred', lwd = 1.25) + + ggplot2::geom_ribbon(fill="royalblue", alpha=0.2) + + ggplot2::geom_line(color = 'darkred', lwd = 1) + ggplot2::theme_light() + ggplot2::ggtitle(title)+ ggplot2::ylab(ylab)+ diff --git a/R/regimes.R b/R/regimes.R index d0a3dc4..a5c261e 100644 --- a/R/regimes.R +++ b/R/regimes.R @@ -35,7 +35,6 @@ #' } #' #' @export - regimes = function( data, # data.frame, matrix, ts, xts, zoo: Endogenous regressors method = 'rf', # string: regime assignment technique ('rf', 'kmeans', 'EM', 'BP') diff --git a/R/sovle_b.R b/R/sovle_b.R new file mode 100644 index 0000000..0ac9276 --- /dev/null +++ b/R/sovle_b.R @@ -0,0 +1,218 @@ + +#------------------------------------------ +# Function to estimate the structural +# error impact matrix +# +# See Cesa-Bianchi's VAR-Toolbox for +# MATLAB implementation +# +# IV-short implementation based on matlab +# code by Nick von Turkovich and Paco +# Vazquez-Grande +#------------------------------------------ +solve_B = function(var, report_iv = TRUE){ + + # function variables + forecast.date = NA + + # set residuals + if(class(var) == 'VAR'){ + residuals = var$residuals[[1]] + }else{ + residuals = var$residuals + } + + if(is.na(var$structure) == TRUE){ + + # retrieve reduced-form residuals + data.residuals = residuals + + # reduced form variance-covariance matrix + cov.matrix = stats::var(stats::na.omit(dplyr::select(data.residuals, -date, -forecast.date))) + + B = cov.matrix + + return(B) + + }else if(var$structure == 'short'){ + + # retrieve reduced-form residuals + data.residuals = residuals + + # reduced form variance-covariance matrix + cov.matrix = stats::var(stats::na.omit(dplyr::select(data.residuals, -date, -forecast.date))) + + # take cholesky decomposition + B = t(chol(cov.matrix)) + + return(B) + + }else if(var$structure == 'IV'){ + + # retrieve reduced-form residuals + data.residuals = residuals + covariates = colnames(dplyr::select(data.residuals, -date, -forecast.date)) + col_to_instrument = var$instrumented + + # retrieve instrument + data.instrument = var$instrument + instrument = colnames(dplyr::select(data.instrument, -date)) + + # combine data + data = + dplyr::inner_join(data.residuals, data.instrument, by = 'date') %>% + stats::na.omit() + + # first stage least squares + model.first_stage = + stats::lm(col_to_instrument ~., + data = data %>% + dplyr::select(col_to_instrument = dplyr::all_of(col_to_instrument), dplyr::all_of(instrument)) %>% + stats::na.omit()) + p_hat = model.first_stage$fitted.values + + if(report_iv == TRUE){ + print(summary(model.first_stage)) + } + + # second stage least squares + # to estimate first column of B + # (automatically scales first entry to 1) + instrumented_shocks = covariates %>% + purrr::map_df(.f = function(X){ + + model.second_stage = stats::lm(data[,X] ~ p_hat) + second_stage_beta = model.second_stage$coefficients[2] + return(second_stage_beta) + + }) %>% + as.vector() + + # scale size of the shock + # see Gertler and Karadi (2015) for background + # see Cesa-Bianchi's VAR-Toolbox for MATLAB implementation + X.demean = + dplyr::select(data, -date, -forecast.date, -dplyr::all_of(instrument)) %>% + dplyr::mutate_all(function(X){return(X - mean(X, na.rm = T))}) %>% + as.matrix() + + sigma_b = 1/(nrow(data) - ncol(var$model$coef) + 1) * t(X.demean) %*% X.demean + + s21s11 = instrumented_shocks[2:length(covariates),] %>% as.matrix() + S11 = sigma_b[1,1] %>% as.numeric() + S21 = sigma_b[2:nrow(sigma_b),1] %>%as.vector() + S22 = sigma_b[2:nrow(sigma_b),2:ncol(sigma_b)] + + Q = (s21s11 * S11) %*% t(s21s11) - (S21 %*% t(s21s11) + as.matrix(s21s11) %*% t(S21)) + S22 + sp = sqrt( S11 - t(S21 - as.matrix(s21s11) %*% S11) %*% as.matrix(solve(Q) %*% (S21 - s21s11 * S11)) ) + + scaled_instrumented_shocks = instrumented_shocks * as.numeric(sp) + + # prepare B matrix + B = matrix(0, ncol = (length(covariates) - 1), nrow = length(covariates)) + B = cbind(scaled_instrumented_shocks, B) + + return(B) + + }else if(var$structure == 'IV-short'){ + + # align instrument and residuals + valid_dates = + dplyr::inner_join( + dplyr::select(stats::na.omit(residuals), date), + dplyr::select(stats::na.omit(var$instrument), date), + by = 'date' + ) + + # set residuals matrix + residuals = residuals %>% + dplyr::inner_join(valid_dates, by = 'date') %>% + dplyr::select(-date, -forecast.date) %>% + as.matrix() + + residuals = -1*residuals + + # set instrument + instrument = var$instrument %>% + data.frame() %>% + dplyr::inner_join(valid_dates, by = 'date') %>% + dplyr::select(-date) + + instrument.mean = instrument %>% + dplyr::mutate_all(mean, na.rm = T) + + # number of observations + n.obs = nrow(residuals) + + # solve for IV implied impact + pi = solve(n.obs^(-1) * t(residuals) %*% residuals) %*% + (n.obs^(-1) * t(residuals) %*% as.matrix(instrument - instrument.mean)) + + phi_sq = + (n.obs^(-1) * t(instrument - instrument.mean) %*% residuals) %*% + solve( n.obs^(-1) * t(residuals) %*% residuals ) %*% + ( n.obs^(-1) * t(residuals) %*% as.matrix(instrument - instrument.mean) ) + + B1 = + n.obs^(-1) * + ( t(residuals) %*% as.matrix(instrument - instrument.mean) ) %*% + (phi_sq)^(-0.5) + + B = matrix(ncol = ncol(residuals), nrow = ncol(residuals)) + B[,1:ncol(instrument)] = B1 + rownames(B) = colnames(B) = colnames(residuals) + + # solve for orthogonalized structural shock + model.first_stage = stats::lm(instrument[,1] ~ residuals) + orthogonal_instrument = instrument - model.first_stage$residuals + orthogonal_instrument = orthogonal_instrument[,] / stats::sd(orthogonal_instrument[,], na.rm = T) + + shock_matrix = matrix(nrow = nrow(residuals), ncol = ncol(residuals)) + shock_matrix[,1] = orthogonal_instrument + + # reduced form variance-covariance matrix + sigma = stats::var(residuals) + + # solve additional entries + # with a lower triangular restriction + order_sequence = c(1:ncol(residuals)) + + for(i in order_sequence){ + + Y = residuals[,i] + X = shock_matrix[,c(1:i)] + model.second_stage = stats::lm(Y~X) + + if(i != utils::tail(order_sequence,1)){ + + B[i,] = + c( model.second_stage$coef[-1], + stats::sd(model.second_stage$residuals), + rep(0, length(order_sequence) - ncol(instrument) - i) ) + + shock_matrix[,i+1] = model.second_stage$residuals / stats::sd(model.second_stage$residuals) + + }else{ + + B[i,] = model.second_stage$coef[-1] + + } + + } + + # make diagonal entries positive + shock_ordering = data.frame(residuals) %>% + dplyr::select(var$instrumented, dplyr::everything()) %>% + colnames() + B.sign = diag(sign(diag(B[shock_ordering,]))) + B = B %*% B.sign + + return(B) + + }else{ + + stop('structure must be set to either NA, "short", "IV", or "IV-short".') + + } + +} diff --git a/R/var.R b/R/var.R index 1ad7879..6829c1c 100644 --- a/R/var.R +++ b/R/var.R @@ -51,7 +51,7 @@ VAR_estimation = function( se = broom::tidy(model) %>% dplyr::select(term, std.error) se$y = target - ll = stats::logLik(model) + ll = stats::logLik(model)[1] # return results return(list(coef = c, se = se, ll = ll)) @@ -73,7 +73,7 @@ VAR_estimation = function( # extract log likelihood ll = purrr::map(models, .f = function(X){return(X$ll)}) %>% - purrr::reduce(dplyr::bind_rows) %>% + purrr::reduce(rbind) %>% sum() # package for return @@ -176,18 +176,18 @@ VAR_estimation = function( #' @param p int: lags #' @param lag.ic string: information criterion to choose the optimal number of lags ('AIC' or 'BIC') #' @param lag.max int: maximum number of lags to test in lag selection -#' @param structure string: type of structural identification strategy to use in model analysis (NULL, 'short', 'IV', or 'IV-short') +#' @param structure string: type of structural identification strategy to use in model analysis (NA, 'short', 'IV', or 'IV-short') #' @param instrument string: name of instrumental variable contained in the data matrix #' @param instrumented string: name of variable to be instrumented in IV and IV-short procedure; default is the first non-date variable in data #' #' @return #' 1. data: data.frame with endogenous variables and 'date' column. -#' 2. model: list with data.frame of model coefficients (in psuedo-companion form), data.frame of coefficient standard errors, integer of lags p, integer of horizons, string of frequency, string of type +#' 2. model: list with data.frame of model coefficients (in psuedo-companion form), data.frame of coefficient standard errors, integer of lags p, integer of horizons, string of frequency, string of deterministic term type, numeric of log-likelihood #' 3. forecasts: list of data.frames per horizon; data.frame with column for date (day the forecast was made), forecast.date (the date being forecasted), target (variable forecasted), and forecast #' 4. residuals: list of data.frames per horizon; data.frame of residuals -#' 5. structure: string denoting which structural identification strategy will be used in analysis (or NULL) +#' 5. structure: string denoting which structural identification strategy will be used in analysis (or NA) #' 6. instrument: data.frame with 'date' column and 'instrument' column (or NULL) -#' 7. instrumented: string denoting which column will be instrumted in 'IV' and 'IV-short' strategies (or NULL) +#' 7. instrumented: string denoting which column will be instrumented in 'IV' and 'IV-short' strategies (or NA) #' #' @seealso [VAR()] #' @seealso [var_irf()] @@ -280,8 +280,8 @@ VAR = function( if(!freq %in% c('day','week','month','quarter','year')){ stop("freq must be one of the following strings: 'day','week','month','quarter','year'") } - if(!structure %in% c('short', 'IV', 'IV-short') & !is.null(structure)){ - stop("strucutre must be one of 'strucutre', 'IV', 'IV-short', or NULL.") + if(!structure %in% c(NA, 'short', 'IV', 'IV-short')){ + stop("structure must be one of 'short, 'IV', 'IV-short', or NA") } if(!is.null(instrument)){ if(!instrument %in% colnames(data)){ @@ -308,7 +308,7 @@ VAR = function( } # detect variable to be instrumented - if(is.null(instrumented)){ + if(is.null(instrumented) & structure %in% c('IV', 'IV-short')){ var_to_instrument = colnames(dplyr::select(data, -date))[1] }else{ var_to_instrument = instrumented @@ -391,28 +391,6 @@ RVAR_estimate = function( type = 'const' # string: type of deterministic terms to add ('none', 'const', 'trend', 'both') ){ - # function warnings - if(!is.matrix(data) & !is.data.frame(data)){ - stop('data must be a matrix or data.frame') - } - if(!is.numeric(p) | p %% 1 != 0){ - stop('p must be an integer') - } - if(!is.numeric(horizon) | horizon %% 1 != 0 | horizon <= 0){ - stop('horizon must be a positive integer') - } - if(!freq %in% c('day','week','month','quarter','year')){ - stop("freq must be one of the following strings: 'day','week','month','quarter','year'") - } - if(!regime %in% colnames(data)){ - stop('regime must be the name of a column in data') - } - - # cast as data frame if ts, xts, or zoo object - if(stats::is.ts(data) | xts::is.xts(data) | zoo::is.zoo(data)){ - data = data.frame(date = zoo::index(date), data) - } - # function variables term = estimate = std.error = model.regime = NULL @@ -458,11 +436,16 @@ RVAR_estimate = function( c = broom::tidy(model) %>% dplyr::select(term, coef = estimate) c$y = target + # standard errors se = broom::tidy(model) %>% dplyr::select(term, std.error) se$y = target + # log likelihood + ll = stats::logLik(model)[1] + # return results - return(list(coef = c, se = se)) + return(list(coef = c, se = se, ll = ll)) + }) # extract coefficients @@ -472,14 +455,20 @@ RVAR_estimate = function( tidyr::pivot_wider(values_from = coef, names_from = term) - # extract coefficients + # extract coefficients' standard errors se = purrr::map(models, .f = function(X){return(X$se)}) %>% purrr::reduce(dplyr::bind_rows) %>% tidyr::pivot_wider(values_from = std.error, names_from = term) + # extract log likelihood + ll = + purrr::map(models, .f = function(X){return(X$ll)}) %>% + purrr::reduce(rbind) %>% + sum() + # package for return - model = list(coef = coef, se = se, p = p, freq = freq, horizon = horizon, regime = regime.val) + model = list(coef = coef, se = se, p = p, freq = freq, type = type, ll = ll, horizon = horizon, regime = regime.val) }) @@ -586,20 +575,20 @@ RVAR_estimate = function( ) ) - }) + }) ### Organize and return output ------- forecasts = as.list(c(1:horizon)) %>% purrr::map(.f = function(h){ - f = fr %>% - purrr::map(.f = function(r){ - return(r$forecast[[paste0('H_',h)]]) - }) %>% - purrr::reduce(dplyr::bind_rows) %>% - dplyr::arrange(date) + f = fr %>% + purrr::map(.f = function(r){ + return(r$forecast[[paste0('H_',h)]]) + }) %>% + purrr::reduce(dplyr::bind_rows) %>% + dplyr::arrange(date) - }) + }) residuals = as.list(c(1:horizon)) %>% @@ -629,18 +618,15 @@ RVAR_estimate = function( ) } -#' Estimate regime-dependent VAR + +#' Estimate regime-dependent VAR, SVAR, or Proxy-SVAR #' #' Estimate a regime-dependent VAR (i.e. a state-dependent VAR), with an exogenous state indicator, of the specification: -#' \deqn{Y_t = X \beta_s + \epsilon_t} -#' where t is the time index, Y is the set of outcome vectors, X the design matrix (of p lagged values of Y), and -#' s is a mutually exclusive state of the world observed at time t-1. When the regime vector is not supplied by the user, then a two-state +#' \deqn{Y_{t+1} = X_t \beta_{s_t} + \epsilon_t} +#' where *t* is the time index, *Y* is the set of outcome vectors, *X* the design matrix (of *p* lagged values of Y), and +#' *s* is a mutually exclusive state of the world observed at time *t*. When the regime vector is not supplied by the user, then a two-state #' regime series is estimated via random forest. #' -#' @details The regime-dependent VAR is a generalization of the popular threshold VAR - which trades off estimating a threshold value for an -#' endogenous variable for accepting an exogenous regime that can be based on information from inside or outside of the system, with or without parametric -#' assumptions, and with or without timing restrictions. -#' #' #' @param data data.frame, matrix, ts, xts, zoo: Endogenous regressors #' @param horizon int: forecast horizons @@ -652,17 +638,25 @@ RVAR_estimate = function( #' @param regime string: name or regime assignment vector in the design matrix (data) #' @param regime.method string: regime assignment technique ('rf', 'kmeans', 'EM', or 'BP') #' @param regime.n int: number of regimes to estimate (applies to kmeans and EM) +#' @param structure string: type of structural identification strategy to use in model analysis (NA, 'short', 'IV', or 'IV-short') +#' @param instrument string: name of instrumental variable contained in the data matrix +#' @param instrumented string: name of variable to be instrumented in IV and IV-short procedure; default is the first non-date variable in data #' -#' @return list of lists, each regime returns its own list with elements `data`, `model`, `forecasts`, `residuals` +#' @return +#' List of lists, where each regime is a list with items: +#' 1. data: data.frame with endogenous variables and 'date' column. +#' 2. model: list with data.frame of model coefficients (in psuedo-companion form), data.frame of coefficient standard errors, integer of lags p, integer of horizons, string of frequency, string of deterministic term type, numeric of log-likelihood, regime indicator +#' 3. forecasts: list of data.frames per horizon; data.frame with column for date (day the forecast was made), forecast.date (the date being forecasted), target (variable forecasted), and forecast +#' 4. residuals: list of data.frames per horizon; data.frame of residuals +#' 5. structure: string denoting which structural identification strategy will be used in analysis (or NA) +#' 6. instrument: data.frame with 'date' column and 'instrument' column (or NULL) +#' 7. instrumented: string denoting which column will be instrumented in 'IV' and 'IV-short' strategies (or NULL) #' #' @seealso [VAR()] -#' @seealso [var_irf()] -#' @seealso [var_fevd()] -#' @seealso [var_hd()] #' @seealso [RVAR()] -#' @seealso [rvar_irf()] -#' @seealso [rvar_fevd()] -#' @seealso [rvar_hd()] +#' @seealso [IRF()] +#' @seealso [FEVD()] +#' @seealso [HD()] #' #' @examples #' \donttest{ @@ -686,17 +680,44 @@ RVAR_estimate = function( #' lag.ic = 'BIC', #' lag.max = 4) #' -#' # impulse response functions -#' rvar.irf = sovereign::rvar_irf(rvar) +#' # impulse response functions +#' rvar.irf = sovereign::rvar_irf(rvar) #' -#' # forecast error variance decomposition -#' rvar.fevd = sovereign::rvar_fevd(rvar) +#' # forecast error variance decomposition +#' rvar.fevd = sovereign::rvar_fevd(rvar) #' -#' # historical shock decomposition -#' rvar.hd = sovereign::rvar_hd(rvar) +#' # historical shock decomposition +#' rvar.hd = sovereign::rvar_hd(rvar) #' #' } #' +#' @details +#' The regime-dependent VAR is a generalization of the popular threshold VAR - which trades off estimating a threshold value for an +#' endogenous variable for accepting an exogenous regime that can be based on information from inside or outside of the system, with or without parametric +#' assumptions, and with or without timing restrictions. Moreover, the RVAR may be extended to include structural shocks, including the use of instrumental +#' variables. +#' +#' **State dependence.** The RVAR augments the traditional VAR by allowing state-dependence in the coefficient matrix. The RVAR differs from the more common threshold VAR (TVAR), due +#' to the fact that states are exogenously determined. As a result, the states (i.e. regimes) do not need to be based on information inside the model, moreover, regimes can be +#' determined by any process the user determines best fits their needs. For example, regimes based on NBER dated recessions and expansions are based on a judgmental process +#' considering hundreds of series, potentially none of which are in the VAR being modeled. Alternatively, a user may use unsupervised machine learning to assign regimes - this is +#' the process the `sovereign::regimes` function facilitates. +#' +#' **Structural shocks.** See Sims (1980) for details regarding the baseline vector-autoregression (VAR) model. The VAR may be augmented to become a structural VAR (SVAR) with one of three different structural identification strategies: +#' 1) short-term impact restrictions via Cholesky decomposition, see Christiano et al (1999) for details **(structure = 'short')** +#' 2) external instrument identification, i.e. a Proxy-SVAR strategy, see Mertens and Ravn (2013) for details **(structure = 'IV')** +#' 3) or a combination of short-term and IV identification via Lunsford (2015) **(structure = 'IV-short')** +#' +#' Note that including structure does not change the estimation of model coefficients or forecasts, but does change impulse response functions, forecast error variance decomposition, +#' and historical decompositions. Historical decompositions will not be available for models using the 'IV' structure. Additionally note that only one instrument may be used in this +#' estimation routine. +#' +#' @references +#' 1. Christiano, Lawrence, Martin Eichenbaum, and Charles Evans "[Monetary policy shocks: What have we learned and to what end?](https://www.sciencedirect.com/science/article/pii/S1574004899010058)" Handbook of Macroeconomics, Vol 1, Part A, 1999. +#' 2. Lunsford, Kurt "[Identifying Structural VARs with a Proxy Variable and a Test for a Weak Proxy](https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2699452#)" 2015. +#' 3. Mertens, Karel and Morten Ravn "[The Dynamic Effects of Personal and Corporate Income Tax Changes in the United States](https://www.aeaweb.org/articles?id=10.1257/aer.103.4.1212)" 2013. +#' 4. Sims, Christopher "[Macroeconomics and Reality](https://www.jstor.org/stable/1912017)" 1980. +#' #' @export # regime-dependent VAR function @@ -710,7 +731,10 @@ RVAR = function( lag.max = NULL, # int: maximum number of lags to test in lag selection regime = NULL, # string: name or regime assignment vector in the design matrix (data) regime.method = 'rf', # string: regime assignment technique ('rf', 'kmeans', 'EM', 'BP') - regime.n = 2 # int: number of regimes to estimate (applies to kmeans and EM) + regime.n = 2, # int: number of regimes to estimate (applies to kmeans and EM) + structure = 'short', # string: type of structural identification strategy to use in model analysis (NULL, 'short', or 'IV') + instrument = NULL, # string: name of instrumental variable contained in the data matrix + instrumented = NULL # string: name of variable to be instrumented in IV and IV-short procedure; default is the first non-date variable in data ){ # function warnings @@ -727,6 +751,52 @@ RVAR = function( stop('lag.max must be an integer if IC-based lag selection is used') } } + if(!structure %in% c(NA, 'short', 'IV', 'IV-short')){ + stop("structure must be one of 'short, 'IV', 'IV-short', or NA") + } + if(!is.null(instrument)){ + if(!instrument %in% colnames(data)){ + stop("instrument must be the name of a variable found in data.") + } + } + if(!is.null(instrumented)){ + if(!instrumented %in% colnames(data)){ + stop("instrumented must be the name of a variable found in data.") + } + } + if(!is.matrix(data) & !is.data.frame(data)){ + stop('data must be a matrix or data.frame') + } + if(!is.numeric(p) | p %% 1 != 0){ + stop('p must be an integer') + } + if(!is.numeric(horizon) | horizon %% 1 != 0 | horizon <= 0){ + stop('horizon must be a positive integer') + } + if(!freq %in% c('day','week','month','quarter','year')){ + stop("freq must be one of the following strings: 'day','week','month','quarter','year'") + } + + # cast as data frame if ts, xts, or zoo object + if(stats::is.ts(data) | xts::is.xts(data) | zoo::is.zoo(data)){ + data = data.frame(date = zoo::index(date), data) + } + + # set aside instruments + if(!is.null(instrument)){ + data.instrument = dplyr::select(data, date, dplyr::all_of(instrument)) + data = dplyr::select(data, -dplyr::all_of(instrument)) + }else{ + data.instrument = NULL + } + + # detect variable to be instrumented + # (defaults to first non-date column if missing) + if(is.null(instrumented) & structure %in% c('IV', 'IV-short')){ + var_to_instrument = colnames(dplyr::select(data, -date, -regime))[1] + }else{ + var_to_instrument = instrumented + } # create regimes if(is.null(regime)){ @@ -780,9 +850,6 @@ RVAR = function( min.ic = which.min(ic.scores) model = models[[min.ic]] - class(model) = 'RVAR' - return(model) - }else{ model = @@ -794,13 +861,15 @@ RVAR = function( freq = freq, type = type ) - - class(model) = 'RVAR' - return(model) - } + # add structure + model$structure = structure + model$instrument = data.instrument + model$instrumented = var_to_instrument -} + class(model) = 'RVAR' + return(model) +} diff --git a/R/var_fevd.R b/R/var_fevd.R index 009e708..2fbc74e 100644 --- a/R/var_fevd.R +++ b/R/var_fevd.R @@ -2,32 +2,35 @@ # Estimate forecast error variance # source code adapted from the MTS package #--------------------------------------------- -fevd = function(Phi, Sig, lag = 4) +fevd = function(Phi, Sig, lag = 4, structure = NA) { if (length(Phi) > 0) { if (!is.matrix(Phi)) Phi = as.matrix(Phi) } - if (!is.matrix(Sig)) - Sig = as.matrix(Sig) if (lag < 1) lag = 1 + p = 0 + if (length(Phi) > 0) { k = nrow(Phi) m = ncol(Phi) p = floor(m/k) } + q = 0 Si = diag(rep(1, k)) m = (lag + 1) * k m1 = (q + 1) * k + if (m > m1) { Si = cbind(Si, matrix(0, k, (m - m1))) } + if (p > 0) { for (i in 1:lag) { if (i <= p) { @@ -52,39 +55,47 @@ fevd = function(Phi, Sig, lag = 4) Si[, (kdx + 1):(kdx + k)] = tmp } } + orSi = NULL - m1 = chol(Sig) + m1 = Sig P = t(m1) orSi = P + for (i in 1:lag) { idx = i * k w1 = Si[, (idx + 1):(idx + k)] w2 = w1 %*% P orSi = cbind(orSi, w2) } + orSi2 = orSi^2 Ome = orSi2[, 1:k] wk = Ome + for (i in 1:lag) { idx = i * k wk = wk + orSi2[, (idx + 1):(idx + k)] Ome = cbind(Ome, wk) } + FeV = NULL OmeRa = Ome[, 1:k] FeV = cbind(FeV, apply(OmeRa, 1, sum)) OmeRa = OmeRa/FeV[, 1] + for (i in 1:lag) { idx = i * k wk = Ome[, (idx + 1):(idx + k)] FeV = cbind(FeV, apply(wk, 1, sum)) OmeRa = cbind(OmeRa, wk/FeV[, (i + 1)]) } + for (i in 1:(lag + 1)) { idx = (i - 1) * k Ratio = OmeRa[, (idx + 1):(idx + k)] } - FEVdec <- list(irf = Si, orthirf = orSi, Omega = Ome, OmegaR = OmeRa) + + FEVdec = list(irf = Si, orthirf = orSi, Omega = Ome, OmegaR = OmeRa) return(FEVdec) } @@ -94,6 +105,9 @@ fevd = function(Phi, Sig, lag = 4) #-------------------------------------------------------- #' Estimate forecast error variance decomposition #' +#' Estimate forecast error variance decomposition for VARs +#' with either short or 'IV-short' structural errors. +#' #' @param var VAR output #' @param horizon int: number of periods #' @param scale boolean: scale variable contribution as percent of total error @@ -133,7 +147,7 @@ fevd = function(Phi, Sig, lag = 4) #' #' # forecast error variance decomposition #' var.fevd = sovereign::var_fevd(var) -#' +#' #' # historical shock decomposition #' var.hd = sovereign::var_hd(var) #' @@ -151,6 +165,12 @@ var_fevd = function( if(!is.numeric(horizon) | horizon %% 1 != 0 | horizon <= 0){ stop('horizon must be a positive integer') } + if(class(var) != 'VAR'){ + stop('var must be a VAR object') + } + if(!var$structure %in% c('short', 'IV-short')){ + stop('FEVD method is only available for VARs with short or IV-short structural errors') + } # function variables forecast.date = y = shock = error = NULL @@ -162,12 +182,14 @@ var_fevd = function( regressors = colnames(dplyr::select(data, -date)) # error covariance matrix - cov.matrix = stats::var(stats::na.omit(dplyr::select(residuals, -date, -forecast.date))) + B = solve_B(var) + B = as.matrix(B) # forecast error variance decomposition - errors = fevd(Phi = dplyr::select(coef, -y, -dplyr::contains('cosnt'), -dplyr::contains('trend')), - Sig = cov.matrix, - lag = horizon) + errors = fevd(Phi = dplyr::select(coef, -y, -dplyr::contains('const'), -dplyr::contains('trend')), + Sig = B, + lag = horizon, + structure = var$structure) # reorganize results response = data.frame(t(errors$OmegaR)) @@ -178,7 +200,15 @@ var_fevd = function( # cast to long-form response = response %>% - tidyr::pivot_longer(cols = regressors, names_to = 'response', values_to = 'error') + tidyr::pivot_longer(!c('shock','horizon'), names_to = 'response', values_to = 'error') + + + if(var$structure %in% c('IV-short')){ + labels = c(var$instrumented, regressors[!regressors %in% var$instrumented]) + labels = rep(labels, length(labels) * (horizon+1)) + response$response = labels + response$shock = labels[order(labels)] + } # scale responses if(scale == TRUE){ @@ -193,6 +223,9 @@ var_fevd = function( #' Estimate regime-dependent forecast error variance decomposition #' +#' Estimate forecast error variance decomposition for RVARs +#' with either short or 'IV-short' structural errors. +#' #' @param rvar RVAR output #' @param horizon int: number of periods #' @param scale boolean: scale variable contribution as percent of total error @@ -235,7 +268,7 @@ var_fevd = function( #' #' # forecast error variance decomposition #' rvar.fevd = sovereign::rvar_fevd(rvar) -#' +#' #' # historical shock decomposition #' rvar.hd = sovereign::rvar_hd(rvar) #' @@ -254,6 +287,12 @@ rvar_fevd = function( if(!is.numeric(horizon) | horizon %% 1 != 0 | horizon <= 0){ stop('horizon must be a positive integer') } + if(class(rvar) != 'RVAR'){ + stop('rvar must be a RVAR object') + } + if(!rvar$structure %in% c('short', 'IV-short')){ + stop('FEVD method is only available for VARs with short or IV-short structural errors') + } # function variables forecast.date = y = shock = error = model.regime = NULL @@ -274,16 +313,38 @@ rvar_fevd = function( # set regime specific data coef = rvar$model[[paste0('regime_',regime.val)]]$coef + residuals = rvar$residuals[[1]] %>% dplyr::filter(model.regime == regime.val) - # error covariance matrix - cov.matrix = stats::var(stats::na.omit(dplyr::select(residuals, -date, -model.regime, -forecast.date))) + # instrument + if(!is.null(rvar$instrument)){ + instrument = dplyr::inner_join(rvar$instrument, dplyr::select(residuals, date), by = 'date') + }else{ + instrument = NULL + } + + # regime-specific model + model = rvar$model[[paste0('regime_',regime.val)]] + + # structural error variance-covariance matrix + B = + solve_B( + var = list( + model = model, + residuals = residuals %>% + dplyr::select(-model.regime), + structure = rvar$structure, + instrument = instrument, + instrumented = rvar$instrumented + ) + ) # forecast error variance decomposition - errors = fevd(Phi = dplyr::select(coef, -y, -dplyr::contains('cosnt'), -dplyr::contains('trend')), - Sig = cov.matrix, - lag = horizon) + errors = fevd(Phi = dplyr::select(coef, -y, -dplyr::contains('const'), -dplyr::contains('trend')), + Sig = B, + lag = horizon, + structure = rvar$structure) # reorganize results response = data.frame(t(errors$OmegaR)) @@ -294,7 +355,15 @@ rvar_fevd = function( # cast to long-form response = response %>% - tidyr::pivot_longer(cols = regressors, names_to = 'response', values_to = 'error') + tidyr::pivot_longer(!c('shock','horizon'), names_to = 'response', values_to = 'error') + + # correct labels for IV-short + if(rvar$structure %in% c('IV-short')){ + labels = c(rvar$instrumented, regressors[!regressors %in% rvar$instrumented]) + labels = rep(labels, length(labels) * (horizon+1)) + response$response = labels + response$shock = labels[order(labels)] + } # scale responses if(scale == TRUE){ diff --git a/R/var_hd.R b/R/var_hd.R index 8be2782..979e2fb 100644 --- a/R/var_hd.R +++ b/R/var_hd.R @@ -4,10 +4,9 @@ #' Estimate historical decomposition #' -#' Estimate historical decomposition with contemporaneous impact -#' restrictions via Cholesky decomposition - taking the variable ordering -#' present in the VAR object. -#' +#' Estimate historical decomposition for VARs with +#' either short or 'IV-short' structural errors. +#' #' @param var VAR output #' #' @return long-from data.frame @@ -57,6 +56,15 @@ var_hd = function( var # VAR object ){ + + # function warnings + if(class(var) != 'VAR'){ + stop('var must be a VAR object') + } + if(!var$structure %in% c('short', 'IV-short')){ + stop('HD method is only available for VARs with short or IV-short structural errors') + } + # set function variables model.regime = target = forecast.date = NULL @@ -73,8 +81,10 @@ var_hd = function( # reduced form error variance-covariance matrix sigma = stats::var(stats::na.omit(residuals)) - # structural error variance-covariance matrix - B = chol(sigma) + # estimate error-covariance matrix or structural impact matrix + B = solve_B(var) + B = as.matrix(B) + # coefficient matrix # (already written in psuedo-companion form) @@ -169,10 +179,8 @@ var_hd = function( #' Estimate regime-dependent historical decomposition #' -#' Estimate historical decomposition with contemporaneous impact -#' restrictions via Cholesky decomposition - taking the variable ordering -#' present in the RVAR object. Estimate one response function per unique -#' state defined by the regime-dependent VAR. +#' Estimate historical decomposition for RVARs with +#' either short or 'IV-short' structural errors. #' #' @param rvar RVAR output #' @@ -226,6 +234,13 @@ rvar_hd = function( rvar # threshold VAR output ){ + # function warnings + if(class(rvar) != 'RVAR'){ + stop('rvar must be a RVAR object') + } + if(!rvar$structure %in% c('short', 'IV-short')){ + stop('HD method is only available for VARs with short or IV-short structural errors') + } # function variables model.regime = target = forecast.date = NULL @@ -243,18 +258,37 @@ rvar_hd = function( purrr::map(.f = function(regime.val){ # recover model dynamics ----------------------------- - # (for now, only use cholesky decomposition) # regime-dependent reduced form residuals residuals = rvar$residuals$H_1 %>% - # dplyr::filter(model.regime == regime.val) %>% dplyr::select(-date, -forecast.date, -model.regime) - # reduced form error variance-covariance matrix - sigma = stats::var(stats::na.omit(residuals)) + is = rvar$data %>% + dplyr::inner_join(dplyr::select(rvar$residuals$H_1, date), by = 'date') %>% + dplyr::select(-regime) + + # instrument + if(!is.null(rvar$instrument)){ + instrument = dplyr::inner_join(rvar$instrument, dplyr::select(is, date), by = 'date') + }else{ + instrument = NULL + } + + # regime-specific model + model = rvar$model[[paste0('regime_',regime.val)]] # structural error variance-covariance matrix - B = chol(sigma) + B = + solve_B( + var = list( + model = model, + residuals = dplyr::inner_join(rvar$residuals$H_1, dplyr::select(is, date), by = 'date') %>% + dplyr::select(-model.regime), + structure = rvar$structure, + instrument = instrument, + instrumented = rvar$instrumented + ) + ) # coefficient matrix # (already written in psuedo-companion form) diff --git a/R/var_irf.R b/R/var_irf.R index b88a446..09960d4 100644 --- a/R/var_irf.R +++ b/R/var_irf.R @@ -1,213 +1,8 @@ - -#------------------------------------------ -# Function to estimate the structural -# error impact matrix -# -# See Cesa-Bianchi's VAR-Toolbox for -# MATLAB implementation -#------------------------------------------ -solve_B = function(var, report_iv = TRUE){ - - if(is.null(var$structure) == TRUE){ - - # retrieve reduced-form residuals - data.residuals = var$residuals[[1]] - - # reduced form variance-covariance matrix - cov.matrix = stats::var(stats::na.omit(dplyr::select(data.residuals, -date, -forecast.date))) - - B = cov.matrix - - return(B) - - }else if(var$structure == 'short'){ - - # retrieve reduced-form residuals - data.residuals = var$residuals[[1]] - - # reduced form variance-covariance matrix - cov.matrix = stats::var(stats::na.omit(dplyr::select(data.residuals, -date, -forecast.date))) - - # take cholesky decomposition - B = t(chol(cov.matrix)) - - return(B) - - }else if(var$structure == 'IV'){ - - # retrieve reduced-form residuals - data.residuals = var$residuals[[1]] - covariates = colnames(dplyr::select(data.residuals, -date, -forecast.date)) - col_to_instrument = var$instrumented - - # retrieve instrument - data.instrument = var$instrument - instrument = colnames(dplyr::select(data.instrument, -date)) - - # combine data - data = - dplyr::inner_join(data.residuals, data.instrument, by = 'date') %>% - na.omit() - - # first stage least squares - model.first_stage = - lm(col_to_instrument ~., - data = data %>% - select(col_to_instrument = dplyr::all_of(col_to_instrument), all_of(instrument)) %>% - na.omit()) - p_hat = model.first_stage$fitted.values - - if(report_iv == TRUE){ - print(summary(model.first_stage)) - } - - # second stage least squares - # to estimate first column of B - # (automatically scales first entry to 1) - instrumented_shocks = covariates %>% - purrr::map_df(.f = function(X){ - - model.second_stage = lm(data[,X] ~ p_hat) - second_stage_beta = model.second_stage$coefficients[2] - return(second_stage_beta) - - }) %>% - as.vector() - - # scale size of the shock - # see Gertler and Karadi (2015) for background - # see Cesa-Bianchi's VAR-Toolbox for MATLAB implementation - X.demean = - select(data, -date, -forecast.date, -dplyr::all_of(instrument)) %>% - mutate_all(function(X){return(X - mean(X, na.rm = T))}) %>% - as.matrix() - - sigma_b = 1/(nrow(data) - ncol(var$model$coef) + 1) * t(X.demean) %*% X.demean - - s21s11 = instrumented_shocks[2:length(covariates),] %>% as.matrix() - S11 = sigma_b[1,1] %>% as.numeric() - S21 = sigma_b[2:nrow(sigma_b),1] %>%as.vector() - S22 = sigma_b[2:nrow(sigma_b),2:ncol(sigma_b)] - - Q = (s21s11 * S11) %*% t(s21s11) - (S21 %*% t(s21s11) + as.matrix(s21s11) %*% t(S21)) + S22 - sp = sqrt( S11 - t(S21 - as.matrix(s21s11) %*% S11) %*% as.matrix(solve(Q) %*% (S21 - s21s11 * S11)) ) - - scaled_instrumented_shocks = instrumented_shocks * as.numeric(sp) - - # prepare B matrix - B = matrix(0, ncol = (length(covariates) - 1), nrow = length(covariates)) - B = cbind(scaled_instrumented_shocks, B) - - return(B) - - }else if(var$structure == 'IV-short'){ - - # align instrument and residuals - valid_dates = - dplyr::inner_join( - dplyr::select(na.omit(var$residuals[[1]]), date), - dplyr::select(na.omit(var$instrument), date), - by = 'date' - ) - - # set residuals matrix - residuals = var$residuals[[1]] %>% - dplyr::inner_join(valid_dates, by = 'date') %>% - dplyr::select(-date, -forecast.date) %>% - as.matrix() - - residuals = -1*residuals - - # set instrument - instrument = var$instrument %>% - data.frame() %>% - dplyr::inner_join(valid_dates, by = 'date') %>% - dplyr::select(-date) - - instrument.mean = instrument %>% - dplyr::mutate_all(mean, na.rm = T) - - # number of observations - n.obs = nrow(var$data) - - # solve for IV implied impact - pi = solve(n.obs^(-1) * t(residuals) %*% residuals) %*% - (n.obs^(-1) * t(residuals) %*% as.matrix(instrument - instrument.mean)) - - phi_sq = - (n.obs^(-1) * t(instrument - instrument.mean) %*% residuals) %*% - solve( n.obs^(-1) * t(residuals) %*% residuals ) %*% - ( n.obs^(-1) * t(residuals) %*% as.matrix(instrument - instrument.mean) ) - - B1 = - n.obs^(-1) * - ( t(residuals) %*% as.matrix(instrument - instrument.mean) ) %*% - (phi_sq)^(-0.5) - - B = matrix(ncol = ncol(residuals), nrow = ncol(residuals)) - B[,1:ncol(instrument)] = B1 - rownames(B) = colnames(B) = colnames(residuals) - - # solve for orthogonalized structural shock - model.first_stage = lm(instrument[,1] ~ residuals) - orthogonal_instrument = instrument - model.first_stage$residuals - orthogonal_instrument = orthogonal_instrument[,] / sd(orthogonal_instrument[,], na.rm = T) - - shock_matrix = matrix(nrow = nrow(residuals), ncol = ncol(residuals)) - shock_matrix[,1] = orthogonal_instrument - - # reduced form variance-covariance matrix - sigma = stats::var(residuals) - - # solve additional entries - # with a lower triangular restriction - order_sequence = c(1:ncol(residuals)) - - for(i in order_sequence){ - - Y = residuals[,i] - X = shock_matrix[,c(1:i)] - model.second_stage = lm(Y~X) - - if(i != tail(order_sequence,1)){ - - B[i,] = - c( model.second_stage$coef[-1], - sd(model.second_stage$residuals), - rep(0, length(order_sequence) - ncol(instrument) - i) ) - - shock_matrix[,i+1] = model.second_stage$residuals / sd(model.second_stage$residuals) - - }else{ - - B[i,] = model.second_stage$coef[-1] - - } - - } - - # make diagonal entries positive - shock_ordering = data.frame(residuals) %>% - select(var$instrumented, dplyr::everything()) %>% - colnames() - B.sign = diag(sign(diag(B[shock_ordering,]))) - B = B %*% B.sign - - return(B) - - }else{ - - stop('structure must be set to either NULL, "short", or "IV".') - - } - -} - #------------------------------------------ # Function to estimate impulse responses # source code adapted from the MTS package #------------------------------------------ -IRF = function (Phi, B, lag, structure = NULL){ +IRF_solve = function(Phi, B, lag, structure = NA){ # cast data as matrices if (!is.matrix(Phi)) @@ -263,7 +58,7 @@ IRF = function (Phi, B, lag, structure = NULL){ irf = Si - if (!is.null(structure)) { + if (!is.na(structure)) { P = B wk1 = cbind(wk1, c(P)) @@ -292,15 +87,15 @@ IRF = function (Phi, B, lag, structure = NULL){ #' Estimate impulse response functions #' -#' Estimate impulse responses with contemporaneous impact restrictions via Cholesky decomposition - -#' taking the variable ordering present in the VAR object. -#' -#' @param var VAR output -#' @param horizon int: number of periods -#' @param bootstraps.num int: number of bootstraps -#' @param CI numeric vector: c(lower ci bound, upper ci bound) +#' @param var VAR output +#' @param horizon int: number of periods +#' @param CI numeric vector: c(lower ci bound, upper ci bound) +#' @param bootstrap.type string: bootstrapping technique to use ('auto', 'standard', or 'wild'); if auto then wild is used for IV or IV-short, else standard is used +#' @param bootstrap.num int: number of bootstraps +#' @param bootstrap.parallel boolean: create IRF draws in parallel +#' @param bootstrap.cores int: number of cores to use in parallel processing; -1 detects and uses half the available cores #' -#' @return list object with elements `irfs`, `ci.lower`, and `ci.upper`; all elements are long-form data.frames +#' @return data.frame with columns `target`, `shock`, `horizon`, `response.lower`, `response`, `response.upper` #' #' @seealso [VAR()] #' @seealso [var_irf()] @@ -330,20 +125,19 @@ IRF = function (Phi, B, lag, structure = NULL){ #' lag.ic = 'BIC', #' lag.max = 4) #' -#' # impulse response functions -#' var.irf = sovereign::var_irf(var) +#' # impulse response functions +#' var.irf = sovereign::var_irf(var) #' -#' # forecast error variance decomposition -#' var.fevd = sovereign::var_fevd(var) +#' # forecast error variance decomposition +#' var.fevd = sovereign::var_fevd(var) #' -#' # historical shock decomposition -#' var.hd = sovereign::var_hd(var) +#' # historical shock decomposition +#' var.hd = sovereign::var_hd(var) #' #' } #' #' @export - var_irf = function( var, # VAR output horizon = 10, # int: number of periods @@ -354,9 +148,6 @@ var_irf = function( bootstrap.cores = -1 # int: number of cores to use in parallel processing; -1 detects and uses half the available cores ){ - # NOTES: scaled wild is used for IV, consistent with the proxy SVAR literature, - # while standard resample is used for others - # function warnings if(!is.numeric(bootstrap.num) | bootstrap.num %% 1 != 0){ stop('bootstrap.num must be an integer') @@ -369,7 +160,7 @@ var_irf = function( } # function variables - y = forecast.date = shock = target = response = response.lower = response.upper = response.med = response.adjust = NULL + y = forecast.date = response.mean = shock = target = response = response.lower = response.upper = response.med = response.adjust = NULL # set data coef = var$model$coef @@ -379,15 +170,11 @@ var_irf = function( # set model variables p = var$model$p freq = var$model$freq - structure = var$structure + type = var$model$type - if('const' %in% colnames(coef) & 'trend' %in% colnames(coef)){ - type = 'both' - }else if('const' %in% colnames(coef)){ - type = 'const' - }else if('trend' %in% colnames(coef)){ - type = 'trend' - } + structure = var$structure + instrument = var$instrument + instrumented = var$instrumented regressors = colnames(dplyr::select(data, -date)) if(!is.null(var$regime)){regressors = regressors[regressors != var$regime$regime]} @@ -404,27 +191,45 @@ var_irf = function( B = as.matrix(B) # estimate IRFs - phi = dplyr::select(coef, -y, -dplyr::contains('cosnt'), -dplyr::contains('trend')) + phi = dplyr::select(coef, -y, -dplyr::contains('const'), -dplyr::contains('trend')) - irf = IRF(Phi = phi, - B = B, - lag = horizon, - structure = structure) + irf = + IRF_solve( + Phi = phi, + B = B, + lag = horizon, + structure = structure + ) # reorganize results irf = data.frame(t(irf)) - if(!is.null(structure)){ - if(structure == 'IV-short'){ - shocks = c(var$instrumented, regressors[!regressors %in% var$instrumented]) - irf$shock = rep(shocks, horizon + 1) - } - }else{ + # assign shocks + if(is.na(structure)){ + irf$shock = rep(regressors, horizon + 1) + irf$horizon = sort(rep(c(0:horizon), length(regressors))) + + }else if(structure == 'IV-short'){ + + shocks = c(instrumented, regressors[!regressors %in% instrumented]) + irf$shock = rep(shocks, horizon + 1) + irf$horizon = sort(rep(c(0:horizon), length(regressors))) + + }else if(structure == 'IV'){ + + irf = irf[(c(1:(horizon))*length(regressors)-length(regressors)+1),] + irf$shock = instrumented + irf$horizon = c(0:(horizon-1)) + + }else if(structure == 'short'){ + + irf$shock = rep(regressors, horizon + 1) + irf$horizon = sort(rep(c(0:horizon), length(regressors))) + } - irf$horizon = sort(rep(c(0:horizon), length(regressors))) - irf = dplyr::arrange(irf, shock, horizon) + # observation identifiers rownames(irf) = NULL colnames(irf) = c(regressors, 'shock', 'horizon') @@ -448,8 +253,7 @@ var_irf = function( # bootstrap residuals if(bootstrap.type == 'wild' | - bootstrap.type == 'auto' & structure == 'IV' | - bootstrap.type == 'auto' & structure == 'IV-short'){ + bootstrap.type == 'auto' & structure %in% c('IV', 'IV-short')){ # 'wild' bootstrap technique for simple distribution # using observed scaled residuals with random sign flip. @@ -458,11 +262,18 @@ var_irf = function( r = sample(c(-1,1), size = nrow(U), replace = T) U = sweep(U, MARGIN = 1, r, `*`) - }else if(bootstrap.type == 'standard' | bootstrap.type == 'auto' & structure != 'IV'){ + # bootstrap instrument + if(!is.null(instrument)){ + instrument.bag = instrument + instrument.bag[,-1] = instrument.bag[,-1] * r + } + + }else if(bootstrap.type == 'standard' | + bootstrap.type == 'auto' & !structure %in% c('IV', 'IV-short')){ - # standard bootstrap technique a al Lutkepohl (2005) + # standard resample bootstrap technique a al Lutkepohl (2005) # draw residuals with replacement - U = na.omit(residuals) + U = stats::na.omit(residuals) U = U[sample(c(1:nrow(U)), size = nrow(residuals), replace = TRUE),] @@ -479,7 +290,7 @@ var_irf = function( for(i in (p+1):nrow(var$data)){ X = Y[(i-p):(i-1),] - X = embed(as.matrix(X), dimension = p) + X = stats::embed(as.matrix(X), dimension = p) X = data.frame(X) if(type %in% c('const', 'both')){X$const = 1} @@ -501,18 +312,12 @@ var_irf = function( p = p, horizon = 1, freq = freq, - type = type, - structure = structure) - - # bootstrap instrument - if(!is.null(structure)){ - if(structure == 'IV' | structure == 'IV-short'){ - var.bag$instrument = var$instrument - var.bag$instrument[,-1] = sweep(var.bag$instrument[,-1], MARGIN = 1, r, `*`) + type = type + ) - var.bag$instrumented = var$instrumented - } - } + var.bag$structure = structure + var.bag$instrumented = instrumented + if(!is.null(instrument)){var.bag$instrument = instrument.bag} # estimate error-covariance matrix or structural impact matrix B.bag = solve_B(var.bag, report_iv = FALSE) @@ -521,60 +326,64 @@ var_irf = function( coef.bag = dplyr::select(var.bag$model$coef, -y, -dplyr::contains('const'), -dplyr::contains('trend')) # estimate IRFs - irf.bag = IRF(Phi = coef.bag, - B = B.bag, - lag = horizon, - structure = structure) + irf.bag = + IRF_solve( + Phi = coef.bag, + B = B.bag, + lag = horizon, + structure = structure + ) # reorganize results irf.bag = data.frame(t(irf.bag)) - if(!is.null(structure)){ - if(structure == 'IV-short'){ - shocks = c(var$instrumented, regressors[!regressors %in% var$instrumented]) - irf.bag$shock = rep(shocks, horizon + 1) - } - }else{ - irf.bag$shock = rep(regressors, horizon + 1) + # assign shocks + if(is.na(structure)){ + + irf.bag$shock = rep(regressors, horizon + 1) + irf.bag$horizon = sort(rep(c(0:horizon), length(regressors))) + + }else if(structure == 'IV-short'){ + + shocks = c(var$instrumented, regressors[!regressors %in% var$instrumented]) + irf.bag$shock = rep(shocks, horizon + 1) + irf.bag$horizon = sort(rep(c(0:horizon), length(regressors))) + + }else if(structure == 'IV'){ + + irf.bag = irf.bag[(c(1:(horizon))*length(regressors)-length(regressors)+1),] + irf.bag$shock = var$instrumented + irf.bag$horizon = c(0:(horizon-1)) + + }else if(structure == 'short'){ + + irf.bag$shock = rep(regressors, horizon + 1) + irf.bag$horizon = sort(rep(c(0:horizon), length(regressors))) + } - irf.bag$horizon = sort(rep(c(0:horizon), length(regressors))) - irf.bag = dplyr::arrange(irf.bag, shock, horizon) + # observation identifiers rownames(irf.bag) = NULL colnames(irf.bag) = c(regressors, 'shock', 'horizon') return(irf.bag) - }) - - # 3. calculate confidence intervals + }, + .options = furrr::furrr_options(seed = TRUE)) + # 2. calculate confidence intervals # collapse to data.frame ci = bagged.irf %>% purrr::reduce(dplyr::bind_rows) - # remove unused shocks in the case of IV - if(!is.null(var$structure)){ - if(var$structure == 'IV'){ - ci = ci %>% - dplyr::filter(shock == regressors[1]) - } - } - - # correct shock names in the case of IV-short - if(!is.null(var$structure)){ - if(var$structure == 'IV-short'){ - - } - } - # estimate bands ci = ci %>% tidyr::pivot_longer(cols = regressors, names_to = 'target', values_to = 'response') %>% - group_by(shock, horizon, target) %>% - summarize(response.lower = quantile(response, probs = p.lower, na.rm = T), - response.upper = quantile(response, probs = p.upper, na.rm = T), - response.mean = mean(response, na.rm = T) ) + dplyr::group_by(shock, horizon, target) %>% + dplyr::summarize( + response.lower = stats::quantile(response, probs = p.lower, na.rm = T), + response.upper = stats::quantile(response, probs = p.upper, na.rm = T), + response.mean = mean(response, na.rm = T) ) # combine point estimates and CI irf = irf %>% @@ -582,7 +391,7 @@ var_irf = function( dplyr::arrange(shock, horizon) irf = - full_join( + dplyr::full_join( irf, ci, by = c('shock', 'target', 'horizon') ) @@ -591,29 +400,29 @@ var_irf = function( # (bias can be introduced in the bootstrapping if residuals are not actually mean zero) irf = irf %>% dplyr::mutate( - response.adjust = response - response.mean, - response.lower = response.lower + response.adjust, - response.upper = response.upper + response.adjust + response.adjust = response - response.mean, + response.lower = response.lower + response.adjust, + response.upper = response.upper + response.adjust ) %>% dplyr::select(target, shock, horizon, response.lower, response, response.upper) %>% dplyr::arrange(target, shock, horizon) ### return output -------------- return(irf) + } #' Estimate regime-dependent impulse response functions #' -#' Estimate impulse responses with contemporaneous impact restrictions via Cholesky decomposition - -#' taking the variable ordering present in the VAR object. Estimate one response function per unique -#' state defined by the regime-dependent VAR. +#' @param rvar RVAR output +#' @param horizon int: number of periods +#' @param CI numeric vector: c(lower ci bound, upper ci bound) +#' @param bootstrap.type string: bootstrapping technique to use ('auto', 'standard', or 'wild'); if auto then wild is used for IV or IV-short, else standard is used +#' @param bootstrap.num int: number of bootstraps +#' @param bootstrap.parallel boolean: create IRF draws in parallel +#' @param bootstrap.cores int: number of cores to use in parallel processing; -1 detects and uses half the available cores #' -#' @param rvar RVAR output -#' @param horizon int: number of periods -#' @param bootstraps.num int: number of bootstraps -#' @param CI numeric vector: c(lower ci bound, upper ci bound) -#' -#' @return list of lists, each regime returns its own list with elements `irfs`, `ci.lower`, and `ci.upper`; all elements are long-form data.frames +#' @return list of regimes, each with data.frame of columns `target`, `shock`, `horizon`, `response.lower`, `response`, `response.upper` #' #' @seealso [VAR()] #' @seealso [var_irf()] @@ -644,29 +453,32 @@ var_irf = function( #' lag.ic = 'BIC', #' lag.max = 4) #' -#' # impulse response functions -#' rvar.irf = sovereign::rvar_irf(rvar) +#' # impulse response functions +#' rvar.irf = sovereign::rvar_irf(rvar) #' -#' # forecast error variance decomposition -#' rvar.fevd = sovereign::rvar_fevd(rvar) +#' # forecast error variance decomposition +#' rvar.fevd = sovereign::rvar_fevd(rvar) #' -#' # historical shock decomposition -#' rvar.hd = sovereign::rvar_hd(rvar) +#' # historical shock decomposition +#' rvar.hd = sovereign::rvar_hd(rvar) #' #' } #' #' @export rvar_irf = function( - rvar, # threshold VAR output - horizon = 10, # int: number of periods - bootstraps.num = 100, # int: number of bootstraps - CI = c(0.1, 0.9) # numeric vector: c(lower ci bound, upper ci bound) + rvar, # threshold VAR output + horizon = 10, # int: number of periods + CI = c(0.1, 0.9), # numeric vector: c(lower ci bound, upper ci bound) + bootstrap.type = 'auto', # string: bootstrapping technique to use ('auto', 'standard', or 'wild'); if auto then wild is used for IV or IV-short, else standard is used + bootstrap.num = 100, # int: number of bootstraps + bootstrap.parallel = FALSE, # boolean: create IRF draws in parallel + bootstrap.cores = -1 # int: number of cores to use in parallel processing; -1 detects and uses half the available cores ){ # function warnings - if(!is.numeric(bootstraps.num) | bootstraps.num %% 1 != 0){ - stop('bootstraps.num must be an integer') + if(!is.numeric(bootstrap.num) | bootstrap.num %% 1 != 0){ + stop('bootstrap.num must be an integer') } if(!is.numeric(CI) | length(CI) != 2 | min(CI) < 0 | max(CI) > 1 | is.na(sum(CI))){ stop('CI must be a two element numeric vector bound [0,1]') @@ -676,25 +488,18 @@ rvar_irf = function( } # function variables - y = forecast.date = shock = target = response = response.lower = response.upper = model.regime = response.med = response.adjust = NULL + y = forecast.date = response.mean = shock = target = response = response.lower = response.upper = model.regime = response.med = response.adjust = NULL # set data data = rvar$data p = rvar$model[[1]]$p freq = rvar$model[[1]]$freq + type = rvar$model[[1]]$type + regime = rvar$regime regressors = colnames(dplyr::select(data, -date, -regime)) - residuals = rvar$residuals[[1]] - - if('const' %in% colnames(rvar$model[[1]]$coef) & - 'trend' %in% colnames(rvar$model[[1]]$coef)){ - type = 'both' - }else if('const' %in% colnames(rvar$model[[1]]$coef)){ - type = 'const' - }else if('trend' %in% colnames(rvar$model[[1]]$coef)){ - type = 'trend' - } + structure = rvar$structure p.lower = CI[1] p.upper = CI[2] @@ -708,135 +513,247 @@ rvar_irf = function( # set regime specific data coef = rvar$model[[paste0('regime_',regime.val)]]$coef - residuals = residuals %>% - dplyr::filter(model.regime == regime.val) + + residuals = rvar$residuals[[1]] %>% + dplyr::filter(model.regime == regime.val) %>% + dplyr::select(-model.regime) + is = data %>% dplyr::inner_join(dplyr::select(residuals, date), by = 'date') %>% - dplyr::rename(regime = regime) - residuals = residuals %>% - dplyr::select(-date, -model.regime, -forecast.date) + dplyr::select(-regime) + + if(!is.null(rvar$instrument)){ + instrument = dplyr::inner_join(rvar$instrument, dplyr::select(is, date), by = 'date') + }else{ + instrument = NULL + } + + model = rvar$model[[paste0('regime_',regime.val)]] ### calculate impulse responses -------------- - # estimate error-covariance matrix - cov.matrix = stats::var(stats::na.omit(residuals)) + + # estimate error-covariance matrix or structural impact matrix + B = + solve_B( + var = list( + model = model, + residuals = residuals, + structure = rvar$structure, + instrument = instrument, + instrumented = rvar$instrumented + ) + ) # estimate IRFs - irf = IRF(Phi = dplyr::select(coef, -y, -dplyr::contains('cosnt'), -dplyr::contains('trend')), - Sig = cov.matrix, - lag = horizon) + phi = dplyr::select(coef, -y, -dplyr::contains('const'), -dplyr::contains('trend')) + + irf = IRF_solve( + Phi = phi, + B = B, + lag = horizon, + structure = rvar$structure + ) # reorganize results - irf = data.frame(t(irf$orthirf)) - irf$shock = rep(regressors, horizon + 1) - irf$horizon = sort(rep(c(0:horizon), length(regressors))) - irf = irf %>% dplyr::arrange(shock, horizon) + irf = data.frame(t(irf)) + + # assign shocks + if(is.na(rvar$structure)){ + + irf$shock = rep(regressors, horizon + 1) + irf$horizon = sort(rep(c(0:horizon), length(regressors))) + + }else if(rvar$structure == 'IV-short'){ + + shocks = c(rvar$instrumented, regressors[!regressors %in% rvar$instrumented]) + irf$shock = rep(shocks, horizon + 1) + irf$horizon = sort(rep(c(0:horizon), length(regressors))) + + }else if(rvar$structure == 'IV'){ + + irf = irf[(c(1:(horizon))*length(regressors)-length(regressors)+1),] + irf$shock = rvar$instrumented + irf$horizon =c(0:(horizon-1)) + + }else if(rvar$structure == 'short'){ + + irf$shock = rep(regressors, horizon + 1) + irf$horizon = sort(rep(c(0:horizon), length(regressors))) + + } + + # observation identifiers rownames(irf) = NULL + colnames(irf) = c(regressors, 'shock', 'horizon') ### bootstrap irf standard errors -------------- - # see Lutkepohl (2005) + + # 0. set up parallel back-end + if(bootstrap.parallel == TRUE){ + if(bootstrap.cores == -1){ + n.cores = floor(future::availableCores() / 2) + }else{ + n.cores = bootstrap.cores + } + future::plan(future::multisession, workers = n.cores) + }else{ + future::plan(future::sequential) + } # 1. create bootstrap time series - bagged.series = as.list(1:bootstraps.num) %>% - purrr::map(.f = function(count){ - - # draw bootstrapped residuals - U = residuals[sample(c(1:nrow(residuals)), - size = nrow(residuals), - replace = TRUE),] - U = U %>% - dplyr::mutate_all(function(X){return(X-mean(X, na.rm = T))}) - - # create lags - X = is %>% - dplyr::select(-regime) %>% - n.lag(lags = p) %>% - dplyr::select(dplyr::contains('.l')) %>% - dplyr::slice(p:nrow(U)) - - if('const' %in% type | 'both' %in% type){X$const = 1} - if('trend' %in% type | 'both' %in% type){X$trend = c(1:nrow(X))} - - # estimate time series - Y = as.matrix(X) %*% as.matrix(t(coef[,-1])) - Y = Y + U - colnames(Y) = regressors - Y = data.frame(Y, date = is$date[1:nrow(Y)]) - - # return synthetic observations - return(Y) - - }) - - # 2. create bootstrapped residuals - bagged.irf = bagged.series %>% - purrr::map(.f = function(synth){ - - # estimate VAR with bootstrapped series + bagged.irf = as.list(1:bootstrap.num) %>% + furrr::future_map(.f = function(count){ + + # bootstrap residuals + if(bootstrap.type == 'wild' | + bootstrap.type == 'auto' & structure %in% c('IV', 'IV-short')){ + + # 'wild' bootstrap technique for simple distribution + # using observed scaled residuals with random sign flip. + # See the Rademacher distribution. + U = residuals[,!colnames(residuals) %in% c('date','forecast.date','model.regime')] + r = sample(c(-1,1), size = nrow(U), replace = T) + U = sweep(U, MARGIN = 1, r, `*`) + + # bootstrap instrument + if(!is.null(instrument)){ + instrument.bag = instrument + instrument.bag[,-1] = instrument.bag[,-1] * r + } + + }else if(bootstrap.type == 'standard' | + bootstrap.type == 'auto' & !structure %in% c('IV', 'IV-short')){ + + # standard bootstrap technique a al Lutkepohl (2005) + # draw residuals with replacement + U = stats::na.omit(residuals) + U = U[sample(c(1:nrow(U)), + size = nrow(residuals), + replace = TRUE),] + U = U %>% + dplyr::select(-date, -forecast.date, -model.regime) %>% + dplyr::mutate_all(function(X){return(X-mean(X, na.rm = T))}) + } + + # recursively build synthetic data + Y = matrix(nrow = nrow(is), ncol = ncol(is)-1) + Y = data.frame(Y); colnames(Y) = regressors + Y[1:p, ] = is[1:p, -1] + + for(i in (p+1):nrow(is)){ + + X = Y[(i-p):(i-1),] + X = stats::embed(as.matrix(X), dimension = p) + X = data.frame(X) + + if(type %in% c('const', 'both')){X$const = 1} + if(type %in% c('trend', 'both')){X$trend = c(1:nrow(X))} + + X.hat = as.matrix(coef[,-1]) %*% t(as.matrix(X)) + X.hat = t(X.hat) + + Y[i, ] = X.hat - U[i,] + + } + + Y$date = is$date + + # estimate VAR with synthetic data var.bag = VAR( - data = synth, + data = Y, p = p, horizon = 1, freq = freq, - type = type) + type = type + ) - # estimate error-covariance matrix - cov.matrix = stats::var(stats::na.omit(dplyr::select(var.bag$residuals[[1]], -date, -forecast.date))) + var.bag$structure = rvar$structure + var.bag$instrumented = rvar$instrumented + if(!is.null(instrument)){var.bag$instrument = instrument.bag} + + # estimate error-covariance matrix or structural impact matrix + B.bag = solve_B(var.bag, report_iv = FALSE) + + # set bagged coef matrix + coef.bag = dplyr::select(var.bag$model$coef, -y, -dplyr::contains('const'), -dplyr::contains('trend')) # estimate IRFs - irf = IRF(Phi = dplyr::select(coef, -y, -dplyr::contains('cosnt'), -dplyr::contains('trend')), - Sig = cov.matrix, - lag = horizon) + irf.bag = + IRF_solve( + Phi = coef.bag, + B = B.bag, + lag = horizon, + structure = rvar$structure + ) # reorganize results - irf = data.frame(t(irf$orthirf)) - irf$shock = rep(regressors, horizon + 1) - irf$horizon = sort(rep(c(0:horizon), length(regressors))) - irf = irf %>% dplyr::arrange(shock, horizon) - rownames(irf) = NULL - - return(irf) - - }) - - # 3. calculate confidence intervals - ci.lower = bagged.irf %>% - purrr::reduce(dplyr::bind_rows) %>% - dplyr::group_by(shock, horizon) %>% - dplyr::summarise_all(stats::quantile, p.lower, na.rm = T) %>% - dplyr::arrange(shock, horizon) %>% - tidyr::pivot_longer(cols = regressors, names_to = 'target', values_to = 'response.lower') %>% - dplyr::arrange(shock, horizon) + irf.bag = data.frame(t(irf.bag)) - ci.upper = bagged.irf %>% - purrr::reduce(dplyr::bind_rows) %>% - dplyr::group_by(shock, horizon) %>% - dplyr::summarise_all(stats::quantile, p.upper, na.rm = T) %>% - dplyr::arrange(shock, horizon) %>% - tidyr::pivot_longer(cols = regressors, names_to = 'target', values_to = 'response.upper') %>% - dplyr::arrange(shock, horizon) + # assign shocks + if(is.na(rvar$structure)){ - ci.med = bagged.irf %>% - purrr::reduce(dplyr::bind_rows) %>% - dplyr::group_by(shock, horizon) %>% - dplyr::summarise_all(stats::quantile, 0.5, na.rm = T) %>% - dplyr::arrange(shock, horizon) %>% - tidyr::pivot_longer(cols = regressors, names_to = 'target', values_to = 'response.med') %>% - dplyr::arrange(shock, horizon) + irf.bag$shock = rep(regressors, horizon + 1) + irf.bag$horizon = sort(rep(c(0:horizon), length(regressors))) + + }else if(rvar$structure == 'IV-short'){ + + shocks = c(rvar$instrumented, regressors[!regressors %in% rvar$instrumented]) + irf.bag$shock = rep(shocks, horizon + 1) + irf.bag$horizon = sort(rep(c(0:horizon), length(regressors))) + + }else if(rvar$structure == 'IV'){ + + irf.bag = irf.bag[(c(1:(horizon))*length(regressors)-length(regressors)+1),] + irf.bag$shock = rvar$instrumented + irf.bag$horizon = c(0:(horizon-1)) + }else if(rvar$structure == 'short'){ + + irf.bag$shock = rep(regressors, horizon + 1) + irf.bag$horizon = sort(rep(c(0:horizon), length(regressors))) + + } + + # observation identifiers + rownames(irf.bag) = NULL + colnames(irf.bag) = c(regressors, 'shock', 'horizon') + + return(irf.bag) + + }, + .options = furrr::furrr_options(seed = TRUE)) + + # 2. calculate confidence intervals + # collapse to data.frame + ci = bagged.irf %>% + purrr::reduce(dplyr::bind_rows) + + # estimate bands + ci = ci %>% + tidyr::pivot_longer(cols = regressors, names_to = 'target', values_to = 'response') %>% + dplyr::group_by(shock, horizon, target) %>% + dplyr::summarize( + response.lower = stats::quantile(response, probs = p.lower, na.rm = T), + response.upper = stats::quantile(response, probs = p.upper, na.rm = T), + response.mean = mean(response, na.rm = T) ) + + # combine point estimates and CI irf = irf %>% tidyr::pivot_longer(cols = regressors, names_to = 'target', values_to = 'response') %>% dplyr::arrange(shock, horizon) irf = - purrr::reduce( - list(irf, ci.lower, ci.upper, ci.med), - dplyr::full_join, - by = c('shock', 'target', 'horizon')) + dplyr::full_join( + irf, ci, + by = c('shock', 'target', 'horizon') + ) + # adjust for bias in CI + # (bias can be introduced in the bootstrapping if residuals are not actually mean zero) irf = irf %>% dplyr::mutate( - response.adjust = response - response.med, + response.adjust = response - response.mean, response.lower = response.lower + response.adjust, response.upper = response.upper + response.adjust ) %>% @@ -855,4 +772,3 @@ rvar_irf = function( } - diff --git a/README.md b/README.md index ac896e5..d9d9dfa 100644 --- a/README.md +++ b/README.md @@ -38,7 +38,8 @@ Local Projections (LP) 1. direct projection forecasting 1. impulse responses -Vector Auto-Regression (VAR) +Vector Auto-Regression (VAR), Structural VAR (SVAR), +and external instrument SVAR (Proxy-SVAR) 1. recursive forecasting 2. impulse responses 3. forecast error variance decomposition @@ -87,7 +88,9 @@ Vector Auto-Regression (VAR) # single-regime var #------------------------------------------ # estimate VAR - # (using IC lag selection) + # (using IC lag selection and short-term + # impact restrictions, i.e. calculate structurals + errors via Cholesky decomposition) var = VAR( data = Data, @@ -110,7 +113,7 @@ Vector Auto-Regression (VAR) # estimate IRF irf = - var_irf( + IRF( var, bootstraps.num = 10, CI = c(0.05,0.95)) @@ -120,13 +123,19 @@ Vector Auto-Regression (VAR) # estimate forecast error variance decomposition fevd = - var_fevd( + FEVD( var, horizon = 10) # plot FEVD plot_fevd(fevd) + # estimate historical decomposition + hd = HD(var) + + # plot HD + plot_hd(hd) + #------------------------------------------- # multi-regime var #------------------------------------------- @@ -141,7 +150,7 @@ Vector Auto-Regression (VAR) # estimate IRF rvar.irf = - rvar_irf( + IRF( rvar, horizon = 10, bootstraps.num = 10, @@ -155,7 +164,7 @@ Vector Auto-Regression (VAR) # estimate forecast error variance decomposition rvar.fevd = - rvar_fevd( + FEVD( rvar, horizon = 10) @@ -165,6 +174,9 @@ Vector Auto-Regression (VAR) # regime 2: high interest rates plot_fevd(rvar.fevd[[2]]) + # estimate HD + rvar.hd = HD(rvar) + #------------------------------------------- # single-regime local projections #------------------------------------------- @@ -178,7 +190,7 @@ Vector Auto-Regression (VAR) freq = 'month') # estimate single-regime IRF - lp.irf = lp_irf(lp) + lp.irf = IRF(lp) # plot IRF plot_irf(lp.irf) @@ -196,7 +208,7 @@ Vector Auto-Regression (VAR) freq = 'month') # estimate multi-regime IRF - rlp.irf = rlp_irf(rlp) + rlp.irf = IRF(rlp) # plot IRF # regime 1: low interest rates @@ -207,4 +219,4 @@ Vector Auto-Regression (VAR) --- ## Contact -If you should have questions, concerns, or wish to collaborate, please contact [Tyler J. Pike](https://tylerjpike.github.io/) +If you should have questions, concerns, or wish to collaborate, please contact [Tyler J. Pike](https://tylerjpike.github.io/) \ No newline at end of file diff --git a/_pkgdown.yml b/_pkgdown.yml index d114765..78cd453 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -27,14 +27,9 @@ reference: - title: Model and Forecast Analysis desc: Functions to anlayze model and forecasting dynamics contents: - - '`var_irf`' - - '`rvar_irf`' - - '`var_fevd`' - - '`rvar_fevd`' - - '`var_hd`' - - '`rvar_hd`' - - '`lp_irf`' - - '`rlp_irf`' + - '`IRF`' + - '`FEVD`' + - '`HD`' - title: Plots desc: Functions to plot model and forecast analysis contents: diff --git a/cran-comments.md b/cran-comments.md new file mode 100644 index 0000000..56cf847 --- /dev/null +++ b/cran-comments.md @@ -0,0 +1,15 @@ +# CRAN Comments + +## Test environments +* local Windows install, R 4.0.5 +* win-builder (devel and release) +* Ubuntu 16.04.6 (on travis-ci), R 4.0.2 +* R-hub Ubuntu Linux 20.04.1 LTS, R-release +* R-hub Fedora Linux, R-devel +* R-hub Windows Server 2008 R2 SP1, R-devel + +## R CMD check results +There were no ERRORs, WARNINGs, or NOTEs. + +## Downstream dependencies +There are currently no downstream dependencies for this package. \ No newline at end of file diff --git a/development_test.R b/development_test.R deleted file mode 100644 index 38cf9da..0000000 --- a/development_test.R +++ /dev/null @@ -1,1047 +0,0 @@ - -TVAR = function( - data, # data.frame, matrix, ts, xts, zoo: Endogenous regressors - threshold, # string: variable to estimate threshold values - horizon = 10, # int: forecast horizons - freq = 'month', # string: frequency of data (day, week, month, quarter, year) - type = 'const', # string: type of deterministic terms to add ('none', 'const', 'trend', 'both') - p = 1, # int: lags - lag.ic = NULL, # string: information criterion to choose the optimal number of lags ('AIC' or 'BIC') - lag.max = NULL, # int: maximum number of lags to test in lag selection - structure = 'short', # string: type of structural identification strategy to use in model analysis (NULL, 'short', or 'IV') - instrument = NULL, # string: name of instrumental variable contained in the data matrix - instrumented = NULL # string: name of variable to be instrumented in IV and IV-short procedure; default is the first non-date variable in data -){ - - # function warnings - if(!is.numeric(p) | p %% 1 != 0){ - stop('p must be an integer') - } - if(!is.null(lag.ic)){ - if(!lag.ic %in% c('BIC','AIC')){ - stop("lag.ic must be either 'BIC', 'AIC', or NULL") - } - } - if(!is.null(lag.max)){ - if(lag.max %% 1 != 0){ - stop('lag.max must be an integer if IC-based lag selection is used') - } - } - if(!is.matrix(data) & !is.data.frame(data)){ - stop('data must be a matrix or data.frame') - } - if(!is.numeric(p) | p %% 1 != 0){ - stop('p must be an integer') - } - if(!is.numeric(horizon) | horizon %% 1 != 0 | horizon <= 0){ - stop('horizon must be a positive integer') - } - if(!freq %in% c('day','week','month','quarter','year')){ - stop("freq must be one of the following strings: 'day','week','month','quarter','year'") - } - if(!structure %in% c('short', 'IV', 'IV-short') & !is.null(structure)){ - stop("strucutre must be one of 'strucutre', 'IV', 'IV-short', or NULL.") - } - if(!is.null(instrument)){ - if(!instrument %in% colnames(data)){ - stop("instrument must be the name of a variable found in data.") - } - } - if(!is.null(instrumented)){ - if(!instrumented %in% colnames(data)){ - stop("instrumented must be the name of a variable found in data.") - } - } - - # cast as data frame if ts, xts, or zoo object - if(stats::is.ts(data) | xts::is.xts(data) | zoo::is.zoo(data)){ - data = data.frame(date = zoo::index(date), data) - } - - # set aside instruments - if(!is.null(instrument)){ - data.instrument = dplyr::select(data, date, dplyr::all_of(instrument)) - data = dplyr::select(data, -dplyr::all_of(instrument)) - }else{ - data.instrument = NULL - } - - # detect variable to be instrumented - if(is.null(instrumented)){ - var_to_instrument = colnames(dplyr::select(data, -date))[1] - }else{ - var_to_instrument = instrumented - } - - - for(t in unique(data[,c(threshold)])){ - - data.threshold = data %>% - mutate(t_regime = if_else(threshold >= t, 1, 0)) - - rvar = - RVAR( - data = data.threshold, - horizon = horizon, - freq = freq, - type = type, - p = p, - lag.ic = lag.ic, - lag.max = lag.max, - regime = t_regime, - ) - - rvar.ll = rvar$ll - - } - - -} -##################################################################### -##################################################################### -##################################################################### - -solve_B = function(var, report_iv = TRUE){ - - # set residuals - if(class(var) == 'VAR'){ - residuals = var$residuals[[1]] - }else{ - residuals = var$residuals - } - - if(is.null(var$structure) == TRUE){ - - # retrieve reduced-form residuals - data.residuals = residuals - - # reduced form variance-covariance matrix - cov.matrix = stats::var(stats::na.omit(dplyr::select(data.residuals, -date, -forecast.date))) - - B = cov.matrix - - return(B) - - }else if(var$structure == 'short'){ - - # retrieve reduced-form residuals - data.residuals = residuals - - # reduced form variance-covariance matrix - cov.matrix = stats::var(stats::na.omit(dplyr::select(data.residuals, -date, -forecast.date))) - - # take cholesky decomposition - B = t(chol(cov.matrix)) - - return(B) - - }else if(var$structure == 'IV'){ - - # retrieve reduced-form residuals - data.residuals = residuals - covariates = colnames(dplyr::select(data.residuals, -date, -forecast.date)) - col_to_instrument = var$instrumented - - # retrieve instrument - data.instrument = var$instrument - instrument = colnames(dplyr::select(data.instrument, -date)) - - # combine data - data = - dplyr::inner_join(data.residuals, data.instrument, by = 'date') %>% - na.omit() - - # first stage least squares - model.first_stage = - lm(col_to_instrument ~., - data = data %>% - select(col_to_instrument = dplyr::all_of(col_to_instrument), all_of(instrument)) %>% - na.omit()) - p_hat = model.first_stage$fitted.values - - if(report_iv == TRUE){ - print(summary(model.first_stage)) - } - - # second stage least squares - # to estimate first column of B - # (automatically scales first entry to 1) - instrumented_shocks = covariates %>% - purrr::map_df(.f = function(X){ - - model.second_stage = lm(data[,X] ~ p_hat) - second_stage_beta = model.second_stage$coefficients[2] - return(second_stage_beta) - - }) %>% - as.vector() - - # scale size of the shock - # see Gertler and Karadi (2015) for background - # see Cesa-Bianchi's VAR-Toolbox for MATLAB implementation - X.demean = - select(data, -date, -forecast.date, -dplyr::all_of(instrument)) %>% - mutate_all(function(X){return(X - mean(X, na.rm = T))}) %>% - as.matrix() - - sigma_b = 1/(nrow(data) - ncol(var$model$coef) + 1) * t(X.demean) %*% X.demean - - s21s11 = instrumented_shocks[2:length(covariates),] %>% as.matrix() - S11 = sigma_b[1,1] %>% as.numeric() - S21 = sigma_b[2:nrow(sigma_b),1] %>%as.vector() - S22 = sigma_b[2:nrow(sigma_b),2:ncol(sigma_b)] - - Q = (s21s11 * S11) %*% t(s21s11) - (S21 %*% t(s21s11) + as.matrix(s21s11) %*% t(S21)) + S22 - sp = sqrt( S11 - t(S21 - as.matrix(s21s11) %*% S11) %*% as.matrix(solve(Q) %*% (S21 - s21s11 * S11)) ) - - scaled_instrumented_shocks = instrumented_shocks * as.numeric(sp) - - # prepare B matrix - B = matrix(0, ncol = (length(covariates) - 1), nrow = length(covariates)) - B = cbind(scaled_instrumented_shocks, B) - - return(B) - - }else if(var$structure == 'IV-short'){ - - # align instrument and residuals - valid_dates = - dplyr::inner_join( - dplyr::select(na.omit(residuals), date), - dplyr::select(na.omit(var$instrument), date), - by = 'date' - ) - - # set residuals matrix - residuals = residuals %>% - dplyr::inner_join(valid_dates, by = 'date') %>% - dplyr::select(-date, -forecast.date) %>% - as.matrix() - - residuals = -1*residuals - - # set instrument - instrument = var$instrument %>% - data.frame() %>% - dplyr::inner_join(valid_dates, by = 'date') %>% - dplyr::select(-date) - - instrument.mean = instrument %>% - dplyr::mutate_all(mean, na.rm = T) - - # number of observations - n.obs = nrow(residuals) - - # solve for IV implied impact - pi = solve(n.obs^(-1) * t(residuals) %*% residuals) %*% - (n.obs^(-1) * t(residuals) %*% as.matrix(instrument - instrument.mean)) - - phi_sq = - (n.obs^(-1) * t(instrument - instrument.mean) %*% residuals) %*% - solve( n.obs^(-1) * t(residuals) %*% residuals ) %*% - ( n.obs^(-1) * t(residuals) %*% as.matrix(instrument - instrument.mean) ) - - B1 = - n.obs^(-1) * - ( t(residuals) %*% as.matrix(instrument - instrument.mean) ) %*% - (phi_sq)^(-0.5) - - B = matrix(ncol = ncol(residuals), nrow = ncol(residuals)) - B[,1:ncol(instrument)] = B1 - rownames(B) = colnames(B) = colnames(residuals) - - # solve for orthogonalized structural shock - model.first_stage = lm(instrument[,1] ~ residuals) - orthogonal_instrument = instrument - model.first_stage$residuals - orthogonal_instrument = orthogonal_instrument[,] / sd(orthogonal_instrument[,], na.rm = T) - - shock_matrix = matrix(nrow = nrow(residuals), ncol = ncol(residuals)) - shock_matrix[,1] = orthogonal_instrument - - # reduced form variance-covariance matrix - sigma = stats::var(residuals) - - # solve additional entries - # with a lower triangular restriction - order_sequence = c(1:ncol(residuals)) - - for(i in order_sequence){ - - Y = residuals[,i] - X = shock_matrix[,c(1:i)] - model.second_stage = lm(Y~X) - - if(i != tail(order_sequence,1)){ - - B[i,] = - c( model.second_stage$coef[-1], - sd(model.second_stage$residuals), - rep(0, length(order_sequence) - ncol(instrument) - i) ) - - shock_matrix[,i+1] = model.second_stage$residuals / sd(model.second_stage$residuals) - - }else{ - - B[i,] = model.second_stage$coef[-1] - - } - - } - - # make diagonal entries positive - shock_ordering = data.frame(residuals) %>% - select(var$instrumented, dplyr::everything()) %>% - colnames() - B.sign = diag(sign(diag(B[shock_ordering,]))) - B = B %*% B.sign - - return(B) - - }else{ - - stop('structure must be set to either NULL, "short", or "IV".') - - } - -} - - -RVAR_estimate = function( - data, # data.frame, matrix, ts, xts, zoo: Endogenous regressors - regime, # string: name or regime assignment vector in the design matrix (data) - p = 1, # int: lags - horizon = 10, # int: forecast horizons - freq = 'month', # string: frequency of data (day, week, month, quarter, year) - type = 'const' # string: type of deterministic terms to add ('none', 'const', 'trend', 'both') -){ - - # function variables - term = estimate = std.error = model.regime = NULL - - # declare regressors - regressors = colnames(dplyr::select(data, -date, -regime)) - - # create regressors - Y = data.frame(data) %>% - dplyr::select(regressors, date) %>% - n.lag(lags = p) %>% - dplyr::full_join( - dplyr::select(data, regime = regime, date), - by = 'date') - - # add deterministic components - if('const' %in% type | 'both' %in% type){Y$const = 1} - if('trend' %in% type | 'both' %in% type){Y$trend = c(1:nrow(Y))} - - # detect regime values - regimes = unique(Y$regime) - - ### estimate coefficients ---------------------- - models = Y %>% - # split by regime - dplyr::group_split(regime) %>% - purrr::map(.f = function(Y){ - - regime.val = unique(Y$regime) - - # calculate equation by equation - models = as.list(regressors) %>% - purrr::map(.f = function(target){ - - X = Y %>% - dplyr::select( - dplyr::contains('.l'), target = target, - dplyr::contains('const'), dplyr::contains('trend')) - - # estimate OLS - model = stats::lm(target ~ . - 1, data = X) - - # coefficients - c = broom::tidy(model) %>% dplyr::select(term, coef = estimate) - c$y = target - - # standard errors - se = broom::tidy(model) %>% dplyr::select(term, std.error) - se$y = target - - # log likelihood - ll = stats::logLik(model)[1] - - # return results - return(list(coef = c, se = se, ll = ll)) - - }) - - # extract coefficients - coef = - purrr::map(models, .f = function(X){return(X$coef)}) %>% - purrr::reduce(dplyr::bind_rows) %>% - tidyr::pivot_wider(values_from = coef, names_from = term) - - - # extract coefficients' standard errors - se = - purrr::map(models, .f = function(X){return(X$se)}) %>% - purrr::reduce(dplyr::bind_rows) %>% - tidyr::pivot_wider(values_from = std.error, names_from = term) - - # extract log likelihood - ll = - purrr::map(models, .f = function(X){return(X$ll)}) %>% - purrr::reduce(rbind) %>% - sum() - - # package for return - model = list(coef = coef, se = se, p = p, freq = freq, type = type, ll = ll, horizon = horizon, regime = regime.val) - - }) - - names(models) = paste0('regime_', unique(Y$regime)) - - - fr = as.list(regimes) %>% - purrr::map(.f = function(regime.val){ - - coef = models[[paste0('regime_', regime.val)]]$coef - - ### estimate forecasts ----------------------- - forecasts = list() - for(i in 1:horizon){ - - # update X - if(i == 1){ - - X = Y %>% - dplyr::select( - dplyr::contains('.l'), - dplyr::contains('const'), - dplyr::contains('trend')) - - X.date = data$date - - }else{ - - X = - dplyr::bind_rows( - dplyr::select(data[1:i,], date, regressors), - dplyr::select(forecast_prev[i+1:nrow(forecast_prev),], date = forecast.date, regressors)) %>% - n.lag(lags = p) - - X.date = X$date - - X = X %>% - dplyr::filter(!is.na(date)) %>% - dplyr::select(dplyr::contains('.l')) - - if('const' %in% type | 'both' %in% type){X$const = 1} - if('trend' %in% type | 'both' %in% type){X$trend = c(1:nrow(X)) + (i-1)} - - } - - # estimate i-step ahead forecast - forecast = as.matrix(X) %*% as.matrix(t(coef[,-1])) - colnames(forecast) = regressors - - # set forecast date - if(i == 1){ - forecast.date = stats::na.omit(X.date) - }else{ - forecast.date = - forecast_date( - forecast.date = stats::na.omit(X.date), - horizon = i-1, - freq = freq - ) - } - - # add in dates - forecast = - data.frame( - date = data$date, - forecast.date = forecast.date, - forecast - ) %>% - dplyr::left_join(dplyr::select(Y, date, model.regime = regime), by = 'date') - - # store forecasts - forecasts[[paste0('H_',i)]] = forecast - forecast_prev = forecast - - } - - ### calculate residuals ----------------------- - residuals = forecasts %>% - # error by forecast horizon - purrr::map(.f = function(forecast){ - - error = data.frame(forecast) - error[,c(regressors)] = forecast[,c(regressors)] - data.frame(data)[, c(regressors)] - - return(error) - - }) - - ### return output ----------------------- - residuals = residuals %>% - purrr::map(.f = function(r){ - return( dplyr::filter(r, regime.val == model.regime) ) - }) - - forecasts = forecasts %>% - purrr::map(.f = function(f){ - return( dplyr::filter(f, regime.val == model.regime) ) - }) - - return( - list( - forecasts = forecasts, - residuals = residuals - ) - ) - - }) - - ### Organize and return output ------- - forecasts = as.list(c(1:horizon)) %>% - purrr::map(.f = function(h){ - - f = fr %>% - purrr::map(.f = function(r){ - return(r$forecast[[paste0('H_',h)]]) - }) %>% - purrr::reduce(dplyr::bind_rows) %>% - dplyr::arrange(date) - - }) - - - residuals = as.list(c(1:horizon)) %>% - purrr::map(.f = function(h){ - - f = fr %>% - purrr::map(.f = function(r){ - return(r$residuals[[paste0('H_',h)]]) - }) %>% - purrr::reduce(dplyr::bind_rows) %>% - dplyr::arrange(date) - - }) - - names(forecasts) = paste0('H_',c(1 : horizon)) - names(residuals) = paste0('H_',c(1 : horizon)) - - ### return output -------------- - return( - list( - model = models, - data = data, - forecasts = forecasts, - residuals = residuals, - regime = regime - ) - ) -} - -RVAR = function( - data, # data.frame, matrix, ts, xts, zoo: Endogenous regressors - horizon = 10, # int: forecast horizons - freq = 'month', # string: frequency of data (day, week, month, quarter, year) - type = 'const', # string: type of deterministic terms to add ('none', 'const', 'trend', 'both') - p = 1, # int: lags - lag.ic = NULL, # string: information criterion to choose the optimal number of lags ('AIC' or 'BIC') - lag.max = NULL, # int: maximum number of lags to test in lag selection - regime = NULL, # string: name or regime assignment vector in the design matrix (data) - regime.method = 'rf', # string: regime assignment technique ('rf', 'kmeans', 'EM', 'BP') - regime.n = 2, # int: number of regimes to estimate (applies to kmeans and EM) - structure = 'short', # string: type of structural identification strategy to use in model analysis (NULL, 'short', or 'IV') - instrument = NULL, # string: name of instrumental variable contained in the data matrix - instrumented = NULL # string: name of variable to be instrumented in IV and IV-short procedure; default is the first non-date variable in data -){ - - # function warnings - if(!is.numeric(p) | p %% 1 != 0){ - stop('p must be an integer') - } - if(!is.null(lag.ic)){ - if(!lag.ic %in% c('BIC','AIC')){ - stop("lag.ic must be either 'BIC', 'AIC', or NULL") - } - } - if(!is.null(lag.max)){ - if(lag.max %% 1 != 0){ - stop('lag.max must be an integer if IC-based lag selection is used') - } - } - if(!structure %in% c('short', 'IV', 'IV-short') & !is.null(structure)){ - stop("strucutre must be one of 'strucutre', 'IV', 'IV-short', or NULL.") - } - if(!is.null(instrument)){ - if(!instrument %in% colnames(data)){ - stop("instrument must be the name of a variable found in data.") - } - } - if(!is.null(instrumented)){ - if(!instrumented %in% colnames(data)){ - stop("instrumented must be the name of a variable found in data.") - } - } - if(!is.matrix(data) & !is.data.frame(data)){ - stop('data must be a matrix or data.frame') - } - if(!is.numeric(p) | p %% 1 != 0){ - stop('p must be an integer') - } - if(!is.numeric(horizon) | horizon %% 1 != 0 | horizon <= 0){ - stop('horizon must be a positive integer') - } - if(!freq %in% c('day','week','month','quarter','year')){ - stop("freq must be one of the following strings: 'day','week','month','quarter','year'") - } - - # cast as data frame if ts, xts, or zoo object - if(stats::is.ts(data) | xts::is.xts(data) | zoo::is.zoo(data)){ - data = data.frame(date = zoo::index(date), data) - } - - # set aside instruments - if(!is.null(instrument)){ - data.instrument = dplyr::select(data, date, dplyr::all_of(instrument)) - data = dplyr::select(data, -dplyr::all_of(instrument)) - }else{ - data.instrument = NULL - } - - # detect variable to be instrumented - # (defaults to first non-date column if missing) - if(is.null(instrumented)){ - var_to_instrument = colnames(dplyr::select(data, -date, -regime))[1] - }else{ - var_to_instrument = instrumented - } - - # create regimes - if(is.null(regime)){ - - data = - regimes( - data, - regime.n = regime.n, - method = regime.method) - - regime = 'regime' - - } - - # create VAR - if(!is.null(lag.ic) & !is.null(lag.max)){ - - ic.scores = vector(length = lag.max+1) - - models = c(1:lag.max) %>% - purrr::map(.f = function(p){ - - # estimate candidate model - model = - RVAR_estimate( - data = data, - p = p, - regime = regime, - horizon = horizon, - freq = freq, - type = type - ) - - # calculate IC - ic.score = - IC( - ic = lag.ic, - errors = model$residuals[[1]], - data = data, - p = p - ) - - ic.scores[p] = ic.score - - # return candidate model - return(model) - - }) - - # return IC minimizing VAR - min.ic = which.min(ic.scores) - model = models[[min.ic]] - - }else{ - - model = - RVAR_estimate( - data = data, - p = p, - regime = regime, - horizon = horizon, - freq = freq, - type = type - ) - } - - # add structure - model$structure = structure - model$instrument = data.instrument - model$instrumented = var_to_instrument - - class(model) = 'RVAR' - return(model) - - -} - -rvar_irf = function( - rvar, # threshold VAR output - horizon = 10, # int: number of periods - CI = c(0.1, 0.9), # numeric vector: c(lower ci bound, upper ci bound) - bootstrap.type = 'auto', # string: bootstrapping technique to use ('auto', 'standard', or 'wild'); if auto then wild is used for IV or IV-short, else standard is used - bootstrap.num = 100, # int: number of bootstraps - bootstrap.parallel = FALSE, # boolean: create IRF draws in parallel - bootstrap.cores = -1 # int: number of cores to use in parallel processing; -1 detects and uses half the available cores -){ - - # function warnings - if(!is.numeric(bootstrap.num) | bootstrap.num %% 1 != 0){ - stop('bootstrap.num must be an integer') - } - if(!is.numeric(CI) | length(CI) != 2 | min(CI) < 0 | max(CI) > 1 | is.na(sum(CI))){ - stop('CI must be a two element numeric vector bound [0,1]') - } - if(!is.numeric(horizon) | horizon %% 1 != 0 | horizon <= 0){ - stop('horizon must be a positive integer') - } - - # function variables - y = forecast.date = shock = target = response = response.lower = response.upper = model.regime = response.med = response.adjust = NULL - - # set data - data = rvar$data - p = rvar$model[[1]]$p - freq = rvar$model[[1]]$freq - regime = rvar$regime - regressors = colnames(dplyr::select(data, -date, -regime)) - structre = rvar$structure - - if('const' %in% colnames(rvar$model[[1]]$coef) & - 'trend' %in% colnames(rvar$model[[1]]$coef)){ - type = 'both' - }else if('const' %in% colnames(rvar$model[[1]]$coef)){ - type = 'const' - }else if('trend' %in% colnames(rvar$model[[1]]$coef)){ - type = 'trend' - } - - p.lower = CI[1] - p.upper = CI[2] - - regimes = dplyr::select(data, regime = regime) - regimes = unique(regimes$regime) - - # estimate impulse responses by regime - results = as.list(regimes) %>% - purrr::map(.f = function(regime.val){ - - # set regime specific data - coef = rvar$model[[paste0('regime_',regime.val)]]$coef - - residuals = rvar$residuals[[1]] %>% - dplyr::filter(model.regime == regime.val) %>% - dplyr::select(-model.regime) - - is = data %>% - dplyr::inner_join(dplyr::select(residuals, date), by = 'date') %>% - dplyr::select(-regime) - - if(!is.null(rvar$instrument)){ - instrument = dplyr::inner_join(rvar$instrument, dplyr::select(is, date), by = 'date') - } - - ### calculate impulse responses -------------- - - # estimate error-covariance matrix or structural impact matrix - B = - solve_B( - var = list( - residuals = residuals, - structure = rvar$structure, - instrument = instrument, - instrumented = rvar$instrumented - ) - ) - - # estimate IRFs - phi = dplyr::select(coef, -y, -dplyr::contains('const'), -dplyr::contains('trend')) - - irf = IRF(Phi = phi, - B = B, - lag = horizon, - structure = rvar$structure) - - # reorganize results - irf = data.frame(t(irf)) - - if(!is.null(rvar$structure)){ - if(rvar$structure == 'IV-short'){ - shocks = c(rvar$instrumented, regressors[!regressors %in% rvar$instrumented]) - irf$shock = rep(shocks, horizon + 1) - } - }else{ - irf$shock = rep(regressors, horizon + 1) - } - - irf$horizon = sort(rep(c(0:horizon), length(regressors))) - irf = dplyr::arrange(irf, shock, horizon) - rownames(irf) = NULL - colnames(irf) = c(regressors, 'shock', 'horizon') - - ### bootstrap irf standard errors -------------- - - # 0. set up parallel back-end - if(bootstrap.parallel == TRUE){ - if(bootstrap.cores == -1){ - n.cores = floor(future::availableCores() / 2) - }else{ - n.cores = bootstrap.parallel - } - future::plan(future::multisession, workers = n.cores) - }else{ - future::plan(future::sequential) - } - - # 1. create bootstrap time series - bagged.irf = as.list(1:bootstrap.num) %>% - furrr::future_map(.f = function(count){ - - # bootstrap residuals - if(bootstrap.type == 'wild' | - bootstrap.type == 'auto' & structure == 'IV' | - bootstrap.type == 'auto' & structure == 'IV-short'){ - - # 'wild' bootstrap technique for simple distribution - # using observed scaled residuals with random sign flip. - # See the Rademacher distribution. - U = residuals[,!colnames(residuals) %in% c('date','forecast.date','model.regime')] - r = sample(c(-1,1), size = nrow(U), replace = T) - U = sweep(U, MARGIN = 1, r, `*`) - - # bootstrap instrument - if(!is.null(instrument)){ - instrument.bag = instrument - instrument.bag[,-1] = sweep(instrument.bag[,-1], MARGIN = 1, r, `*`) - } - - }else if(bootstrap.type == 'standard' | - bootstrap.type == 'auto' & structure != 'IV' & structure != 'IV-short'){ - - # standard bootstrap technique a al Lutkepohl (2005) - # draw residuals with replacement - U = na.omit(residuals) - U = U[sample(c(1:nrow(U)), - size = nrow(residuals), - replace = TRUE),] - U = U %>% - dplyr::select(-date, -forecast.date, -model.regime) %>% - dplyr::mutate_all(function(X){return(X-mean(X, na.rm = T))}) - } - - # recursively build synthetic data - Y = matrix(nrow = nrow(is), ncol = ncol(is)-1) - Y = data.frame(Y); colnames(Y) = regressors - Y[1:p, ] = is[1:p, -1] - - for(i in (p+1):nrow(is)){ - - X = Y[(i-p):(i-1),] - X = embed(as.matrix(X), dimension = p) - X = data.frame(X) - - if(type %in% c('const', 'both')){X$const = 1} - if(type %in% c('trend', 'both')){X$trend = c(1:nrow(X))} - - X.hat = as.matrix(coef[,-1]) %*% t(as.matrix(X)) - X.hat = t(X.hat) - - Y[i, ] = X.hat - U[i,] - - } - - Y$date = is$date - - # estimate VAR with synthetic data - var.bag = - VAR( - data = Y, - p = p, - horizon = 1, - freq = freq, - type = type - ) - - var.bag$structure = rvar$structure - var.bag$instrumented = rvar$instrumented - if(!is.null(instrument)){var.bag$instrument = instrument.bag} - - # estimate error-covariance matrix or structural impact matrix - B.bag = solve_B(var.bag, report_iv = FALSE) - - # set bagged coef matrix - coef.bag = dplyr::select(var.bag$model$coef, -y, -dplyr::contains('const'), -dplyr::contains('trend')) - - # estimate IRFs - irf.bag = IRF(Phi = coef.bag, - B = B.bag, - lag = horizon, - structure = rvar$structure) - - # reorganize results - irf.bag = data.frame(t(irf.bag)) - - if(!is.null(rvar$structure)){ - if(rvar$structure == 'IV-short'){ - shocks = c(rvar$instrumented, regressors[!regressors %in% rvar$instrumented]) - irf.bag$shock = rep(shocks, horizon + 1) - } - }else{ - irf.bag$shock = rep(regressors, horizon + 1) - } - - irf.bag$horizon = sort(rep(c(0:horizon), length(regressors))) - irf.bag = dplyr::arrange(irf.bag, shock, horizon) - rownames(irf.bag) = NULL - colnames(irf.bag) = c(regressors, 'shock', 'horizon') - - return(irf.bag) - - }) - - # 2. calculate confidence intervals - # collapse to data.frame - ci = bagged.irf %>% - purrr::reduce(dplyr::bind_rows) - - # remove unused shocks in the case of IV - if(!is.null(rvar$structure)){ - if(rvar$structure == 'IV'){ - ci = dplyr::filter(ci, shock == rvar$instrumented) - } - } - - # estimate bands - ci = ci %>% - tidyr::pivot_longer(cols = regressors, names_to = 'target', values_to = 'response') %>% - group_by(shock, horizon, target) %>% - summarize(response.lower = quantile(response, probs = p.lower, na.rm = T), - response.upper = quantile(response, probs = p.upper, na.rm = T), - response.mean = mean(response, na.rm = T) ) - - # combine point estimates and CI - irf = irf %>% - tidyr::pivot_longer(cols = regressors, names_to = 'target', values_to = 'response') %>% - dplyr::arrange(shock, horizon) - - irf = - full_join( - irf, ci, - by = c('shock', 'target', 'horizon') - ) - - # adjust for bias in CI - # (bias can be introduced in the bootstrapping if residuals are not actually mean zero) - irf = irf %>% - dplyr::mutate( - response.adjust = response - response.mean, - response.lower = response.lower + response.adjust, - response.upper = response.upper + response.adjust - ) %>% - dplyr::select(target, shock, horizon, response.lower, response, response.upper) %>% - dplyr::arrange(target, shock, horizon) - - ### return output -------------- - return(irf) - - }) - - names(results) = paste0('regime_', regimes) - - ### return output -------------- - return(results) - -} - -plot_individual_irf = function( - irf, # irf object - shock.var, # string: name of variable to treat as the shock - response.var, # string: name of variable to treat as the response - title, # string: title of the chart - ylab # string: y-axis label -){ - - # function variables - response = error = horizon = shock = target = response.lower = response.upper = NULL - - # filter for one shock and one target - plotdata = irf %>% - dplyr::filter(shock == shock.var & target == response.var) %>% - dplyr::select(horizon, response, response.lower, response.upper) - - # plot - irf.plot <- plotdata %>% - ggplot2::ggplot(ggplot2::aes(x=horizon, y=response, ymin=response.lower, ymax=response.upper)) + - ggplot2::geom_hline(yintercept = 0, color="black") + - ggplot2::geom_ribbon(fill="royalblue", alpha=0.2) + - ggplot2::geom_line(color = 'darkred', lwd = 1) + - ggplot2::theme_light() + - ggplot2::ggtitle(title)+ - ggplot2::ylab(ylab)+ - ggplot2::xlab("") + - ggplot2::theme(plot.title = ggplot2::element_text(size = 11, hjust=0.5), - axis.title.y = ggplot2::element_text(size=11)) - - irf.plot -} - -##################################################################### -# TEST - -library(tidyverse) -library(lubridate) - -# Data -data.macro = read_csv('/scratch/m1tjp01/Andrea/FinancialStability/Data/Intermediate/model_monthly_macro_covariates.csv') -data.macro = data.macro %>% filter(year(date) >= 1991) -data.macro = data.macro %>% select(date, PCE, UNEMP, PR, EBP, mp_shock) - -# MODEL -# var = -# VAR( -# data = data.macro, -# horizon = 1, -# freq = 'month', -# p = 1, -# structure = 'IV-short', -# instrument = 'mp_shock', -# instrumented = 'PR' -# ) -# -# irf = var_irf(var, horizon = 20, bootstrap.type = 'wild', bootstrap.num = 10, bootstrap.parallel = TRUE) -# plot_irf(irf, shocks = 'PR') - -rvar = - RVAR( - data = data.macro, - horizon = 1, - freq = 'month', - p = 1, - structure = 'IV-short', - instrument = 'mp_shock', - instrumented = 'PR' - ) - -irf = rvar_irf(rvar, horizon = 20, bootstrap.type = 'wild', bootstrap.num = 10, bootstrap.parallel = TRUE) -plot_irf(irf[[1]], shocks = 'PR') - diff --git a/docs/LICENSE.html b/docs/LICENSE.html index 7d15b60..5c16da1 100644 --- a/docs/LICENSE.html +++ b/docs/LICENSE.html @@ -87,7 +87,7 @@ sovereign - 1.1.0 + 1.2.0 diff --git a/docs/authors.html b/docs/authors.html index 3c04416..86c09b5 100644 --- a/docs/authors.html +++ b/docs/authors.html @@ -87,7 +87,7 @@ sovereign - 1.1.0 + 1.2.0 diff --git a/docs/index.html b/docs/index.html index 5fc025c..b62cdf2 100644 --- a/docs/index.html +++ b/docs/index.html @@ -17,7 +17,7 @@ + historical decompositions, and forecast error variance decompositions. "> + + + + + + + +Estimate forecast error variance decomposition — FEVD • sovereign + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ + + + +
+ +
+
+
+

Estimate forecast error variance decomposition

+ Source: R/analysis_wrappers.R + +
+ +
+

Estimate the forecast error variance decomposition for VARs with +either short or 'IV-short' structural errors. See VAR +and RVAR documentation for details regarding structural errors.

+
+ + 
FEVD(model, horizon = 10, scale = TRUE)
+ +

Arguments

+ + + + + + + + + + + + + + +
model

VAR or RVAR class object

horizon

int: number of periods

scale

boolean: scale variable contribution as percent of total error

+ +

Value

+ +

long-form data.frame

+

See also

+ +

VAR()

+

var_fevd()

+

RVAR()

+

rvar_fevd()

+ +

Examples

+ 
# \donttest{
+
+ # simple time series
+ AA = c(1:100) + rnorm(100)
+ BB = c(1:100) + rnorm(100)
+ CC = AA + BB + rnorm(100)
+ date = seq.Date(from = as.Date('2000-01-01'), by = 'month', length.out = 100)
+ Data = data.frame(date = date, AA, BB, CC)
+
+ # estimate VAR
+  var =
+    sovereign::VAR(
+      data = Data,
+      horizon = 10,
+      freq = 'month',
+      lag.ic = 'BIC',
+      lag.max = 4)
+
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion

+# impulse response functions
+var.irf = sovereign::IRF(var)
+
#> 
+#> Attaching package: ‘purrr’
#> The following object is masked from ‘package:testthat’:
+#> 
+#>     is_null

+# forecast error variance decomposition
+var.fevd = sovereign::FEVD(var)
+
+# historical shock decomposition
+var.hd = sovereign::HD(var)
+
+# }
+
+
+
+ +
+ + + +
+ + + + + + + + + + + diff --git a/docs/reference/HD.html b/docs/reference/HD.html new file mode 100644 index 0000000..ce9df2b --- /dev/null +++ b/docs/reference/HD.html @@ -0,0 +1,256 @@ + + + + + + + + +Estimate historical decomposition — HD • sovereign + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ + + + +
+ +
+
+
+

Estimate historical decomposition

+ Source: R/analysis_wrappers.R + +
+ +
+

Estimate the historical decomposition for VARs with +either 'short' or 'IV-short' structural errors. See VAR +and RVAR documentation for details regarding structural errors.

+
+ + 
HD(model)
+ +

Arguments

+ + + + + + +
model

VAR or RVAR class object

+ +

Value

+ +

long-from data.frame

+

See also

+ +

VAR()

+

var_hd()

+

RVAR()

+

rvar_hd()

+ +

Examples

+ 
# \donttest{
+
+ # simple time series
+ AA = c(1:100) + rnorm(100)
+ BB = c(1:100) + rnorm(100)
+ CC = AA + BB + rnorm(100)
+ date = seq.Date(from = as.Date('2000-01-01'), by = 'month', length.out = 100)
+ Data = data.frame(date = date, AA, BB, CC)
+
+ # estimate VAR
+  var =
+    sovereign::VAR(
+      data = Data,
+      horizon = 10,
+      freq = 'month',
+      lag.ic = 'BIC',
+      lag.max = 4)
+
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion

+# impulse response functions
+var.irf = sovereign::IRF(var)
+
#> 
+#> Attaching package: ‘purrr’
#> The following object is masked from ‘package:testthat’:
+#> 
+#>     is_null

+# forecast error variance decomposition
+var.fevd = sovereign::FEVD(var)
+
+# historical shock decomposition
+var.hd = sovereign::HD(var)
+
+# }
+
+
+
+ +
+ + + +
+ + + + + + + + + + + diff --git a/docs/reference/LP.html b/docs/reference/LP.html index 539a088..17c0725 100644 --- a/docs/reference/LP.html +++ b/docs/reference/LP.html @@ -88,7 +88,7 @@ sovereign - 1.0.1 + 1.2.0 @@ -211,6 +211,13 @@

Arg

Value

list object with elements data, model, forecasts, residuals; if there is more than one forecast horizon estimated, then model, forecasts, residuals will each be a list where each element corresponds to a single horizon

+

References

+ + +
    +

  1. Jorda, Oscar "Estimation and Inference of Impulse Responses by Local Projections" 2005.

    2. +
+

See also

LP()

diff --git a/docs/reference/RLP.html b/docs/reference/RLP.html index 76f2029..3fedb1e 100644 --- a/docs/reference/RLP.html +++ b/docs/reference/RLP.html @@ -46,7 +46,10 @@ - + @@ -88,7 +91,7 @@ sovereign - 1.0.1 + 1.2.0
@@ -147,7 +150,10 @@

Estimate regime-dependent local projections

-

Estimate regime-dependent local projections

+

Estimate a regime-dependent local projection (i.e. a state-dependent LP), with an exogenous state indicator, of the specification: +$$Y_{t+h} = X_t \beta_{s_t} + \epsilon_t$$ +where t is the time index, and s is a mutually exclusive state of the world observed at time t. When the regime vector is +not supplied by the user, then a two-state regime series is estimated via random forest.

RLP(
@@ -225,7 +231,16 @@ 
Arg
 
     
Value

 
-    
list object with elements data, model, forecasts, residuals; if there is more than one forecast horizon estimated, then model, forecasts, residuals will each be a list where each element corresponds to a single horizon

+    
list of lists, one list per regime, each regime with objects with elements data, model, forecasts, residuals;
+if there is more than one forecast horizon estimated, then model, forecasts, residuals
+will each be a list where each element corresponds to a single horizon

+    
References

+
+    
+
    
+
  1. Jorda, Oscar "Estimation and Inference of Impulse Responses by Local Projections" 2005.
    2. 
+

+
     
See also

 
     
LP()

diff --git a/docs/reference/RVAR.html b/docs/reference/RVAR.html
index d59fa5b..01d25ec 100644
--- a/docs/reference/RVAR.html
+++ b/docs/reference/RVAR.html
@@ -6,7 +6,7 @@
 
 
 
-Estimate regime-dependent VAR — RVAR • sovereign
+Estimate regime-dependent VAR, SVAR, or Proxy-SVAR — RVAR • sovereign
 
 
 
@@ -45,11 +45,11 @@
 
 
 
-
+
 
 
 
@@ -92,7 +92,7 @@
       
       
         sovereign
-        1.1.0
+        1.2.0
       
     

 
@@ -145,16 +145,16 @@
 

   

     

-    
Estimate regime-dependent VAR

+    
Estimate regime-dependent VAR, SVAR, or Proxy-SVAR

     Source: R/var.R
     
     

 
     

     
Estimate a regime-dependent VAR (i.e. a state-dependent VAR), with an exogenous state indicator, of the specification:
-$$Y_t = X \beta_s + \epsilon_t$$
-where t is the time index, Y is the set of outcome vectors, X the design matrix (of p lagged values of Y), and
-s is a mutually exclusive state of the world observed at time t-1. When the regime vector is not supplied by the user, then a two-state
+$$Y_{t+1} = X_t \beta_{s_t} + \epsilon_t$$
+where t is the time index, Y is the set of outcome vectors, X the design matrix (of p lagged values of Y), and
+s is a mutually exclusive state of the world observed at time t. When the regime vector is not supplied by the user, then a two-state
 regime series is estimated via random forest.

     

 
@@ -168,7 +168,10 @@ 
Estimate regime-dependent VAR

   lag.max = NULL,
   regime = NULL,
   regime.method = "rf",
-  regime.n = 2
+  regime.n = 2,
+  structure = "short",
+  instrument = NULL,
+  instrumented = NULL
 )

Arguments

@@ -214,26 +217,69 @@

Arg regime.n

int: number of regimes to estimate (applies to kmeans and EM)

+ + structure +

string: type of structural identification strategy to use in model analysis (NA, 'short', 'IV', or 'IV-short')

+ + + instrument +

string: name of instrumental variable contained in the data matrix

+ + + instrumented +

string: name of variable to be instrumented in IV and IV-short procedure; default is the first non-date variable in data

+

Value

-

list of lists, each regime returns its own list with elements data, model, forecasts, residuals

+

List of lists, where each regime is a list with items:

    +

  1. data: data.frame with endogenous variables and 'date' column.

    2. +

  2. model: list with data.frame of model coefficients (in psuedo-companion form), data.frame of coefficient standard errors, integer of lags p, integer of horizons, string of frequency, string of deterministic term type, numeric of log-likelihood, regime indicator

    3. +

  3. forecasts: list of data.frames per horizon; data.frame with column for date (day the forecast was made), forecast.date (the date being forecasted), target (variable forecasted), and forecast

    4. +

  4. residuals: list of data.frames per horizon; data.frame of residuals

    5. +

  5. structure: string denoting which structural identification strategy will be used in analysis (or NA)

    6. +

  6. instrument: data.frame with 'date' column and 'instrument' column (or NULL)

    7. +

  7. instrumented: string denoting which column will be instrumted in 'IV' and 'IV-short' strategies (or NULL)

    8. +
+

Details

The regime-dependent VAR is a generalization of the popular threshold VAR - which trades off estimating a threshold value for an endogenous variable for accepting an exogenous regime that can be based on information from inside or outside of the system, with or without parametric -assumptions, and with or without timing restrictions.

+assumptions, and with or without timing restrictions. Moreover, the RVAR may be extended to include structural shocks, including the use of instrumental +variables.

+

State dependence. The RVAR augments the traditional VAR by allowing state-dependence in the coefficient matrix. The RVAR differs from the more common threshold VAR (TVAR), due +to the fact that states are exegonesouly determined. As a result, the states (i.e. regimes) do not need to be based on information inside the model, moreover, regimes can be +determined by any process the user determines best fits their needs. For example, regimes based on NBER dated recessions and expansions are based on a judgmental process +considering hundreds of series, potentially none of which are in the VAR being modeled. Alternatively, a user may use unsupervised machine learning to assign regimes - this is +the process the sovereign::regimes function facilitates.

+

Structural shocks. See Sims (1980) for details regarding the baseline vector-autoregression (VAR) model. The VAR may be augmented to become a structural VAR (SVAR) with one of three different structural identification strategies:

    +

  1. short-term impact restrictions via Cholesky decomposition, see Christiano et al (1999) for details (structure = 'short')

    2. +

  2. external instrument identification, i.e. a Proxy-SVAR strategy, see Mertens and Ravn (2013) for details (structure = 'IV')

    3. +

  3. or a combination of short-term and IV identification via Lunsford (2015) (structure = 'IV-short')

    4. +
+ +

Note that including structure does not change the estimation of model coefficients or forecasts, but does change impulse response functions, forecast error variance decomposition, +and historical decompositions. Historical decompositions will not be available for models using the 'IV' structure. Additionally note that only one instrument may be used in this +estimation routine.

+

References

+ + +
    +

  1. Christiano, Lawrence, Martin Eichenbaum, and Charles Evans "Monetary policy shocks: What have we learned and to what end?" Handbook of Macroeconomics, Vol 1, Part A, 1999.

    2. +

  2. Lunsford, Kurt "Identifying Structural VARs with a Proxy Variable and a Test for a Weak Proxy" 2015.

    3. +

  3. Mertens, Karel and Morten Ravn "The Dynamic Effects of Personal and Corporate Income Tax Changes in the United States" 2013.

    4. +

  4. Sims, Christopher "Macroeconomics and Reality" 1980.

    5. +
+

See also

VAR()

-

var_irf()

-

var_fevd()

-

var_hd()

RVAR()

-

rvar_irf()

-

rvar_fevd()

-

rvar_hd()

+

IRF()

+

FEVD()

+

HD()

Examples

# \donttest{
@@ -256,16 +302,16 @@ 
Examp
       regime.n = 2,
       lag.ic = 'BIC',
       lag.max = 4)
-
#> Error in n.lag(., lags = p): could not find function "n.lag"

-# impulse response functions
-rvar.irf = sovereign::rvar_irf(rvar)
-
#> Error in sovereign::rvar_irf(rvar): object 'rvar' not found

-# forecast error variance decomposition
-rvar.fevd = sovereign::rvar_fevd(rvar)
-
#> Error in sovereign::rvar_fevd(rvar): object 'rvar' not found

-# historical shock decomposition
-rvar.hd = sovereign::rvar_hd(rvar)
-
#> Error in sovereign::rvar_hd(rvar): object 'rvar' not found

+
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion

+ # impulse response functions
+ rvar.irf = sovereign::rvar_irf(rvar)
+
+ # forecast error variance decomposition
+ rvar.fevd = sovereign::rvar_fevd(rvar)
+
+ # historical shock decomposition
+ rvar.hd = sovereign::rvar_hd(rvar)
+
 # }
diff --git a/docs/reference/VAR.html b/docs/reference/VAR.html index 1a0b816..50b8ab0 100644 --- a/docs/reference/VAR.html +++ b/docs/reference/VAR.html @@ -88,7 +88,7 @@ sovereign - 1.1.0 + 1.2.0 @@ -196,7 +196,7 @@

Arg structure -

string: type of structural identification strategy to use in model analysis (NULL, 'short', 'IV', or 'IV-short')

+

string: type of structural identification strategy to use in model analysis (NA, 'short', 'IV', or 'IV-short')

instrument @@ -213,12 +213,12 @@

Value

  1. data: data.frame with endogenous variables and 'date' column.

    2. -

  2. model: list with data.frame of model coefficients (in psuedo-companion form), data.frame of coefficient standard errors, integer of lags p, integer of horizons, string of frequency, string of type

    3. +

  3. model: list with data.frame of model coefficients (in psuedo-companion form), data.frame of coefficient standard errors, integer of lags p, integer of horizons, string of frequency, string of deterministic term type, numeric of log-likelihood

  4. forecasts: list of data.frames per horizon; data.frame with column for date (day the forecast was made), forecast.date (the date being forecasted), target (variable forecasted), and forecast

  5. residuals: list of data.frames per horizon; data.frame of residuals

    6. -

  6. structure: string denoting which structural identification strategy will be used in analysis (or NULL)

    7. +

  7. structure: string denoting which structural identification strategy will be used in analysis (or NA)

  8. instrument: data.frame with 'date' column and 'instrument' column (or NULL)

    9. -

  9. instrumented: string denoting which column will be instrumted in 'IV' and 'IV-short' strategies (or NULL)

    10. +

  10. instrumented: string denoting which column will be instrumted in 'IV' and 'IV-short' strategies (or NA)

Details

@@ -266,16 +266,16 @@

Examp freq = 'month', lag.ic = 'BIC', lag.max = 4) -
#> Error in n.lag(., lags = p): could not find function "n.lag"
+
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
# impulse response functions var.irf = sovereign::var_irf(var) -
#> Error in var$model: object of type 'closure' is not subsettable
+ # forecast error variance decomposition var.fevd = sovereign::var_fevd(var) -
#> Error in var$model: object of type 'closure' is not subsettable
+ # historical shock decomposition var.hd = sovereign::var_hd(var) -
#> Error in var$model: object of type 'closure' is not subsettable
# } +# }
diff --git a/docs/reference/index.html b/docs/reference/index.html index 476b74f..7382046 100644 --- a/docs/reference/index.html +++ b/docs/reference/index.html @@ -87,7 +87,7 @@ sovereign - 1.1.0 + 1.2.0 @@ -193,7 +193,7 @@

RVAR()

-

Estimate regime-dependent VAR

+

Estimate regime-dependent VAR, SVAR, or Proxy-SVAR

@@ -228,51 +228,21 @@

var_irf()

+

IRF()

Estimate impulse response functions

-

rvar_irf()

- -

Estimate regime-dependent impulse response functions

- - - -

var_fevd()

+

FEVD()

Estimate forecast error variance decomposition

-

rvar_fevd()

- -

Estimate regime-dependent forecast error variance decomposition

- - - -

var_hd()

+

HD()

Estimate historical decomposition

- - - -

rvar_hd()

- -

Estimate regime-dependent historical decomposition

- - - -

lp_irf()

- -

Estimate impulse response functions

- - - -

rlp_irf()

- -

Estimate regime-dependent impulse response functions

diff --git a/docs/reference/irf.html b/docs/reference/irf.html new file mode 100644 index 0000000..713184d --- /dev/null +++ b/docs/reference/irf.html @@ -0,0 +1,291 @@ + + + + + + + + +Estimate impulse response functions — IRF • sovereign + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ + + + +
+ +
+
+
+

Estimate impulse response functions

+ Source: R/analysis_wrappers.R + +
+ +
+

See VAR, RVAR, LP, and RLP documentation for details +regarding models and structural errors.

+
+ + 
IRF(
+  model,
+  horizon = 10,
+  CI = c(0.1, 0.9),
+  bootstrap.type = "auto",
+  bootstrap.num = 100,
+  bootstrap.parallel = FALSE,
+  bootstrap.cores = -1
+)
+ +

Arguments

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
model

VAR, RVAR, LP, or RLP class object

horizon

int: number of periods

CI

numeric vector: c(lower ci bound, upper ci bound)

bootstrap.type

string: bootstrapping technique to use ('auto', 'standard', or 'wild'); if auto then wild is used for IV or IV-short, else standard is used

bootstrap.num

int: number of bootstraps

bootstrap.parallel

boolean: create IRF draws in parallel

bootstrap.cores

int: number of cores to use in parallel processing; -1 detects and uses half the available cores

+ +

Value

+ +

data frame with columns target, shock, horizon, response.lower, response, response.upper; regime-based models return a list with a data frame per regime.

+

See also

+ +

var_irf()

+

rvar_irf()

+

lp_irf()

+

rlp_irf()

+ +

Examples

+ 
# \donttest{
+
+ # simple time series
+ AA = c(1:100) + rnorm(100)
+ BB = c(1:100) + rnorm(100)
+ CC = AA + BB + rnorm(100)
+ date = seq.Date(from = as.Date('2000-01-01'), by = 'month', length.out = 100)
+ Data = data.frame(date = date, AA, BB, CC)
+
+ # estimate VAR
+  var =
+    sovereign::VAR(
+      data = Data,
+      horizon = 10,
+      freq = 'month',
+      lag.ic = 'BIC',
+      lag.max = 4
+    )
+
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion

+ # impulse response function
+ var.irf = sovereign::IRF(var)
+
+  # local projection forecasts
+  lp =
+    sovereign::LP(
+      data = Data,
+      horizon = c(1:10),
+      lag.ic = 'AIC',
+      lag.max = 4,
+      type =  'both',
+      freq = 'month')
+
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion

+  # LP impulse response function
+  lp.irf = sovereign::IRF(lp)
+
+# }
+
+
+
+ +
+ + +
+
+

Developed by Tyler J. Pike.

+
+ +
+

Site built with pkgdown 1.6.1.

+
+ +
+
+ + + + + + + + + + + diff --git a/docs/reference/lp_irf.html b/docs/reference/lp_irf.html index 68a784d..e60f512 100644 --- a/docs/reference/lp_irf.html +++ b/docs/reference/lp_irf.html @@ -88,7 +88,7 @@ sovereign - 1.0.1 + 1.2.0 diff --git a/docs/reference/rvar_fevd.html b/docs/reference/rvar_fevd.html index bd8a295..9b3a822 100644 --- a/docs/reference/rvar_fevd.html +++ b/docs/reference/rvar_fevd.html @@ -46,7 +46,8 @@ - + @@ -88,7 +89,7 @@ sovereign - 1.1.0 + 1.2.0 @@ -147,7 +148,8 @@

Estimate regime-dependent forecast error variance decomposition

-

Estimate regime-dependent forecast error variance decomposition

+

Estimate forecast error variance decomposition for RVARs +with either short or 'IV-short' structural errors.

rvar_fevd(rvar, horizon = 10, scale = TRUE)
@@ -204,16 +206,16 @@

Examp regime.n = 2, lag.ic = 'BIC', lag.max = 4) -
#> Error in n.lag(., lags = p): could not find function "n.lag"
+
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
# impulse response functions rvar.irf = sovereign::rvar_irf(rvar) -
#> Error in sovereign::rvar_irf(rvar): object 'rvar' not found
+ # forecast error variance decomposition rvar.fevd = sovereign::rvar_fevd(rvar) -
#> Error in sovereign::rvar_fevd(rvar): object 'rvar' not found
+ # historical shock decomposition rvar.hd = sovereign::rvar_hd(rvar) -
#> Error in sovereign::rvar_hd(rvar): object 'rvar' not found
+ # }
diff --git a/docs/reference/rvar_hd.html b/docs/reference/rvar_hd.html index 7095685..f6d1b5c 100644 --- a/docs/reference/rvar_hd.html +++ b/docs/reference/rvar_hd.html @@ -46,10 +46,8 @@ - + @@ -91,7 +89,7 @@ sovereign - 1.1.0 + 1.2.0 @@ -150,10 +148,8 @@

Estimate regime-dependent historical decomposition

-

Estimate historical decomposition with contemporaneous impact -restrictions via Cholesky decomposition - taking the variable ordering -present in the RVAR object. Estimate one response function per unique -state defined by the regime-dependent VAR.

+

Estimate historical decomposition for RVARs with +either short or 'IV-short' structural errors.

rvar_hd(rvar)
@@ -202,16 +198,16 @@

Examp regime.n = 2, lag.ic = 'BIC', lag.max = 4) -
#> Error in n.lag(., lags = p): could not find function "n.lag"
+
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
# impulse response functions rvar.irf = sovereign::rvar_irf(rvar) -
#> Error in sovereign::rvar_irf(rvar): object 'rvar' not found
+ # forecast error variance decomposition rvar.fevd = sovereign::rvar_fevd(rvar) -
#> Error in sovereign::rvar_fevd(rvar): object 'rvar' not found
+ # historical shock decomposition rvar.hd = sovereign::rvar_hd(rvar) -
#> Error in sovereign::rvar_hd(rvar): object 'rvar' not found
+ # }
diff --git a/docs/reference/rvar_irf.html b/docs/reference/rvar_irf.html index 73442ac..fc12e16 100644 --- a/docs/reference/rvar_irf.html +++ b/docs/reference/rvar_irf.html @@ -46,9 +46,7 @@ - + @@ -90,7 +88,7 @@ sovereign - 1.1.0 + 1.2.0 @@ -149,12 +147,18 @@

Estimate regime-dependent impulse response functions

-

Estimate impulse responses with contemporaneous impact restrictions via Cholesky decomposition - -taking the variable ordering present in the VAR object. Estimate one response function per unique -state defined by the regime-dependent VAR.

+

Estimate regime-dependent impulse response functions

- 
rvar_irf(rvar, horizon = 10, bootstraps.num = 100, CI = c(0.1, 0.9))
+ 
rvar_irf(
+  rvar,
+  horizon = 10,
+  CI = c(0.1, 0.9),
+  bootstrap.type = "auto",
+  bootstrap.num = 100,
+  bootstrap.parallel = FALSE,
+  bootstrap.cores = -1
+)

Arguments

@@ -168,18 +172,30 @@

Arg

- + + + + + + + + + - - + + + + + +

int: number of periods
bootstraps.numCI

numeric vector: c(lower ci bound, upper ci bound)

bootstrap.type

string: bootstrapping technique to use ('auto', 'standard', or 'wild'); if auto then wild is used for IV or IV-short, else standard is used

bootstrap.num

int: number of bootstraps
CI

numeric vector: c(lower ci bound, upper ci bound)

bootstrap.parallel

boolean: create IRF draws in parallel

bootstrap.cores

int: number of cores to use in parallel processing; -1 detects and uses half the available cores

Value

-

list of lists, each regime returns its own list with elements irfs, ci.lower, and ci.upper; all elements are long-form data.frames

+

list of regimes, each with data.frame of columns target, shock, horizon, response.lower, response, response.upper

See also

VAR()

@@ -210,16 +226,16 @@

Examp regime.n = 2, lag.ic = 'BIC', lag.max = 4) -

#> Error in n.lag(., lags = p): could not find function "n.lag"
-# impulse response functions -rvar.irf = sovereign::rvar_irf(rvar) -
#> Error in sovereign::rvar_irf(rvar): object 'rvar' not found
-# forecast error variance decomposition -rvar.fevd = sovereign::rvar_fevd(rvar) -
#> Error in sovereign::rvar_fevd(rvar): object 'rvar' not found
-# historical shock decomposition -rvar.hd = sovereign::rvar_hd(rvar) -
#> Error in sovereign::rvar_hd(rvar): object 'rvar' not found
+
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
+ # impulse response functions + rvar.irf = sovereign::rvar_irf(rvar) + + # forecast error variance decomposition + rvar.fevd = sovereign::rvar_fevd(rvar) + + # historical shock decomposition + rvar.hd = sovereign::rvar_hd(rvar) + # }
diff --git a/docs/reference/var_fevd.html b/docs/reference/var_fevd.html index cad8f58..2e5e3b6 100644 --- a/docs/reference/var_fevd.html +++ b/docs/reference/var_fevd.html @@ -46,7 +46,8 @@ - + @@ -88,7 +89,7 @@ sovereign - 1.1.0 + 1.2.0 @@ -147,7 +148,8 @@

Estimate forecast error variance decomposition

-

Estimate forecast error variance decomposition

+

Estimate forecast error variance decomposition for VARs +with either short or 'IV-short' structural errors.

var_fevd(var, horizon = 10, scale = TRUE)
@@ -201,16 +203,16 @@

Examp freq = 'month', lag.ic = 'BIC', lag.max = 4) -
#> Error in n.lag(., lags = p): could not find function "n.lag"
+
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
# impulse response functions var.irf = sovereign::var_irf(var) -
#> Error: object of type 'closure' is not subsettable
+ # forecast error variance decomposition var.fevd = sovereign::var_fevd(var) -
#> Error: object of type 'closure' is not subsettable
+ # historical shock decomposition var.hd = sovereign::var_hd(var) -
#> Error: object of type 'closure' is not subsettable
+ # }
diff --git a/docs/reference/var_hd.html b/docs/reference/var_hd.html index e2dd842..469bc47 100644 --- a/docs/reference/var_hd.html +++ b/docs/reference/var_hd.html @@ -46,9 +46,8 @@ - + @@ -90,7 +89,7 @@ sovereign - 1.1.0 + 1.2.0 @@ -149,9 +148,8 @@

Estimate historical decomposition

-

Estimate historical decomposition with contemporaneous impact -restrictions via Cholesky decomposition - taking the variable ordering -present in the VAR object.

+

Estimate historical decomposition for VARs with +either short or 'IV-short' structural errors.

var_hd(var)
@@ -197,16 +195,16 @@

Examp freq = 'month', lag.ic = 'BIC', lag.max = 4) -
#> Error in n.lag(., lags = p): could not find function "n.lag"
+
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
# impulse response functions var.irf = sovereign::var_irf(var) -
#> Error: object of type 'closure' is not subsettable
+ # forecast error variance decomposition var.fevd = sovereign::var_fevd(var) -
#> Error: object of type 'closure' is not subsettable
+ # historical shock decomposition var.hd = sovereign::var_hd(var) -
#> Error: object of type 'closure' is not subsettable
+ # }
diff --git a/docs/reference/var_irf.html b/docs/reference/var_irf.html index c836f2d..bb01fad 100644 --- a/docs/reference/var_irf.html +++ b/docs/reference/var_irf.html @@ -46,8 +46,7 @@ - + @@ -89,7 +88,7 @@ sovereign - 1.1.0 + 1.2.0 @@ -148,8 +147,7 @@

Estimate impulse response functions

-

Estimate impulse responses with contemporaneous impact restrictions via Cholesky decomposition - -taking the variable ordering present in the VAR object.

+

Estimate impulse response functions

var_irf(
@@ -178,14 +176,26 @@ 
Arg
       
numeric vector: c(lower ci bound, upper ci bound)

     
     
-      bootstraps.num
+      bootstrap.type
+      
string: bootstrapping technique to use ('auto', 'standard', or 'wild'); if auto then wild is used for IV or IV-short, else standard is used

+    
+    
+      bootstrap.num
       
int: number of bootstraps

     
+    
+      bootstrap.parallel
+      
boolean: create IRF draws in parallel

+    
+    
+      bootstrap.cores
+      
int: number of cores to use in parallel processing; -1 detects and uses half the available cores

+    
     
 
     
Value

 
-    
list object with elements irfs, ci.lower, and ci.upper; all elements are long-form data.frames

+    
data.frame with columns target, shock, horizon, response.lower, response, response.upper

     
See also

 
     
VAR()

@@ -215,16 +225,16 @@ 
Examp
       freq = 'month',
       lag.ic = 'BIC',
       lag.max = 4)
-
#> Error in n.lag(., lags = p): could not find function "n.lag"

-# impulse response functions
-var.irf = sovereign::var_irf(var)
-
#> Error: object of type 'closure' is not subsettable

-# forecast error variance decomposition
-var.fevd = sovereign::var_fevd(var)
-
#> Error: object of type 'closure' is not subsettable

-# historical shock decomposition
-var.hd = sovereign::var_hd(var)
-
#> Error: object of type 'closure' is not subsettable

+
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion
#> Warning: NAs introduced by coercion

+ # impulse response functions
+ var.irf = sovereign::var_irf(var)
+
+ # forecast error variance decomposition
+ var.fevd = sovereign::var_fevd(var)
+
+ # historical shock decomposition
+ var.hd = sovereign::var_hd(var)
+
 # }
diff --git a/docs/to-do.html b/docs/to-do.html index 5627cf6..c1a3383 100644 --- a/docs/to-do.html +++ b/docs/to-do.html @@ -87,7 +87,7 @@ sovereign - 1.1.0 + 1.2.0 @@ -145,100 +145,50 @@

To-do

-
+

-High Priority (Version 1.0.1)

+High Priority (Version 1.2.0)
    -
  1. code -
      -
    1. -historical decomposition +
    2. Code
        -
      1. code
        2. -
      2. clean code
        3. -
      3. test
        4. +
      4. threshold VAR
      5. -documentation
        +more structural options
        6. -
      -
      3. -
    3. -LP volatility correction +
    4. Forecast performance
        -
      1. code
        2. -
      2. clean and add warnings
        3. -
      3. test
        4. -
      4. -documentation
        -
        5. -
      -
      5. -
    5. Add classes (e.g. class(var))
      6. -
    6. -Make n.lag not require a date
      +
    7. add BIC and AIC
      8. +
    8. add RMSE, MSE, MAE
    2. -
  2. documentation -
      -
    1. refine discussion of shocks
      2. -
    2. -define equations
      -
    3. -
  3. website +
  4. Documentation
      -
    1. -update home page
      -
      2. -
    2. -update references
      -
      3. -
    3. -update vignette
      +
    4. check that new VAR values section is correct
      5. +
    5. document proxy VAR code and methodology
      6. +
    6. link to helpful external resources (e.g. Bianchi’s slides)
      7. +
    7. add academic references to methods
      8. -
    8. update news
    5. -
-
-
-

-Medium Priority (Version 1.0.2)

+
  • Website
      -
    1. Code -
        -
      1. endeogenously estimated threshold VAR
        2. -
      2. more structural options -
          -
        1. no structure
          2. -
        2. short term (cholesky)
          3. -
        3. long term
          4. -
        4. instrumental variable (proxy var)
          5. -
        +
      3. update home page
        4. -
      4. Forecast performance -
          -
        1. add BIC and AIC
          2. -
        2. add RMSE, MSE, MAE
          +
        3. update references
          4. -
        -
        5. -
      +
    2. update vignette
      3. -
    3. Documentation -
        -
      1. regimes function methods
        2. -
      2. link to helpful external resources (e.g. Bianchi’s slides)
        3. -
      3. add academic references to methods
        4. +
      4. update news
    • -
    +

    -Low Priority (> Version 1.0.2)

    +Low Priority (> Version 1.2.0)
    1. Code
        @@ -246,6 +196,7 @@

      1. State-dependent ML
      2. Bayesian Treshold VAR
        3. +
      3. Data cleaning (e.g. hamilton filter)

    2. Documentation @@ -255,6 +206,17 @@

    +
    +

    +FOMC data

    +

    write a set of functions, perhaps its own package, with nicely cleaned FOMC data (set the functions to pull from a google drive, github, or some public location)

    +
      +
    • FOMC meeting and announcement dates
      • +
    • policy rates
      • +
    • QE actions
      • +
    • and derived MP shocks
      • +
    +
    diff --git a/man/FEVD.Rd b/man/FEVD.Rd new file mode 100644 index 0000000..a7e1d07 --- /dev/null +++ b/man/FEVD.Rd @@ -0,0 +1,63 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/analysis_wrappers.R +\name{FEVD} +\alias{FEVD} +\title{Estimate forecast error variance decomposition} +\usage{ +FEVD(model, horizon = 10, scale = TRUE) +} +\arguments{ +\item{model}{VAR or RVAR class object} + +\item{horizon}{int: number of periods} + +\item{scale}{boolean: scale variable contribution as percent of total error} +} +\value{ +long-form data.frame +} +\description{ +Estimate the forecast error variance decomposition for VARs with +either short or 'IV-short' structural errors. See VAR +and RVAR documentation for details regarding structural errors. +} +\examples{ +\donttest{ + + # simple time series + AA = c(1:100) + rnorm(100) + BB = c(1:100) + rnorm(100) + CC = AA + BB + rnorm(100) + date = seq.Date(from = as.Date('2000-01-01'), by = 'month', length.out = 100) + Data = data.frame(date = date, AA, BB, CC) + + # estimate VAR + var = + sovereign::VAR( + data = Data, + horizon = 10, + freq = 'month', + lag.ic = 'BIC', + lag.max = 4) + +# impulse response functions +var.irf = sovereign::IRF(var) + +# forecast error variance decomposition +var.fevd = sovereign::FEVD(var) + +# historical shock decomposition +var.hd = sovereign::HD(var) + +} + +} +\seealso{ +\code{\link[=VAR]{VAR()}} + +\code{\link[=var_fevd]{var_fevd()}} + +\code{\link[=RVAR]{RVAR()}} + +\code{\link[=rvar_fevd]{rvar_fevd()}} +} diff --git a/man/HD.Rd b/man/HD.Rd new file mode 100644 index 0000000..d4d6e76 --- /dev/null +++ b/man/HD.Rd @@ -0,0 +1,59 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/analysis_wrappers.R +\name{HD} +\alias{HD} +\title{Estimate historical decomposition} +\usage{ +HD(model) +} +\arguments{ +\item{model}{VAR or RVAR class object} +} +\value{ +long-from data.frame +} +\description{ +Estimate the historical decomposition for VARs with +either 'short' or 'IV-short' structural errors. See VAR +and RVAR documentation for details regarding structural errors. +} +\examples{ +\donttest{ + + # simple time series + AA = c(1:100) + rnorm(100) + BB = c(1:100) + rnorm(100) + CC = AA + BB + rnorm(100) + date = seq.Date(from = as.Date('2000-01-01'), by = 'month', length.out = 100) + Data = data.frame(date = date, AA, BB, CC) + + # estimate VAR + var = + sovereign::VAR( + data = Data, + horizon = 10, + freq = 'month', + lag.ic = 'BIC', + lag.max = 4) + +# impulse response functions +var.irf = sovereign::IRF(var) + +# forecast error variance decomposition +var.fevd = sovereign::FEVD(var) + +# historical shock decomposition +var.hd = sovereign::HD(var) + +} + +} +\seealso{ +\code{\link[=VAR]{VAR()}} + +\code{\link[=var_hd]{var_hd()}} + +\code{\link[=RVAR]{RVAR()}} + +\code{\link[=rvar_hd]{rvar_hd()}} +} diff --git a/man/IRF.Rd b/man/IRF.Rd new file mode 100644 index 0000000..3a5e8f4 --- /dev/null +++ b/man/IRF.Rd @@ -0,0 +1,86 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/analysis_wrappers.R +\name{IRF} +\alias{IRF} +\title{Estimate impulse response functions} +\usage{ +IRF( + model, + horizon = 10, + CI = c(0.1, 0.9), + bootstrap.type = "auto", + bootstrap.num = 100, + bootstrap.parallel = FALSE, + bootstrap.cores = -1 +) +} +\arguments{ +\item{model}{VAR, RVAR, LP, or RLP class object} + +\item{horizon}{int: number of periods} + +\item{CI}{numeric vector: c(lower ci bound, upper ci bound)} + +\item{bootstrap.type}{string: bootstrapping technique to use ('auto', 'standard', or 'wild'); if auto then wild is used for IV or IV-short, else standard is used} + +\item{bootstrap.num}{int: number of bootstraps} + +\item{bootstrap.parallel}{boolean: create IRF draws in parallel} + +\item{bootstrap.cores}{int: number of cores to use in parallel processing; -1 detects and uses half the available cores} +} +\value{ +data frame with columns \code{target}, \code{shock}, \code{horizon}, \code{response.lower}, \code{response}, \code{response.upper}; regime-based models return a list with a data frame per regime. +} +\description{ +See VAR, RVAR, LP, and RLP documentation for details +regarding models and structural errors. +} +\examples{ +\donttest{ + + # simple time series + AA = c(1:100) + rnorm(100) + BB = c(1:100) + rnorm(100) + CC = AA + BB + rnorm(100) + date = seq.Date(from = as.Date('2000-01-01'), by = 'month', length.out = 100) + Data = data.frame(date = date, AA, BB, CC) + + # estimate VAR + var = + sovereign::VAR( + data = Data, + horizon = 10, + freq = 'month', + lag.ic = 'BIC', + lag.max = 4 + ) + + # impulse response function + var.irf = sovereign::IRF(var) + + # local projection forecasts + lp = + sovereign::LP( + data = Data, + horizon = c(1:10), + lag.ic = 'AIC', + lag.max = 4, + type = 'both', + freq = 'month') + + # LP impulse response function + lp.irf = sovereign::IRF(lp) + +} + +} +\seealso{ +\code{\link[=var_irf]{var_irf()}} + +\code{\link[=rvar_irf]{rvar_irf()}} + +\code{\link[=lp_irf]{lp_irf()}} + +\code{\link[=rlp_irf]{rlp_irf()}} +} diff --git a/man/LP.Rd b/man/LP.Rd index 1282125..4245675 100644 --- a/man/LP.Rd +++ b/man/LP.Rd @@ -69,6 +69,11 @@ Estimate local projections } +} +\references{ +\enumerate{ +\item Jorda, Oscar "\href{https://www.aeaweb.org/articles?id=10.1257/0002828053828518}{Estimation and Inference of Impulse Responses by Local Projections}" 2005. +} } \seealso{ \code{\link[=LP]{LP()}} diff --git a/man/RLP.Rd b/man/RLP.Rd index 2fbeac3..e5be863 100644 --- a/man/RLP.Rd +++ b/man/RLP.Rd @@ -48,10 +48,15 @@ RLP( \item{regime.n}{int: number of regimes to estimate (applies to kmeans and EM)} } \value{ -list object with elements \code{data}, \code{model}, \code{forecasts}, \code{residuals}; if there is more than one forecast horizon estimated, then \code{model}, \code{forecasts}, \code{residuals} will each be a list where each element corresponds to a single horizon +list of lists, one list per regime, each regime with objects with elements \code{data}, \code{model}, \code{forecasts}, \code{residuals}; +if there is more than one forecast horizon estimated, then \code{model}, \code{forecasts}, \code{residuals} +will each be a list where each element corresponds to a single horizon } \description{ -Estimate regime-dependent local projections +Estimate a regime-dependent local projection (i.e. a state-dependent LP), with an exogenous state indicator, of the specification: +\deqn{Y_{t+h} = X_t \beta_{s_t} + \epsilon_t} +where \emph{t} is the time index, and \emph{s} is a mutually exclusive state of the world observed at time \emph{t}. When the regime vector is +not supplied by the user, then a two-state regime series is estimated via random forest. } \examples{ \donttest{ @@ -83,6 +88,11 @@ Estimate regime-dependent local projections } +} +\references{ +\enumerate{ +\item Jorda, Oscar "\href{https://www.aeaweb.org/articles?id=10.1257/0002828053828518}{Estimation and Inference of Impulse Responses by Local Projections}" 2005. +} } \seealso{ \code{\link[=LP]{LP()}} diff --git a/man/RVAR.Rd b/man/RVAR.Rd index 8f30977..f93475c 100644 --- a/man/RVAR.Rd +++ b/man/RVAR.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/var.R \name{RVAR} \alias{RVAR} -\title{Estimate regime-dependent VAR} +\title{Estimate regime-dependent VAR, SVAR, or Proxy-SVAR} \usage{ RVAR( data, @@ -14,7 +14,10 @@ RVAR( lag.max = NULL, regime = NULL, regime.method = "rf", - regime.n = 2 + regime.n = 2, + structure = "short", + instrument = NULL, + instrumented = NULL ) } \arguments{ @@ -37,21 +40,54 @@ RVAR( \item{regime.method}{string: regime assignment technique ('rf', 'kmeans', 'EM', or 'BP')} \item{regime.n}{int: number of regimes to estimate (applies to kmeans and EM)} + +\item{structure}{string: type of structural identification strategy to use in model analysis (NA, 'short', 'IV', or 'IV-short')} + +\item{instrument}{string: name of instrumental variable contained in the data matrix} + +\item{instrumented}{string: name of variable to be instrumented in IV and IV-short procedure; default is the first non-date variable in data} } \value{ -list of lists, each regime returns its own list with elements \code{data}, \code{model}, \code{forecasts}, \code{residuals} +List of lists, where each regime is a list with items: +\enumerate{ +\item data: data.frame with endogenous variables and 'date' column. +\item model: list with data.frame of model coefficients (in psuedo-companion form), data.frame of coefficient standard errors, integer of lags p, integer of horizons, string of frequency, string of deterministic term type, numeric of log-likelihood, regime indicator +\item forecasts: list of data.frames per horizon; data.frame with column for date (day the forecast was made), forecast.date (the date being forecasted), target (variable forecasted), and forecast +\item residuals: list of data.frames per horizon; data.frame of residuals +\item structure: string denoting which structural identification strategy will be used in analysis (or NA) +\item instrument: data.frame with 'date' column and 'instrument' column (or NULL) +\item instrumented: string denoting which column will be instrumted in 'IV' and 'IV-short' strategies (or NULL) +} } \description{ Estimate a regime-dependent VAR (i.e. a state-dependent VAR), with an exogenous state indicator, of the specification: -\deqn{Y_t = X \beta_s + \epsilon_t} -where t is the time index, Y is the set of outcome vectors, X the design matrix (of p lagged values of Y), and -s is a mutually exclusive state of the world observed at time t-1. When the regime vector is not supplied by the user, then a two-state +\deqn{Y_{t+1} = X_t \beta_{s_t} + \epsilon_t} +where \emph{t} is the time index, \emph{Y} is the set of outcome vectors, \emph{X} the design matrix (of \emph{p} lagged values of Y), and +\emph{s} is a mutually exclusive state of the world observed at time \emph{t}. When the regime vector is not supplied by the user, then a two-state regime series is estimated via random forest. } \details{ The regime-dependent VAR is a generalization of the popular threshold VAR - which trades off estimating a threshold value for an endogenous variable for accepting an exogenous regime that can be based on information from inside or outside of the system, with or without parametric -assumptions, and with or without timing restrictions. +assumptions, and with or without timing restrictions. Moreover, the RVAR may be extended to include structural shocks, including the use of instrumental +variables. + +\strong{State dependence.} The RVAR augments the traditional VAR by allowing state-dependence in the coefficient matrix. The RVAR differs from the more common threshold VAR (TVAR), due +to the fact that states are exegonesouly determined. As a result, the states (i.e. regimes) do not need to be based on information inside the model, moreover, regimes can be +determined by any process the user determines best fits their needs. For example, regimes based on NBER dated recessions and expansions are based on a judgmental process +considering hundreds of series, potentially none of which are in the VAR being modeled. Alternatively, a user may use unsupervised machine learning to assign regimes - this is +the process the \code{sovereign::regimes} function facilitates. + +\strong{Structural shocks.} See Sims (1980) for details regarding the baseline vector-autoregression (VAR) model. The VAR may be augmented to become a structural VAR (SVAR) with one of three different structural identification strategies: +\enumerate{ +\item short-term impact restrictions via Cholesky decomposition, see Christiano et al (1999) for details \strong{(structure = 'short')} +\item external instrument identification, i.e. a Proxy-SVAR strategy, see Mertens and Ravn (2013) for details \strong{(structure = 'IV')} +\item or a combination of short-term and IV identification via Lunsford (2015) \strong{(structure = 'IV-short')} +} + +Note that including structure does not change the estimation of model coefficients or forecasts, but does change impulse response functions, forecast error variance decomposition, +and historical decompositions. Historical decompositions will not be available for models using the 'IV' structure. Additionally note that only one instrument may be used in this +estimation routine. } \examples{ \donttest{ @@ -75,32 +111,34 @@ assumptions, and with or without timing restrictions. lag.ic = 'BIC', lag.max = 4) -# impulse response functions -rvar.irf = sovereign::rvar_irf(rvar) + # impulse response functions + rvar.irf = sovereign::rvar_irf(rvar) -# forecast error variance decomposition -rvar.fevd = sovereign::rvar_fevd(rvar) + # forecast error variance decomposition + rvar.fevd = sovereign::rvar_fevd(rvar) -# historical shock decomposition -rvar.hd = sovereign::rvar_hd(rvar) + # historical shock decomposition + rvar.hd = sovereign::rvar_hd(rvar) } +} +\references{ +\enumerate{ +\item Christiano, Lawrence, Martin Eichenbaum, and Charles Evans "\href{https://www.sciencedirect.com/science/article/pii/S1574004899010058}{Monetary policy shocks: What have we learned and to what end?}" Handbook of Macroeconomics, Vol 1, Part A, 1999. +\item Lunsford, Kurt "\href{https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2699452#}{Identifying Structural VARs with a Proxy Variable and a Test for a Weak Proxy}" 2015. +\item Mertens, Karel and Morten Ravn "\href{https://www.aeaweb.org/articles?id=10.1257/aer.103.4.1212}{The Dynamic Effects of Personal and Corporate Income Tax Changes in the United States}" 2013. +\item Sims, Christopher "\href{https://www.jstor.org/stable/1912017}{Macroeconomics and Reality}" 1980. +} } \seealso{ \code{\link[=VAR]{VAR()}} -\code{\link[=var_irf]{var_irf()}} - -\code{\link[=var_fevd]{var_fevd()}} - -\code{\link[=var_hd]{var_hd()}} - \code{\link[=RVAR]{RVAR()}} -\code{\link[=rvar_irf]{rvar_irf()}} +\code{\link[=IRF]{IRF()}} -\code{\link[=rvar_fevd]{rvar_fevd()}} +\code{\link[=FEVD]{FEVD()}} -\code{\link[=rvar_hd]{rvar_hd()}} +\code{\link[=HD]{HD()}} } diff --git a/man/VAR.Rd b/man/VAR.Rd index b79f225..d4fd756 100644 --- a/man/VAR.Rd +++ b/man/VAR.Rd @@ -32,7 +32,7 @@ VAR( \item{lag.max}{int: maximum number of lags to test in lag selection} -\item{structure}{string: type of structural identification strategy to use in model analysis (NULL, 'short', 'IV', or 'IV-short')} +\item{structure}{string: type of structural identification strategy to use in model analysis (NA, 'short', 'IV', or 'IV-short')} \item{instrument}{string: name of instrumental variable contained in the data matrix} @@ -41,12 +41,12 @@ VAR( \value{ \enumerate{ \item data: data.frame with endogenous variables and 'date' column. -\item model: list with data.frame of model coefficients (in psuedo-companion form), data.frame of coefficient standard errors, integer of lags p, integer of horizons, string of frequency, string of type +\item model: list with data.frame of model coefficients (in psuedo-companion form), data.frame of coefficient standard errors, integer of lags p, integer of horizons, string of frequency, string of deterministic term type, numeric of log-likelihood \item forecasts: list of data.frames per horizon; data.frame with column for date (day the forecast was made), forecast.date (the date being forecasted), target (variable forecasted), and forecast \item residuals: list of data.frames per horizon; data.frame of residuals -\item structure: string denoting which structural identification strategy will be used in analysis (or NULL) +\item structure: string denoting which structural identification strategy will be used in analysis (or NA) \item instrument: data.frame with 'date' column and 'instrument' column (or NULL) -\item instrumented: string denoting which column will be instrumted in 'IV' and 'IV-short' strategies (or NULL) +\item instrumented: string denoting which column will be instrumted in 'IV' and 'IV-short' strategies (or NA) } } \description{ diff --git a/man/rvar_fevd.Rd b/man/rvar_fevd.Rd index 8f2143d..23113a7 100644 --- a/man/rvar_fevd.Rd +++ b/man/rvar_fevd.Rd @@ -17,7 +17,8 @@ rvar_fevd(rvar, horizon = 10, scale = TRUE) list, each regime returns its own long-form data.frame } \description{ -Estimate regime-dependent forecast error variance decomposition +Estimate forecast error variance decomposition for RVARs +with either short or 'IV-short' structural errors. } \examples{ \donttest{ diff --git a/man/rvar_hd.Rd b/man/rvar_hd.Rd index 1f12fec..1e06fa3 100644 --- a/man/rvar_hd.Rd +++ b/man/rvar_hd.Rd @@ -13,10 +13,8 @@ rvar_hd(rvar) long form data.frames } \description{ -Estimate historical decomposition with contemporaneous impact -restrictions via Cholesky decomposition - taking the variable ordering -present in the RVAR object. Estimate one response function per unique -state defined by the regime-dependent VAR. +Estimate historical decomposition for RVARs with +either short or 'IV-short' structural errors. } \examples{ \donttest{ diff --git a/man/rvar_irf.Rd b/man/rvar_irf.Rd index fae85c2..b0fcc60 100644 --- a/man/rvar_irf.Rd +++ b/man/rvar_irf.Rd @@ -4,24 +4,36 @@ \alias{rvar_irf} \title{Estimate regime-dependent impulse response functions} \usage{ -rvar_irf(rvar, horizon = 10, bootstraps.num = 100, CI = c(0.1, 0.9)) +rvar_irf( + rvar, + horizon = 10, + CI = c(0.1, 0.9), + bootstrap.type = "auto", + bootstrap.num = 100, + bootstrap.parallel = FALSE, + bootstrap.cores = -1 +) } \arguments{ \item{rvar}{RVAR output} \item{horizon}{int: number of periods} -\item{bootstraps.num}{int: number of bootstraps} - \item{CI}{numeric vector: c(lower ci bound, upper ci bound)} + +\item{bootstrap.type}{string: bootstrapping technique to use ('auto', 'standard', or 'wild'); if auto then wild is used for IV or IV-short, else standard is used} + +\item{bootstrap.num}{int: number of bootstraps} + +\item{bootstrap.parallel}{boolean: create IRF draws in parallel} + +\item{bootstrap.cores}{int: number of cores to use in parallel processing; -1 detects and uses half the available cores} } \value{ -list of lists, each regime returns its own list with elements \code{irfs}, \code{ci.lower}, and \code{ci.upper}; all elements are long-form data.frames +list of regimes, each with data.frame of columns \code{target}, \code{shock}, \code{horizon}, \code{response.lower}, \code{response}, \code{response.upper} } \description{ -Estimate impulse responses with contemporaneous impact restrictions via Cholesky decomposition - -taking the variable ordering present in the VAR object. Estimate one response function per unique -state defined by the regime-dependent VAR. +Estimate regime-dependent impulse response functions } \examples{ \donttest{ @@ -45,14 +57,14 @@ state defined by the regime-dependent VAR. lag.ic = 'BIC', lag.max = 4) -# impulse response functions -rvar.irf = sovereign::rvar_irf(rvar) + # impulse response functions + rvar.irf = sovereign::rvar_irf(rvar) -# forecast error variance decomposition -rvar.fevd = sovereign::rvar_fevd(rvar) + # forecast error variance decomposition + rvar.fevd = sovereign::rvar_fevd(rvar) -# historical shock decomposition -rvar.hd = sovereign::rvar_hd(rvar) + # historical shock decomposition + rvar.hd = sovereign::rvar_hd(rvar) } diff --git a/man/var_fevd.Rd b/man/var_fevd.Rd index 4d1b31d..be29b07 100644 --- a/man/var_fevd.Rd +++ b/man/var_fevd.Rd @@ -17,7 +17,8 @@ var_fevd(var, horizon = 10, scale = TRUE) long-form data.frame } \description{ -Estimate forecast error variance decomposition +Estimate forecast error variance decomposition for VARs +with either short or 'IV-short' structural errors. } \examples{ \donttest{ diff --git a/man/var_hd.Rd b/man/var_hd.Rd index 89e448b..53ffde7 100644 --- a/man/var_hd.Rd +++ b/man/var_hd.Rd @@ -13,9 +13,8 @@ var_hd(var) long-from data.frame } \description{ -Estimate historical decomposition with contemporaneous impact -restrictions via Cholesky decomposition - taking the variable ordering -present in the VAR object. +Estimate historical decomposition for VARs with +either short or 'IV-short' structural errors. } \examples{ \donttest{ diff --git a/man/var_irf.Rd b/man/var_irf.Rd index a820c06..c27fa20 100644 --- a/man/var_irf.Rd +++ b/man/var_irf.Rd @@ -21,14 +21,19 @@ var_irf( \item{CI}{numeric vector: c(lower ci bound, upper ci bound)} -\item{bootstraps.num}{int: number of bootstraps} +\item{bootstrap.type}{string: bootstrapping technique to use ('auto', 'standard', or 'wild'); if auto then wild is used for IV or IV-short, else standard is used} + +\item{bootstrap.num}{int: number of bootstraps} + +\item{bootstrap.parallel}{boolean: create IRF draws in parallel} + +\item{bootstrap.cores}{int: number of cores to use in parallel processing; -1 detects and uses half the available cores} } \value{ -list object with elements \code{irfs}, \code{ci.lower}, and \code{ci.upper}; all elements are long-form data.frames +data.frame with columns \code{target}, \code{shock}, \code{horizon}, \code{response.lower}, \code{response}, \code{response.upper} } \description{ -Estimate impulse responses with contemporaneous impact restrictions via Cholesky decomposition - -taking the variable ordering present in the VAR object. +Estimate impulse response functions } \examples{ \donttest{ @@ -49,14 +54,14 @@ taking the variable ordering present in the VAR object. lag.ic = 'BIC', lag.max = 4) -# impulse response functions -var.irf = sovereign::var_irf(var) + # impulse response functions + var.irf = sovereign::var_irf(var) -# forecast error variance decomposition -var.fevd = sovereign::var_fevd(var) + # forecast error variance decomposition + var.fevd = sovereign::var_fevd(var) -# historical shock decomposition -var.hd = sovereign::var_hd(var) + # historical shock decomposition + var.hd = sovereign::var_hd(var) } diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf index 8267677..160ffd8 100644 Binary files a/tests/testthat/Rplots.pdf and b/tests/testthat/Rplots.pdf differ diff --git a/tests/testthat/test-lp.R b/tests/testthat/test-lp.R index 68306a6..ed7d69b 100644 --- a/tests/testthat/test-lp.R +++ b/tests/testthat/test-lp.R @@ -63,6 +63,7 @@ test_that("Local projection workflow", { rlp = RLP( data = Data, + regime = 'reg', p = 1, horizon = c(1:10), NW = FALSE, diff --git a/tests/testthat/test-threvhold_var.R b/tests/testthat/test-regime_var.R similarity index 96% rename from tests/testthat/test-threvhold_var.R rename to tests/testthat/test-regime_var.R index 233fea7..5ae46f9 100644 --- a/tests/testthat/test-threvhold_var.R +++ b/tests/testthat/test-regime_var.R @@ -3,7 +3,7 @@ test_that("threshold VAR workflow", { # simple time series AA = c(1:100) + rnorm(100) BB = c(1:100) + rnorm(100) - CC = AA + BB + rnorm(100) + CC = AA + BB^2 + rnorm(100) date = seq.Date(from = as.Date('2000-01-01'), by = 'month', length.out = 100) Data = data.frame(date = date, AA, BB, CC) @@ -43,7 +43,7 @@ test_that("threshold VAR workflow", { irf = rvar_irf( rvar, - bootstraps.num = 10, + bootstrap.num = 10, CI = c(0.05,0.95)) expect_true(is.data.frame(irf[[1]])) diff --git a/tests/testthat/test-structure.R b/tests/testthat/test-structure.R new file mode 100644 index 0000000..3262324 --- /dev/null +++ b/tests/testthat/test-structure.R @@ -0,0 +1,197 @@ +test_that("VAR workflow", { + + # simple time series + AA = c(1:100) + rnorm(100) + BB = c(1:100) + rnorm(100) + CC = AA + BB + rnorm(100) + inst = rbinom(100, size = 1, prob = .5) + date = seq.Date(from = as.Date('2015-01-01'), by = 'month', length.out = 100) + Data = data.frame(date, inst, AA, BB, CC) + + # estimate VAR without structure + var = + VAR( + data = Data, + p = 2, + horizon = 10, + freq = 'month', + type = 'both', + structure = NA) + + expect_true(is.list(var)) + expect_true(is.list(var$model)) + expect_true(is.list(var$forecasts)) + expect_true(is.list(var$residuals)) + + var.irf = IRF(var, bootstrap.num = 10) + + expect_true(is.data.frame(var.irf)) + + # estimate VAR with short-run restrictions + var.short = + VAR( + data = Data, + p = 2, + horizon = 10, + freq = 'month', + type = 'both', + structure = 'short') + + expect_true(is.list(var.short)) + expect_true(is.list(var.short$model)) + expect_true(is.list(var.short$forecasts)) + expect_true(is.list(var.short$residuals)) + + var.irf = IRF(var.short, bootstrap.num = 10) + var.hd = HD(var.short) + var.fevd = FEVD(var.short) + + expect_true(is.data.frame(var.irf)) + expect_true(is.data.frame(var.hd)) + expect_true(is.data.frame(var.fevd)) + + # estimate VAR with external instrument + var.iv = + VAR( + data = Data, + p = 2, + horizon = 10, + freq = 'month', + type = 'both', + structure = 'IV', + instrument = 'inst', + instrumented = 'AA') + + expect_true(is.list(var.iv)) + expect_true(is.list(var.iv$model)) + expect_true(is.list(var.iv$forecasts)) + expect_true(is.list(var.iv$residuals)) + + var.irf = IRF(var.iv, bootstrap.num = 10) + + expect_true(is.data.frame(var.irf)) + + + # estimate VAR with external instrument and short-run restrictions + var.iv_short = + VAR( + data = Data, + p = 2, + horizon = 10, + freq = 'month', + type = 'both', + structure = 'IV-short', + instrument = 'inst', + instrumented = 'AA') + + expect_true(is.list(var.iv_short)) + expect_true(is.list(var.iv_short$model)) + expect_true(is.list(var.iv_short$forecasts)) + expect_true(is.list(var.iv_short$residuals)) + + var.irf = IRF(var.iv_short, bootstrap.num = 10) + var.hd = HD(var.iv_short) + var.fevd = FEVD(var.iv_short) + + expect_true(is.data.frame(var.irf)) + expect_true(is.data.frame(var.hd)) + expect_true(is.data.frame(var.fevd)) + + # estimate RVAR without structure + + Data = dplyr::mutate(Data, reg = dplyr::if_else(AA > median(AA), 1, 0)) + + var = + RVAR( + data = Data, + p = 2, + regime = 'reg', + horizon = 10, + freq = 'month', + type = 'both', + structure = NA) + + expect_true(is.list(var)) + expect_true(is.list(var$model)) + expect_true(is.list(var$forecasts)) + expect_true(is.list(var$residuals)) + + var.irf = IRF(var, bootstrap.num = 10) + + expect_true(is.list(var.irf)) + + # estimate VAR with short-run restrictions + var.short = + RVAR( + data = Data, + p = 2, + regime = 'reg', + horizon = 10, + freq = 'month', + type = 'both', + structure = 'short') + + expect_true(is.list(var.short)) + expect_true(is.list(var.short$model)) + expect_true(is.list(var.short$forecasts)) + expect_true(is.list(var.short$residuals)) + + var.irf = IRF(var.short, bootstrap.num = 10) + var.hd = HD(var.short) + var.fevd = FEVD(var.short) + + expect_true(is.list(var.irf)) + expect_true(is.data.frame(var.hd)) + expect_true(is.list(var.fevd)) + + # estimate VAR with external instrument + var.iv = + RVAR( + data = Data, + p = 2, + regime = 'reg', + horizon = 10, + freq = 'month', + type = 'both', + structure = 'IV', + instrument = 'inst', + instrumented = 'AA') + + expect_true(is.list(var.iv)) + expect_true(is.list(var.iv$model)) + expect_true(is.list(var.iv$forecasts)) + expect_true(is.list(var.iv$residuals)) + + var.irf = IRF(var.iv, bootstrap.num = 10) + + expect_true(is.list(var.irf)) + + + # estimate VAR with external instrument and short-run restrictions + var.iv_short = + RVAR( + data = Data, + p = 2, + regime = 'reg', + horizon = 10, + freq = 'month', + type = 'both', + structure = 'IV-short', + instrument = 'inst', + instrumented = 'AA') + + expect_true(is.list(var.iv_short)) + expect_true(is.list(var.iv_short$model)) + expect_true(is.list(var.iv_short$forecasts)) + expect_true(is.list(var.iv_short$residuals)) + + var.irf = IRF(var.iv_short, bootstrap.num = 10) + var.hd = HD(var.iv_short) + var.fevd = FEVD(var.iv_short) + + expect_true(is.list(var.irf)) + expect_true(is.data.frame(var.hd)) + expect_true(is.list(var.fevd)) + + +}) diff --git a/tests/testthat/test-var.R b/tests/testthat/test-var.R index d1e149e..47e467c 100644 --- a/tests/testthat/test-var.R +++ b/tests/testthat/test-var.R @@ -8,7 +8,7 @@ test_that("VAR workflow", { Data = data.frame(date = date, AA, BB, CC) - # estimate VAR (with lag selection) + # estimate VAR (without lag selection) var = VAR( data = Data, @@ -51,7 +51,7 @@ test_that("VAR workflow", { irf = var_irf( var, - bootstraps.num = 10, + bootstrap.num = 10, CI = c(0.05,0.95)) expect_true(is.data.frame(irf)) diff --git a/vignettes/getting_started.Rmd b/vignettes/getting_started.Rmd index 8dd02bc..b7c99e3 100644 --- a/vignettes/getting_started.Rmd +++ b/vignettes/getting_started.Rmd @@ -160,7 +160,7 @@ Second we estimate the associate impulse response functions, bootstrapped confid var.irf = sovereign::var_irf( var, - bootstraps.num = 10, + bootstrap.num = 10, CI = c(0.05,0.95)) # plot IRF @@ -225,7 +225,7 @@ rvar.irf = sovereign::rvar_irf( rvar, horizon = 10, - bootstraps.num = 10, + bootstrap.num = 10, CI = c(0.05,0.95)) # plot IRF