Nothing
R/lp.R
In sovereign: State-Dependent Empirical Analysis
Defines functions
RLP
RLP_estimate
LP
LP_estimate
Documented in
LP
RLP
#------------------------------------------
# Function to estimate LP
#------------------------------------------
# LP estimate LP models, forecasts and residuals
LP_estimate = function(
data, # data.frame, matrix, ts, xts, zoo: Endogenous regressors
p = 1, # int: lags
horizons = 1, # int: forecast horizons, can be a numeric vector with multiple horizons
freq = 'month', # string: frequency of data ('day', 'week', 'month', 'quarter', or 'year')
type = 'const', # string: type of deterministic terms to add ('none', 'const', 'trend', or'both')
# OLS-based IRF parameters
NW = FALSE, # Newey-West correction on variance-covariance matrix
NW_lags = NULL, # number of lags to use in Newey-West correction
NW_prewhite = NULL # prewhite option for Newey-West correction (see sandwich::NeweyWest function)
){
# 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 = NULL
# declare regressors
regressors = colnames(dplyr::select(data, -date))
# create regressors
Y = data.frame(data) %>%
n.lag(lags = p)
# remove date
Y = Y %>% dplyr::select(-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))}
# operate by horizon
outputs = as.list(horizons) %>%
purrr::map(.f = function(horizon){
### estimate coefficients ----------------------
models = as.list(regressors) %>%
purrr::map(.f = function(target){
# lead target
if(horizon > 1){
X = Y %>%
dplyr::select(dplyr::contains('.l'), target = target,
dplyr::contains('const'), dplyr::contains('trend')) %>%
dplyr::mutate(target = dplyr::lead(target, horizon-1))
}else{
X = Y %>% dplyr::select(dplyr::contains('.l'), target = target,
dplyr::contains('const'), dplyr::contains('trend'))
}
# estimate OLS
model = stats::lm(target ~ . -1, data = X)
# correct standard errors
if(NW == TRUE){
model = lmtest::coeftest(model, vcov = sandwich::NeweyWest(model, lag = NW_lags, prewhite = NW_prewhite))
}
# coefficients
c = broom::tidy(model) %>% dplyr::select(term, coef = estimate)
c$y = target
se = broom::tidy(model) %>% dplyr::select(term, std.error)
se$y = target
# return results
return(list(coef = c, se = se))
})
# 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
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)
# package for return
model = list(coef = coef, se = se, p = p, freq = freq, horizon = horizon)
### estimate forecasts -----------------------
# set design matrix
X = Y %>% dplyr::select(dplyr::contains('.l'),
dplyr::contains('const'), dplyr::contains('trend'))
# estimate i-step ahead forecast
forecast = as.matrix(X) %*% as.matrix(t(coef[,-1]))
colnames(forecast) = regressors
# add in dates
forecasts =
data.frame(
date = data$date,
forecast.date =
forecast_date(
forecast.date = data$date,
horizon = horizon-1,
freq = freq),
forecast)
### calculate residuals -----------------------
residuals = data.frame(forecasts)
residuals[,c(regressors)] = forecast[,c(regressors)] - data.frame(data)[, c(regressors)]
### return output --------------
return(
list(
model = model,
forecasts = forecasts,
residuals = residuals
)
)
})
# reorganize output
if(length(horizons) > 1){
model = purrr::map(outputs, .f = function(X){return(X$model)}); names(model) = paste0('H_',horizons)
forecasts = purrr::map(outputs, .f = function(X){return(X$forecasts)}); names(forecasts) = paste0('H_',horizons)
residuals = purrr::map(outputs, .f = function(X){return(X$residuals)}); names(residuals) = paste0('H_',horizons)
}else{
model = purrr::map(outputs, .f = function(X){return(X$model)})[[1]]
forecasts = outputs[[1]]$forecasts
residuals = outputs[[1]]$residuals
}
return(
list(
model = model,
data = data,
forecasts = forecasts,
residuals = residuals
)
)
}
#' Estimate local projections
#'
#' @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')
#' @param type string: type of deterministic terms to add ('none', 'const', 'trend', or 'both')
#' @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 NW boolean: Newey-West correction on variance-covariance matrix
#' @param NW_lags int: number of lags to use in Newey-West correction
#' @param NW_prewhite boolean: TRUE prewhite option for Newey-West correction (see sandwich::NeweyWest)
#'
#' @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
#'
#' @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{
#'
#' # 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)
#'
#' # local projection forecasts
#' lp =
#' sovereign::LP(
#' data = Data,
#' horizon = c(1:10),
#' lag.ic = 'AIC',
#' lag.max = 4,
#' type = 'both',
#' freq = 'month')
#'
#' # impulse response function
#' irf = sovereign::lp_irf(lp)
#'
#' }
#'
#' @export
LP = function(
data, # data.frame, matrix, ts, xts, zoo: Endogenous regressors
horizons = 1, # int: forecast horizons, can be a numeric vector with multiple 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')
# lag selection
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
# OLS-based IRF parameters
NW = FALSE, # Newey-West correction on variance-covariance matrix
NW_lags = NULL, # number of lags to use in Newey-West correction
NW_prewhite = NULL # prewhite option for Newey-West correction (see sandwich::NeweyWest function)
){
# 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(!type %in% c('none', 'const', 'trend', 'both')){
stop('type must be one of the following strings: "none", "const", "trend", "both"')
}
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(!freq %in% c('day','week','month','quarter','year')){
stop("freq must be one of the following strings: 'day','week','month','quarter','year'")
}
# estimate LP
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 =
LP_estimate(
data,
p = p,
horizons = horizons,
freq = freq,
type = type,
NW = NW,
NW_lags = NW_lags,
NW_prewhite = NW_prewhite
)
# calculate IC
if(length(horizons) > 1){
ic.score =
IC(
ic = lag.ic,
errors = model$residuals[[1]],
data = data,
p = p
)
}else{
ic.score =
IC(
ic = lag.ic,
errors = model$residuals,
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]]
class(model) = 'LP'
return(model)
}else{
model =
LP_estimate(
data,
p = p,
horizons = horizons,
freq = freq,
type = type,
NW = NW,
NW_lags = NW_lags,
NW_prewhite = NW_prewhite
)
class(model) = 'LP'
return(model)
}
}
#------------------------------------------
# Function to estimate threshold LP
#------------------------------------------
# estimate multi-regime LP models, forecasts, and residuals
RLP_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
horizons = 1, # int: forecast horizons, can be a numeric vector with multiple 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')
# OLS-based IRF parameters
NW = FALSE, # Newey-West correction on variance-covariance matrix
NW_lags = NULL, # number of lags to use in Newey-West correction
NW_prewhite = TRUE # prewhite option for Newey-West correction (see sandwich::NeweyWest function)
){
# 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
# 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 ----------------------
# iterate by regime
models = as.list(regimes) %>%
purrr::map(.f = function(regime.val){
# operate by horizon
outputs = as.list(horizons) %>%
purrr::map(.f = function(horizon){
models = as.list(regressors) %>%
purrr::map(.f = function(target){
# set and lead target
if(horizon > 1){
X = Y %>%
dplyr::select(dplyr::contains('.l'), target = target, regime = regime,
dplyr::contains('const'), dplyr::contains('trend')) %>%
dplyr::mutate(target = dplyr::lead(target, horizon-1)) %>%
dplyr::filter(regime == regime.val) %>%
dplyr::select(-regime)
}else{
X = Y %>%
dplyr::select(dplyr::contains('.l'), target = target, regime = regime,
dplyr::contains('const'), dplyr::contains('trend')) %>%
dplyr::filter(regime == regime.val) %>%
dplyr::select(-regime)
}
# estimate OLS
model = stats::lm(target ~ .-1, data = X)
# correct standard errors
if(NW == TRUE){
model = lmtest::coeftest(model, vcov = sandwich::NeweyWest(model, lag = NW_lags, prewhite = NW_prewhite))
}
# coefficients
c = broom::tidy(model) %>% dplyr::select(term, coef = estimate)
c$y = target
se = broom::tidy(model) %>% dplyr::select(term, std.error)
se$y = target
# return results
return(list(coef = c, se = se))
})
# 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
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)
# package for return
model = list(coef = coef, se = se, p = p, freq = freq, horizon = horizon)
return(model)
})
names(outputs) = paste0('H_',horizons)
return(outputs)
})
names(models) = paste0('regime_',regimes)
### estimate forecasts and residuals --------
# iterate by horizon
fr = as.list(horizons) %>%
purrr::map(.f = function(horizon){
# operate by regime
outputs = as.list(regimes) %>%
purrr::map(.f = function(regime.val){
# calculate forecasts
coef = models[[paste0('regime_', regime.val)]][[paste0('H_',horizon)]]$coef
X = Y %>%
dplyr::select(
dplyr::contains('.l'),
dplyr::contains('const'),
dplyr::contains('trend'))
forecast = as.matrix(X) %*% as.matrix(t(coef[,-1]))
colnames(forecast) = regressors
forecasts =
data.frame(
date = Y$date,
forecast.date = forecast_date(
forecast.date = Y$date,
horizon = horizon-1,
freq = freq),
forecast) %>%
dplyr::left_join(dplyr::select(Y, date, model.regime = regime), by = 'date')
# calculate residuals
residuals = data.frame(forecasts)
residuals[,c(regressors)] = forecast[,c(regressors)] - data.frame(data)[, c(regressors)]
# return output
residuals = residuals %>%
dplyr::filter(regime.val == model.regime)
forecasts = forecasts %>%
dplyr::filter(regime.val == model.regime)
return(
list(
forecasts = forecasts,
residuals = residuals
)
)
})
names(outputs) = paste0('regime_',regimes)
forecasts = purrr::map(outputs, .f = function(X){return(X$forecasts)}) %>%
purrr::reduce( dplyr::bind_rows) %>%
dplyr::arrange(date)
residuals = purrr::map(outputs, .f = function(X){return(X$residuals)})%>%
purrr::reduce( dplyr::bind_rows) %>%
dplyr::arrange(date)
return(list(forecasts = forecasts, residuals = residuals))
})
### Organize and return output -------
forecasts = purrr::map(fr, .f = function(X){return(X$forecasts)})
names(forecasts) = paste0('H_',horizons)
residuals = purrr::map(fr, .f = function(X){return(X$residuals)})
names(residuals) = paste0('H_',horizons)
if(length(horizons) == 1){
models = purrr::map(models, .f = function(X){return(X$H_1)})
forecasts = forecasts[[1]]
residuals = residuals[[1]]
}
return(
list(
models = models,
data = data,
forecasts = forecasts,
residuals = residuals,
regime = regime
)
)
}
#' 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')
#' @param type string: type of deterministic terms to add ('none', 'const', 'trend', or 'both')
#' @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 NW boolean: Newey-West correction on variance-covariance matrix
#' @param NW_lags int: number of lags to use in Newey-West correction
#' @param NW_prewhite boolean: TRUE prewhite option for Newey-West correction (see sandwich::NeweyWest)
#' @param regime string: name or regime assignment vector in the design matrix (data)
#' @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 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{
#'
#' # 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)
#' # add regime
#' Data = dplyr::mutate(Data, reg = dplyr::if_else(AA > median(AA), 1, 0))
#'
#' # local projection forecasts
#' rlp =
#' sovereign::RLP(
#' data = Data,
#' regime = 'reg',
#' horizon = c(1:10),
#' freq = 'month',
#' p = 1,
#' type = 'const',
#' NW = TRUE,
#' NW_lags = 1,
#' NW_prewhite = FALSE)
#'
#' # impulse response function
#' rirf = sovereign::rlp_irf(rlp)
#'
#' }
#'
#' @export
RLP = function(
data, # data.frame, matrix, ts, xts, zoo: Endogenous regressors
horizons = 1, # int: forecast horizons, can be a numeric vector with multiple 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')
# lag selection
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
# OLS-based IRF parameters
NW = FALSE, # Newey-West correction on variance-covariance matrix
NW_lags = NULL, # number of lags to use in Newey-West correction
NW_prewhite = NULL, # prewhite option for Newey-West correction (see sandwich::NeweyWest function)
# regimes
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)
){
# 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(!freq %in% c('day','week','month','quarter','year')){
stop("freq must be one of the following strings: 'day','week','month','quarter','year'")
}
if(!type %in% c('none', 'const', 'trend', 'both')){
stop('type must be one of the following strings: "none", "const", "trend", "both"')
}
# create regimes
if(is.null(regime)){
data =
regimes(
data,
regime.n = regime.n,
method = regime.method)
regime = 'regime'
}
# estimate LP
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 =
RLP_estimate(
data,
p = p,
regime = regime,
horizons = horizons,
freq = freq,
type = type,
NW = NW,
NW_lags = NW_lags,
NW_prewhite = NW_prewhite
)
# calculate IC
if(length(horizons) > 1){
ic.score =
IC(
ic = lag.ic,
errors = model$residuals[[1]],
data = data,
p = p
)
}else{
ic.score =
IC(
ic = lag.ic,
errors = model$residuals,
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]]
class(model) = 'RLP'
return(model)
}else{
model =
RLP_estimate(
data,
p = p,
regime = regime,
horizons = horizons,
freq = freq,
type = type,
NW = NW,
NW_lags = NW_lags,
NW_prewhite = NW_prewhite
)
class(model) = 'RLP'
return(model)
}
}
Try the sovereign package in your browser
Any scripts or data that you put into this service are public.
sovereign documentation built on Jan. 5, 2022, 1:08 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions RLP RLP_estimate LP LP_estimate
Documented in LP RLP
#------------------------------------------
# Function to estimate LP
#------------------------------------------
# LP estimate LP models, forecasts and residuals
LP_estimate = function(
data, # data.frame, matrix, ts, xts, zoo: Endogenous regressors
p = 1, # int: lags
horizons = 1, # int: forecast horizons, can be a numeric vector with multiple horizons
freq = 'month', # string: frequency of data ('day', 'week', 'month', 'quarter', or 'year')
type = 'const', # string: type of deterministic terms to add ('none', 'const', 'trend', or'both')
# OLS-based IRF parameters
NW = FALSE, # Newey-West correction on variance-covariance matrix
NW_lags = NULL, # number of lags to use in Newey-West correction
NW_prewhite = NULL # prewhite option for Newey-West correction (see sandwich::NeweyWest function)
){
# 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 = NULL
# declare regressors
regressors = colnames(dplyr::select(data, -date))
# create regressors
Y = data.frame(data) %>%
n.lag(lags = p)
# remove date
Y = Y %>% dplyr::select(-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))}
# operate by horizon
outputs = as.list(horizons) %>%
purrr::map(.f = function(horizon){
### estimate coefficients ----------------------
models = as.list(regressors) %>%
purrr::map(.f = function(target){
# lead target
if(horizon > 1){
X = Y %>%
dplyr::select(dplyr::contains('.l'), target = target,
dplyr::contains('const'), dplyr::contains('trend')) %>%
dplyr::mutate(target = dplyr::lead(target, horizon-1))
}else{
X = Y %>% dplyr::select(dplyr::contains('.l'), target = target,
dplyr::contains('const'), dplyr::contains('trend'))
}
# estimate OLS
model = stats::lm(target ~ . -1, data = X)
# correct standard errors
if(NW == TRUE){
model = lmtest::coeftest(model, vcov = sandwich::NeweyWest(model, lag = NW_lags, prewhite = NW_prewhite))
}
# coefficients
c = broom::tidy(model) %>% dplyr::select(term, coef = estimate)
c$y = target
se = broom::tidy(model) %>% dplyr::select(term, std.error)
se$y = target
# return results
return(list(coef = c, se = se))
})
# 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
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)
# package for return
model = list(coef = coef, se = se, p = p, freq = freq, horizon = horizon)
### estimate forecasts -----------------------
# set design matrix
X = Y %>% dplyr::select(dplyr::contains('.l'),
dplyr::contains('const'), dplyr::contains('trend'))
# estimate i-step ahead forecast
forecast = as.matrix(X) %*% as.matrix(t(coef[,-1]))
colnames(forecast) = regressors
# add in dates
forecasts =
data.frame(
date = data$date,
forecast.date =
forecast_date(
forecast.date = data$date,
horizon = horizon-1,
freq = freq),
forecast)
### calculate residuals -----------------------
residuals = data.frame(forecasts)
residuals[,c(regressors)] = forecast[,c(regressors)] - data.frame(data)[, c(regressors)]
### return output --------------
return(
list(
model = model,
forecasts = forecasts,
residuals = residuals
)
)
})
# reorganize output
if(length(horizons) > 1){
model = purrr::map(outputs, .f = function(X){return(X$model)}); names(model) = paste0('H_',horizons)
forecasts = purrr::map(outputs, .f = function(X){return(X$forecasts)}); names(forecasts) = paste0('H_',horizons)
residuals = purrr::map(outputs, .f = function(X){return(X$residuals)}); names(residuals) = paste0('H_',horizons)
}else{
model = purrr::map(outputs, .f = function(X){return(X$model)})[[1]]
forecasts = outputs[[1]]$forecasts
residuals = outputs[[1]]$residuals
}
return(
list(
model = model,
data = data,
forecasts = forecasts,
residuals = residuals
)
)
}
#' Estimate local projections
#'
#' @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')
#' @param type string: type of deterministic terms to add ('none', 'const', 'trend', or 'both')
#' @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 NW boolean: Newey-West correction on variance-covariance matrix
#' @param NW_lags int: number of lags to use in Newey-West correction
#' @param NW_prewhite boolean: TRUE prewhite option for Newey-West correction (see sandwich::NeweyWest)
#'
#' @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
#'
#' @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{
#'
#' # 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)
#'
#' # local projection forecasts
#' lp =
#' sovereign::LP(
#' data = Data,
#' horizon = c(1:10),
#' lag.ic = 'AIC',
#' lag.max = 4,
#' type = 'both',
#' freq = 'month')
#'
#' # impulse response function
#' irf = sovereign::lp_irf(lp)
#'
#' }
#'
#' @export
LP = function(
data, # data.frame, matrix, ts, xts, zoo: Endogenous regressors
horizons = 1, # int: forecast horizons, can be a numeric vector with multiple 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')
# lag selection
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
# OLS-based IRF parameters
NW = FALSE, # Newey-West correction on variance-covariance matrix
NW_lags = NULL, # number of lags to use in Newey-West correction
NW_prewhite = NULL # prewhite option for Newey-West correction (see sandwich::NeweyWest function)
){
# 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(!type %in% c('none', 'const', 'trend', 'both')){
stop('type must be one of the following strings: "none", "const", "trend", "both"')
}
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(!freq %in% c('day','week','month','quarter','year')){
stop("freq must be one of the following strings: 'day','week','month','quarter','year'")
}
# estimate LP
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 =
LP_estimate(
data,
p = p,
horizons = horizons,
freq = freq,
type = type,
NW = NW,
NW_lags = NW_lags,
NW_prewhite = NW_prewhite
)
# calculate IC
if(length(horizons) > 1){
ic.score =
IC(
ic = lag.ic,
errors = model$residuals[[1]],
data = data,
p = p
)
}else{
ic.score =
IC(
ic = lag.ic,
errors = model$residuals,
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]]
class(model) = 'LP'
return(model)
}else{
model =
LP_estimate(
data,
p = p,
horizons = horizons,
freq = freq,
type = type,
NW = NW,
NW_lags = NW_lags,
NW_prewhite = NW_prewhite
)
class(model) = 'LP'
return(model)
}
}
#------------------------------------------
# Function to estimate threshold LP
#------------------------------------------
# estimate multi-regime LP models, forecasts, and residuals
RLP_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
horizons = 1, # int: forecast horizons, can be a numeric vector with multiple 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')
# OLS-based IRF parameters
NW = FALSE, # Newey-West correction on variance-covariance matrix
NW_lags = NULL, # number of lags to use in Newey-West correction
NW_prewhite = TRUE # prewhite option for Newey-West correction (see sandwich::NeweyWest function)
){
# 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
# 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 ----------------------
# iterate by regime
models = as.list(regimes) %>%
purrr::map(.f = function(regime.val){
# operate by horizon
outputs = as.list(horizons) %>%
purrr::map(.f = function(horizon){
models = as.list(regressors) %>%
purrr::map(.f = function(target){
# set and lead target
if(horizon > 1){
X = Y %>%
dplyr::select(dplyr::contains('.l'), target = target, regime = regime,
dplyr::contains('const'), dplyr::contains('trend')) %>%
dplyr::mutate(target = dplyr::lead(target, horizon-1)) %>%
dplyr::filter(regime == regime.val) %>%
dplyr::select(-regime)
}else{
X = Y %>%
dplyr::select(dplyr::contains('.l'), target = target, regime = regime,
dplyr::contains('const'), dplyr::contains('trend')) %>%
dplyr::filter(regime == regime.val) %>%
dplyr::select(-regime)
}
# estimate OLS
model = stats::lm(target ~ .-1, data = X)
# correct standard errors
if(NW == TRUE){
model = lmtest::coeftest(model, vcov = sandwich::NeweyWest(model, lag = NW_lags, prewhite = NW_prewhite))
}
# coefficients
c = broom::tidy(model) %>% dplyr::select(term, coef = estimate)
c$y = target
se = broom::tidy(model) %>% dplyr::select(term, std.error)
se$y = target
# return results
return(list(coef = c, se = se))
})
# 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
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)
# package for return
model = list(coef = coef, se = se, p = p, freq = freq, horizon = horizon)
return(model)
})
names(outputs) = paste0('H_',horizons)
return(outputs)
})
names(models) = paste0('regime_',regimes)
### estimate forecasts and residuals --------
# iterate by horizon
fr = as.list(horizons) %>%
purrr::map(.f = function(horizon){
# operate by regime
outputs = as.list(regimes) %>%
purrr::map(.f = function(regime.val){
# calculate forecasts
coef = models[[paste0('regime_', regime.val)]][[paste0('H_',horizon)]]$coef
X = Y %>%
dplyr::select(
dplyr::contains('.l'),
dplyr::contains('const'),
dplyr::contains('trend'))
forecast = as.matrix(X) %*% as.matrix(t(coef[,-1]))
colnames(forecast) = regressors
forecasts =
data.frame(
date = Y$date,
forecast.date = forecast_date(
forecast.date = Y$date,
horizon = horizon-1,
freq = freq),
forecast) %>%
dplyr::left_join(dplyr::select(Y, date, model.regime = regime), by = 'date')
# calculate residuals
residuals = data.frame(forecasts)
residuals[,c(regressors)] = forecast[,c(regressors)] - data.frame(data)[, c(regressors)]
# return output
residuals = residuals %>%
dplyr::filter(regime.val == model.regime)
forecasts = forecasts %>%
dplyr::filter(regime.val == model.regime)
return(
list(
forecasts = forecasts,
residuals = residuals
)
)
})
names(outputs) = paste0('regime_',regimes)
forecasts = purrr::map(outputs, .f = function(X){return(X$forecasts)}) %>%
purrr::reduce( dplyr::bind_rows) %>%
dplyr::arrange(date)
residuals = purrr::map(outputs, .f = function(X){return(X$residuals)})%>%
purrr::reduce( dplyr::bind_rows) %>%
dplyr::arrange(date)
return(list(forecasts = forecasts, residuals = residuals))
})
### Organize and return output -------
forecasts = purrr::map(fr, .f = function(X){return(X$forecasts)})
names(forecasts) = paste0('H_',horizons)
residuals = purrr::map(fr, .f = function(X){return(X$residuals)})
names(residuals) = paste0('H_',horizons)
if(length(horizons) == 1){
models = purrr::map(models, .f = function(X){return(X$H_1)})
forecasts = forecasts[[1]]
residuals = residuals[[1]]
}
return(
list(
models = models,
data = data,
forecasts = forecasts,
residuals = residuals,
regime = regime
)
)
}
#' 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')
#' @param type string: type of deterministic terms to add ('none', 'const', 'trend', or 'both')
#' @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 NW boolean: Newey-West correction on variance-covariance matrix
#' @param NW_lags int: number of lags to use in Newey-West correction
#' @param NW_prewhite boolean: TRUE prewhite option for Newey-West correction (see sandwich::NeweyWest)
#' @param regime string: name or regime assignment vector in the design matrix (data)
#' @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 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{
#'
#' # 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)
#' # add regime
#' Data = dplyr::mutate(Data, reg = dplyr::if_else(AA > median(AA), 1, 0))
#'
#' # local projection forecasts
#' rlp =
#' sovereign::RLP(
#' data = Data,
#' regime = 'reg',
#' horizon = c(1:10),
#' freq = 'month',
#' p = 1,
#' type = 'const',
#' NW = TRUE,
#' NW_lags = 1,
#' NW_prewhite = FALSE)
#'
#' # impulse response function
#' rirf = sovereign::rlp_irf(rlp)
#'
#' }
#'
#' @export
RLP = function(
data, # data.frame, matrix, ts, xts, zoo: Endogenous regressors
horizons = 1, # int: forecast horizons, can be a numeric vector with multiple 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')
# lag selection
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
# OLS-based IRF parameters
NW = FALSE, # Newey-West correction on variance-covariance matrix
NW_lags = NULL, # number of lags to use in Newey-West correction
NW_prewhite = NULL, # prewhite option for Newey-West correction (see sandwich::NeweyWest function)
# regimes
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)
){
# 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(!freq %in% c('day','week','month','quarter','year')){
stop("freq must be one of the following strings: 'day','week','month','quarter','year'")
}
if(!type %in% c('none', 'const', 'trend', 'both')){
stop('type must be one of the following strings: "none", "const", "trend", "both"')
}
# create regimes
if(is.null(regime)){
data =
regimes(
data,
regime.n = regime.n,
method = regime.method)
regime = 'regime'
}
# estimate LP
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 =
RLP_estimate(
data,
p = p,
regime = regime,
horizons = horizons,
freq = freq,
type = type,
NW = NW,
NW_lags = NW_lags,
NW_prewhite = NW_prewhite
)
# calculate IC
if(length(horizons) > 1){
ic.score =
IC(
ic = lag.ic,
errors = model$residuals[[1]],
data = data,
p = p
)
}else{
ic.score =
IC(
ic = lag.ic,
errors = model$residuals,
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]]
class(model) = 'RLP'
return(model)
}else{
model =
RLP_estimate(
data,
p = p,
regime = regime,
horizons = horizons,
freq = freq,
type = type,
NW = NW,
NW_lags = NW_lags,
NW_prewhite = NW_prewhite
)
class(model) = 'RLP'
return(model)
}
}
Try the sovereign package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.