devtools::install_github('tylerJPike/sovereign')

install.packages("sovereign")

# load packages
library(sovereign)         # analysis
library(dplyr)             # general cleaning
library(lubridate)         # date functions

#-------------------------------------------
# create data
#-------------------------------------------
# pull and prepare data from FRED
quantmod::getSymbols.FRED(
    c('UNRATE','INDPRO','GS10'), 
    env = globalenv())

Data = cbind(UNRATE, INDPRO, GS10)

Data = data.frame(Data, date = zoo::index(Data)) %>%
    filter(lubridate::year(date) >= 1990,
           lubridate::year(date) <  2020) %>% 
    na.omit()

# create a regime explicitly   
Data.threshold = Data %>%
    mutate(mp = if_else(GS10 > median(GS10), 1, 0))

#------------------------------------------
# learn regimes
#------------------------------------------
# assign regimes based on unsurpervised kmeans
#  (will not be used further)
regimes = 
    regimes(
        data = Data, 
        regime.n = 3, 
        engine = 'kmeans')

#------------------------------------------
# single-regime var
#------------------------------------------
# estimate VAR
# (using IC lag selection and short-term 
#  impact restrictions, i.e. calculate structurals
   errors via Cholesky decomposition)
var =
    VAR(
        data = Data,
        horizon = 10,
        freq = 'month',
        lag.ic = 'BIC',
        lag.max = 4)

# if impacts of the COVID shock are present,
#  in the user's data, then the user can 
#  implement the Lenza-Primiceri COVID volatility
#  correction for VARs: 
# var = covid_volatility_correction(var)

# plot forecasts
plot_forecast(var$forecasts[[1]])

# plot residuals
plot_error(var$residuals[[1]])

# estimate IRF
irf =
    IRF(
        var,
        bootstraps.num = 10,
        CI = c(0.05,0.95))

# plot IRF
plot_irf(irf)

# estimate forecast error variance decomposition
fevd =
    FEVD(
        var,
        horizon = 10)

# plot FEVD
plot_fevd(fevd)

# estimate historical decomposition
hd = HD(var)

# plot HD
plot_hd(hd)

#-------------------------------------------
# multi-regime var
#-------------------------------------------
# estimate multi-regime VAR
var =
    RVAR(
        data = Data.threshold,
        regime = 'mp',
        p = 1,
        horizon = 1,
        freq = 'month')

# estimate IRF
rvar.irf =
    IRF(
        rvar,
        horizon = 10,
        bootstraps.num = 10,
        CI = c(0.05,0.95))

# plot IRF
# regime 1: low interest rates
plot_irf(rvar.irf[[1]])
# regime 2: high interest rates
plot_irf(rvar.irf[[2]])

# estimate forecast error variance decomposition
rvar.fevd =
    FEVD(
        rvar,
        horizon = 10)

# plot FEVD
# regime 1: low interest rates
plot_fevd(rvar.fevd[[1]])
# regime 2: high interest rates
plot_fevd(rvar.fevd[[2]])

# estimate HD
rvar.hd = HD(rvar)

#-------------------------------------------
# single-regime local projections
#-------------------------------------------
# estimate single-regime forecasts 
#  (one or multiple horizons may be estimated)
lp = 
    LP(
        data = Data,
        p = 1,
        horizon = 1,
        freq = 'month')

# estimate single-regime IRF
lp.irf = IRF(lp)

# plot IRF
plot_irf(lp.irf)

#-------------------------------------------
# multi-regime local projections
#-------------------------------------------
# estimate multi-regime IRF
rlp = 
    RLP(
        data = Data,
        regime = 'mp',
        p = 1,
        horizon = 1,
        freq = 'month')

# estimate multi-regime IRF
rlp.irf = IRF(rlp)

# plot IRF
# regime 1: low interest rates
plot_irf(rlp.irf[[1]])
# regime 2: high interest rates
plot_irf(rlp.irf[[2]])